Quantitative assessment method for risk of pollutant release from sediment based on sediment incipient motion mode

Abstract

The application discloses a sediment pollution release risk quantitative evaluation method based on a sediment start-up mode, relates to the technical field of water environment monitoring, and obtains water power parameters, bed surface structure parameters and pollutant activity parameters at a sediment interface; extracts a turbulent burst event sequence, determines a surface microstructure yield threshold; establishes a three-dimensional parameter coupling mapping relationship, modulates the yield threshold and calculates a micro-disturbance erosion thickness, and corrects the effective thickness of the interface diffusion boundary layer; calculates a continuous sub-start disturbance intensity index, substitutes the index into an interface effective diffusion coefficient dynamic equation, and obtains a pollutant instantaneous release flux; reversely substitutes the release flux into a pore water concentration attenuation equation, reversely calculates a micro-disturbance duration and compares the micro-disturbance duration with a measured value, corrects the yield threshold and the effective thickness when a deviation exceeds a threshold value, and iterates until the deviation falls within a threshold value range, and outputs a release risk grade, so that the method can quantitatively evaluate the release risk of sediment pollutants under real water power conditions.

Description

Technical Field

[0001] This invention relates to the field of water environment monitoring technology, specifically to a method for quantitatively assessing the risk of pollutant release from sediment based on sediment initiation modes. Background Technology

[0002] In existing technologies, the risk assessment of sediment pollutant release mainly relies on two approaches: one is based on chemical equilibrium distribution or speciation experiments, calculating static release potential through solid content and pore water concentration; the other borrows the critical initiation criterion from river sediment dynamics, activating the release module only when the shear stress of the flow exceeds a threshold, assuming that interface exchange is negligible if initiation is not reached. However, under real hydrodynamic conditions, before macroscopic collective initiation occurs, there are microscopic interface disturbances caused by bottom turbulence in the sediment. These disturbances are sufficient to disrupt the surface microstructure and accelerate the actual release of pollutants into the water body. Traditional binary initiation models precisely ignore this intermediate state, leading to a systematic underestimation of release fluxes under moderate hydrological events, and the assessment parameters become severely ineffective when transferred between different sedimentary environments. The fundamental reason is that existing methods lack the ability to physically couple the near-bed surface motion states of sediment particles of different orders of magnitude with the pollutant desorption dynamics, and cannot quantify the impact of micro-disturbances on the diffusion coefficient and interface update frequency. The conventional improvement approach in this field is to improve the fitting accuracy by introducing empirical correction coefficients or adding field correction factors based on the existing critical initiation formula. This type of forward optimization method always retains the starting threshold as a logical node, thus failing to fundamentally avoid missing the contribution of intermediate state releases. In other words, regardless of how the correction coefficient is adjusted, as long as the logical switch structure is retained, there will inevitably be a critical boundary point for determining whether to start or not, and the existence of this boundary point itself constitutes a forced truncation of the continuous physical process. Current practices have not yet recognized the need to eliminate this logical switch at the evaluation architecture level.

[0003] To address the aforementioned shortcomings, this invention proposes a quantitative evaluation method based on the decomposition of sediment initiation modes and the coupling of release responses. This method uses the microscopic motion state of sediment particles as the true triggering factor for release, thereby solving the problem of risk misjudgment caused by neglecting sub-initiation interface disturbances in traditional methods. Summary of the Invention

[0004] The purpose of this invention is to provide a method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation modes, in order to solve the problems mentioned in the background art.

[0005] To address the aforementioned technical problems, this invention provides the following technical solution: a method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation modes, comprising the following steps:

[0006] Step 1: Obtain hydrodynamic parameters, bed structure parameters, and pollutant activity parameters at the sediment interface in the assessment area;

[0007] Step 2: Extract the time series of turbulent burst events based on the hydrodynamic parameters, determine the yield threshold of the surface microstructure based on the bed structure parameters, and determine the initial concentration gradient at the interface based on the pollutant activity parameters;

[0008] Step 3: Establish the coupling mapping relationship between the hydrodynamic parameters, the bed structure parameters, and the pollutant activity parameters. Using the turbulent burst event as the excitation boundary, dynamically modulate the yield threshold, and calculate the micro-disturbance erosion thickness using the modulated yield threshold. Then, use the micro-disturbance erosion thickness to correct the effective thickness of the diffusion boundary layer. The physical basis of the above coupling mapping is that the turbulent pulsation energy does not directly act on the pollutant molecules, but sequentially changes the stress state of the bed microstructure, thereby affecting the range of particles that can be eroded, and finally, the degree of unobstructed diffusion path at the interface is determined by the amount of erosion. The three links are in a causal transmission relationship, and the parameter change in any link will be transmitted downstream along this causal chain. If the parameters of each dimension are calculated independently and then simply superimposed, the causal transmission path will be cut off, and the amplification effect of micro-disturbance on interface diffusion will not be reflected. Therefore, the three need to be included in the same coupling framework for synchronous calculation.

