A multi-process coupling-based pollutant production and transport calculation method

Abstract

The application provides a pollution production and transport calculation method based on multi-process coupling, comprising: calculating rainfall driving data, underlying surface spatial attribute parameters and total mass of migratable pollutants of a calculation unit; performing hydrological simulation based on the parameters to obtain multiple hydrological transport paths, determine corresponding hydrodynamic parameters and sediment concentration; calculating a dynamic distribution coefficient based on the parameters and the concentration, distributing the total mass to each hydrological transport path to obtain initial pollutant load; performing transport calculation on each hydrological transport path respectively, performing secondary redistribution of the phase state of the pollutants across the paths based on the spatial variation of the sediment concentration in the transport to obtain output flux of each path; and integrating to obtain total pollutant output flux. The application can improve the accuracy and adaptability of pollutant simulation in a complex hydrological scenario.

Description

Technical Field

[0001] This invention belongs to the field of watershed hydrology and environmental simulation technology, specifically a pollutant generation and transport calculation method based on multi-process coupling. Background Technology

[0002] The generation and transport of pollutants in a watershed are influenced by a combination of factors, including rainfall runoff, soil erosion, and multi-path transport through complex media. Establishing computational models of physical processes can reveal the spatiotemporal distribution patterns of pollutants and quantify pollutant load characteristics under different hydraulic conditions. High-resolution simulations of physical mechanisms can provide data support for understanding the hydrodynamic evolution and material cycling patterns of a watershed.

[0003] Current mainstream watershed pollutant simulation schemes mostly employ empirical load functions or physical models based on steady-state assumptions. In existing technical schemes, the physicochemical partitioning of pollutants between the solid and liquid phases is usually set as a one-time partitioning operation based on fixed empirical constants; in the spatial transport stage, the multi-layered surface and subsurface water flow paths are often generalized into a single hybrid calculation channel, and in this process, the attachment state of pollutants is assumed to be an irreversible and constant binding process.

[0004] Furthermore, existing methods struggle to balance the accuracy of interface response and the fidelity of multi-media heterogeneous transport under complex hydrodynamic conditions when handling coupled simulations of watershed water, sediment, and materials at multiple temporal and spatial scales. Therefore, it is necessary to investigate a new computational method. Summary of the Invention

[0005] The purpose of this invention is to provide a pollutant production and transport calculation method based on multi-process coupling, in order to solve the above-mentioned problems in the prior art.

[0006] The technical solution, on the one hand, is a method for calculating pollutant generation and transport based on multi-process coupling, including:

[0007] Acquire rainfall-driven data, underlying surface spatial property parameters, and total mass of mobile pollutants for the computing unit;

[0008] Hydrological simulations were performed based on rainfall-driven data and underlying surface spatial attribute parameters to separate multiple hydrological transport paths and determine the hydrodynamic parameters and sediment concentrations corresponding to each hydrological transport path.

[0009] The dynamic distribution coefficient is calculated based on hydrodynamic parameters and sediment concentration, and the total mass of mobile pollutants is distributed to each hydrological transport path using the dynamic distribution coefficient to obtain the initial pollutant load corresponding to each hydrological transport path.

[0010] Transport calculations were performed on the initial pollutant loads corresponding to each hydrological transport path, and a secondary redistribution of pollutant phases across the path was performed based on the spatial changes in sediment concentration during the transport process to obtain the output flux corresponding to each hydrological transport path.

[0011] The output fluxes corresponding to each hydrological transport pathway are integrated to obtain the total pollutant output flux.

[0012] Beneficial effects: This invention can improve the accuracy and adaptability of pollutant simulation in complex hydrological scenarios. Specific effects will be described in detail later. Attached Figure Description

[0013] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings:

[0014] Figure 1 This paper presents a technical roadmap for constructing a pollutant generation and transport model in the upstream watershed of a reservoir and integrating it with the reservoir area model for optimization.

[0015] Figure 2 This is a flowchart of a pollutant generation and transport calculation method based on multi-process coupling according to an embodiment of this application;

[0016] Figure 3 This is a flowchart illustrating the calculation of dynamic distribution coefficients based on hydrodynamic parameters and sediment concentration, according to an embodiment of this application.

[0017] Figure 4 This is a flowchart illustrating how to determine the sediment concentration corresponding to each hydrological transport path, according to an embodiment of this application.

[0018] Figure 5 This is a flowchart illustrating an embodiment of the present application for obtaining rainfall-driven data and underlying surface spatial attribute parameters of a computing unit;

[0019] Figure 6 This is a flowchart illustrating the calculation of dynamic allocation coefficients according to an embodiment of this application; Detailed Implementation

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

[0021] It should be noted that the terms "first," "second," etc., in the specification and accompanying drawings of this invention are used to distinguish similar objects and are not necessarily used to describe a predetermined order or sequence. It should be understood that such data can be interchanged where appropriate so that embodiments of the invention described herein can be implemented in sequences other than those illustrated or described herein. Furthermore, the terms "including" and "having," and any variations thereof, are intended to cover non-exclusive inclusion. For example, a process, method, system, product, or apparatus that includes a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0022] In light of the above issues, the applicant conducted a search and analysis and found that:

[0023] When faced with high-intensity rainstorms or violent hydrodynamic field fluctuations, calculation methods based on static parameter solidification and single-channel hybrid generalization are insufficient to characterize the physicochemical evolution process between complex media interfaces, resulting in significant calculation deviations in the output mass flux in terms of temporal peak and spatial distribution.

[0024] Combination Figures 1 to 6 To address the above problems, the present invention will be specifically described through the following embodiments.

[0025] Furthermore, some nouns or terms that appear in the description of the embodiments of this application shall be interpreted as follows:

[0026] In this application, the target pollutant is a type of migratory pollutant.

[0027] The quadratic phase redistribution network can also be called the phase quadratic redistribution network.

[0028] Hydraulic characteristic parameters, also known as hydrodynamic parameters;

[0029] The computational unit represents the basic spatial simulation grid or sub-basin structure divided during the spatial discretization process of the watershed.

[0030] Rainfall-driven data represent the time series of rainfall in the study area during the target simulation period.

[0031] The underlying surface spatial attribute parameters include slope and slope length extracted from the digital elevation model, saturated permeability coefficient, field water holding capacity and equivalent soil thickness extracted from soil spatial distribution attributes, and Manning roughness coefficient determined according to land use type.

[0032] Total mass of mobile pollutants characterizes the total physical mass of target pollutants mobilized from the soil surface into the soil and water system under the scouring action of the current rainfall event.

[0033] Multiple hydrological transport pathways are classified into surface runoff pathways located above the surface and subsurface runoff pathways and underground runoff pathways located below the surface.

[0034] Hydrodynamic parameters, including the flow velocity and equivalent water depth for each path.

[0035] The dynamic allocation coefficient is not a single static empirical constant, but a complex function parameter that evolves continuously with sediment concentration and hydrodynamic parameters in the spatiotemporal grid.

[0036] This embodiment provides a pollutant generation and transport calculation method based on multi-process coupling, the method comprising:

[0037] Step 101: Obtain rainfall-driven data, underlying surface spatial attribute parameters, and total mass of mobile pollutants from the computing unit;

[0038] When extracting rainfall-driven data at each computation time step, the system strictly enforces historical backtracking time window constraints. Accordingly, only meteorological observation time series within the current time step and the historical backtracking window are extracted to avoid computational distortion caused by referencing data from time periods that have not yet occurred during time iterations.

[0039] As an optional implementation method, the spatial attribute parameters of the underlying surface can also be dynamically extracted and updated by combining multi-source remote sensing images to reflect the spatiotemporal variation characteristics of vegetation cover in different seasons.

[0040] Step 102: Based on rainfall-driven data and underlying surface spatial attribute parameters, perform hydrological simulation to separate multiple hydrological transport paths and determine the hydrodynamic parameters and sediment concentrations corresponding to each hydrological transport path;

