method for simulating mutual feedback between topographic erosion and habitat succession

Abstract

This invention discloses a method for simulating the interaction between topographic erosion and sedimentation and habitat succession. The method includes: acquiring initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass features of a target water area; converting the vegetation biomass features into hydrodynamic resistance terms, and inputting them, along with the initial topography and hydrodynamic conditions, into a pre-constructed three-dimensional dynamic geomorphological numerical model to calculate updated topographic features; inputting the updated topographic features into a pre-constructed staged vegetation growth dynamics model to calculate updated vegetation biomass features; and feeding back the updated vegetation biomass features in real time to the three-dimensional dynamic geomorphological numerical model to update the hydrodynamic resistance terms, driving the calculation of the next time step. This invention establishes a two-way feedback loop between hydrodynamics, geomorphology, and habitat, improving the accuracy and numerical computation stability of long-term, multi-scale ecological geomorphological evolution simulations.

Description

Technical Field

[0001] This invention belongs to the field of estuarine and coastal dynamic geomorphology and numerical simulation technology, and in particular, it is a method for simulating the interaction between topographic erosion and sedimentation and habitat succession. Background Technology

[0002] The habitat systems of estuaries and other land-sea confluence areas are directly constrained by water flow shear and sediment transport, while surface cover morphology also constrains local flow field structure and material deposition. Revealing and reproducing the coupling mechanisms of non-steady hydrodynamics, sediment transport, topographic evolution, and ecological processes across multiple timescales can improve the accuracy of estuarine evolution prediction and optimize the long-term hydrological response assessment of water-related projects.

[0003] Current dynamic geomorphological numerical models, when dealing with water-related ecological areas, typically generalize the bottom surface environment as static boundary conditions, often using globally constant roughness parameters or fixed artificial empirical resistance terms to approximate the impediment of water flow by the bottom surface. At the level of sediment dynamics calculation, existing studies generally rely on changes in bottom stress to indirectly deduce erosion and deposition evolution, and in terms of time stepping strategies, they mainly follow the traditional globally unified acceleration factor to drive long-term topographic extrapolation.

[0004] However, existing technologies suffer from computational challenges due to the potential for spatiotemporal distortion and numerical divergence during long-cycle simulations. Therefore, further research and innovation are needed to alleviate these problems. Summary of the Invention

[0005] The objective of this invention is:

[0006] This paper proposes a simulation method for the interaction between topographic erosion and sedimentation and habitat succession, which can at least alleviate the computational problems caused by long-term simulations in existing technologies, such as spatiotemporal distortion and numerical divergence.

[0007] The technical solution of this application is:

[0008] A method for simulating the interaction between topographic erosion and sedimentation and habitat succession includes:

[0009] Acquire the initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass characteristics of the target water area;

[0010] The initial vegetation biomass characteristics are transformed into hydrodynamic resistance terms. The hydrodynamic resistance terms, along with the initial topographic features and initial hydrodynamic conditions, are input into a pre-constructed three-dimensional dynamic geomorphological numerical model to calculate the updated topographic features.

[0011] The updated terrain features are input into a pre-constructed staged vegetation growth dynamics model to calculate the updated vegetation biomass features.

[0012] The updated vegetation biomass characteristics are fed back in real time to the three-dimensional dynamic geomorphological numerical model to update the hydrodynamic resistance term, driving the evolution calculation of the next time step, forming a two-way feedback closed loop of topographic erosion and sedimentation and habitat succession, and outputting dynamic simulation results of topography and habitat.

[0013] The beneficial effects of this application are:

[0014] This method enables high-frequency capture of ecological mutation processes without divergence at multi-decade simulation scales, aiming to alleviate the computational burden of spatiotemporal distortion and numerical divergence caused by long-term simulations. Attached Figure Description

[0015] 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:

[0016] Figure 1 This is an execution flowchart of a method for simulating the interaction between topographic erosion and sedimentation and habitat succession according to an embodiment of this application;

[0017] Figure 2 This is an execution flowchart of a colonization stage model according to an embodiment of this application;

[0018] Figure 3 This is an execution flowchart of a root development stage model according to an embodiment of this application;

[0019] Figure 4 This is an execution flowchart of a dynamic growth stage model according to an embodiment of this application;

[0020] Figure 5 This is a flowchart illustrating the calculation of environmental carrying capacity according to an embodiment of this application;

[0021] Figure 6 This is an execution flowchart of converting initial vegetation biomass characteristics or updated vegetation biomass characteristics into hydrodynamic resistance terms according to an embodiment of this application.

[0022] Figure 7 This is an execution flowchart illustrating how a graphics processor can be used to accelerate the solution process of a feedback closed loop in parallel, according to an embodiment of this application. Detailed Implementation

[0023] 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.

[0024] 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.

[0025] To overcome the aforementioned problems in the existing technology, the inventors, through literature review and technical analysis, verified that existing related solutions generally have the following shortcomings:

[0026] Conventional numerical systems often struggle to capture local flow structure variations during succession when dealing with high-frequency, large-gradient changing conditions, and also fail to accurately map the nonlinear effects of boundary deformation on the overall sediment transport rate.

[0027] Existing numerical models have limitations in terms of dynamic feedback mechanisms at multiphase interfaces and long-term computational stability. Especially when facing complex conditions such as high-frequency hydrodynamic stress and severe geomorphic deformation, traditional algorithms have low process mapping fidelity and are prone to numerical divergence and spatiotemporal distortion in long-term simulations.

[0028] To overcome the shortcomings of existing technologies, referring to Figures 1-7 The implementation of this application is described in detail using the following embodiments.

[0029] The following is a description of some of the parameters involved in this application:

[0030] Among them, the characteristic parameter of eddy current oscillation force is the KC number.

[0031] For example, a method for simulating the interaction between topographic erosion and sedimentation and habitat succession specifically includes the following steps:

[0032] Step 101: Obtain the initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass characteristics of the target water area.

[0033] In practice, the environmental baseline data of the target water area is read. This target water area usually refers to a key ecological area of ​​land-sea interaction, which may be a estuary shoal or a salt marsh wetland.

[0034] The initial terrain features can be specifically represented as the spatially discretized grid elevation field matrix.

[0035] Initial hydrodynamic conditions include external wave boundary field parameters such as significant wave height, spectral peak period, and wave direction, as well as tidal level and runoff boundary conditions.

[0036] The initial vegetation biomass characteristics reflect the spatial distribution density and absolute mass value of vegetation at each grid node at the start of the simulation.

[0037] Correspondingly, the data loading and initialization process can be deployed in the video memory of a graphics processor that supports a unified computing architecture, thereby providing a hardware parallel acceleration foundation for multi-scale dynamic evolution.

[0038] Step 102: The initial vegetation biomass characteristics are converted into hydrodynamic resistance terms, and the hydrodynamic resistance terms, along with the initial topographic features and initial hydrodynamic conditions, are input into the pre-constructed three-dimensional dynamic geomorphological numerical model to calculate the updated topographic features.

[0039] In other cases, this step can also yield the corresponding hydrodynamic state parameters;

[0040] Among them, hydrodynamic state parameters can be state quantities that characterize the dynamic changes of water bodies in the study area, including but not limited to water level, water depth, water flow velocity, inundation duration, inundation frequency, and water level rise and fall rate.