[0009] Step 4: Calculate the sub-initiation disturbance intensity index based on the time series of the turbulent burst event and the modulated yield threshold. This sub-initiation disturbance intensity index is a physical quantity with a value within a continuous range, used to quantify the real-time intensity level of micro-detachment and local transition of sediment surface particles under current hydrodynamic conditions. The index output in this step does not have any fixed critical boundary value, and its value responds continuously with the change in turbulent fluctuation intensity. This makes the evaluation result no longer dependent on a pre-set logic switch, but directly reflects the actual degree of interface micro-disturbance, providing a continuous driving variable that matches the physical process for subsequent flux release calculations.

[0010] Step 5: Substitute the sub-start disturbance intensity index into the dynamic equation of the effective diffusion coefficient of the interface to calculate the instantaneous release flux of pollutants;

[0011] Step Six: Substitute the instantaneous release flux of the pollutant back into the bed surface pore water concentration decay equation to deduce the required micro-disturbance duration. Compare the deduce micro-disturbance duration with the measured turbulent burst duration. If the deviation exceeds a preset threshold, correct the yield threshold and the effective thickness, and repeat the calculation of the sub-start disturbance intensity index and subsequent steps until the deviation falls within the preset threshold range. Output the final release risk level.

[0012] Preferably, the coupling mapping relationship is established as follows: the pulsation intensity of the turbulent burst event is used as the first driving variable to dynamically reduce the yield threshold of the surface microstructure, and the reduction magnitude is positively correlated with the pulsation intensity; at the same time, the amount of fine particles released by micro-disturbance erosion is used as the weighting coefficient to compress the effective thickness of the diffusion boundary layer, and the compression magnitude is positively correlated with the amount of fine particles released.

[0013] Preferably, the micro-disturbance erosion thickness is calculated as follows: the upper limit of the particle size that can be eroded is determined according to the modulated yield threshold, and then the volume ratio of each particle size component with a particle size smaller than the upper limit is accumulated and the corresponding layer thickness is multiplied by the particle size distribution curve of the surface micro-aggregates to obtain the micro-disturbance erosion thickness.

[0014] Preferably, the sub-start disturbance intensity index is calculated as follows: the cumulative probability density of turbulent burst energy exceeding the yield threshold is statistically analyzed per unit time, and a positive feedback coefficient of the bed pore water pressure on microcrack expansion is introduced. The cumulative probability density is multiplied by the positive feedback coefficient and then normalized to obtain the sub-start disturbance intensity index.

[0015] Preferably, the preset threshold is a combination of an energy conservation constraint boundary and a mass conservation verification threshold. The energy conservation constraint boundary is used to limit the range of the difference between the micro-disturbance input energy and the energy dissipated by the bed deformation, and the mass conservation verification threshold is used to limit the upper limit of the allowable deviation between the disturbance duration calculated by the release flux and the measured disturbance duration.

[0016] Preferably, when repeatedly performing the calculation of the sub-start disturbance intensity index and subsequent steps, the yield threshold is corrected by increasing or decreasing the value of the yield threshold according to the deviation direction between the previously back-calculated disturbance duration and the measured disturbance duration; the effective thickness is corrected by decreasing or increasing the value of the effective thickness according to the magnitude by which the previously output instantaneous release flux exceeds the preset range.

[0017] Preferably, the hydrodynamic parameters include at least the instantaneous shear stress spectrum of bottom turbulence and the near-bed velocity gradient; the bed structure parameters include at least the particle size distribution of surface micro-aggregates and the extracellular polymer cementation index; and the pollutant activity parameters include at least the ratio of dissolved pollutant concentration to particulate pollutant concentration in the surface sediments.

[0018] Preferably, in the dynamic equation of the effective diffusion coefficient of the interface, the effective diffusion coefficient increases continuously with the increase of the sub-start disturbance intensity index, and the effective diffusion coefficient is obtained by interpolation and table lookup using the sub-start disturbance intensity index as the independent variable through a pre-calibrated sediment type correlation curve.

[0019] Preferably, the hydrodynamic parameters further include the bottom turbulence burst frequency and the pulsating pressure amplitude, wherein the bottom turbulence burst frequency is used to extract the time series of the turbulence burst event, and the pulsating pressure amplitude is used to calculate the turbulence burst energy.

[0020] Preferably, the time series of the turbulent burst event is divided into four levels according to the burst intensity, the sub-start disturbance intensity component under each level is calculated, and the components are weighted and summed according to the time proportion to obtain the comprehensive sub-start disturbance intensity index. The comprehensive sub-start disturbance intensity index is used to replace the sub-start disturbance intensity index and is substituted into the dynamic equation of the effective diffusion coefficient of the interface.