[0041] Accordingly, a three-source runoff generation mechanism combined with water balance is adopted to decompose the runoff;

[0042] For each separated hydrological transport path, the corresponding hydrodynamic parameters are calculated by combining the surface friction resistance and subsurface infiltration characteristics of the corresponding path. For example, for the slope runoff path, the Manning formula is used to calculate the slope runoff velocity and slope runoff depth; for the interflow path, the interflow velocity is estimated based on Darcy's law combined with the soil saturated permeability coefficient and hydraulic gradient.

[0043] Furthermore, after determining the hydrodynamic parameters, the mass of eroded material is calculated using a soil erosion model based on the basic principles of hydrodynamics. Then, the sediment concentration of the slope runoff path at the corresponding time step is calculated based on the water flow volume, which can provide environmental driving variables for the phase evolution of pollutants.

[0044] Step 103: Calculate the dynamic distribution coefficient based on hydrodynamic parameters and sediment concentration, and use the dynamic distribution coefficient to distribute the total mass of mobile pollutants to each hydrological transport path to obtain the initial pollutant load corresponding to each hydrological transport path.

[0045] Furthermore, after obtaining the dynamic allocation coefficient, based on the physical adsorption potential energy characterized by the allocation coefficient, the total mass of migratable pollutants is decomposed into the mass of solid pollutants in the adsorption state and the mass of liquid pollutants in the dissolved state according to the corresponding proportion.

[0046] Next, based on the different physical and dynamic capabilities of each hydrological transport path, the initial space station deployment and allocation of mass payloads will be carried out.

[0047] As an example, the mass of dissolved liquid pollutants enters each hydrological transport pathway according to the proportion of water flow in each pathway.

[0048] The mass of solid pollutants in the adsorbed state is dynamically weighted and allocated according to the ability of each path to transport sediment particles.

[0049] After the allocation calculation is completed, the sum of the allocation results constitutes the initial pollutant load of each hydrological transport path at the starting point of the current calculation time step.

[0050] Step 104: Perform transport calculations on the initial pollutant loads corresponding to each hydrological transport path, and perform secondary redistribution of pollutant phases across the path based on the spatial changes in sediment concentration during the transport process to obtain the output flux corresponding to each hydrological transport path.

[0051] For each hydrological transport pathway, perform independent and differentiated transport calculations.

[0052] As an example, in the subsurface flow path and underground runoff path, pollutants travel through the pores of the soil matrix and are subject to physical adsorption delay and biochemical degradation by matrix particles. Their mass output exhibits a smooth and lagging fluctuation characteristic in the time series.

[0053] For surface runoff pathways, there exists a cross-physical interface co-evolution that differs from traditional steady distribution theories. Correspondingly, as sediment particles are transported downstream along the slope with the water flow, some particles detach from the main flow due to gravity settling, resulting in a physical decrease in sediment concentration along the spatial migration distance. With the decrease in sediment concentration, the total adsorption capacity of the water-sediment mixing system for pollutants experiences a net decrease.

[0054] At this point, the phase equilibrium established at the original location is disrupted. Some adsorbed pollutants originally attached to the solid sediment detach from the sediment carrier and enter the dissolved water phase, resulting in desorption and release. The pollutants released through desorption are no longer carried by the surface sediment, but are transferred across the spatial interface to the interflow path through infiltration mechanisms, completing a secondary phase redistribution.

[0055] After undergoing differentiated evolution along the path and cross-path transfer, each path independently outputs the output flux to the exit section of the target computing node.

[0056] Step 105: Integrate the output fluxes corresponding to each hydrological transport path to obtain the total pollutant output flux.

[0057] Accordingly, the output fluxes of the slope runoff path, interflow path, and underground runoff path at the same time point are superimposed and merged to obtain the total pollutant output flux representing the absolute pollution intensity of the target calculation unit.

[0058] Optionally, after obtaining the total pollutant output flux, a time definite integral is performed on the total pollutant output flux over the entire duration of the target rainfall event.

[0059] Alternatively, by performing an integration operation in the time dimension, the instantaneous flux physical quantity, measured by the rate of change over time, can be converted into a physical quantity representing the total pollutant load of the cumulative environmental output of a single event. This can be used as an input parameter and directly connected to the downstream water quality early warning platform or the water environment carrying capacity safety assessment calculation network.

[0060] In conjunction with the preceding embodiments, multiple hydrological transport pathways include:

[0061] Slope runoff pathways, interflow pathways, and subsurface runoff pathways;

[0062] Accordingly, hydrodynamic parameters include the slope runoff velocity along the slope runoff path.

[0063] In this step, the boundaries of the transport space of the hydrological model are defined.

[0064] That is, slope runoff paths above the ground surface have the physical characteristics of high flow velocity and strong sediment carrying capacity.

[0065] The flow path of soil in shallow soil is constrained by soil porosity, and its flow velocity exhibits sluggish characteristics and weak sediment carrying capacity.

[0066] The underground runoff path located in deep aquifers has the longest residence time and does not have the capacity to carry sediment due to filtration.

[0067] Based on this, the hydrodynamic parameters were specified as slope runoff velocity, providing a fundamental variable for characterizing the dynamic adsorption equilibrium time of the water-sediment mixing system.

[0068] In some embodiments, the dynamic distribution coefficient is calculated based on hydrodynamic parameters and sediment concentration, including:

[0069] Optionally, in each calculation time step, the rate of change of sediment concentration at the current time step relative to the previous time step is calculated;

[0070] In the actual implementation process, since the sediment concentration changes slowly during the non-peak period of rainfall events, if the entire set of thermodynamic allocation calculations is performed in every small time step, a large amount of ineffective computing power will be consumed.

[0071] Therefore, a time-triggered mechanism based on the concentration change rate is constructed to calculate this change rate, and the corresponding calculation formula is as follows:

[0072] r _c (t)=abs(C _s (t)-C _s (t-1)) / C _s (t-1);

[0073] Where, r _c (t) represents the rate of change of concentration at the current time step t, abs() represents the operation of taking the absolute value, C _s (t) represents the sediment concentration at the current time step, C _s (t-1) represents the sediment concentration at the previous time step t-1.

[0074] Optionally, when the concentration change rate is greater than the preset time trigger threshold, the dynamic distribution coefficient is calculated based on the pre-configured reference sediment concentration, reference flow velocity and benchmark distribution coefficient.

[0075] When the concentration change rate is not greater than the time trigger threshold, the dynamic allocation coefficient of the previous time step is used as the dynamic allocation coefficient of the current time step.

[0076] Among them, when calculating the dynamic distribution coefficient, a particle concentration effect correction is introduced based on the negative power function relationship of the ratio of sediment concentration to reference sediment concentration, so as to characterize the decay characteristics of the distribution coefficient as the sediment concentration increases.

[0077] A hydrodynamic correction is introduced based on the saturation inhibition function relationship of the ratio of slope runoff velocity to reference velocity to characterize the attenuation characteristics of the distribution coefficient as the velocity increases.

[0078] Accordingly, the preset time trigger threshold can be set to a specific value between 0.05 and 0.15.

[0079] When the calculated concentration change rate exceeds the set value, the coefficient update calculation for this step is triggered; otherwise, the allocation coefficient from the previous step is directly used. During the calculation, based on the improved isothermal adsorption theory, a fundamental dynamic allocation coefficient incorporating particle concentration effects and hydrodynamic corrections is constructed. Its calculation formula is as follows:

[0080] K _d_base (t)=K _d0 ×(C _s (t) / C _ref ) -β ×[1 / (1+α×v _s (t) / v _ref )];

[0081] Among them, K _d_base (t) is the basic dynamic allocation coefficient, K _d0 For the pre-configured baseline allocation coefficient, C _s (t) represents the sediment concentration at the current time step, C _ref The pre-configured reference sediment concentration is given by β, which is the particle concentration effect index, and α is the hydrodynamic correction factor. _s (t) represents the slope runoff velocity at the current time step, v _ref This is a pre-configured reference flow rate.