[0041] Accordingly, the initial vegetation biomass characteristics are extracted, converted into spatially distributed additional drag source terms, and directly loaded into the bottom stress term of the tidal momentum equation.

[0042] Among them, the three-dimensional dynamic geomorphological numerical model is responsible for simulating the dynamic interaction between the unsteady flow field and the bed boundary. The wave field can be handled by a combination of parameterized input and analytical attenuation within the model.

[0043] Furthermore, the attenuation of waves within the vegetated area can be calculated using the exponential decay model, and the bed track velocity can be calculated using linear wave theory.

[0044] By combining the effects of water flow and waves, the non-equilibrium sediment transport mode is driven to perform empirical calculations of sediment initiation velocity and sediment carrying force, thereby solving the topographic scouring and deposition equations and outputting the scalar matrix of grid node elevation changes after the end of the current hydrodynamic calculation cycle, which can be defined as updating topographic features.

[0045] Step 103: Input the updated terrain features into the pre-constructed phased vegetation growth dynamics model to calculate the updated vegetation biomass features.

[0046] In other cases, the inputs in this step may also include hydrodynamic state parameters;

[0047] In this step, after the local topography undergoes erosion and deposition deformation, the uplift or downcutting of the bed elevation will directly change the actual water depth and inundation duration in the local area, thereby reshaping the microhabitat conditions of the ecological entity.

[0048] Among them, the phased vegetation growth dynamics model is used to receive updated terrain features and simultaneously refresh the upper limit of spatial environmental carrying capacity and hydrodynamic stripping stress index of each grid unit.

[0049] After undergoing screening of planting time windows, judgment of root erosion resistance critical thresholds, and dynamic growth model deduction that integrates multiple stress factors, the survival quality and spatial density of benthic vegetation undergo nonlinear evolution.

[0050] Next, the execution unit calculates the biomass stock after experiencing topographic constraints and natural metabolism, and outputs updated vegetation biomass characteristics.

[0051] Step 104: The updated vegetation biomass characteristics are fed back in real time to the three-dimensional dynamic geomorphological numerical model to update the hydrodynamic resistance term, so as to drive the evolution calculation of the next time step, forming a two-way feedback closed loop of topographic erosion and sedimentation and habitat succession, and outputting the dynamic simulation results of topography and habitat.

[0052] Accordingly, the characteristics of vegetation biomass are extracted and updated, and the bottom stress distribution matrix is ​​refreshed, thereby changing the local velocity distribution and wave propagation energy dissipation rate in the hydrodynamic flow field, which leads to the spatial reconstruction of suspended and bedload sediment transport rates in the next stage.

[0053] It should be understood that the evolution of hydrodynamics leads to changes in topography, which in turn constrains vegetation growth and vegetation succession, thus altering water flow resistance and causing secondary adjustments to the hydrodynamic field, forming a continuous process of two-way feedback loop in dynamics.

[0054] As the time step continues to advance in a loop, under the parallel scheduling strategy of the graphics processor, it continuously outputs grid elevation time series data and biomass spatial distribution maps covering a time scale of decades, namely the dynamic simulation results of terrain and habitat.

[0055] Another example includes an offline construction step to obtain the basic drag coefficients required for a three-dimensional dynamic geomorphological numerical model, before acquiring the initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass characteristics of the target water area.

[0056] Step 201: In a pre-constructed wave and water flow channel, different levels of wave conditions and unidirectional flow conditions are combined and set to determine the water flow resistance pattern of the vegetation area under different layout forms.

[0057] In practice, before starting the online dynamic geomorphological simulation, a wind and wave flow tank needs to be built in the experimental environment.

[0058] In some scenarios, the specific dimensions of the wave channel can be set as follows: 175m long, 1.2m wide, and 1.6m high.

[0059] Other size options may also be selected; as long as the wave propagation requirements are met, no single size is required.

[0060] Furthermore, vegetation models in non-submerged or semi-submerged states were placed in the experimental section at the bottom of the water tank.

[0061] For example, vegetation arrangement can be divided into neat arrangement and staggered arrangement. The neat arrangement specifically adopts a square grid arrangement, while the staggered arrangement specifically adopts a triangular grid arrangement with adjacent rows staggered by half a distance.

[0062] During the testing phase, regular waves with preset wave height and period are input by controlling the wave-generating equipment at the front end of the water tank.

[0063] For example, the wave height can be between 0.03 and 0.25 m; the period can be between 0.8 and 2.5 seconds.

[0064] At the same time, the system applies unidirectional flow velocities of different directions and magnitudes.

[0065] Among them, the hydrodynamic combination conditions include pure wave conditions, conditions where waves and unidirectional currents in the same direction are superimposed, and conditions where waves and unidirectional currents in the opposite direction are superimposed.

[0066] Furthermore, by deploying three-dimensional acoustic Doppler current meters and wave height meters in the vegetated area, the three-dimensional velocity time series and wave attenuation amplitude of the flow field can be collected.

[0067] Based on the collected data, the distribution of stress state along the path of vegetation communities under various hydrodynamic combinations is extracted, and the corresponding water flow resistance curve is generated.

[0068] As an alternative, under predetermined flow patterns, the spatial distribution of output water flow resistance differs. Specifically, when waves and unidirectional flow are superimposed, the maximum resistance experienced by vegetation on the wave-facing side is greater than that experienced by vegetation downstream; when waves and reverse unidirectional flow are superimposed, the maximum resistance experienced by vegetation decreases monotonically along the wave propagation direction.

[0069] The above provides offline support for the spatial gradation assignment of drag coefficients in numerical simulations.

[0070] Step 202: Extract the Reynolds number and KC number that are correlated with the water flow resistance law to a preset standard;

[0071] It is necessary to screen out key dimensionless control parameters from the water flow resistance law, so as to convert the experimental results of the water tank model into empirical formulas that can be implanted into the numerical model.

[0072] Among them, the Morrison equation can be used to separate the wave resultant force on a single plant into a drag term and an inertial term. After removing the inertial force interference, the corresponding measured vegetation drag force coefficient can be calculated in reverse.

[0073] After obtaining the measured vegetation drag force coefficient sequence, a multidimensional regression matrix of the coefficient and the physical state parameters of the flow field is constructed.

[0074] Multiple dimensionless parameters, including relative wave height, relative wavelength, Froude number, Reynolds number, and KC number, are included in the correlation determination criteria.

[0075] Statistical tests revealed that the correlation coefficient between relative wave height and Froude number and resistance fluctuation was less than 0.1, so these waves were excluded.

[0076] The Reynolds number, which reflects the viscous force of the fluid, and the KC number, which reflects the periodic characteristics of the waves, have a correlation test index that exceeds the set lower limit threshold. They are retained as core feature quantities, that is, the Reynolds number and the KC number are extracted.

[0077] Step 203: Based on the Reynolds number and KC number fitting, an empirical formula is established to obtain the pre-constructed basic drag force coefficient, and the basic drag force coefficient is configured into the three-dimensional dynamic geomorphological numerical model.

[0078] Accordingly, after locking the core dimensionless parameters, two-dimensional scatter plots of Reynolds number and measured vegetation drag coefficient, and KC number and measured vegetation drag coefficient were established respectively.