[0021] A sediment pollutant release risk quantification assessment system based on sediment initiation modes includes:

[0022] The data acquisition unit is used to acquire hydrodynamic parameters, bed structure parameters, and pollutant activity parameters at the sediment interface in the assessment area.

[0023] The coupling mapping unit is used to establish the coupling mapping relationship between the hydrodynamic parameters, the bed structure parameters and the pollutant activity parameters, use the turbulent burst event as the excitation boundary to dynamically modulate the yield threshold, calculate the micro-disturbance erosion thickness, and correct the effective thickness of the diffusion boundary layer.

[0024] The intensity calculation unit is used to calculate the sub-start disturbance intensity index based on the time series of the turbulent burst event and the modulated yield threshold.

[0025] The flux calculation unit is used to substitute the sub-start disturbance intensity index into the dynamic equation of the effective diffusion coefficient of the interface to calculate the instantaneous release flux of pollutants.

[0026] The closed-loop verification unit is used to back-substitute the instantaneous release flux of the pollutant into the bed surface pore water concentration decay equation, back-calculate the duration of micro-disturbance, compare the back-calculated time with the measured duration of turbulent burst, trigger iterative correction when the deviation exceeds a preset threshold, until the deviation falls within the threshold range, and output the final release risk level.

[0027] This invention provides a method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation modes. It has the following beneficial effects:

[0028] By employing a continuous sub-initiation disturbance intensity index to replace the traditional binary initiation discrimination logic, this method can capture the microscopic stripping and local transition behavior of the sediment surface before macroscopic initiation. This microscopic behavior truly exists under the sudden action of water flow turbulence. Traditional methods ignore this stage because they rely on fixed critical thresholds. The evaluation path of this invention no longer sets logic switches, and the evaluation results can cover a wider range of hydrodynamic effects. This method helps to reduce the risk of underestimating release flux under moderate hydrological conditions. At the same time, this method establishes a coupled mapping relationship between hydrodynamic parameters, bed structure parameters, and pollutant activity parameters. The parameters of the three dimensions are mutually modulated during the calculation process, rather than being simply superimposed after independent calculation. This coupling mechanism allows the effective diffusion coefficient to be smoothly adjusted with the change of disturbance intensity, and the interface update frequency no longer relies on empirical constants for ex-post correction. This approach helps to improve the model's adaptability to different sedimentary environments.

[0029] Through a closed-loop iterative verification mechanism, the calculated release flux is back-substituted into the pore water concentration decay equation. By back-calculating the required disturbance time and comparing it with the measured disturbance time, a closed computational link is formed. This link can automatically correct the initial values ​​of the yield threshold and the effective thickness of the diffusion boundary layer. This process does not rely on external measured data for manual calibration. This mechanism helps reduce the impact of input parameter uncertainty on the final evaluation result. At the same time, this method further introduces a graded weighting method and an auxiliary parameter verification method in the weighting. The graded weighting method helps smooth the interference of time-series fluctuations on the output flux, and the auxiliary parameter verification method can maintain the integrity of event extraction when the quality of hydraulic monitoring data fluctuates. The above weighting methods are independent and can be used in combination. This method as a whole provides a complete computational link from micro-disturbance identification to risk level output. This link does not rely on prior judgment of the starting state, and the output risk level can reflect the release tendency of sediment pollutants under real hydrodynamic conditions.

[0030] This method supports batch processing of long-term series monitoring data. By statistically analyzing the calculation results of multiple turbulent events, it can extract regionally representative release risk patterns and characteristic thresholds, thus taking into account both the needs of instantaneous response assessment and long-term trend analysis in prototype observation. It can provide an operable quantitative tool for water environment management and ecological risk assessment. Attached Figure Description

[0031] Figure 1 This is a flowchart illustrating the method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode, as described in this invention.

[0032] Figure 2 This is the timing diagram for the closed-loop iterative verification of the present invention. Detailed Implementation

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

[0034] Example 1: Please refer to Figure 1 and Figure 2 This invention provides a method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation modes. This method is applicable to water bodies such as rivers, lakes, reservoirs, and nearshore sea areas with fine-particle sediments and pollutant enrichment. In practical applications, the monitoring device can operate continuously within the assessment area, collecting long-term hydrodynamic and pollutant concentration data. Each turbulent burst event corresponds to an instantaneous release flux calculation result. By statistically summarizing the calculation results of multiple events, an assessment conclusion reflecting the sediment release pattern in the area can be formed.