[0082] In other words, this step introduces the particle concentration effect from the field of environmental chemistry into hydrological calculations. Correspondingly, as the concentration of suspended sediment in water increases, the competitive adsorption of dissolved pollutants by particles intensifies, while the particle surface may be obscured by colloidal microparticles, leading to a decrease in the effective adsorption area per unit mass of sediment. The negative power function term maps this physical phenomenon.

[0083] Meanwhile, increased flow velocity shortens the contact time at the water-sand interface, making it difficult for the adsorption reaction to reach static equilibrium. The saturation inhibition function term in the formula ensures that the correction term approaches 0 when the flow velocity approaches infinity.

[0084] In some scenarios, such as extreme rainstorms with sufficient computing power or unstable sediment concentrations, a fixed-time-interval forced update mechanism can be used instead of the concentration change rate trigger mechanism. For example, a pre-defined number of time steps can be used to force a recalculation of the dynamically allocated coefficients to prevent error accumulation.

[0085] Furthermore, the steps for calculating the dynamic allocation coefficients can also be as follows:

[0086] To obtain real-time water temperature;

[0087] Based on the water temperature and a pre-configured reference temperature, a temperature correction term is constructed using the Arrhenius equation;

[0088] The dynamic allocation coefficient is obtained by using a temperature correction term to perform thermodynamic attenuation or enhancement correction on the allocation coefficient results that are superimposed with particle concentration effect correction and hydrodynamic correction.

[0089] Correspondingly, the physical adsorption process of most pollutants is an exothermic reaction; based on the laws of thermodynamics, an increase in the temperature of the water environment will cause the adsorption equilibrium to shift towards desorption.

[0090] By acquiring real-time water temperature data and substituting it into the temperature correction term derived from the Arrhenius equation, the basic dynamic distribution coefficients are multiplied to obtain the correction formula:

[0091] K _d (t)=K _d_base (t)×exp[(E _a / R)×(1 / T(t)-1 / T _ref )];

[0092] Among them, K _d (t) represents the final dynamic allocation coefficient, K _d_base (t) is the basic dynamic allocation coefficient, exp is an exponential function with the natural constant as the base, E _a Let T(t) be the apparent activation energy of the adsorption reaction, R be the universal gas constant, and T(t) be the measured water temperature at the current time step. _ref This is the pre-configured reference temperature.

[0093] With this correction, when the actual water temperature is higher than the reference temperature, the calculated result in the exponential term is less than 0, which in turn makes the exponential function output a decay coefficient less than 1, thus achieving decay correction in the thermodynamic direction.

[0094] In a further embodiment, the total mass of migratable pollutants is allocated to each hydrological transport path using a dynamic allocation coefficient to obtain the initial pollutant load corresponding to each hydrological transport path, including:

[0095] Accordingly, the dynamic distribution coefficient and sediment concentration are combined for dimensional unification, and the adsorbed state mass ratio and dissolved state mass ratio are calculated to decompose the total mass of migratable pollutants into the adsorbed pollutant mass and the dissolved pollutant mass.

[0096] Furthermore, by introducing a unit conversion constant, the proportion of the target pollutant in the solid phase and the proportion in the liquid phase are calculated using the following formula:

[0097] f _ads (t)=[K _d (t)×C _s (t)×μ] / [1+K _d (t)×C _s [(t)×μ];

[0098] f _dis (t)=1-f _ads (t);

[0099] Among them, f _ads (t) represents the mass ratio of adsorbed states, K _d (t) is the dynamic allocation coefficient, C _s (t) represents the sediment concentration, μ is the unit conversion constant used to match the dynamic distribution coefficient with the volume unit of sediment concentration, and f _dis (t) represents the mass ratio of dissolved substances.

[0100] Next, the total mass of the acquired migratory pollutants will be multiplied by two proportional parameters to calculate the mass of adsorbed pollutants and dissolved pollutants entering the water-soil mixed environment at the current time step, thereby completing the initial isolation of the chemical phase.

[0101] Accordingly, based on hydrodynamic parameters, the sediment carrying capacity index and runoff of each hydrological transport path are determined, and the runoff of each hydrological transport path separated in the hydrological simulation is extracted.

[0102] In this step, the dynamic characteristic data of each path are extracted;

[0103] Among them, runoff can be derived from the decomposition results of water balance;

[0104] The sediment carrying capacity index needs to be further quantified based on the actual water depth and flow velocity characteristics of each path.

[0105] Among them, the sediment carrying capacity indicators for each hydrological transport pathway are determined, including:

[0106] Furthermore, the sediment carrying capacity index of the underground runoff path is set to zero;

[0107] Based on the power function relationship between flow velocity and water depth for both slope runoff and interflow paths, the corresponding sediment carrying capacity index is calculated.

[0108] Among them, the sediment carrying capacity index of the interflow path is also reduced and corrected using the measured effective porosity of the soil.

[0109] Alternatively, the effective porosity of soil can be derived from the spatial properties of the underlying surface.

[0110] Correspondingly, since underground aquifers lack sediment transport capabilities, their sediment carrying capacity is forcibly set at 0. The carrying capacity of surface and shallow soils is calculated using the following formula:

[0111] ω _s (t)=v _s (t) a1×h _s (t) a2 ;

[0112] ω _i (t)=v _i (t) a1 ×h _i (t) a2 ×θ _e ;

[0113] In the formula, ω _s (t) represents the sediment carrying capacity index of the slope runoff path, v _s (t) represents the surface runoff velocity, h _s (t) represents the slope runoff depth, and a1 and a2 are preset empirical indices representing the contribution weights of flow velocity and water depth, respectively. ω _i (t) represents the sediment carrying capacity index of the interflow path, v _i (t) represents the interflow velocity, h _i (t) represents the equivalent water depth of the interflow, θ _e The effective porosity of the soil.

[0114] It should be understood that introducing the effective porosity of the soil for multiplicative reduction can reflect the pore filtration constraint encountered by the interflow when it penetrates the soil matrix.

[0115] Accordingly, the mass of adsorbed pollutants is allocated according to the normalized weight of the sediment carrying capacity index of each hydrological transport path, and the mass of dissolved pollutants is allocated according to the proportion of the runoff of each hydrological transport path. The allocation results are integrated to obtain the initial pollutant load corresponding to each hydrological transport path.

[0116] In this step, the adsorbed state weights assigned to the slope runoff path are calculated, specifically by calculating ω. _s The ratio of (t) to the sum of the total carrying capacity index of the slope and the soil flow is obtained. Multiplying this weight by the mass of adsorbed pollutants yields the adsorbed load obtained from the slope path.

[0117] Similarly, the ratio of the runoff from each path to the total production flow is used as the water volume weight to allocate the mass of dissolved pollutants.

[0118] Based on this, the adsorbed pollutant mass and dissolved pollutant mass allocated to each path are added together to construct the initial pollutant load for the hydrological transport path to enter the next stage of spatial evolution.

[0119] Optionally, if the accuracy of the sensing equipment is insufficient to support the continuous calculation of flow velocity and water depth for each path, the unit width flow rate can be used instead of the composite power function of velocity and water depth. That is, the unit width flow rate ratio can be directly used as a simplified estimation basis for the sediment carrying capacity index.

[0120] As one possible implementation, a secondary redistribution of pollutant phases across the transport pathway is performed based on the spatial variation of sediment concentration during transport, specifically:

[0121] The slope runoff path is discretized into multiple spatial nodes along the direction of water flow. The decrease in sediment concentration along the path due to sedimentation is calculated to obtain the local sediment concentration at each spatial node.

[0122] At each spatial node, the spatial concentration change rate between adjacent spatial nodes is calculated based on the local sediment concentration at adjacent spatial nodes.

[0123] When the rate of change of spatial concentration is greater than the preset spatial trigger threshold, a secondary phase redistribution operation is triggered at the current spatial node.

[0124] When the rate of change of spatial concentration is not greater than the spatial trigger threshold, the pollutant load of the current spatial node is directly transferred to the downstream spatial node along the slope runoff path.

[0125] Accordingly, this step involves spatial gridding of surface water and sediment transport;

[0126] That is, the slope runoff path with continuous physical length is divided into multiple discrete spatial nodes along the direction of water flow, and the distance between adjacent nodes is determined by dividing the total slope length by the preset total number of nodes.

[0127] As sediment particles move towards the foot of the slope with the water flow, larger particles are deposited first due to gravity settling and the interception effect of surface vegetation roughness, resulting in a monotonically decreasing concentration of suspended sediment in the water flow along the way.

[0128] To address this sediment attenuation process along the route, sediment concentration values ​​at two adjacent spatial nodes are extracted, and the spatial concentration change rate is calculated. The specific calculation formula is as follows:

[0129] r _x (x _k ,t)=abs(C _s (x _k ,t)-C _s (x _k-1 ,t)) / C _s (x _k-1 ,t);

[0130] Where, r _x (x _k ,t) represents the current time step at spatial node x. _k With the previous spatial node x _k-1 The spatial concentration change rate between them, abs() characterizes the absolute value operation, C_s (x _k ,t) represents the spatial node x reached at the current time step. _k The local sediment concentration at C _s (x _k-1 ,t) represents the current time step at the previous spatial node x. _k-1 Local sediment concentration at the location.

[0131] In this case, after obtaining the spatial concentration change rate, it is compared with the preset spatial trigger threshold for determination;

[0132] When the value of the rate of change is greater than the preset spatial trigger threshold, it indicates that the concentration decay has broken the original chemical balance of the water-sediment system, thereby triggering the secondary redistribution of phases.

[0133] It should be understood that this triggering mechanism can eliminate redundant iterative calculations caused by minor fluctuations in concentration.

[0134] In some scenarios, phase redistribution operations include:

[0135] Based on the local sediment concentration at the current spatial node, the local dynamic distribution coefficient is recalculated, and the local adsorbed mass ratio is updated using the local dynamic distribution coefficient.

[0136] When the updated local adsorbed state mass ratio is lower than the adsorbed state mass ratio of the previous spatial node, it indicates a decrease in the environmental adsorption capacity. Based on the difference in the decrease in the adsorbed state mass ratio and the total mass of remaining pollutants that have arrived at the current spatial node, the amount of adsorbed state desorption and release is calculated.

[0137] The dissolved pollutants corresponding to the adsorbed desorption release amount are deducted from the slope runoff path and added to the pollutant load of the interflow path across the path, thus outputting the output flux of the slope runoff path.

[0138] When the triggering condition is met, the local sediment concentration after the decay of the current spatial node is substituted into the negative power function relationship to update the local dynamic distribution coefficient, and then the new local adsorbed state mass ratio is recalculated.

[0139] The new ratio value is compared by subtracting it from the original ratio value of the upstream node.

[0140] In another embodiment, the phase redistribution operation can also be performed by the following steps:

[0141] Based on the local sediment concentration of the current spatial node, the local dynamic distribution coefficient is recalculated, and the local adsorbed state mass ratio of the current spatial node is calculated using the local dynamic distribution coefficient and the local sediment concentration.

[0142] When the mass ratio of locally adsorbed states is lower than that of the adsorbed states at the previous spatial node, it indicates a decrease in the adsorption capacity of the environment.

[0143] The amount of adsorbed desorption and release is calculated based on the difference in the decrease of the adsorbed mass ratio and the total mass of remaining pollutants that have reached the current space node.

[0144] The adsorbed state mass ratio of the previous spatial node is given by the initial ratio determined by the dynamic allocation coefficient when it is allocated in the source region at the first spatial node, and by the local adsorbed state mass ratio calculated at the previous spatial node at subsequent spatial nodes.

[0145] The total mass of remaining pollutants is the remaining mass of the initial pollutant load after being gradually transported and deducted from each space node and arriving at the current space node;

[0146] When the local adsorbed state mass ratio is not lower than the adsorbed state mass ratio of the previous spatial node, no desorption and release operation is performed at the current spatial node, and the current pollutant load continues to be transferred to the downstream spatial node.

[0147] The dissolved pollutants corresponding to the adsorbed desorption release amount are deducted from the slope runoff path and added to the initial pollutant load of the interflow path across the path, thus outputting the output flux of the slope runoff path.

[0148] In some scenarios, the reference sediment concentration, baseline distribution coefficient, and particle concentration effect index are set to preset values, and no other correction conditions are introduced for the time being.

[0149] Among them, the upstream nodes have higher sediment concentrations and corresponding lower dynamic distribution coefficients, and the overall adsorption capacity of the water-sediment system remains at a high level.

[0150] During the downstream transport of water, sediment settles under the influence of gravity, the sediment concentration at nodes decreases, and the local dynamic distribution coefficient increases accordingly.

[0151] The increase in the distribution coefficient cannot offset the decrease in sediment concentration, resulting in a net decrease in the overall adsorption capacity of downstream nodes compared to upstream nodes. Consequently, the proportion of adsorbed pollutants decreases, which in turn triggers the desorption and release of pollutants from the sediment surface.

[0152] Furthermore, due to the net decrease in adsorption capacity, excess physical mass must detach from the surface of the solid particles and enter the liquid phase. This detachment process is quantified by calculating the percentage decrease and multiplying it by the current node's mass base; the formula can be expressed as:

[0153] ΔM _des (x _k ,t)=[f _ads_old (x _k-1 ,t)-f _ads_new (x_k ,t)]×M _s_remain (x _k ,t);

[0154] Where, ΔM _des (x _k ,t) represents the space node x _k The amount of adsorbed desorption released at the site, f _ads_old (x _k-1 ,t) is the previous spatial node x _k-1 The mass ratio of adsorbed states at f _ads_new (x _k ,t) represents the current spatial node x _k The updated local adsorbed state mass ratio, M _s_remain (x _k ,t) represents the current reached spatial node x _k The total mass of remaining pollutants at the location.

[0155] Accordingly, after calculating the amount of desorption and release in the adsorbed state, its physical properties are changed from the adsorbed state to the dissolved state;

[0156] This portion of dissolved mass, along with surface infiltration water, crosses the surface physical interface and enters the subsurface soil flow path.

[0157] Therefore, it is necessary to subtract the current mass load from the slope runoff path and then add it to the interflow path. After recursive subtraction and attenuation calculations at all spatial nodes, the remaining mass reaching the terminal node at the toe of the slope is calculated. Combined with the time step size, the final flux of the slope runoff path is converted and output.

[0158] In other scenarios, in areas lacking high-resolution spatial elevation data, the continuous one-dimensional discrete recursive process can be ignored. Instead, the overall single desorption amount assessment can be performed directly based on the initial sediment concentration at the top of the slope and the terminal sediment concentration observed at the bottom of the slope. This improves the computational convergence speed while sacrificing some spatial evolution accuracy.

[0159] In some embodiments, transport calculations are performed on the initial pollutant loads corresponding to each hydrological transport path. For interflow paths, this includes:

[0160] The initial pollutant load corresponding to the interflow path is combined with the additional adsorbed desorption release to obtain the total pollutant input load of the interflow path.

[0161] The transport of total pollutant input load is calculated using a pre-built transfer function model, which includes a hysteresis factor to characterize adsorption and desorption in the soil matrix and a first-order decay coefficient to characterize degradation.

[0162] The output flux of the interflow path after hysteresis and concentration decay is output based on the transfer function model.

[0163] Correspondingly, the interflow path not only receives the initial pollutant load from the source area, but also receives the adsorbed desorption and release from the infiltration along the way. The sum of the two constitutes the total physical mass entering the soil transport channel.

[0164] Unlike the instantaneous characteristics of surface flow, interflow in porous media exhibits a strong reactive delay. The transfer function model is used to characterize this delay and loss process, and its actual calculation formula is as follows:

[0165] F _i (t)=[M _i_total (t-τ _i ) / Δt]×exp[-λ×D _soil / v _i (t)];

[0166] Among them, F _i (t) represents the output flux of the midstream path at the current time step, M _i_total (t-τ _i ) represents the time elapsed after the lag τ _i The total pollutant input load after delay, Δt represents the time step size, exp represents the natural constant exponential function, λ represents the preset first-order decay coefficient characterizing the biodegradation intensity, and D _soil Indicates the equivalent thickness of the soil layer, v _i (t) represents the current velocity of the interflow in the soil at the current time step.

[0167] Alternatively, D _soil / v _i (t) characterizes the equivalent contact time between the interflow and the soil matrix during its passage through the soil layer, and is used to measure the cumulative intensity of degradation.

[0168] Furthermore, the hysteresis delay time τ _i It is directly determined by the hysteresis factor. The specific value of the hysteresis factor can be calculated by extracting the soil dry bulk density and effective porosity; the calculation formula is as follows:

[0169] R _f =1+(ρ _b ×K _D_soil ) / θ _e ;

[0170] Among them, R _f ρ represents the hysteresis factor. _b This indicates the obtained dry bulk density of the soil, K. _D_soil This represents the static distribution coefficient of the target pollutant on the target soil matrix.

[0171] The above calculations correspond to the following pattern: the slower the flow velocity of the interflow, the longer the pollutants remain underground, and the more they are degraded and lost by soil microorganisms.

[0172] In other embodiments, transport calculations are performed on the initial pollutant loads corresponding to each hydrological transport path, including for subsurface runoff paths:

[0173] Obtain the pre-configured background groundwater concentration and groundwater path attenuation coefficient;

[0174] Based on the groundwater path attenuation coefficient, long-path attenuation calculations are performed on the initial pollutant load corresponding to the groundwater runoff path. The attenuation calculation results are then superimposed with the background groundwater concentration to output the output flux of the groundwater runoff path.

[0175] In this stage, for the underground runoff path at the bottom layer, because the aquifer path is long and the residence time usually spans multiple hydrological cycles, the initial pollutant load infiltrated by the current rainfall event has already undergone a high degree of exponential decay by the time it reaches the outlet section.

[0176] Therefore, the background concentration of groundwater, established in advance through observations during the dry season, is extracted. This background concentration is then multiplied by the groundwater runoff to obtain the basement background flux.

[0177] Next, the residual load of the current event after long-path attenuation will be added to the background flux of the base to smoothly output the output flux of the underground runoff path.

[0178] As an optional implementation, transport calculations are performed on the initial pollutant loads corresponding to each hydrological transport path. For groundwater runoff paths, the calculations may also include:

[0179] Obtain the pre-configured background groundwater concentration and groundwater path attenuation coefficient;

[0180] Based on the groundwater path attenuation coefficient, long-path attenuation calculation is performed on the initial pollutant load corresponding to the groundwater runoff path to obtain the attenuated event pollutant flux.

[0181] The background flux is obtained by multiplying the background concentration of groundwater by the runoff volume of the groundwater path.

[0182] The attenuated event pollutant flux is superimposed with the background flux to output the output flux of the subsurface runoff path.

[0183] In conjunction with the foregoing embodiments, the total mass of migratable pollutants is further decomposed into the mass of adsorbed pollutants and the mass of dissolved pollutants, including:

[0184] Based on pre-configured physicochemical properties of pollutants, the ionic speciation characteristics of target pollutants are determined.

[0185] If the ionic form is determined to be cationic, then the dynamic partition coefficient and sediment concentration are combined to perform dimensional unification processing, and the adsorbed mass ratio and dissolved mass ratio are calculated.

[0186] If the ionic state is determined to be anionic, the dynamic allocation calculation is skipped, the adsorbed state mass ratio is set to 0, and the dissolved state mass ratio is set to 100%, so that the total mass of the migratable pollutants in the anionic state is allocated to the dissolved pollutant mass.

[0187] Accordingly, different treatment architectures are provided for pollutants with different physicochemical properties.

[0188] In nature, soil colloids typically carry a negative charge on their surface, and their electrochemical properties determine the significant differences in the adsorption capacity of soil particles for different ions. Therefore, after obtaining data on the target pollutant, its ionic property parameters are extracted.

[0189] For example, in the simulation of nitrogen pollutants, if the input target pollutant is ammonium nitrogen, it is determined that it carries a positive charge, which is consistent with the characteristics of a cationic state.

[0190] Because positive and negative charges attract each other, ammonium nitrogen is easily adsorbed by soil colloids. Therefore, it is introduced into the dynamic distribution coefficient calculation and dimensionless processing process to dynamically separate its adsorbed and dissolved mass.

[0191] Conversely, if the target pollutant is nitrate nitrogen, it is determined to be negatively charged, which is consistent with the characteristics of anionic state;

[0192] Due to the repulsion of like charges, nitrate nitrogen is difficult to adhere to the surface of soil colloids, and its migration in the soil and water environment almost entirely depends on the aqueous solution. Therefore, the blocking process is directly triggered, bypassing negative power functions and exponential calculations, forcing its adsorbed mass ratio to be equal to 0 and its dissolved mass ratio to be equal to 100%.

[0193] It should be understood that by using forced truncation to determine the transport patterns of anionic pollutants, redundant iterations can be accurately reflected, and convergence efficiency can be improved.

[0194] In a further embodiment, prior to hydrological simulation and pollutant transport calculation, the method further includes constructing offline physical parameters for obtaining pre-configured reference sediment concentration, reference flow velocity, and baseline distribution coefficients, specifically including:

[0195] Soil and sediment samples were obtained from the target watershed, and batch adsorption experiments were performed on the soil and sediment samples under set hydrodynamic conditions of still water or slow flow and at set standard sediment concentrations.

[0196] The static partition coefficient at which the batch adsorption experiment reaches the physicochemical adsorption equilibrium state is obtained and used as the baseline partition coefficient.

[0197] The established standard sediment concentration and baseline hydrodynamic conditions were respectively calibrated as reference sediment concentration and reference flow velocity.

[0198] Used to establish basic data support through offline physical experiments before the model starts online computation.

[0199] Since the subsequent dynamic allocation coefficients are based on a benchmark value for correction, the accuracy of this benchmark value directly affects the physical fidelity of the entire process. Accordingly, operators collect undisturbed soil samples in the actual environment of the target computing unit and prepare them into standard sediment particles. Artificial water tanks or shaking reactors are constructed in the laboratory, and predetermined flow rates and solid-liquid ratios are set to execute isothermal batch adsorption experiments.

[0200] Furthermore, the ratio of solid phase concentration to liquid phase concentration when the experimental system reaches adsorption-desorption dynamic equilibrium is extracted, and this ratio is established as the static partition coefficient, which is then mapped to the baseline partition coefficient in online calculation.

[0201] Based on this, the sediment concentration values ​​maintained by the environment during the experiment were recorded and mapped to reference sediment concentrations, and the experimental flow rates were recorded and mapped to reference flow rates.

[0202] In some embodiments, determining the sediment concentration corresponding to each hydrological transport pathway includes:

[0203] Soil erosion is calculated based on parameters such as path runoff, topographic slope, and slope length.

[0204] The calculated soil erosion amount is divided by the path runoff at the corresponding time step to determine the sediment concentration used to drive dynamic phase distribution.

[0205] In this step, the mechanism for obtaining sediment concentration is clarified to originate from physical erosion calculations rather than empirical assumptions, namely:

[0206] The system extracts the topographic slope and slope length parameters output by the digital elevation model, combines them with the path runoff output by the hydrodynamic model, and substitutes them into the soil erosion calculation logic to calculate the total mass of soil stripping and transport within a predetermined time step.

[0207] In one example, the formula for calculating this erosion amount is as follows:

[0208] E(t) = k _e ×Q _s (t) γ1 ×(S _0 ) γ2 ×(L _s )γ3 ×C _v ×P _v ;

[0209] Where E(t) is the soil erosion at the current time step, and k _e Q is the comprehensive empirical erosion coefficient used to balance the dimensions of both sides of the equation. _s (t) represents the path flow rate at the current time step, S _0 For the terrain slope, L _s Let γ1, γ2, and γ3 be the corresponding dimensionless hydraulic geometric indices, and C be the slope length parameter. _v P is a vegetation cover management factor. _v Factors related to soil and water conservation measures.

[0210] Based on this, after calculating the soil erosion amount, the step time for that time step is extracted, the total water flow volume flowing through the spatial grid at the current time step is calculated, and then the sediment concentration is calculated using the logic of mass divided by volume. The conversion formula is as follows:

[0211] C _s (t)=E(t) / (Q _s (t)×Δt);

[0212] Among them, C _s Q(t) represents the sediment concentration at the current time step, E(t) represents the soil erosion at the current time step, and Q(t) represents the soil erosion at the current time step. _s (t) represents the path flow rate at the current time step, and Δt represents the time step size.

[0213] Furthermore, the calculated sediment concentration is then sent to the dynamic phase allocation module as an input variable for determining the concentration change rate and correcting the particle concentration effect.

[0214] In another embodiment, determining the sediment concentration corresponding to each hydrological transport path may further include:

[0215] Soil erosion is calculated based on the path runoff, topographic slope and slope length parameters in the underlying surface spatial attribute parameters output by hydrological simulation.

[0216] The calculated soil erosion amount is divided by the product of the path runoff and the step length at the corresponding time step to determine the sediment concentration used to drive dynamic phase distribution.

[0217] In another embodiment, for sub-basins with flat terrain and low runoff, the exponential influence of terrain slope and slope length parameters in the erosion formula is reduced. In this case, a critical shear force model can be used. By determining whether the actual shear force of the water flow is greater than the soil's critical shear force, Boolean variables are used to control whether quantitative calculation of erosion is triggered, further optimizing the allocation of computing power during drought or low rainfall events.

[0218] In one possible implementation, acquiring rainfall-driven data and underlying surface spatial attribute parameters of the computing unit specifically includes:

[0219] Real-time collection of rainfall time series data is achieved through automatic weather monitoring stations and hydrological sensing nodes deployed in the target watershed;

[0220] The pre-stored digital elevation model and remote sensing land use type data are called to parse and obtain the corresponding underlying surface spatial attribute parameters.

[0221] For example, the automatic weather monitoring station is equipped with a tipping bucket rain gauge and a Doppler weather radar to collect rainfall intensity data at a preset time frequency and generate rainfall time series data;

[0222] The preset time frequency can be set comprehensively according to the rainfall monitoring accuracy requirements.

[0223] For example, a hydrological sensing node can be deployed at the outlet section of the target computing unit to acquire measured water level and flow data as boundary conditions for the model flow generation module.

[0224] For example, the extraction of spatial attributes can be achieved by connecting to a geographic information database and calling a pre-stored digital elevation model.

[0225] Accordingly, an eight-direction flow algorithm based on grid flow direction is used to extract terrain elevation differences and calculate the terrain slope and maximum confluence path length of each computation grid, i.e., slope and slope length parameters.

[0226] Furthermore, multispectral remote sensing images are acquired, and supervised classification algorithms are used to analyze land cover features and identify the corresponding remote sensing land use type data.

[0227] One example is mapping the physical grid of the land surface into woodland, grassland, arable land, and construction land.

[0228] Next, different land use types are mapped to a preset physical attribute lookup table to match the Manning roughness coefficient and vegetation cover management factor of the corresponding grid.

[0229] In another embodiment, in remote watersheds where the ground meteorological monitoring network is sparsely distributed, satellite precipitation inversion products provided by global precipitation measurement missions can also be obtained. Using spatial downscaling algorithms, these products can be fused with data from sparse ground monitoring points to generate continuous rainfall time series data with high spatial resolution.

[0230] In some embodiments, after obtaining the total pollutant output flux, the method further includes:

[0231] The total pollutant output flux is input into the preset water environment carrying capacity assessment model to determine the environmental capacity status of the target calculation unit;

[0232] When the environmental capacity status characterization is overloaded, the quality contribution information of the output flux corresponding to each hydrological transport path is extracted as multipath load source tracing information, and a watershed water quality early warning strategy containing multipath load source tracing information is generated based on the multipath load source tracing information.

[0233] When the environmental capacity status indicator is not overloaded, the output result is a safe environmental capacity status.

[0234] The preset water environment carrying capacity assessment model is configured with the maximum allowable daily load of the target water body under predetermined hydrological conditions.

[0235] Accordingly, after obtaining the total pollutant load calculated through time integration, it is numerically compared with the maximum permissible daily load. When the calculated total pollutant load is greater than the maximum permissible daily load, the environmental capacity is determined to be in an overload state; otherwise, it is determined to be in a safe state.

[0236] Furthermore, when an overload state is triggered, intermediate calculation variables from previous steps are extracted to construct a watershed water quality early warning strategy that includes multi-path load source tracing information. Correspondingly, the data stream of the output node is traced back to extract the percentage of mass contribution of the slope runoff path, interflow path, and groundwater runoff path to the total output flux, and the proportion of physical mass in the adsorbed state and the proportion of physical mass in the dissolved state are calculated separately.

[0237] For example, if the source tracing information indicates that the pollutant load exceeds a predetermined proportion and originates from the slope runoff path and is mainly in the adsorbed state, the watershed water quality early warning strategy will generate an intervention instruction to increase the riverside vegetation buffer zone and soil and water conservation measures in the source area.

[0238] In another example, if the source tracing information indicates that the mass contribution rate of dissolved pollutants in the soil flow path exceeds the predetermined mass contribution rate, then the early warning strategy generates an instruction to deploy subsurface flow artificial wetlands to prolong the underground transport time and enhance degradation loss.

[0239] In another embodiment, for downstream receiving water bodies with three-dimensional hydrodynamic and water quality coupling calculation conditions, the total pollutant output flux can be obtained as the upper boundary inflow condition to drive a two-dimensional or three-dimensional river network hydrodynamic model, outputting the spatiotemporal distribution field of pollutant concentration at the downstream control section, thereby supporting dynamic water quality early warning based on gridded spatial prediction.

[0240] Furthermore, the total load of a single rainfall event is quantitatively extracted. The total pollutant output flux, whose physical dimension is mass over time, needs to be integrated over time. Accordingly, a discrete-time step summation algorithm is used, and its calculation formula is as follows:

[0241] L _event =Σ[F _total [(t)×Δt];

[0242] Among them, L _event The total pollutant output load for a single rainfall event is represented by Σ, which indicates the summation operation performed over all computational time steps throughout the entire duration of the target rainfall event. _total (t) represents the total pollutant output flux calculated at the current time step, and Δt is the discrete time step. This integral operation transforms the variable from a rate variable to an absolute mass variable.

[0243] Furthermore, it is necessary to convert the dimensionless hysteresis factor into a specific physical time delay parameter, that is:

[0244] τ _i =(L _i / v _i_avg )×R _f ;

[0245] Where, τ _i L represents the lag time of the interflow path. _i v represents the average confluence path length of the interflow within the target space grid. _i_avg R is the average water flow velocity along this path during the corresponding time period. _f This is the hysteresis factor obtained from the aforementioned calculation. Through this multiplicative transformation, the chemical adsorption resistance of soil pores is mapped to a computational delay in the time dimension.

[0246] In an exemplary verification scenario, a typical subtropical hilly watershed with moderate rainfall intensity is selected as the target calculation unit. The same rainfall event is simulated and compared using the method of this invention and the traditional single-channel method with fixed allocation coefficients.

[0247] The results show that the total pollutant output flux time series output by the method of the present invention is closer to the measured value in terms of peak time and peak amplitude, and the flux decay curve shape in the receding stage is better than that of the measured process line than the traditional method.

[0248] It is understandable that the specific improvement in accuracy may vary depending on the characteristics of the watershed, rainfall conditions, and the type of pollutants.

[0249] According to one aspect of this application, a pollutant generation and transport calculation system based on multi-process coupling is provided, which can be used to execute the method provided in the embodiments of this application. The calculation system includes:

[0250] The system includes a data acquisition module, a hydrological simulation and feature extraction module, a dynamic phase allocation module, a differential transport module, and a flux integration and output module.

[0251] Correspondingly, the data acquisition module is used to acquire rainfall-driven data, underlying surface spatial attribute parameters, and total mass of mobile pollutants from the computing unit. Specifically, this module is configured to connect to an external meteorological station interface and a geographic information database, extract time-series data at a preset sampling frequency, and convert the extracted data into a standardized matrix data structure for use by downstream computing nodes.

[0252] The hydrological simulation and feature extraction module communicates with the data acquisition module and is used to perform hydrological simulation based on rainfall-driven data and underlying surface spatial attribute parameters. It separates multiple hydrological transport paths and determines the hydrodynamic parameters and sediment concentrations corresponding to each hydrological transport path. The module integrates a three-source runoff physical algorithm and a soil erosion algorithm and is responsible for outputting the water flow volume and dynamic characteristic variables of the surface and subsurface paths in parallel.

[0253] The dynamic phase allocation module communicates with the hydrological simulation and feature extraction module. It is used to calculate the dynamic allocation coefficient based on hydrodynamic parameters and sediment concentration. The dynamic allocation coefficient is used to allocate the total mass of mobile pollutants to each hydrological transport path to obtain the initial pollutant load corresponding to each hydrological transport path. The module has a built-in mathematical operation unit, which is responsible for performing floating-point multiplication and addition operations including the negative power function of particle concentration effect and hydrodynamic saturation inhibition function, so as to output the allocation weights that are dynamically updated over time.

[0254] The differentiated transport module, which communicates with the dynamic phase allocation module, is used to perform transport calculations on the initial pollutant loads corresponding to each hydrological transport path. During the transport process, it performs secondary redistribution of pollutant phases across paths based on the spatial changes in sediment concentration to obtain the output flux corresponding to each hydrological transport path. This module adopts a spatial grid unidirectional recursive computational architecture and is equipped with a state monitoring interface to capture the difference in adsorption capacity in real time and trigger calculation instructions for desorption mass transfer.

[0255] The flux integration output module communicates with the differentiated transport module and is used to integrate the output fluxes corresponding to each hydrological transport path to obtain the total pollutant output flux.

[0256] According to another aspect of this application, an electronic device is provided.

[0257] Accordingly, the electronic device includes a processor and a memory communicatively connected to the processor. The memory stores computer instructions that can be executed by the processor;

[0258] When the processor executes computer instructions, it implements the method provided in any of the embodiments of this application.

[0259] Based on this, the processor can be implemented using physical hardware such as a central processing unit, microprocessor, digital signal processor or application-specific integrated circuit;

[0260] Memory, such as random access memory, read-only memory, or flash memory, is a non-volatile storage medium. When the electronic device is in operation, the processor calls computer instructions from the memory through the system data bus to sequentially complete a series of calculation tasks, including data extraction, hydrodynamic coupling separation, phase dynamic update, and secondary desorption calculation, and outputs quantitative flux data that can be used for water environment carrying capacity assessment.

[0261] Accordingly, this embodiment also provides a computer-readable storage medium;

[0262] A computer program is stored on a computer-readable storage medium. When the computer program is executed by a processor of an electronic device, it can implement the steps and procedures provided in any of the method embodiments.

[0263] Computer-readable storage media can be implemented using physical carriers such as optical discs, magnetic hard disks, or solid-state drives.

[0264] In summary, this invention discloses a method for calculating pollutant generation and transport based on multi-process coupling, including:

[0265] Acquire rainfall-driven data, underlying surface spatial property parameters, and total mass of mobile pollutants for the target computational unit;

[0266] Runoff generation is calculated based on rainfall-driven data and underlying surface spatial attribute parameters, and the runoff volume and hydraulic characteristic parameters of slope runoff path, interflow path and groundwater runoff path are obtained.

[0267] Erosion simulation was performed based on the runoff volume and hydraulic characteristic parameters of the slope runoff path to obtain the sediment concentration of the slope runoff.

[0268] The dynamic allocation coefficient is calculated by using the slope runoff sediment concentration and hydraulic characteristic parameters as driving variables. Based on the dynamic allocation coefficient, the total mass of mobile pollutants is decomposed into the mass of adsorbed pollutants and the mass of dissolved pollutants. Based on the dynamic conditions of each path, the mass of adsorbed pollutants and the mass of dissolved pollutants are weighted and allocated to each path to obtain the initial pollutant load of each path.

[0269] Transport calculations are performed separately for the initial pollutant loads along each path;

[0270] In the path of slope runoff transport, the amount of adsorbed desorption and release is determined based on the change in environmental adsorption capacity caused by the spatial decay of sediment concentration. The amount of adsorbed desorption and release is then transferred to the interflow path to obtain the output flux of each path to the confluence node.

[0271] The output fluxes of each path are superimposed to obtain the total pollutant output flux at the outlet section of the target calculation unit.

[0272] In other embodiments, parts of the method of the present invention may also include:

[0273] Analyze the changes in key non-point source pollution parameters in the catchment area, small watershed, and tributaries flowing into the reservoir, and optimize the pollution parameters;

[0274] Improve model computation speed by utilizing second-order algorithms and GPU parallel technology;

[0275] Construct a pollutant generation and transport model based on the coupling of multiple processes including watershed rainfall, runoff and sediment production, and pollutant flux.

[0276] Furthermore, a watershed-reservoir transition unit for finite element structural simulation is constructed, a multidimensional hydrodynamic and water quality model of the watershed and reservoir is coupled, the relationship of water quality parameters at the connection section between the two is supplemented, and it is integrated into a digital twin model of water quality of a certain reservoir.

[0277] By integrating multidimensional hydrodynamics and hydrological and water quality simulation, the accuracy of the integrated model prediction is optimized and improved. By combining data assimilation and machine learning methods, the computational efficiency of the model is improved, enabling the simulation of pollutant flux and real-time prediction of water quality.

[0278] According to one aspect of this application, for the process of sediment attenuation along the slope, a first-order spatial attenuation model is adopted, which combines the sediment particle settling characteristics and the slope hydrodynamic conditions to recursively solve the local sediment concentration along the runoff direction, node by node.

[0279] Based on the sediment concentration at the upstream nodes, the sedimentation loss effect of sediment during runoff transport is characterized by taking into account sediment particle settling velocity, slope runoff velocity, water depth and node spatial step size, and the sediment concentration distribution at each downstream node is obtained in turn.

[0280] At the first spatial node, the overall sediment concentration is used as the initial value.

[0281] The overall settling velocity of sediment particles can be determined by using the well-known Stokes settling formula or empirical settling velocity table, based on the sediment particle size distribution characteristics of the study area.

[0282] According to another aspect of this application, the output flux of the underground runoff path is determined by the superposition relationship between pollutant load, attenuation loss and background concentration within the underground runoff path;

[0283] By combining the initial pollutant load of groundwater runoff, the calculation time step, the groundwater attenuation coefficient, and the average residence time, the natural attenuation effect of pollutants during underground transport is characterized. The output flux of the groundwater path is obtained by superimposing the contributions of the background concentration of groundwater and the real-time groundwater runoff.

[0284] The average residence time of groundwater can be determined based on regional hydrogeological exploration data or tracer test results.

[0285] For example, the criteria for determining the decomposition of three water sources can be expressed as follows:

[0286] Slope runoff occurs when rainfall intensity exceeds the soil saturation permeability coefficient.

[0287] Rainfall infiltration exceeds field capacity, forming interflow; deep infiltration replenishes aquifers, forming groundwater runoff.

[0288] In other embodiments, the total mass of migratable contaminants can be obtained in the following ways:

[0289] The determination is based on the product relationship between the pollutant concentration in the soil surface layer within the target calculation unit, the thickness of the active layer affected by rainfall erosion, and the area of ​​the calculation unit.

[0290] In other embodiments, the particle concentration effect index and hydrodynamic correction coefficient can be obtained by performing batch adsorption experiments on soil samples from the target watershed under different solid-liquid ratios and different water flow velocities, and then using the least squares fitting method to perform parameter inversion based on experimental data and the aforementioned formula.

[0291] In one alternative implementation, the spatial trigger threshold can be adaptively adjusted according to the sediment particle size distribution characteristics of the target watershed and the required simulation accuracy.

[0292] The parameters and thresholds for which no setting conditions or numerical examples are mentioned in this application can be set and adjusted by engineers based on experience, actual working conditions, etc.

[0293] This invention introduces particle concentration effects and hydrodynamic correction terms to construct a phase allocation coefficient calculation mechanism that adapts to the spatiotemporal dynamic evolution of the water and sediment environment. This mechanism helps to suppress the interference of competitive adsorption of suspended particles and water scouring on chemical equilibrium in transient scenarios, and reduces static assessment errors caused by periods of drastic concentration fluctuations.

[0294] Furthermore, this application not only performs multi-path dynamic weighted allocation based on physical carrying capacity in the source region, but also constructs a secondary phase redistribution network along the route. By quantifying the net decrease in adsorption capacity caused by sediment deposition, the desorption and release during transport and their cross-interface transport to underground space are reconstructed, overcoming the limitation of the traditional model that the phase state is transported from one point to the end.

[0295] Furthermore, this application, in conjunction with a pore transport model incorporating physical hysteresis and degradation attenuation, and a differential blocking determination of anions and cations based on electrochemical characteristics, can be used to decouple and restore the heterogeneous evolution characteristics of multiple physical, chemical, and hydrological processes within a natural watershed, thereby elevating the reliability and spatiotemporal adaptability of watershed multi-process coupled simulation to a new level.

[0296] It should be noted that the various specific technical features described in the above embodiments can be combined in any suitable manner without contradiction. To avoid unnecessary repetition, the present invention will not describe the various possible combinations separately.

[0297] The parameters and thresholds for which no setting conditions or numerical examples are mentioned in this application can be set and adjusted by engineers based on experience, actual working conditions, etc.

Claims

1. A method for calculating pollutant generation and transport based on multi-process coupling, characterized in that, include: Acquire rainfall-driven data, underlying surface spatial property parameters, and total mass of mobile pollutants for the computing unit; Hydrological simulations were performed based on rainfall-driven data and spatial property parameters of the underlying surface to separate multiple hydrological transport paths and determine the hydrodynamic parameters and sediment concentrations corresponding to each hydrological transport path. Based on hydrodynamic parameters and sediment concentration, a dynamic allocation coefficient is calculated. The total mass of mobile pollutants is then allocated to each hydrological transport path using the dynamic allocation coefficient, thus obtaining the initial pollutant load corresponding to each hydrological transport path. Transport calculations are performed on the initial pollutant loads corresponding to each hydrological transport path. During the transport process, a secondary redistribution of pollutant phases across the path is performed based on the spatial changes in sediment concentration to obtain the output flux corresponding to each hydrological transport path. The output fluxes corresponding to each hydrological transport pathway are integrated to obtain the total pollutant output flux; Among them, multiple hydrological transport pathways include: slope runoff pathways, interflow pathways, and groundwater runoff pathways; hydrodynamic parameters include slope runoff velocity of the slope runoff pathway. The calculation of the dynamic allocation coefficient based on hydrodynamic parameters and sediment concentration includes: calculating the rate of change of sediment concentration at the current time step relative to the previous time step; when the rate of change of concentration exceeds a preset time trigger threshold, calculating the dynamic allocation coefficient based on a pre-configured reference sediment concentration, reference flow velocity, and baseline allocation coefficient; wherein, in calculating the dynamic allocation coefficient, a particle concentration effect correction is introduced based on the negative power function relationship of the ratio of sediment concentration to reference sediment concentration; and a hydrodynamic correction is introduced based on the saturation inhibition function relationship of the ratio of slope runoff velocity to reference velocity. The calculation of the dynamic allocation coefficient also includes: obtaining the real-time water temperature; constructing a temperature correction term using the Arrhenius equation based on the water temperature and a pre-configured reference temperature; and using the temperature correction term to perform thermodynamic attenuation or enhancement correction on the allocation coefficient result that superimposed the particle concentration effect correction and the hydrodynamic correction to obtain the dynamic allocation coefficient. The determination of sediment concentration corresponding to each hydrological transport path includes: calculating soil erosion based on path runoff, topographic slope, and slope length parameters; dividing the calculated soil erosion by the path runoff at the corresponding time step to determine the sediment concentration used to drive dynamic phase distribution. The secondary phase redistribution includes: recalculating the local dynamic allocation coefficient based on the local sediment concentration at the current spatial node, updating the local adsorbed mass ratio using the local dynamic allocation coefficient; when the updated local adsorbed mass ratio is lower than the adsorbed mass ratio at the previous spatial node, it indicates a decrease in environmental adsorption capacity; based on the difference in the decrease in the adsorbed mass ratio and the total remaining pollutant mass arriving at the current spatial node, calculating the adsorbed desorption release amount; deducting the dissolved pollutants corresponding to the adsorbed desorption release amount from the slope runoff path, adding them across the path to the pollutant load of the interflow path, and outputting the output flux of the slope runoff path. The process involves performing transport calculations on the initial pollutant loads corresponding to each hydrological transport path. For the interflow path, this includes: combining the initial pollutant loads corresponding to the interflow path with the additional adsorbed desorption and release amounts to obtain the total pollutant input load of the interflow path; calling a transfer function model to perform transport calculations on the total pollutant input load, wherein the transfer function model includes a hysteresis factor to characterize the adsorption and desorption of the soil matrix and a first-order decay coefficient to characterize the degradation process; and outputting the output flux of the interflow path after hysteresis delay and concentration decay based on the transfer function model.

2. The method according to claim 1, characterized in that, Obtain rainfall-driven data and underlying surface spatial attribute parameters for the computational unit, specifically including: Real-time collection of rainfall time series data is achieved through automatic weather monitoring stations and hydrological sensing nodes deployed in the target watershed; By calling the digital elevation model and remote sensing land use type data, the corresponding underlying surface spatial attribute parameters are obtained through parsing.

3. The method according to claim 1, characterized in that, After obtaining the total pollutant output flux, the following is also included: The total pollutant output flux is input into the preset water environment carrying capacity assessment model to determine the environmental capacity status of the target calculation unit; When the environmental capacity status characterization is overloaded, a watershed water quality early warning strategy containing multi-path load source tracing information is generated.

4. The method according to claim 1, characterized in that, Transport calculations were performed on the initial pollutant loads corresponding to each hydrological transport pathway, including those for groundwater runoff pathways: Obtain the background concentration of groundwater and the groundwater path attenuation coefficient; Based on the groundwater path attenuation coefficient, long-path attenuation calculations are performed on the initial pollutant load corresponding to the groundwater runoff path. The attenuation calculation results are then superimposed with the background concentration of groundwater to output the output flux of the groundwater runoff path.

Images