[0079] Based on the convergence trend of the scatter distribution, a family of nonlinear functions is selected for least-squares optimization. For the Reynolds number, an inverse proportional function model that is negatively correlated with it is established.

[0080] For the KC number, a power function model in negative exponential form is established. In one case, after parameter optimization and convergence, the corresponding empirical formula is generated as follows:

[0081] C _d,v =28.32 / Re+0.6717;

[0082] In the above formula, C _d,vRe represents the vegetation base drag coefficient, and Re represents the Reynolds number.

[0083] In another case, C _d,v =36.4×KC -0.7245 ;

[0084] In the above formula, C _d,v KC represents the drag coefficient of the vegetation base.

[0085] By applying the above empirical formulas, a database of offline parameters for calculating the basic drag force coefficient is established.

[0086] In actual use, when the online simulation system can obtain the instantaneous Reynolds number of the flow field, the basic drag force coefficient can be calculated using an empirical formula based on the Reynolds number.

[0087] When the simulated working condition is dominated by wave motion and the KC number can be directly determined by the wave period and amplitude, an empirical formula based on the KC number is used for calculation; an appropriate empirical formula can be selected according to the characteristics of the actual hydrodynamic boundary conditions.

[0088] After offline construction is completed, the computational library is directly deployed into the three-dimensional dynamic geomorphological numerical model. Under the baseline condition that the vegetation is assumed to be a rigid or semi-rigid body, it provides the initial assignment basis for the fluid dynamic resistance loading process inside the model.

[0089] According to an embodiment of this application, an optional implementation method is described in which updated terrain features are input into a pre-constructed staged vegetation growth dynamics model to calculate updated vegetation biomass features.

[0090] The phased vegetation growth dynamics model includes a planting stage model, a root development stage model, and a dynamic growth stage model.

[0091] The life cycle of vegetation in tidal flat environments exhibits heterogeneous dynamic characteristics.

[0092] Optionally, a model module containing three timing succession phases may be invoked;

[0093] Based on the real-time growth time and survival status of vegetation in the grid cells, this module dynamically distributes the data to the corresponding computation branches to perform habitat evolution simulation.

[0094] The colonization phase model is specifically used to execute the following steps:

[0095] Step 301: Based on the water level time series results output by the three-dimensional dynamic geomorphology numerical model, extract the duration of no flooding in the current grid cell, and compare the duration of no flooding with the preset planting opportunity time threshold to determine the planting opportunity window.

[0096] Optionally, time-series water level data within the hydrodynamic calculation cycle can be analyzed to statistically determine the continuous duration of the target grid exposed above the water surface after low tide. This statistical value can be compared with the minimum unsubmerged time threshold required for species germination. If the actual exposure time exceeds the set threshold, the calculation channel is opened to allow the ecological parameters to enter the probability determination stage; if it is less than or equal to the set threshold, it is determined that the conditions for ecological germination are not met in this cycle.

[0097] Step 302: Under the condition of satisfying the planting opportunity window, obtain the flooding frequency of the current grid cell, and determine the planting probability of the current grid cell based on the pre-configured correspondence between flooding frequency and planting probability.

[0098] After the planting window is opened, the tidal inundation frequency index accumulated over the long-term evolution of the grid is extracted.

[0099] Among them, the probability of establishment and the frequency of flooding show a monotonically negative correlation.

[0100] Optionally, a nonlinear sigmoid function is used to fit the sampling characteristic curve, and the specific calculation formula is as follows:

[0101] P _col =1 / (1+exp(12×(f _inun -0.5)));

[0102] In the above formula, P _col f represents the colonization probability. _inun This represents the flooding frequency, the constant 12 represents the curve steepness adjustment parameter, and 0.5 represents the logical inflection point for determining the flooding frequency.

[0103] Accordingly, the output of this function is extracted, and when the flooding frequency is close to 1, the calculated probability converges to 0.

[0104] When the frequency of flooding approaches 0, the probability converges to 1.

[0105] Step 303: Determine the vegetation establishment result based on the establishment probability, and output the initial seedling density as the updated vegetation biomass characteristic of the grid cell when the establishment is successful.

[0106] Optionally, the system generates uniformly distributed pseudo-random numbers within the interval of 0 and 1;

[0107] The pseudo-random number is compared with the obtained planting probability.

[0108] When the random number is less than or equal to the planting probability, the state of the grid is changed to planting survival, and the pre-configured planting stock value is assigned to the grid. The updated vegetation biomass characteristics are output for subsequent time steps.

[0109] When the random number is greater than the planting probability, the grid remains unvegetated.

[0110] The root development stage model is specifically used to perform the following steps:

[0111] Optionally, the growth time of the vegetation is obtained, and the critical impact shear stress at the current moment is calculated based on a preset rule that the critical impact shear stress of the seedling increases linearly with the growth time.

[0112] The growth period of vegetation can be recorded cumulatively from the time of successful planting;

[0113] During the period between successful vegetation establishment and the maturity stage, the seedlings enter the root development model branch. During this stage, the vegetation's pull-out resistance is dynamically increasing.

[0114] Furthermore, a linearly increasing function can be invoked once to update the impact shear stress threshold at the current moment, i.e.:

[0115] τ _b,crit,ph2 =a _τ ×(T _t -T _ph1 )+b _τ ;

[0116] In the formula, τ _b,crit,ph2 That is, the critical impact shear stress, a _τ That is, the linear growth rate of the root system's resistance to erosion, T _t That is, the current actual growth time of the vegetation, T _ph1 That is, the time of establishment and development onset, b _τ This refers to the basic pull-out resistance value during the early stages of root development.

[0117] Optionally, obtain the actual bed shear stress of the current mesh element;

[0118] Correspondingly, the scalar sum of the bottom stress tensor transmitted by the underlying hydrodynamic interaction system is read synchronously and used as the input of the destructive force load applied to the seedling surface by the current instantaneous water flow.

[0119] Optionally, the actual bed shear stress is compared with the critical impact shear stress, and the survival status or scour-peeling status of the seedlings is determined based on the comparison results. If the seedlings are determined to be in a survival status, their growth is maintained.

[0120] Furthermore, the difference between the real-time extracted destructive force and the dynamically calculated biological pull-out force is determined;

[0121] When the actual destructive force exceeds the shear resistance limit, the seedling density data of the grid is forcibly reset to zero, and the scour and stripping event is recorded.

[0122] When the actual destructive force fails to exceed the upper limit, the ecological entity data is retained, and the model is advanced to the next iteration step after accumulating the growth time.

[0123] The dynamic growth phase model is specifically used to perform the following steps:

[0124] Construct a logistic growth term based on the maximum growth rate of vegetation;

[0125] Construct a natural mortality item based on a preset natural mortality rate;

[0126] Obtain the current actual water depth output from the three-dimensional dynamic geomorphological numerical model, and construct an inundation loss term based on the difference between the maximum tolerable water depth of vegetation and the current actual water depth.

[0127] The hydrodynamic scour loss term is constructed based on the difference between the maximum bed shear stress that the vegetation can withstand and the current actual bed shear stress.

[0128] Based on logistic growth, natural mortality, inundation loss, and hydrodynamic scour loss terms, a biomass spatiotemporal variation control equation is constructed. Combined with environmental carrying capacity constrained by both elevation and sediment thickness, the updated vegetation biomass characteristics are calculated and output.