[0035] After selecting the assessment area, a sediment interface monitoring device is first deployed. This device includes at least an underwater acoustic Doppler current meter installed 30 to 50 cm above the sediment surface, a miniature pressure sensor installed at the bed surface (level with the sediment-water interface), and a pore water sampling probe buried 2 to 5 cm below the surface. All three are arranged along the same vertical line, with the current meter and pressure sensor positioned to ensure they are within the sensing range of the same turbulent vortex structure. Once monitoring is initiated, instantaneous flow velocity data of the bottom water, bed surface pulsating pressure data, and dissolved pollutant concentration data in the interstitial water of the surface sediment are continuously collected. Simultaneously, undisturbed sediment samples collected in the field are tested in the laboratory to obtain the particle size distribution curve of surface microaggregates, the relative content index of extracellular polymers, and the initial concentration ratio of particulate pollutants to dissolved pollutants. These collected and tested results together constitute the input parameter set for subsequent calculations.

[0036] Next, the collected instantaneous velocity data underwent time-frequency transformation processing to extract the turbulent burst event sequence per unit time. Each burst event in this sequence recorded its start time, duration, and peak pulsation energy. Furthermore, the particle size corresponding to a cumulative volume percentage of 65% in the surface micro-aggregate particle size distribution curve was used as the initial yield threshold. This percentage was selected based on the following criteria: when the cumulative volume percentage is below 65%, the corresponding particle size is too small, which would overestimate the erosion resistance of the surface microstructure; when the percentage is above 75%, the corresponding particle size is too large, which would introduce large particles that do not provide structural support into the yield judgment range. Indoor comparative tests on different types of sediments verified that within the range of 65 to 75%, the yield threshold has a relatively stable impact on the subsequent flux calculation results. In this embodiment, 65 was used as the initial value, which will be automatically adjusted through a closed-loop iteration process.

[0037] Then, a multi-dimensional parameter coupling mapping relationship is established: the peak pulse energy of each turbulent burst event is used as the driving variable. According to a pre-calibrated energy-weakening function, the instantaneous reduction of the yield threshold due to this event is calculated. The calibration process of this function is as follows: in a laboratory tank, using the same type of sediment sample, different levels of turbulent pulse energy are applied, and the actual yield stress value when the surface micro-aggregates begin to displace under this energy is recorded. The experimental data for each group are plotted as a scatter plot with pulse energy as the horizontal axis and the yield stress reduction ratio as the vertical axis. A monotonically increasing calibration curve is obtained by fitting a piecewise linear interpolation method, i.e., connecting adjacent experimental data points with a straight line and applying monotonicity constraints to ensure that the vertical axis value corresponding to any horizontal axis on the curve is unique and does not decrease. In actual calculations, the peak pulse energy of the current turbulent burst event is input as the horizontal axis into the calibration curve to find the corresponding yield threshold reduction ratio. Multiplying this reduction ratio by the initial yield threshold yields the effective yield threshold at the current moment.

[0038] The energy-weakening function exhibits a monotonically increasing relationship, meaning that the higher the pulsation energy, the greater the reduction in the yield threshold. The reduced yield threshold is taken as the effective yield threshold at the current moment. Simultaneously, based on this effective yield threshold and combined with the surface micro-aggregate particle size distribution curve, all particle sizes smaller than the effective yield threshold are screened out. The product of the volume percentage of these components and their corresponding layer thickness is accumulated to obtain the micro-perturbation erosion thickness caused by the sudden event.

[0039] The erosion thickness is determined according to the following calculation rules:

[0040] Let the total number of particle size groups obtained by screening be n, and let each particle size group be numbered i, where i = 1, 2, 3, ..., n; let V be the volume percentage of each particle size group. i, which is a dimensionless percentage; the layer thickness corresponding to each particle size group is denoted as d. i .

[0041] Accumulate V of all granular groups i With d i The product of these values ​​is the accumulated value, which is the micro-disturbance erosion thickness D.

[0042] The above calculation rules are executed independently in each burst event. The specific value of the effective yield threshold determines the size of the number of grain groups n participating in the accumulation. The lower the effective yield threshold, the more grain groups participate in the accumulation and the greater the erosion thickness.

[0043] The erosion thickness is further used to correct the effective thickness of the diffusion boundary layer. The correction method is to subtract a reduction amount proportional to the erosion thickness from the initial effective thickness, and the reduced value is taken as the effective thickness of the diffusion boundary layer at the current moment.

[0044] Based on the above coupling mapping, the continuous sub-start disturbance intensity index is calculated; the calculation process of this index is as follows:

[0045] The peak energy of all turbulent burst events within a unit time window is statistically analyzed. Each peak energy is compared with the effective yield threshold at the current moment, and all events whose peak energy exceeds the threshold are recorded.

[0046] The peak energy values ​​of these transcendental events are extracted, and their cumulative proportion in the total energy of all events is calculated. This is then multiplied by the positive feedback coefficient of pore water pressure on microcrack expansion. This positive feedback coefficient characterizes the promoting effect of each unit of microcrack opening under pore water pressure on the further instability of surrounding particles, thus obtaining the initial strength value. The value of this positive feedback coefficient is determined in advance as follows: In laboratory tests, a known water pressure is applied to the sediment sample, and the change in the opening area of ​​the surface microcracks is observed. Then, the water pressure is gradually increased, and the new opening area is recorded. The slope between the water pressure increment and the opening area increment is calculated. This slope is then compared with the pure mechanical erosion under the same conditions without considering water pressure; the resulting ratio is the positive feedback coefficient. For sediments with different degrees of compaction, this coefficient varies between 1.2 and 3.0, with lower values ​​for higher compaction levels.

[0047] The initial intensity value is then divided by the maximum peak energy of all events within the time window to complete the normalization process. The normalized value is the sub-initiation disturbance intensity index corresponding to that time window. This index, with a value within a continuous interval, characterizes the real-time intensity level of microscopic peeling and local transition of sediment surface particles. The complete calculation rule for this index is as follows: set the total number of turbulent burst events within the current time window as N, and number each event in the order of occurrence as j, where j = 1, 2, 3, ..., N; the peak energy of each event is denoted as E. j The effective yield threshold at the current moment is denoted as τ; the set of events exceeding this threshold is denoted as S, and the number of events in the set is M.

[0048] Calculate the sum of peak energies of all events within set S, divide it by the sum of peak energies of all events within the window to obtain the excess energy accumulation ratio R, which is a dimensionless value. Multiply R by the positive feedback coefficient α to obtain the initial intensity value I0, which is also a dimensionless value. Divide I0 by the maximum peak energy E among all events within the window. max The sub-start disturbance intensity index I is obtained for this time window. This index is a dimensionless value and its value range is [0,1].

[0049] The normalized value is the sub-initiation disturbance intensity index corresponding to that time window. This index takes values ​​within a continuous range, characterizing the real-time intensity level of microscopic detachment and local transition of sediment surface particles.

[0050] Subsequently, the calculated sub-initiation disturbance intensity index was substituted into the dynamic equation for the effective diffusion coefficient of the interface. This dynamic equation was constructed as follows: For different types of sediment samples, the effective diffusion coefficients corresponding to different disturbance intensity levels were measured in a laboratory flume, forming sediment type correlation curves. During actual calculation, the corresponding correlation curve was selected based on the current sediment type, and the sub-initiation disturbance intensity index was used as the lookup table input. The effective diffusion coefficient value at the current moment was obtained through interpolation. The interpolation lookup table rule is as follows: the pre-calibrated sediment type correlation curve is set to contain K discrete data points, each data point numbered k, where k = 1, 2, 3, ..., K; the intensity index value corresponding to each data point is I. k The effective diffusion coefficient value corresponding to each data point is D. k The intensity index calculated so far is denoted as I. We then find two adjacent data points in the correlation curve such that their intensity indices satisfy I. k ≤I≤I k+1 The current effective diffusion coefficient D is determined using linear interpolation, and its value lies between D0 and D1. k With D k+1The interpolation method ensures that the effective diffusion coefficient continuously increases with the increase of the sub-initiation disturbance intensity index. The interpolation result is the effective diffusion coefficient value at the current moment. Multiplying the effective diffusion coefficient by the initial concentration gradient at the interface yields the instantaneous pollutant release flux.

[0051] Next, closed-loop verification and iterative correction are performed. The instantaneous release flux is substituted backwards into the bed pore water concentration decay equation. This equation starts with the measured initial pore water concentration and uses the release flux as the consumption term, calculating backwards along the time axis. The calculation continues from the current moment until the backward-calculated pore water concentration decrease equals the measured concentration decrease. The elapsed time is recorded; this time is the duration of the micro-perturbation required for the backward calculation. The calculation rule for this backward calculation process is: set the initial bed pore water concentration as C0; the measured pore water concentration at the current moment is C... t Concentration decrease ΔC = C0 - C t The flux released is F; the duration of the back-calculation is T. b Starting from C0 and using F as the concentration consumption rate, we trace backward along the time axis to reach T. b At time t, the total decrease in pore water concentration equals ΔC; record the time T for this reverse calculation. b This duration is related to the measured total duration T of turbulent bursts. m The deviation between the two is δ=|T b -T m |

[0052] Simultaneously, the durations of all events in the measured turbulent burst event sequence from the first step are statistically analyzed and summed to obtain the total measured duration of the turbulent burst. The time obtained by reverse calculation is then compared with this total measured duration, and the deviation between the two is calculated.