[0129] Accordingly, once the seedlings have crossed the root development stage and entered the maturity stage, the process is transferred to the dynamic growth stage model. The absolute mass increment of the ecological entity can be iteratively calculated using ordinary differential equations, with the corresponding formula as follows:

[0130] dB / dt=r _max ×B×(1-B / K)-d×B-max(hh _crit ,0)×PE _h -max(τ-τ _crit ,0)×PE _τ ;

[0131] Within the formula:

[0132] dB / dt represents the rate of change of biomass per unit time step;

[0133] r _max Indicates the maximum growth rate of vegetation;

[0134] B represents the current vegetation biomass density, i.e., the vegetation biomass characteristics;

[0135] K represents environmental carrying capacity, including r _max The combination of B and K represents the logistic growth term;

[0136] d represents the natural mortality rate, corresponding to the natural death item;

[0137] h represents the current actual water depth;

[0138] h _crit Indicates the maximum water depth that can be tolerated;

[0139] PE _h This represents the inundation stress coefficient, and the corresponding product term represents the inundation loss term.

[0140] τ represents the current actual shear stress on the bed surface;

[0141] τ _crit This indicates the maximum bed surface shear stress that it can withstand;

[0142] PE _τ This represents the scouring stress coefficient, and the corresponding product term represents the hydrodynamic scouring loss term.

[0143] Furthermore, the judgment logic of the extraction function max is read. When the actual water depth or shear force does not exceed the limit, the corresponding loss term returns to zero. Once the tolerance threshold is exceeded, the biomass is increased proportionally to the magnitude of the excess to reduce the load.

[0144] In some scenarios, a non-negative lower bound constraint can be added, i.e. (B=max(B,0)), to avoid abnormal results of negative biomass caused by numerical errors.

[0145] The environmental carrying capacity is calculated in the following way:

[0146] Obtain the pre-configured ideal environmental carrying capacity, the actual elevation of the current grid cell, and the sediment thickness;

[0147] Based on the relative relationship between the actual elevation and the preset optimal vegetation growth elevation, the Gaussian distribution response term is calculated.

[0148] The nonlinear saturated response term is calculated based on the ratio of sediment thickness to a preset critical sediment thickness.

[0149] The environmental carrying capacity is calculated by multiplying the ideal environmental carrying capacity, the Gaussian distributed response term, and the nonlinear saturated response term.

[0150] Sediment thickness can be output from a three-dimensional dynamic geomorphological numerical model. To quantitatively describe the limitation of topographic morphology on biomass, a carrying capacity calculation logic containing multiple constraint functions was constructed, and the calculation formula is as follows:

[0151] K=K _0 ×exp(-β×(zz _opt ) 2 )×(h _s / (h _s +h _c ));

[0152] In the formula, K represents the environmental carrying capacity, K _0 Represents the ideal environmental carrying capacity, exp(-β×(zz) _opt ) 2 ) represents the Gaussian distributed response term, β represents the elevation stress coefficient, and z represents the actual elevation. _opt Indicates the optimal growth elevation, h _s / (h _s +h _c ) represents the nonlinear saturated response term, h _s h represents the thickness of the sediment. _c Indicates the critical thickness of the sediment.

[0153] In the above, when the elevation of the bottom grid deviates from the suitable survival range, the decay weight between 0 and 1 is calculated using an exponential function.

[0154] Simultaneously assess sediment deposition, limiting the upper limit parameter before the deposition thickness exceeds the critical point; after exceeding the critical point, this parameter tends to a constant of 1.

[0155] According to an embodiment of this application, the current actual bed shear stress in the hydrodynamic scour loss term is calculated in the following manner:

[0156] The instantaneous wave velocity amplitude, water flow velocity, and wave-current interaction correction coefficient of the current grid cell are obtained from the three-dimensional dynamic geomorphological numerical model.

[0157] The wave-current nonlinear superposition calculation is performed based on the instantaneous wave velocity amplitude and the water flow velocity to obtain the superimposed velocity term. The wave-current interaction correction coefficient is used to characterize the wave-current coupling strength during the superposition calculation.

[0158] Obtain the vegetation water-facing surface area that is positively correlated with the current vegetation biomass density; or, in other words, calculate the vegetation water-facing area projection per unit bed area based on the current vegetation biomass density; the vegetation water-facing area projection is the product of vegetation stem density, stem diameter and effective stem height, which characterizes the proportion of vegetation frontal shading of incoming flow per unit bed area.

[0159] The actual bed shear stress is calculated based on the product of water density, preset drag coefficient, vegetation surface area facing water, and superimposed flow velocity term.

[0160] In the shallow sea dynamic domain, to objectively reconstruct the hydrodynamic stripping efficiency, the velocity scalars of the wave field and current field are read concurrently. By introducing a correction coefficient to reconstruct the quadratic expansion formula, the engine's calculation formula is expressed as:

[0161] τ=ρ _w ×C _D ×A×((U _max ) 2 +u 2+2×γ×U _max ×u);

[0162] In the formula, τ represents the actual shear stress on the bed surface, and ρ _w C represents the density of water. _D U represents the drag coefficient. _max The instantaneous wave velocity amplitude is represented by γ, the water flow velocity is represented by u, the wave-current interaction correction coefficient is represented by γ, the polynomial combination inside the parentheses represents the superimposed velocity term, and A represents the water-facing projection area of ​​vegetation per unit bed area, which is specifically the product of vegetation stem density, stem diameter and effective interception height.

[0163] Furthermore, when extracting the data sequence, the sign of the water flow velocity directly reflects the flow direction. When encountering wave-current reversal conditions, the internal product term turns negative, thereby deducting the calculation error of artificially inflated energy caused by dynamic offsetting.

[0164] In one alternative implementation, the submergence stress coefficient and the scour stress coefficient can be determined by least squares fitting or regression analysis based on measured data of biomass decline of the target species under different water depths and shear stress conditions.

[0165] The elevation stress coefficient can be determined by fitting a Gaussian function to the measured elevation distribution bandwidth of the target species; the specific value can be adaptively adjusted according to the actual species characteristics and environmental conditions of the study area.

[0166] Among them, the wave-current interaction correction coefficient can be used to characterize the nonlinear coupling strength between wave trajectory motion and steady-state water flow.

[0167] In other words, γ takes a positive value when the wave and the water flow propagate in the same direction;

[0168] When waves and water flow propagate in opposite directions, γ takes a negative value; its absolute value can be determined by comparing and calibrating with the theoretical superposition model based on the measured bottom velocity data under wave-current coexistence conditions.

[0169] The above reconstructs the phased growth equation controlled by wave-current extreme stripping and dual environmental carrying capacity, capturing the dynamic response of vegetation from the seedling vulnerable stage to the mature stage, and extending the simulation dimension of ecological evolution to the level of complex systems with multiple species and multiple stresses.

[0170] In some embodiments, the species parameter configuration branches in the ecological evolution process are further explained. Since estuarine wetlands typically exhibit multi-species competition and spatial zonation, providing parallel parameter set branches supports differentiated response calculations for various ecological entities.

[0171] According to the embodiments of this application, the phased vegetation growth dynamics model is configured with multiple parallel parameter sets for estuarine salt marsh vegetation types. The estuarine salt marsh vegetation types include at least one of Suaeda salsa, Reed, or Spartina alterniflora. The parameter sets for different vegetation types respectively contain different flood tolerance elevation, drought tolerance elevation, and ideal environmental carrying capacity parameters.