[0053] This deviation value must simultaneously satisfy two preset constraints: the first constraint is an energy conservation constraint, meaning the deviation value cannot exceed the difference between the total energy input by the micro-disturbance and the total energy dissipated by the bed deformation. The second constraint is a mass conservation constraint, meaning the deviation value cannot exceed the upper limit of the allowable deviation between the flux back-calculation time and the measured disturbance time. The above two constraints are determined according to the following rules:

[0054] The first constraint is an energy conservation constraint: the total energy of the micro-perturbation input is set to E. in The total energy dissipated by bed deformation is E. diss The upper limit of the allowable difference range is ΔE. max This constraint requires: E in -E diss ≤ΔE max ;

[0055] The second constraint is the mass conservation constraint: the time for back-calculating the flux release is set to T. b The measured total duration of turbulent sudden onset is T. m The upper limit of the allowable deviation is ΔT. max This constraint requires: |T b -T m |≤ΔT max The deviation value is considered to fall within the preset threshold range only when both of the above constraints are met, the iteration is terminated and the risk level is output; if either constraint is not met, the iteration correction is started.

[0056] If the deviation value does not meet any of the above constraints, it is judged as out of tolerance, and iterative correction is initiated. The correction rules are as follows: if the back-calculated duration is greater than the measured duration, the effective yield threshold is appropriately increased, and the effective thickness of the diffusion boundary layer is correspondingly enlarged to reduce the released flux in subsequent calculations; if the back-calculated duration is less than the measured duration, the effective yield threshold is appropriately decreased, and the effective thickness of the diffusion boundary layer is correspondingly reduced to increase the released flux in subsequent calculations. After correction, the sub-start disturbance intensity index calculation, released flux calculation, and reverse back-substitution comparison steps are re-executed.

[0057] This process is repeated until the deviation value simultaneously meets two preset constraints. At this point, the iteration terminates, and the final calculated instantaneous pollutant release flux is used as the output value. Based on the preset numerical range within which this output value falls, the release risk level of the assessment area is determined. This preset numerical range is pre-defined as follows: several sediment samples with known risk levels are collected under different hydrodynamic conditions; the instantaneous release flux of each sample is calculated using this method; the calculation results are compared with the known risk levels to determine the flux numerical boundaries corresponding to each risk level. In this embodiment, the flux is divided into three risk levels—low risk, medium risk, and high risk—from low to high. Overlapping buffer zones are set between adjacent intervals to avoid level jumps caused by edge fluctuations.

[0058] Example 2, based on Example 1, uses a graded weighting method instead of a single statistical method when calculating the sub-start disturbance intensity index. The process is as follows:

[0059] After obtaining the sequence of turbulent burst events within a unit time window, the peak energy value corresponding to each event in the sequence is first extracted. These values ​​are then arranged in ascending order, and the minimum and maximum peak energies in the sequence are determined. The numerical span between the minimum and maximum peak energies is evenly divided into four consecutive energy level intervals. The number of levels should not be too large, otherwise the sample size within each level will be too small, leading to statistical instability; nor should it be too small, otherwise it will be unable to reflect the differentiated contributions of disturbances of different intensities. Sensitivity tests under different hydrodynamic conditions show that the four levels can achieve a balance between computational stability and distinguishing granularity. Each level interval has the same energy width, ensuring that each burst event is uniquely assigned to one of the level intervals based on its peak energy magnitude.

[0060] After classifying the events into different intensity levels, for each level interval, all burst events with peak energies falling within that interval are selected, forming a subset of events for that level. For this subset, following the same procedure as the main process, the cumulative proportion of events whose peak energies exceed the current effective yield threshold is calculated. Simultaneously, a positive feedback coefficient for pore water pressure on microfracture propagation is introduced, and this cumulative proportion is multiplied by the positive feedback coefficient to obtain the initial intensity value for that level. This initial intensity value is then divided by the largest peak energy among all events in the subset to complete the normalization process within that level range, yielding the intensity component corresponding to that level.

[0061] Subsequently, the time weight corresponding to each level interval is calculated. By summing the durations of all events in the subset of events for that level, the total duration of that level is obtained. This total duration is then divided by the total length of the entire time window, and the resulting ratio is used as the time proportion weight for that level. This weight reflects the time share of the disturbance of that level intensity throughout the entire observation period.

[0062] After obtaining the intensity components and time-based weights for each level, the intensity component of each level is multiplied by its corresponding time-based weight to obtain the weighted contribution value for that level. The weighted contribution values ​​for all levels are then summed to obtain the comprehensive sub-initiation disturbance intensity index. This comprehensive index replaces the previously calculated single index and is substituted into the dynamic equation of the effective diffusion coefficient at the interface. By adopting this level-based weighting method, the contributions of turbulent bursts of different magnitudes to pollutant release are treated differently based on their actual duration, resulting in a smoother calculated response for release fluxes and avoiding excessive fluctuations in release fluxes caused by an excessively high proportion of individual extreme burst events in the calculation.