[0172] When encountering an ecological entity botanically classified as Suaeda salsa, the corresponding pre-existing physiological boundary constraint matrix is ​​extracted.

[0173] In some scenarios, the flood-resistant elevation limit for this type can be set to 0.8m, the drought-resistant elevation limit to 2.0m, and the optimal growth elevation constraint range can be set to 1.2-0.6m.

[0174] In some application scenarios, when calculating the bearing capacity parameters during the dynamic growth stage, the ideal environmental bearing capacity can be configured between 0.8-1.5 kg / m². 2 Within the numerical range, the corresponding annual natural mortality rate is configured as 8-10%.

[0175] Furthermore, based on the specific shape of the Gaussian distribution attenuation curve constrained by the interval parameters, the corresponding environmental carrying capacity index is output.

[0176] When dealing with ecological entities botanically classified as reeds, another set of parallel parameter subsets is invoked across the grid. In some cases, the flood tolerance elevation limit for this type is set at 0.5m, the drought tolerance elevation limit is set at 1.8m, and the optimal growth elevation constraint range can be between 1.0 and 1.4m.

[0177] In some application scenarios, the ideal environmental load-bearing capacity can be adjusted to 2.5-4.0 kg / m². 2 The corresponding annual natural mortality rate is configured as 5-7%. The logistic equation operation under this branch is constrained by a wider elevation survival window.

[0178] When dealing with ecological entities botanically classified as Spartina alterniflora, constraint parameters with stronger hydrodynamic tolerance are read.

[0179] In some application scenarios, the flood-resistant elevation limit for this type can be set to 0.3m, the drought-resistant elevation limit can be set to 1.5m, and the optimal growth elevation constraint range is between 0.8-1.2m.

[0180] Its ideal environmental bearing capacity is 3.0-5.0 kg / m³. 2 The annual natural mortality rate is set at 1-3%.

[0181] Based on this, the dynamic model can quantitatively calculate the biomass stock of this species under high-frequency flooding conditions in the low-tide shoal area.

[0182] In this embodiment, based on the actual elevation section to which the grid cell belongs and the initial species allocation status, predetermined flood tolerance and bearing capacity parameters can be adaptively matched and issued to complete the quantitative simulation of the evolution of multi-species distribution patterns.

[0183] The values ​​and ranges given above are not the only limitations and can be adjusted comprehensively based on vegetation type and tidal conditions through field monitoring and literature research. Other values ​​and ranges in this application can be set and adjusted by technical personnel according to actual needs, industry standards and literature research. Specific calculations can be performed using well-known methods, which will not be elaborated on here.

[0184] According to the embodiments of this application, a cross-scale morphological transformation and hydrodynamic flexible avoidance correction mechanism is provided for real flexible vegetation groups in nature.

[0185] The specific process of converting initial vegetation biomass characteristics or regenerated vegetation biomass characteristics into hydrodynamic resistance terms includes:

[0186] Step 501: Obtain the pre-configured allometric growth mapping relationship and the basic drag coefficient of vegetation;

[0187] When processing real estuarine salt marsh vegetation data, the biomass characteristics output by the ecological model are two-dimensional scalars and cannot be directly used in the calculation of three-dimensional wave flow resistance.

[0188] Optionally, a predefined allometric growth mapping relationship can be read from memory. This relationship is used to describe the nonlinear geometric allometric distribution law of plant individual morphology with the accumulation of population biomass.

[0189] The dimensionless parameter matrix output from the previous water tank test calibration is retrieved simultaneously to extract the basic drag force coefficient.

[0190] Step 502: Through the allometric growth mapping relationship, the initial vegetation biomass characteristics or updated vegetation biomass characteristics, which are scalar data, are nonlinearly mapped to the three-dimensional morphological parameters of vegetation.

[0191] Accordingly, based on environmental carrying capacity indicators, the input biomass density scalar can be deconstructed into three-dimensional morphological parameters including stem density, stem height, and stem diameter.

[0192] The operation process performs nonlinear transformation of independent channels according to the allometric growth law.

[0193] When calculating spatial stem density, the ratio of biomass characteristics to environmental carrying capacity is extracted, and a power function including the density allometric index is used for calculation. The corresponding conversion formula is as follows:

[0194] N _s =N _s,max×(B / K) α_N ;

[0195] In the formula, N _s This represents the calculated stem density;

[0196] N _s,max This indicates the preset maximum stem density;

[0197] B represents vegetation biomass characteristics;

[0198] K represents environmental carrying capacity;

[0199] α _N The density allotropy represents the self-thinning effect of vegetation during its growth process.

[0200] When calculating the stem height of vegetation, an exponential saturation function can be used to reflect the growth distribution strategy of plants that first elongate longitudinally and then thicken laterally. The calculation formula is as follows:

[0201] h _v =h _v,max ×(1-exp(-λ _h ×B / K));

[0202] In the formula, h _v The calculated stem height;

[0203] h _v,max Indicates the preset maximum stem height;

[0204] B represents vegetation biomass characteristics;

[0205] K represents environmental carrying capacity;

[0206] λ _h This represents the height allocation constraint constant.

[0207] When calculating the stem diameter of vegetation, a power function scaling is applied based on a given extreme value range. The calculation formula is as follows:

[0208] d _s =d _s,min +(d _s,max -d _s,min )×(B / K) α_d ;

[0209] In the formula, d _s This represents the calculated stem diameter;

[0210] d _s,min Indicates the preset minimum stem diameter;

[0211] d _s,max Indicates the preset maximum stem diameter;

[0212] B represents vegetation biomass characteristics;

[0213] K represents environmental carrying capacity;

[0214] α _d This indicates the stem diameter allotropy index.

[0215] Step 503: Calculate the flexible reduction coefficient based on the three-dimensional morphological parameters and flow field parameters using the pre-constructed flexible reduction function;

[0216] Among them, the three-dimensional morphological parameters are obtained from the three-dimensional dynamic geomorphological numerical model;

[0217] Natural vegetation bends and deforms when subjected to strong hydrodynamic impacts to reduce water resistance on the upstream side. Furthermore, dimensionless numbers can be introduced based on fluid mechanics principles to construct a reduction discrimination logic.

[0218] For example, the stem height and stem diameter are retrieved from the three-dimensional morphological parameters, and the Cauchy number is calculated by combining the basic drag coefficient and the real-time bottom flow velocity characteristic parameters. The calculation formula can be expressed as:

[0219] Ca=(ρ×C _D ×d _s ×(u _rms ) 2 ×(h _v ) 3 ) / EI;

[0220] In the formula, Ca represents the calculated Cauchy number, which reflects the ratio of the fluid action torque to the vegetation elastic restoring torque;

[0221] ρ represents the fluid density;

[0222] C _D Indicates the drag coefficient;

[0223] d _s Represents the stem diameter in the three-dimensional morphological parameters;

[0224] u _rms Represents flow velocity characteristic parameters;

[0225] h _v Represents the stem height in the three-dimensional morphological parameters;

[0226] EI represents the preset bending stiffness of vegetation.

[0227] For example, after obtaining the Cauchy number of the current grid, the determination unit compares it numerically with a preset critical Cauchy number.

[0228] When the Cauchy number is determined to be less than or equal to the critical Cauchy number, the vegetation is determined to exhibit a rigid response, and the flexibility reduction factor is assigned to 1.

[0229] When the Cauchy number is determined to be greater than the critical Cauchy number, it is determined that the vegetation is undergoing flexible bending to avoid the obstacle, and the reduction equation is activated. The flexibility reduction coefficient is calculated using the ratio of the critical Cauchy number to the Cauchy number. The calculation formula can be expressed as:

[0230] φ=(Ca _cr / Ca) 1 / 3 ;

[0231] In the formula, φ is the flexibility reduction coefficient, and Ca _cr Here, Ca is the preset critical Cauchy number, and Ca is the Cauchy number calculated for the current grid.

[0232] The output value of this function is between 0 and 1, and is used to characterize the magnitude of drag reduction caused by bending deformation.

[0233] Step 504: Multiply the basic drag coefficient by the flexibility reduction coefficient to obtain the effective drag coefficient, and use the effective drag coefficient as the hydrodynamic resistance term.

[0234] Optionally, the basic drag force coefficient read from the offline library is extracted and multiplied with the flexibility reduction coefficient dynamically obtained according to the current flow state and three-dimensional morphology.

[0235] The resulting product eliminates the overestimation of energy dissipation under the rigidity assumption, and the output is the effective drag force coefficient. This coefficient matrix is ​​formatted and written into the underlying data bus, serving as the updated hydrodynamic drag term sent to the bottom stress calculation interface of the 3D dynamic geomorphological numerical model, supporting the subsequent elevation deformation extrapolation of the mesh.

[0236] The density allotropy index, height distribution constraint constant, and stem diameter allotropy index can be determined by nonlinear regression fitting based on field-measured biomass-morphology data of the target vegetation species. The critical Cauchy number can be determined by back-calculation based on the hydrodynamic conditions at which vegetation begins to bend and deform under fluid action.

[0237] In some embodiments, the three-dimensional morphological parameters include stem density, stem height, and stem diameter;

[0238] The process of performing nonlinear mapping through allometric growth mapping specifically includes:

[0239] The stem density is calculated based on the preset maximum stem density and the ratio of biomass characteristics to environmental carrying capacity;

[0240] Based on the preset maximum stem height and ratio relationship, the stem height is calculated using an exponential saturation function;

[0241] The stem diameter is calculated based on the preset minimum and maximum stem diameters and their ratio.

[0242] In this embodiment, allometric growth mapping and fluid Cauchy number reduction mechanisms are introduced to refine the characterization of water flow resistance from empirical mean to the physical scale. This is used to alleviate the problem of low mapping fidelity in existing models.

[0243] On the one hand, traditional models often rely on the roughness of the substrate to indirectly regulate sediment deposition. On the other hand, this application quantifies the interception effect of three-dimensional vegetation on suspended particles in water by reconstructing the mass conservation equation of the underlying layer.

[0244] Furthermore, the three-dimensional dynamic geomorphological numerical model includes suspended sediment transport equations for calculating sediment transport patterns and topographic erosion and deposition equations for calculating topographic evolution. Specifically, it involves performing the following steps:

[0245] Step 601: Obtain the direct interception efficiency, vegetation stem density, vegetation stem diameter, effective interception depth, flow velocity, and suspended sediment content of the current target water area.

[0246] Optionally, vegetation morphology geometric parameters of the target grid are extracted from the spatial allotropic mapping matrix to obtain vegetation stem density and vegetation stem diameter;

[0247] The height of the vegetation structure can be compared with the actual water depth at the current tide level, and the minimum value between the two can be extracted as the effective interception depth.

[0248] Simultaneously, the hydrodynamic solver outputs the three-dimensional velocity vector of the grid node at the current moment and the concentration state of suspended sediment in the water body, i.e., suspended sediment content.

[0249] Step 602: Based on the multiplicative relationship between direct interception efficiency, vegetation stem density, vegetation stem diameter, effective interception water depth, flow velocity and suspended sediment content, calculate the vegetation interception sedimentation flux that reflects the stem and leaf interception effect.

[0250] This step can be used to deduce the mass of sediment that is stripped from the water and forcibly settled by vegetation stems and leaves per unit time and per unit area.

[0251] Among them, the total amount of sediment particles intercepted by vegetation communities, as porous media, is positively correlated with the geometric scale of the cross-section and the intensity of water transport. The interception effect can be quantified by performing the following multiplication operation:

[0252] F _trap =η _trap ×N _s ×d _s ×min(h _v ,h)×|u|×S;

[0253] In the formula, F _trap η represents the vegetation interception of sedimentation flux. _trap N represents the efficiency of direct interception._s Indicates the density of vegetation stems, d _s Min(h) represents the stem diameter of vegetation. _v ,h) represents the effective interception depth of the minimum value between the extracted stem height and the actual water depth, |u| represents the absolute value of the flow velocity, and S represents the suspended sediment content.

[0254] Step 603: Introduce the vegetation interception sedimentation flux as an additional sink term into the suspended sediment transport equation to update the sediment transport pattern, and superimpose the sediment mass corresponding to the vegetation interception sedimentation flux into the topographic erosion and deposition equation to calculate and output the updated topographic features.

[0255] Within the framework of non-equilibrium sediment transport, the spatiotemporal variation of suspended sediment concentration is often controlled by advection diffusion and bed scouring and deposition fluxes. A negative sedimentation flux factor can be added to the right-hand side of the two-dimensional suspended sediment transport equation to subtract the mass of suspended particles trapped by the stem-leaf network.

[0256] Correspondingly, the intercepted material settles on the surface of the substrate after leaving the water body.

[0257] Optionally, the vegetation interception sedimentation flux is divided by the dry density of sediment to convert it into volumetric flux, and then added to the sedimentation term of the bottom topographic erosion and deposition equation.

[0258] After a discrete integral solution over one hydrodynamic time step, a matrix parameter containing the seabed elevation deformation is generated, and the updated terrain features are output.

[0259] In actual execution, the direct interception efficiency involved in the above steps is obtained through the following process:

[0260] For example, the pre-configured sediment particle feature scale is obtained, and the vegetation stem diameter of the current grid cell is extracted;

[0261] Read the median or average particle size of the sediment in the subsoil pre-mixed by geological survey and set it as the benchmark parameter characterizing the fineness of suspended matter;

[0262] Since the estuary area is mostly composed of fine-grained cohesive sediment, its movement trajectory in the water is easily affected by the boundary layer viscosity effect of the rigid cylindrical surface of plants. The extraction of this parameter can provide a basis for subsequent calculations.

[0263] For example, the direct interception efficiency is calculated based on the square of the ratio of the characteristic scale of sediment particles to the diameter of vegetation stems, in order to characterize the interception effect of vegetation stems and leaves on fine-particle sticky sediment.

[0264] The magnitude of the interception efficiency depends on the collision probability when suspended particulate matter penetrates the gaps in plant stems.

[0265] Optionally, the extracted sediment particle size parameters and plant geometric morphology parameters are divided to achieve dimensionless transformation, and a nonlinear amplification strategy is adopted. The specific calculation formula is as follows:

[0266] η _trap =C _trap ×(d _p / d _s ) 2 ;