[0063] Example 3, based on Example 1, introduces bottom turbulence burst frequency and pulsating pressure amplitude as auxiliary parameters in addition to the basic hydrodynamic parameters. These parameters are used for integrity verification of event extraction and as backup inputs for energy calculation, respectively. The process is as follows:

[0064] Regarding the auxiliary application of burst frequency, during continuous monitoring, the system continuously counts the total number of turbulent burst events extracted within the current time window. This total number divided by the window duration gives the current burst frequency. Simultaneously, the system maintains a sliding historical database, storing burst frequency values ​​for several consecutive time windows prior to the current moment. The average and standard deviation are calculated from these historical values ​​to serve as the historical baseline and its fluctuation range.

[0065] When the current burst frequency extracted in real time deviates from the historical baseline average, the system compares the deviation to the standard deviation. If the deviation does not exceed 2.5 times the standard deviation, the current event extraction result is considered complete and valid. If the deviation exceeds 2.5 times the standard deviation, it is determined that the event extraction process may have missed or false detections. The basis for this 2.5 times value is: for a random process that approximately follows a normal distribution, the probability of an observation exceeding the range of the average plus or minus 2.5 times the standard deviation is less than one percent. When the current burst frequency falls outside this range, it indicates a high probability of systematic bias in the event extraction process, necessitating supplementary extraction or removal operations.

[0066] At this point, the system automatically lowers the initial amplitude detection threshold used to identify burst events. The lowering method involves progressively decreasing the threshold by 5% of the current value. After each reduction, the system rescans the original flow velocity pulsation signal within the current time window, adding low-amplitude events previously missed due to an excessively high threshold to the sequence, while removing high-frequency noise points mistakenly included due to an excessively low threshold. After completing these addition and removal operations, the burst frequency of the current window is recalculated and compared with the historical baseline again. This adjustment process is repeated until the recalculated burst frequency falls within the range of the historical baseline average plus or minus 2.5 standard deviations. At this point, the burst event sequence of the current window is confirmed to meet the integrity requirements, and the sequence is then used for subsequent coupling mapping and intensity index calculations.

[0067] Regarding the auxiliary application of pulsating pressure amplitude: When deploying the monitoring device, pressure sensors are simultaneously installed at the bed surface to collect real-time pulsating pressure time-history data. This time-history data is squared and then multiplied by the water density value to obtain an equivalent estimate of turbulent pulsating energy. This estimate serves as a backup input for the conventional energy value calculated based on flow velocity data.

[0068] Simultaneously, the system performs real-time evaluation of the flow velocity data quality. The evaluation method involves monitoring the echo intensity signal from the acoustic Doppler current meter. When the echo intensity falls below a set threshold, the current flow velocity data quality is deemed insufficient to support conventional energy calculations. In this case, the equivalent energy estimate obtained from the aforementioned calculation based on pulsating pressure amplitude is directly fed into the energy-weakening function, replacing the original flow velocity-based energy value, as the driving variable for the current moment to calculate the instantaneous reduction in the yield threshold. If the flow velocity data quality recovers, the system switches back to the conventional flow velocity-based energy calculation method.

[0069] As a further refinement of the pressure amplitude application, when the flow velocity data quality is at a critical level, a complete switch can be avoided, and a fusion approach can be adopted instead. In this case, a weighted average is taken between the conventional energy value calculated based on flow velocity and the equivalent estimate calculated based on pressure amplitude. The weighting coefficients of the two are dynamically determined based on the specific value of the current echo intensity. The method is as follows: the echo intensity is normalized to the range of zero to one, the weight of the conventional energy value based on flow velocity is equal to the normalized echo intensity value, and the weight of the equivalent pressure amplitude is equal to one minus the normalized value; the lower the echo intensity, the greater the weight of the equivalent pressure amplitude, and the sum of the two weights is always one. The fused weighted average is then fed into the energy-weakening function. This fusion method can avoid data jumps during the switching process and further ensure the continuity of the weakening function input value over time. The above two refined application methods of auxiliary parameters can be used independently in conjunction with the main process, or they can be used simultaneously to enhance the adaptability of the entire evaluation method under different hydrodynamic monitoring conditions.