[0267] In the formula, η _trap Indicates the efficiency of direct interception, C _trap d represents the empirical interception coefficient that comprehensively considers both particle surface adhesion and biofilm adhesion effects. _p d represents the characteristic scale of sediment particles. _s Indicates the diameter of the plant stem.

[0268] The above-mentioned empirical interception coefficient can be determined by technicians through indoor flume experiments on a vegetation model group, measuring the difference in suspended sediment concentration between upstream and downstream, back-calculating the interception efficiency, and then fitting the results.

[0269] The above-mentioned method of incorporating stem and leaf interception sedimentation flux into the sediment transport evolution supplements the calculation method of fine-particle cohesive sediment deposition in porous vegetation media, which is helpful for simulating the positive feedback process of topographic sedimentation.

[0270] In one possible implementation, a graphics processing unit (GPU) is used to accelerate the solution process of the mutual feedback closed loop in parallel. This parallel acceleration specifically includes:

[0271] Step 701: Obtain the vegetation activity of each grid cell at the current mutual feedback evolution moment;

[0272] Optionally, the current stock data and evolution rate of the target grid are read via the memory bus to assess the degree of change in the vegetation life cycle state of the area. Due to the spatial heterogeneity of habitat conditions in different grids, this comprehensive index is quantified separately for each discrete computing node.

[0273] The vegetation activity level is calculated as follows:

[0274] Step 7011: Calculate the biomass saturation based on the current vegetation biomass characteristics and environmental carrying capacity of the grid cells;

[0275] Correspondingly, the current vegetation biomass characteristics and the environmental carrying capacity of the corresponding grid can be extracted, and a division operation can be performed. This ratio reflects the relative crowding degree of the current vegetation stock from the upper limit of the spatial ecological capacity.

[0276] Step 7012: Calculate the relative rate of change of biomass based on the changes in the current vegetation biomass characteristics in adjacent time steps;

[0277] Optionally, the biomass data of the previous vegetation update cycle can be retrieved from the cache, the absolute difference can be calculated by combining it with the current biomass data, and then normalized by dividing it by the product of the time step and the maximum possible growth increment.

[0278] Accordingly, the absolute value of the obtained rate of change is taken to account for both the rapid expansion of vegetation and the large-scale erosion and stripping.

[0279] Step 7013: Compare the biomass saturation with the relative rate of biomass change, and select the maximum value of the two as the vegetation activity.

[0280] Accordingly, the two extracted normalized indicators are compared, and the one with the larger value is selected as the final evaluation parameter. The calculation formula is as follows:

[0281] I _veg =max(B / K,|ΔB / Δt _veg | / (r _max ×K));

[0282] In the formula, I _veg B represents vegetation activity, K represents environmental carrying capacity, B / K represents biomass saturation, ΔB represents biomass change between adjacent time steps, and Δt represents biomass change between adjacent time steps. _veg Indicates the vegetation regeneration time step, r _max Represents the maximum growth rate of vegetation, |ΔB / Δt _veg | / (r _max ×K) represents the relative rate of change of biomass after extracting the absolute value.

[0283] Suppose that in a certain computational iteration step, the biomass saturation of the predetermined grid cells is calculated to be 0.8 after normalization, indicating that the vegetation in this area is close to the spatial carrying capacity limit;

[0284] Meanwhile, because this grid is in a period of succession stagnation, its relative rate of biomass change is calculated to be 0.1. After comparison and determination, 0.8 was selected as the vegetation activity level of this grid unit.

[0285] Step 702: Construct a spatially adaptive morphological acceleration field based on vegetation activity, and dynamically adjust the morphological acceleration factor of each grid unit according to the spatially adaptive morphological acceleration field.

[0286] After obtaining the activity scalar of the global grid, a morphological acceleration distribution matrix with the same dimension as the spatial computing grid is constructed.

[0287] Among them, the morphological acceleration factor decreased as vegetation activity increased.

[0288] Optionally, a linear decay feedback equation controlled by the activity index is introduced, and the calculation formula is as follows:

[0289] f _mor =f _mor,max ×(1-κ _veg ×I _veg );

[0290] In the formula, f _mor f represents the dynamically adjusted morphological acceleration factor. _mor,max κ represents the preset maximum morphological acceleration factor for the bare beach area. _veg I represents the vegetation constraint coefficient. _veg Indicates vegetation activity.

[0291] In summary, when vegetation is in a period of growth or degradation, the activity index is close to the upper limit, and the morphological acceleration factor is suppressed to a smaller baseline, thereby improving the temporal resolution to restore the biological evolution process.

[0292] When in a bare beach area without vegetation or during a period of evolutionary stagnation, the activity index tends to the lower limit, and the morphological acceleration factor recovers to a larger baseline, which is used to maintain the propulsion efficiency of bottom hydrodynamic calculations.

[0293] Step 703: Use the adjusted morphological acceleration factor to control the time step of geomorphic evolution in the three-dimensional dynamic geomorphic numerical model in order to achieve long-period dynamic simulation calculation.

[0294] In the process of performing multi-threaded concurrent computations in the computing architecture, a layered time-stepping strategy can be adopted to synchronize the underlying non-constant hydrodynamic calculations and vegetation ecological updates.

[0295] For example, hydrodynamics and sediment transport operate in smaller time steps, while vegetation ecological regeneration operates in daily time steps.

[0296] Furthermore, the morphological acceleration factor of each grid is extracted to locally amplify the bed scouring and sedimentation output from the single-step hydrodynamic calculation.

[0297] In one alternative implementation, when processing sediment flux exchange between adjacent grid cells, the harmonic mean of the morphological acceleration factor of two adjacent cells is extracted and weighted for allocation.

[0298] By controlling the frequency of data exchange between the central processing unit and the graphics processing unit, and performing a write-back of the video memory containing vegetation biomass parameters after each preset number of hydrodynamic bottom-level cycles, the system outputs evolutionary simulation results that support long-term timescales while ensuring the authenticity of ecological mutation periods.

[0299] In another optional implementation, the value of the vegetation constraint coefficient can be determined by comparing the simulation stability of the globally unified acceleration scheme and the adaptive acceleration scheme in the verification condition, so as to ensure that the morphological acceleration factor is not lower than the preset minimum safety limit when the vegetation activity reaches its maximum.

[0300] In this method, a spatially adaptive morphological acceleration field is constructed based on the local evolutionary activity of vegetation. This mechanism can maintain the efficiency of global long-term geomorphological evolution and adaptively reduce the time step in areas of dramatic ecological change, thus promoting the high-frequency capture of ecological mutation processes at cross-decadal simulation scales and suppressing numerical divergence in long-term simulations.

[0301] In some scenarios, such as for complex estuarine beds with mixed distributions of multiple particle sizes, parallel calculation branches for sediment gradation and mixing layer adjustment can be added when updating topographic features in the calculation output.

[0302] Optionally, a mixing layer region is defined in the vertical direction of the bed surface, and the thickness of this region fluctuates dynamically with the hydrodynamic intensity.

[0303] When the topographic scouring and deposition equation determines that the current state is continuous scouring, the lower interface of the mixed layer is controlled to cut downwards towards the subsoil, and the sediment particle size distribution data of the deep original subsoil is continuously read and replenished into the mixed layer.

[0304] It should be understood that this method can avoid the calculation divergence problem caused by the erosion of surface sediment in long-cycle scour simulations.

[0305] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A method for simulating the interaction between topographic erosion and sedimentation and habitat succession, characterized in that, include: Acquire the initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass characteristics of the target water area; The initial vegetation biomass characteristics are transformed into hydrodynamic resistance terms. The hydrodynamic resistance terms, along with the initial topographic features and initial hydrodynamic conditions, are input into a pre-constructed three-dimensional dynamic geomorphological numerical model to calculate the updated topographic features. The updated terrain features are input into a pre-constructed staged vegetation growth dynamics model to calculate the updated vegetation biomass features. The updated vegetation biomass characteristics are fed back in real time to the three-dimensional dynamic geomorphological numerical model to update the hydrodynamic resistance term, driving the evolution calculation of the next time step, forming a two-way feedback closed loop of topographic erosion and sedimentation and habitat succession, and outputting the dynamic simulation results of topography and habitat. In the process of parallel acceleration of solving the mutual feedback closed loop, the specific steps include: obtaining the vegetation activity of each grid cell at the current mutual feedback evolution moment; A spatially adaptive morphological acceleration field is constructed based on vegetation activity, and the morphological acceleration factor of each grid unit is dynamically adjusted according to the spatially adaptive morphological acceleration field. The morphological acceleration factor decreases as vegetation activity increases; the adjusted morphological acceleration factor is used to control the time step of geomorphological evolution in the three-dimensional dynamic geomorphological numerical model in order to achieve long-term dynamic simulation calculation. The vegetation activity level is calculated as follows: Biomass saturation is calculated based on the current vegetation biomass characteristics and environmental carrying capacity of the grid cells; the relative rate of biomass change is calculated based on the changes in current vegetation biomass characteristics over adjacent time steps; the biomass saturation level and the relative rate of biomass change are compared, and the maximum value between the two is selected as the vegetation activity level I. _veg ;I _veg =max(B / K,|ΔB / Δt _veg | / (r _max ×K); where B represents the current vegetation biomass characteristics, K represents the environmental carrying capacity, B / K represents the biomass saturation, ΔB represents the biomass change between adjacent time steps, and Δt _veg Indicates the vegetation regeneration time step, r _max Represents the maximum growth rate of vegetation, |ΔB / Δt _veg | / (r _max ×K) represents the relative rate of change of biomass after extracting the absolute value.

2. The method according to claim 1, characterized in that, By inputting updated terrain features into a pre-constructed phased vegetation growth dynamics model, the updated vegetation biomass characteristics are calculated. The phased vegetation growth dynamics model includes a planting stage model, a root development stage model, and a dynamic growth stage model.

3. The method according to claim 2, characterized in that, The colonization phase model is specifically used to perform the following steps: Extract the duration of no flooding in the current grid cell, compare the duration of no flooding with the preset planting opportunity time threshold, and determine the planting opportunity window; Under the condition of satisfying the planting opportunity window, the planting probability of the current grid cell is determined based on the pre-configured correspondence between flooding frequency and planting probability. The planting result of vegetation is determined based on the planting probability. When the planting is successful, the initial seedling density is output as the updated vegetation biomass feature of the grid cell.

4. The method according to claim 2, characterized in that, The root development stage model is specifically used to perform the following steps: The growth time of vegetation was obtained, and the critical impact shear stress of seedlings was set to increase linearly with the growth time. Obtain the actual bed shear stress of the current mesh element; The actual bed shear stress is compared with the critical impact shear stress. The survival status or scour and peeling status of the seedlings is determined based on the comparison results. If the seedlings are determined to be in a survival status, their growth is maintained.

5. The method according to claim 2, characterized in that, The dynamic growth phase model is specifically used to perform the following steps: Construct a logistic growth term based on the maximum growth rate of vegetation; Construct a natural mortality item based on a preset natural mortality rate; The inundation loss term is constructed based on the difference between the maximum tolerable water depth of vegetation and the current actual water depth. The hydrodynamic scour loss term is constructed based on the difference between the maximum bed shear stress that the vegetation can withstand and the current actual bed shear stress. Based on logistic growth, natural mortality, inundation loss, and hydrodynamic scour loss terms, a biomass spatiotemporal variation control equation is constructed. Combined with the environmental carrying capacity constrained by both elevation and sediment thickness, the updated vegetation biomass characteristics are calculated and output.

6. The method according to claim 5, characterized in that, Environmental carrying capacity is calculated in the following way: Obtain the pre-configured ideal environmental carrying capacity, the actual elevation of the current grid cell, and the sediment thickness; Based on the relative relationship between the actual elevation and the preset optimal vegetation growth elevation, the Gaussian distribution response term is calculated. Based on the ratio of sediment thickness to a preset critical sediment thickness, the nonlinear saturated response term is calculated. The environmental carrying capacity is calculated by multiplying the ideal environmental carrying capacity, the Gaussian distributed response term, and the nonlinear saturated response term.

7. The method according to claim 1, characterized in that, The specific process of converting initial vegetation biomass characteristics or regenerating vegetation biomass characteristics into hydrodynamic resistance terms includes: Obtain the pre-configured allometric growth mapping relationship and the basic drag coefficient of vegetation; Through the allometric growth mapping relationship, the initial vegetation biomass characteristics or updated vegetation biomass characteristics, which are scalar data, are nonlinearly mapped to the three-dimensional morphological parameters of vegetation. The flexible reduction coefficient is calculated based on three-dimensional morphological parameters and flow field parameters using a pre-constructed flexible reduction function. The effective drag coefficient is obtained by multiplying the basic drag coefficient by the flexibility reduction coefficient, and then used as the hydrodynamic resistance term.

8. The method according to claim 1, characterized in that, The three-dimensional dynamic geomorphological numerical model contains suspended sediment transport equations for calculating sediment transport patterns and topographic erosion and deposition equations for calculating topographic evolution. The process of calculating updated topographic features by inputting a pre-constructed three-dimensional dynamic geomorphological numerical model specifically includes the following steps: To obtain the direct interception efficiency, vegetation stem density, vegetation stem diameter, effective interception depth, flow velocity, and suspended sediment content of the target water area; Based on the multiplicative relationship between direct interception efficiency, vegetation stem density, vegetation stem diameter, effective interception water depth, flow velocity, and suspended sediment content, the vegetation interception sedimentation flux reflecting the stem and leaf interception effect is calculated. The vegetation interception sedimentation flux is introduced as an additional sink term into the suspended sediment transport equation to update the sediment transport pattern. The sediment mass corresponding to the vegetation interception sedimentation flux is superimposed into the topographic erosion and deposition equation, and the updated topographic features are calculated and output.

9. The method according to claim 1, characterized in that, Before obtaining the initial topographic features, initial hydrodynamic conditions, and initial vegetation biomass characteristics of the target water area, the offline construction step is also included to obtain the basic drag coefficients required for the three-dimensional dynamic geomorphological numerical model: Different levels of wave conditions and unidirectional flow conditions were combined and set in a pre-constructed wind and wave channel to determine the water flow resistance pattern of the vegetation area under different layout forms. Extract Reynolds number and eddy current oscillation force characteristic parameters that are correlated with the water flow resistance law to a preset standard; An empirical formula is established based on the Reynolds number and the characteristic parameters of eddy current oscillation force to obtain the pre-constructed basic drag force coefficient.

Images