[0070] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for quantitatively assessing the risk of pollutant release from sediment based on sediment initiation modes, characterized in that, Includes the following steps: Step 1: Obtain hydrodynamic parameters, bed structure parameters, and pollutant activity parameters at the sediment interface in the assessment area; Step 2: Extract the time series of turbulent burst events based on the hydrodynamic parameters, determine the yield threshold of the surface microstructure based on the bed structure parameters, and determine the initial concentration gradient at the interface based on the pollutant activity parameters; Step 3: Establish the coupling mapping relationship between the hydrodynamic parameters, the bed structure parameters and the pollutant activity parameters. Using the turbulent burst event as the excitation boundary, dynamically modulate the yield threshold, calculate the micro-disturbance erosion thickness with the modulated yield threshold, and then correct the effective thickness of the diffusion boundary layer with the micro-disturbance erosion thickness. Step 4: Calculate the sub-initiation disturbance intensity index based on the time series of the turbulent burst event and the modulated yield threshold; Step 5: Substitute the sub-start disturbance intensity index into the dynamic equation of the effective diffusion coefficient of the interface to calculate the instantaneous release flux of pollutants; Step Six: Substitute the instantaneous release flux of the pollutant back into the bed surface pore water concentration decay equation to deduce the required micro-disturbance duration. Compare the deduce micro-disturbance duration with the measured turbulent burst duration. If the deviation exceeds a preset threshold, correct the yield threshold and the effective thickness, and repeat the calculation of the sub-start disturbance intensity index and subsequent steps until the deviation falls within the preset threshold range. Output the final release risk level.

2. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 1, characterized in that: The coupling mapping relationship is established as follows: the pulsation intensity of the turbulent burst event is used as the first driving variable to dynamically reduce the yield threshold of the surface microstructure, and the reduction magnitude is positively correlated with the pulsation intensity; at the same time, the amount of fine particles released by micro-disturbance erosion is used as the weighting coefficient to compress the effective thickness of the diffusion boundary layer, and the compression magnitude is positively correlated with the amount of fine particles released.

3. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 2, characterized in that: The micro-disturbance erosion thickness is calculated as follows: the upper limit of the particle size that can be eroded is determined based on the modulated yield threshold, and then the particle size distribution curve of the surface micro-aggregates is combined with the product of the volume ratio of each particle size component with a particle size smaller than the upper limit and the corresponding layer thickness to obtain the micro-disturbance erosion thickness.

4. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 3, characterized in that: The sub-start disturbance intensity index is calculated as follows: the cumulative probability density of turbulent burst energy exceeding the yield threshold is statistically analyzed per unit time, and the positive feedback coefficient of bed pore water pressure on microcrack expansion is introduced. The cumulative probability density is multiplied by the positive feedback coefficient and then normalized to obtain the sub-start disturbance intensity index.

5. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 1, characterized in that: The preset threshold is a combination of the energy conservation constraint boundary and the mass conservation verification threshold. The energy conservation constraint boundary is used to limit the range of the difference between the micro-disturbance input energy and the energy dissipated by the bed deformation. The mass conservation verification threshold is used to limit the upper limit of the allowable deviation between the disturbance duration calculated by the release flux and the measured disturbance duration.

6. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 5, characterized in that: When repeatedly performing the calculation of the sub-start disturbance intensity index and subsequent steps, the yield threshold is corrected by increasing or decreasing the value of the yield threshold according to the deviation direction between the previously back-calculated disturbance duration and the measured disturbance duration; the effective thickness is corrected by decreasing or increasing the value of the effective thickness according to the magnitude by which the previously output instantaneous release flux exceeds the preset range.

7. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 1, characterized in that: The hydrodynamic parameters include at least the instantaneous shear stress spectrum of bottom turbulence and the near-bed velocity gradient; the bed structure parameters include at least the particle size distribution of surface micro-aggregates and the degree of extracellular polymer cementation index; and the pollutant activity parameters include at least the ratio of dissolved pollutant concentration to particulate pollutant concentration in the surface sediments.

8. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 7, characterized in that: In the dynamic equation of the effective diffusion coefficient of the interface, the effective diffusion coefficient increases continuously with the increase of the sub-start disturbance intensity index, and the effective diffusion coefficient is obtained by interpolation and table lookup using the sub-start disturbance intensity index as the independent variable through a pre-calibrated sediment type correlation curve.

9. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 1, characterized in that: The hydrodynamic parameters also include the bottom turbulence burst frequency and the pulsating pressure amplitude. The bottom turbulence burst frequency is used to extract the time series of the turbulence burst event, and the pulsating pressure amplitude is used to calculate the turbulence burst energy.

10. The method for quantitatively assessing the risk of sediment pollutant release based on sediment initiation mode according to claim 1, characterized in that: The time series of the turbulent burst events is divided into four levels according to the burst intensity. The sub-start disturbance intensity component under each level is calculated. Then, the components are weighted and summed according to their time proportion to obtain the comprehensive sub-start disturbance intensity index. The comprehensive sub-start disturbance intensity index is used to replace the sub-start disturbance intensity index and is substituted into the dynamic equation of the effective diffusion coefficient of the interface.

Images