A simulation calculation method for the coupled evolution of planar two-dimensional coastal aeolian landforms and vegetation growth.

Abstract

The present invention relates to a simulation calculation method for the coupled evolution of two-dimensional coastal aeolian landforms and vegetation growth, solves the control equations such as aeolian convection and the growth and extinction of vegetation communities in the coastal direction, considers the influence of oblique driving force and uneven coastal terrain on the response of the ecological geomorphic system and the evolution of the plane pattern, and realizes the simulation of the vegetation growth process and the landform evolution under the mutual restriction of the backshore aeolian sand movement within a certain period of time.

Description

Technical Field

[0001] The invention relates to a simulation calculation method for the coupled evolution of a two-dimensional coastal sand landform and vegetation growth, and belongs to the field of coastal sand dune protection and restoration and sand control. Background Art

[0002] With global climate change and rising sea levels, the risk of overtopping and flooding caused by storm surges has increased, leading to a sudden increase in the pressure on coastal disaster prevention and reduction. Coastal sand dunes are the last line of defense against coastal floods and play an important role in protecting the lives and property of people in coastal areas. However, the continuous impact of sea breezes and waves will continue to impact the sand dunes, causing changes in the shape and structure of the sand dunes, and even causing the collapse and loss of the sand dunes. Therefore, the protection and restoration of sandy coastal dunes has become one of the most important construction contents of my country's coastal protection and restoration projects.

[0003] Coastal aeolian sand movement refers to the movement of near-bottom sediment mainly in a saltation manner driven by wind, and is the main process of sediment exchange in coastal dunes. Under different dynamic conditions, aeolian sand movement causes rapid (or slow) changes in coastal dune topography. Although during storm surges, coastal sediment movement is dominated by water and sand movement, during normal waves (or storm recovery intervals), aeolian sand movement plays a dominant role, which is fully reflected on Yuanzhao Beach in Pingtan, Fujian Province and Guanyinshan Beach in Xiamen, my country. Accurate simulation and calculation of the evolution of backshore-dune landforms can provide an effective method for the protection and restoration of coastal dunes. However, the movement of wind-blown sand on the coast is complex. On a temporal scale, due to the high randomness of wind power, changes are difficult to accurately simulate and quantify, resulting in a high degree of dynamics in the movement of wind-blown sand. On a spatial scale, due to the complex dune terrain and rich geomorphic elements on the actual coast, and the influence of vegetation growth and distribution, the growth of vegetation is restricted by many factors, such as soil moisture, salinity, wind conditions, etc., and the changes in its distribution and density have an important impact on the evolution of the landform, but it is very difficult to accurately simulate the growth and changes of vegetation. The quantitative relationship between vegetation and wind and sand is not clear enough, and the coupling between sand and wind field and vegetation makes the calculation of the evolution of beach backshore-dune landform more cumbersome.

[0004] At present, the bottleneck of the simulation of the coupled evolution of beach sand landforms and vegetation growth is that under the action of two-dimensional oblique wind power, the simulation of vegetation growth and the interaction of coastal sand landforms has temporal and spatial complexity. For example, during the growth cycle, the wind-cutting and sand-fixing effect of vegetation affects beach erosion and accumulation, and the change of beach elevation will react on the growth of vegetation. This process will also be affected by the humidity of the beach and the distribution of sand sources in time and space. Therefore, the simulation of the evolution of backshore-sand dunes and vegetation is a difficult point, and a new simulation calculation method is urgently needed to address this simulation bottleneck. Summary of the invention

[0005] The present invention provides a two-dimensional coastal sand landform and vegetation growth coupled evolution simulation calculation method, which can simulate the vegetation growth process and the landform evolution under the mutual restriction of the backshore sand movement within a certain period of time.

[0006] The technical solution adopted by the present invention to solve its technical problem is:

[0007] A simulation calculation method for the coupled evolution of two-dimensional coastal sand landform and vegetation growth includes the following steps:

[0008] Step S1: Obtain key information, including coastal topography data, representative coastal sediment particle size, vegetation growth cycle, tide level, and boundary wind speed;

[0009] Step S2: The wind speed at each location is obtained through the boundary wind speed and the coastal terrain data. The wind speed at each location includes the established three-dimensional coordinate system, in which the offshore direction is set as the x-axis direction, the shoreward direction is set as the positive direction, the along-shore direction is set as the y-axis direction, the downstream direction is set as the positive direction, and the ground elevation direction in the coastal terrain data is set as the z-axis direction, and the vertically upward direction is set as the positive direction;

[0010] Step S3: simulating the evolution of backshore vegetation based on the vegetation growth cycle, generalizing the growth elements of the vegetation growth cycle into two parameters, namely, vegetation growth height and vegetation coverage;

[0011] Step S4: Based on the tide level and representative sediment particle size information, the beach surface humidity of the shoreline is calculated, and the beach surface moisture content is further calculated;

[0012] Step S5: Based on the vegetation growth height and vegetation coverage obtained in step S3, the beach surface moisture content obtained in step S4, and the wind speed at each location, the equilibrium sediment transport rate at each location on the backshore beach is calculated;

[0013] Step S6: Calculate the actual sediment transport rate at each location on the backshore beach surface;

[0014] Step S7: Based on the actual sediment transport rate obtained in step S6, calculate the actual sediment transport rate in the x direction and the y direction;

[0015] Step S8: Calculate the change of the elevation of the backshore within a time step, and then perform iterative loop calculation from step S1 until the simulation time specified by the user is completed;

[0016] Furthermore, in step S2, it is assumed that the wind at the boundary is at an angle The wind speed components in the x and y directions at a distance z from the ground and the total wind speed are:

[0017]

[0018] In formula (1), is the angle between the incident wind and the positive direction of the x-axis, is the surface shear stress disturbance value in the x direction, is the surface shear stress disturbance value in the y direction. Specifically,

[0019]

[0020] In formulas (2) and (3), is the partial derivative of ground elevation in the x direction, is the partial derivative of the ground elevation in the y direction, α and β are both functions of L / z′, L is the characteristic length of the ground elevation change, z′ is the effective roughness length, α and β are set to 3 and 0.2 on the windward slope of the coast, and α and β are set to 2 and 1 on the leeward slope of the coast;

[0021] Furthermore, in step S3, the control equations of the vegetation growth height and the vegetation coverage are converted into linear equations by the finite element method for subsequent solution;

[0022] Specifically, the control equations for vegetation growth height and vegetation coverage are:

[0023]

[0024] In formula (4), h v is the vegetation growth height, H v,max is the maximum growth height of vegetation, C v is the vegetation coverage, C v,max is the maximum vegetation coverage, T g For the growth cycle;

[0025] Among them, the vegetation coverage rate C v At the same time, affected by the accumulation of wind and sand, the influence formula is:

[0026]

[0027] In formula (5), is the vegetation coverage rate at the previous moment, H v,max is the maximum growth height of vegetation, z is the ground elevation at that moment, and z pre is the ground elevation at the previous moment;

[0028] Multiplying formula (4) with the weight function Φ and integrating it over the calculation area, we get

[0029]

[0030] In formula (6), is the vegetation growth height at the previous moment, and the weight function Φ is:

[0031]

[0032] make Deriving the unit equation

[0033]

[0034] Based on the unit matrix multiplication form of the linear equation system formula (8), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100;

[0035] Furthermore, in step S4, calculating the water content of the shoreline beach surface includes the following steps:

[0036] Step S41, obtaining the height of the capillary water action on the beach, the calculation formula is:

[0037]

[0038] In formula (9), h c is the height of the capillary water action on the beach; γ is the surface tension of water, which is 0.0728N / m; θ is the angle between the pore water and the pore wall, which is 0; ρ w is the water density; r m is the average pore radius, taking the value d m / 5;d m is the average radius of the representative particle size of sediment;

[0039] Step S42, assuming that the water content at the shoreline is 25%:

[0040]

[0041] In formula (10), w 20mm is the water content of the beach surface layer at 20 mm, h c is the height of capillary water action on the beach, h is the elevation of the groundwater level on the beach, which matches the tide level;

[0042] Furthermore, in step S5, based on the vegetation growth height and vegetation coverage obtained in step S3, the beach surface moisture content obtained in step S4, and the wind speed at each location, the equilibrium sediment transport rate at each location on the backshore beach is calculated.

[0043] The calculation formula of equilibrium sediment transport rate is:

[0044]

[0045] In formula (11), Q sat is the equilibrium sediment transport rate; α p is a constant used to convert the height z above the ground bThe wind speed is converted into the near-bottom friction velocity u * ; C is the parameter representing the particle size distribution width of sediment; ρ a is the air density, in kg / m 3 ; g is the acceleration due to gravity, in m / s 2 ;d n is the representative particle size of sediment, in m; D n is the reference particle size, in m; u zb is the height above ground z b The wind speed, u tb is the height above ground z b Critical starting wind speed;

[0046] α p According to the Prandtl-VonKarman law, the calculation formula is:

[0047]

[0048] In formula (12), is the Karman constant, z b is the height above the ground at that moment, z' b is the effective roughness length;

[0049] u tb The calculation formula is:

[0050]

[0051] In formula (13), is the critical friction velocity near the bottom;

[0052] u zb The calculation formula is:

[0053] u zb =u * / α P (14)

[0054] In formula (14), u * is the friction velocity near the bottom;

[0055] Furthermore, in formula (13), the critical friction velocity near the bottom is The impacts include:

[0056]

[0057] In formula (15), A is the empirical coefficient; g is the acceleration due to gravity, in m / s 2 ρ P is the density of sediment particles, in kg / m 3 ρa is the air density, in kg / m 3 ;d n is the representative particle size of sediment, in m;

[0058] Critical friction velocity near the bottom The influence of the beach surface water content is also included, namely:

[0059] u * th,moist =α w u * (16)

[0060] In formula (16), is the critical friction velocity near the bottom under the set beach water content, u * is the critical friction velocity near the bottom when the beach water content is low, a w is the conversion coefficient, α w =1+0.1(d 50,ref / d 50 ) 20mm , d 50,ref is the reference particle size, w 20mm is the water content of the beach surface layer at 20 mm, d 50 is the median particle size of sediment particles;

[0061] Furthermore, in formula (14), the effect of vegetation on the ground elevation z b Wind speed u zb The impact is the vegetation growth height or vegetation coverage.

[0062] If affected by the height of vegetation growth,

[0063] Effect of vegetation growth height on the friction velocity u near the bottom * The impact is:

[0064]

[0065] In formula (17), ρ v is the vegetation factor, ρ v =(h v / H v,max ) 2 ,h v is the vegetation growth height, H v,max is the maximum growth height of vegetation; u *,0 is the friction velocity near the bottom without the influence of vegetation; Γ is the parameter representing the resistance of plant morphology, with a value of 4;

[0066] If affected by vegetation coverage,

[0067] The effect of vegetation coverage on wind speed at 10m from the surface is:

[0068] u 10,veg =u 10 (1-kC ν ) (18)

[0069] In formula (18), u 10,veg is the wind speed at 10 m above the surface; k is a parameter that characterizes the effect of vegetation coverage on wind speed at 10 m above the surface, and takes different values ​​for different plant species. For typical small upright or flat herbaceous dune plants, k = 1.8, and for small round stemless plants, k = 4.6; C v is the vegetation coverage rate;

[0070] Furthermore, in step S6, the actual sediment transport rate control equation at each location on the backshore beach is:

[0071]

[0072] In formula (19), C c is the sediment concentration per unit area, in kg / m 2 ξ is the conversion coefficient between wind speed and sediment movement speed, which is 0.8; T is the time required for erosion or accumulation to complete, which is 1s; C u is the saturated sediment concentration, in kg / m 2 ; p is the reduction coefficient of the source term when the sand source is insufficient;

[0073] Among them, C u =Q sat / u z ; The reduction factor p is S ε is the amount of erodible sediment, in kg / m 2 ;

[0074] Solve formula (19) and combine formula (19) with Multiply and integrate over the computational area, and we get

[0075]

[0076] make The time term is derived by backward difference to derive the unit equation

[0077]

[0078] In formula (20), τ s is the stabilization parameter, u is the velocity vector, is the gradient of the interpolation basis function;

[0079]

[0080] In formula (22), Δt is the time step, Δl is the grid length, and the average value is taken here; u x is the wind speed component in the x direction; u y is the wind speed component in the y direction;

[0081] Based on the unit matrix multiplication form of the linear equation system formula (21), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100;

[0082] Further, based on the actual sediment transport rate obtained in step S6, the actual sediment transport rates in the x direction and the y direction are obtained, that is,

[0083]

[0084] In formula (23), Q b,x is the actual sediment transport rate in the x direction, Q b,y is the actual sediment transport rate in the y direction; the current coastal elevation change is calculated using the mass conservation method, that is:

[0085]

[0086] In formula (24), ρ s is the sediment density, λ is the porosity, and z is the ground elevation;

[0087] Similarly, the mass conservation equation for simulating landform evolution is obtained as follows:

[0088]

[0089] Based on the unit matrix multiplication form of the linear equation system formula (25), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100.

[0090] Through the above technical solution, compared with the prior art, the present invention has the following beneficial effects:

[0091] The simulation calculation method for the coupled evolution of two-dimensional coastal aeolian landforms and vegetation growth provided by the present invention fully considers the ambiguity of the interaction mechanism between vegetation and wind and sand, and realizes the simulation of the vegetation growth process and the landform evolution under the mutual constraint of the backshore aeolian sand movement within a certain period of time. BRIEF DESCRIPTION OF THE DRAWINGS

[0092] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0093] Figure 1It is a schematic flow chart of a simulation calculation method for the coupled evolution of two-dimensional coastal wind-sand landform and vegetation growth provided by the present invention;

[0094] Figure 2 The first embodiment provided by the present invention is about the changes in vegetation height, wind speed and sand transport rate when the simulated terrain is a two-dimensional ideal sand dune;

[0095] Figure 3 The first embodiment provided by the present invention is about simulating the change and deflection of wind speed on a two-dimensional plane during the growth of vegetation when the terrain is a two-dimensional ideal sand dune;

[0096] Figure 4 The first embodiment provided by the present invention is about simulating the change of the sediment transport rate on a two-dimensional plane during the vegetation growth process when the terrain is a two-dimensional ideal sand dune;

[0097] Figure 5 is a schematic diagram of a simulated beach profile, tide level and incident conditions according to a second embodiment of the present invention;

[0098] Figure 6 The second embodiment of the present invention provides a second embodiment of the present invention, which is about the changes in vegetation height, wind speed, sediment transport rate and beach profile at the initial and final moments of the simulation of a selected representative profile without considering the effect of vegetation on wind speed attenuation;

[0099] Figure 7 The second embodiment provided by the present invention is about the changes in vegetation height, wind speed, sediment transport rate and beach profile at the initial and final moments of the simulation of a selected representative profile taking into account the effect of vegetation on wind speed attenuation. DETAILED DESCRIPTION

[0100] The present invention will now be described in further detail with reference to the accompanying drawings. In the description of the present application, it should be understood that the orientation or positional relationship indicated by the terms "left side", "right side", "upper part", "lower part", etc. is based on the orientation or positional relationship shown in the accompanying drawings, and is only for the convenience of describing the present invention and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and "first", "second", etc. do not indicate the importance of the components, and therefore cannot be understood as a limitation on the present invention. The specific dimensions used in this embodiment are only for illustrating the technical solution by example, and do not limit the scope of protection of the present invention.

[0101] As described in the background technology, coastal dunes are the last line of defense against coastal floods and play an important role in protecting the lives and property of people in coastal areas. Accurate simulation and calculation of the evolution of shore-dune landforms can provide an effective method for the protection and restoration of coastal dunes. At present, most of the models provided for simulation calculations make some simple assumptions to simplify the design process, resulting in deviations between the calculation results and the actual situation; and many parameters need to be input when the model is formed, and the measurement and determination of these parameters have certain errors and uncertainties.

[0102] Therefore, in order to fully consider the complexity of coastal wind and sand movement, the present application provides a two-dimensional simulation and calculation method for the coupled evolution of coastal wind and sand landforms and vegetation growth, solves the control equations such as wind and sand convection and the growth and extinction of vegetation communities in the coastal direction, considers the influence of oblique driving force and uneven coastal terrain on the response of the ecological geomorphological system and the evolution of the plane pattern, and realizes the simulation of the coupled evolution of beach-backshore ecological geomorphology.

[0103] Figure 1 As shown, it is a flow chart of the entire simulation calculation method provided by this application, which specifically includes the following steps:

[0104] Step S1: Read terrain, boundary, grid and keyword files to obtain key information such as coastal terrain data, representative coastal sediment particle size, vegetation growth cycle, tide level and boundary wind speed.

[0105] Step S2: The wind speed at each location is obtained through the boundary wind speed and the coastal terrain data. In each location, the three-dimensional coordinate system is established, and the offshore direction is set as the x-axis direction (the shoreward direction is the positive direction), the coastal direction is set as the y-axis direction (the downstream direction is the positive direction), and the ground elevation direction in the coastal terrain data is set as the z-axis direction, and the vertically upward direction is the positive direction;

[0106] Assume the wind at the boundary is at an angle The wind speed components in the x and y directions at a distance z from the ground and the total wind speed are:

[0107]

[0108]

[0109] In formula (1), is the angle between the incident wind and the positive direction of the x-axis, is the surface shear stress disturbance value in the x direction, is the surface shear stress disturbance value in the y direction. Specifically, the formula that characterizes the distribution of wind shear stress on the terrain profile is:

[0110]

[0111] In formulas (2) and (3), is the partial derivative of ground elevation in the x direction, is the partial derivative of the ground elevation in the y direction, α and β are both functions of L / z′, L is the characteristic length of the ground elevation change, z′ is the effective roughness length, and α and β are set to 3 and 0.2 on the windward slope of the coast, and α and β are set to 2 and 1 on the leeward slope of the coast.

[0112] Step S3: simulating the evolution of backshore vegetation based on the vegetation growth cycle, generalizing the growth elements of the vegetation growth cycle into two parameters, namely, vegetation growth height and vegetation coverage;

[0113] Among them, the control equations of vegetation growth height and vegetation coverage are:

[0114]

[0115] In formula (4), h v is the vegetation growth height, H v,max is the maximum growth height of vegetation, C v is the vegetation coverage, C v,max is the maximum vegetation coverage, T g For the growth cycle;

[0116] In formula (4), the maximum growth height of vegetation H v,max and the maximum vegetation coverage C v,max It is assumed that it is related to the height difference of the groundwater surface line, that is, it is believed that the rising sea water due to capillary action has an adverse effect on the growth of vegetation. The details are shown in Table 1:

[0117] Table 1 Relationship between the maximum height and maximum coverage of backshore vegetation and the height difference of groundwater surface

[0118]

[0119] According to Table 1, the maximum growth height H of vegetation is determined by interpolation. v,max and the maximum vegetation coverage C v,max .

[0120] In addition, vegetation growth is also affected by wind and sand movement. In the backshore vegetation evolution numerical module, the vegetation growth height is related to the erosion or accumulation height at that location, that is, the height of the vegetation increases / decreases with the increase of erosion / accumulation height. However, when the erosion exceeds a certain level, the vegetation is considered dead and the growth height is set to 0. The maximum erosion height of vegetation death set in the module is 0.2m. Vegetation coverage C v At the same time, affected by the accumulation of wind and sand, the influence formula is:

[0121]

[0122] In formula (5), is the vegetation coverage rate at the previous moment, H v,max is the maximum growth height of vegetation, z is the ground elevation at that moment, and z pre is the ground elevation at the previous moment.

[0123] Next, the control equations of vegetation growth height and vegetation coverage (Formula (4)) are transformed into linear equations through the finite element method for subsequent solution; in order to better achieve the fit of the simulation area boundary, the computational grid division of this application adopts triangular irregular grid, the finite element method is used in the numerical format, and the time term is discretized by backward difference; specifically, the linear interpolation basis function is used, and according to the Galerkin finite element method, the weight function is equal to the interpolation function.

[0124] Multiplying formula (4) with the weight function Φ and integrating it over the calculation area, we get

[0125]

[0126] In formula (6), is the vegetation growth height at the previous moment, is the vegetation coverage rate at the previous moment. The weight function Φ is:

[0127]

[0128] make Deriving the unit equation

[0129]

[0130] Based on the unit matrix multiplication form of the linear equation system formula (8), the overall equation is assembled and then solved by the generalized minimum residual method (GMRES), with the maximum number of iterations set to 100.

[0131] Step S4: Based on the tide level and representative particle size information of sediment, the beach surface humidity of the shoreline is calculated, taking into account the influence of wave climbing and underground capillary water on beach surface humidity and wind-blown sand movement near the beach shoreline. Due to the capillary effect between beach sand particles, the groundwater depth of the beach can continue to rise for a certain height. The calculation formula for obtaining the height of the capillary water effect on the beach is:

[0132]

[0133] In formula (9), h c is the height of the capillary water action on the beach; γ is the surface tension of water, which is 0.0728N / m; θ is the angle between the pore water and the pore wall, which is 0; ρ w is the water density; r mis the average pore radius, taking the value d m / 5;d m is the average radius of the representative particle size of sediment;

[0134] h calculated in formula (9) c The starting position is the groundwater surface line. Assuming that the height of the groundwater surface line is the same as the static water level line and the height along the way is constant, the beach surface humidity changes linearly with the elevation. It is assumed that the moisture content at the shoreline is 25%, that is,

[0135] w 20mm =25%(h c -h),h=[0,h c ] (10)

[0136] In formula (10), w 20mm is the water content of the beach surface layer at 20 mm, h c is the height of capillary water action on the beach, and h is the elevation of the groundwater level on the beach, which matches the tide level and is expressed as the tide level.

[0137] The above step S3 is to calculate the vegetation growth height and vegetation coverage, and step S4 is to find the beach surface moisture content to quantify the critical friction velocity near the bottom. That is, step S5, combined with the wind speed at each position, calculates the equilibrium sediment transport rate at each position of the backshore beach surface.

[0138] The calculation formula of equilibrium sediment transport rate is:

[0139]

[0140] In formula (11), Q sat is the equilibrium sediment transport rate; α p is a constant used to convert the height z above the ground b The wind speed is converted into the near-bottom friction velocity u * ; C is the parameter representing the particle size distribution width of sediment; ρ a is the air density, in kg / m 3 ; g is the acceleration due to gravity, in m / s 2 ;d n is the representative particle size of sediment, in m; D n is the reference particle size, in m; u zb is the height above ground z b The wind speed, u tb is the height above ground z b Critical starting wind speed;

[0141] α p According to the Prandtl-VonKarman law, the calculation formula is:

[0142]

[0143] In formula (12), is the Karman constant, z b is the height above the ground at that moment, z' b is the effective roughness length;

[0144] u tb The calculation formula is:

[0145]

[0146] In formula (13), is the critical friction velocity near the bottom;

[0147] u zb The calculation formula is:

[0148] u zb =u * / α P (14)

[0149] In formula (14), u * is the friction velocity near the bottom.

[0150] Formula (13) and formula (14) represent the factors affecting the near-bottom critical friction velocity and the near-bottom friction velocity. For formula (13), the influence of the bottom critical friction velocity includes:

[0151]

[0152] In formula (15), A is the empirical coefficient; g is the acceleration due to gravity, in m / s 2 ρ P is the density of sediment particles, in kg / m 3 ρ a is the air density, in kg / m 3 ;d n is the representative particle size of sediment, in m;

[0153] Critical friction velocity near the bottom The influence of the beach surface water content is also included, namely:

[0154] u * th,moist =α w u * (16)

[0155] In formula (16), is the critical friction velocity near the bottom under the set beach water content, u * is the critical friction velocity near the bottom when the beach water content is low, aw is the conversion coefficient, α w =1+0.1(d 50,ref / d 50 ) 20mm , d 50,ref is the reference particle size, w 20mm is the water content of the beach surface layer at 20 mm, d 50 is the median particle size of sediment particles.

[0156] For formula (14), the vegetation effect on the ground elevation z b Wind speed u zb The impact is the vegetation growth height or vegetation coverage, which is reflected in the impact on the near-bottom friction flow velocity.

[0157] If affected by the height of vegetation growth,

[0158] Effect of vegetation growth height on the friction velocity u near the bottom * The impact is:

[0159]

[0160] In formula (17), ρ v is the vegetation factor, reflecting the ability of vegetation to cut wind and fix sand, ρ v =(h v / H v,max ) 2 ,h v is the vegetation growth height, H v,max is the maximum growth height of vegetation; u *,0 is the friction velocity near the bottom without the influence of vegetation; Γ is the parameter representing the resistance of plant morphology, with a value of 4;

[0161] If affected by vegetation coverage,

[0162] The effect of vegetation coverage on wind speed at 10m from the surface is:

[0163] u 10,veg =u 10 (1-kC ν ) (18)

[0164] In formula (18), u 10,veg is the wind speed at 10 m above the surface; k is a parameter that characterizes the effect of vegetation coverage on wind speed at 10 m above the surface, and takes different values ​​for different plant species. For typical small upright or flat herbaceous dune plants, k = 1.8, and for small round stemless plants, k = 4.6; C v is the vegetation coverage rate.

[0165] It should be noted that in actual calculations, the wind speed u under the influence of vegetation is zb The user can choose to use formula (17) or formula (18). The results calculated by the two formulas are the friction flow velocity near the bottom and the wind speed 10m above the ground. Therefore, the results are substituted into formula (14) to convert the wind speed at a specific height.

[0166] Step S6: Calculate the actual sediment transport rate at each location on the backshore beach surface;

[0167] The actual sediment transport rate control equation is:

[0168]

[0169] In formula (19), C c is the sediment concentration per unit area, in kg / m 2 ξ is the conversion coefficient between wind speed and sediment movement speed, which is 0.8; T is the time required for erosion or accumulation to complete, which is 1s; C u is the saturated sediment concentration, in kg / m 2 ; p is the reduction coefficient of the source term when the sand source is insufficient;

[0170] Among them, C u =Q sat / u z ; The reduction factor p is S ε is the amount of erodible sediment, in kg / m 2 Considering the insufficient effect of sediment source in erosion, we first assume that the maximum erosion depth d max , then the amount of eroded sediment S ε for:

[0171] Since the governing equation (19) contains only convection terms and is a strong convection equation, this type of problem is prone to numerical oscillations in traditional finite element methods. Therefore, the Streamline Upwind / Petrov-Galerkin (SUPG) method is used to solve the equation, which can effectively deal with problems dominated by convection and reduce possible oscillations in the numerical solution. The SUPG method enhances stability by adding an additional term to the interpolation basis function, which is aligned with the streamline direction. This means that the interpolation basis function Φ is replaced by τ s is the stabilization parameter, u is the velocity vector, is the gradient of the interpolation basis function; the stabilization term helps to reduce numerical oscillations caused by the dominance of convection.

[0172] Specifically, formula (19) and Multiply and integrate over the computational area, and we get

[0173]

[0174] make The time term is derived by backward difference to derive the unit equation

[0175]

[0176] By adjusting the interpolation basis function, the SUPG method can better capture the physical process dominated by convection, improve the stability and accuracy of the numerical solution, and improve the stability of the stabilization parameter τ. s The selection of is preliminarily determined using the following formula:

[0177]

[0178] In formula (22), Δt is the time step, Δl is the grid length, and the average value is taken here; u x is the wind speed component in the x direction; u y is the wind speed component in the y direction;

[0179] Based on the unit matrix multiplication form of the linear equation system formula (21), the overall equation is assembled and the generalized minimum residual method (GMRES) is used to solve the overall equation, and the maximum number of iterations is set to 100.

[0180] Then, in step S7, based on the actual sediment transport rate obtained in step S6, the actual sediment transport rate in the x direction and the y direction is obtained:

[0181]

[0182] In formula (23), Q b,x is the actual sediment transport rate in the x direction, O b,y is the actual sediment transport rate in the y direction; the current coastal elevation change is calculated using the mass conservation method, that is:

[0183]

[0184] In formula (24), ρ s is the sediment density, λ is the porosity, and z is the ground elevation;

[0185] Similarly, the mass conservation equation for simulating landform evolution is also assembled using the finite element method.

[0186]

[0187] Based on the unit matrix multiplication form of the linear equation system formula (25), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100.

[0188] Step S8: Calculate the change of the elevation of the coastal backshore within a time step, and then perform iterative loop calculation from step S1 until the simulation time specified by the user is completed.

[0189] In order to verify the superiority of the simulation calculation method provided above, this application provides two embodiments of using the above method to deduce the vegetation growth process and the coupled evolution with the coastal backshore-dune landform. The two embodiments respectively simulate the deflection of the oblique wind incident from the boundary in the simulation area and the dynamic process of the coupling of vegetation growth and landform evolution.

[0190] Figure 2-Figure 4 This is the first embodiment. The simulated terrain is a two-dimensional ideal sand dune. The representative sections selected are 3100S, 9200S and 18300S, the vegetation growth period is set to 12000S, and the simulation time is 33h. The boundary conditions are that the incident wind direction is 30-60 degrees to the positive direction of the x-axis, and the incident wind speed is defined as the wind speed at a height of 10m above the ground, with a magnitude of 10m / s, and incident from the position x=-83m. Figure 2 Figures 2a-2c show the changes in sediment transport rate, wind speed, and vegetation height at 3100S, 9200S, and 18300S, respectively; Figure 3 3a-3c in the middle show the change and deflection of wind speed on a two-dimensional plane during vegetation growth at 3100S, 9200S and 18300S, respectively; Figure 4 Figures 4a-4c show the changes in sediment transport rate on a two-dimensional plane during vegetation growth at 3100S, 9200S, and 18300S, respectively.

[0191] From the first embodiment, it can be seen that the simulation calculation method provided by the present application can accurately simulate the spatial distribution of oblique wind under two-dimensional complex terrain, and can better restore the acceleration and deceleration process of wind passing through sand dunes in local areas due to changes in slope.

[0192] Figure 5-Figure 7 This is the second embodiment, which adds idealized sand dunes to the measured beach profile of Pingtan, Fujian. The vegetation growth period is set to 12000S, and the simulation time is 33h. The boundary condition is that the wind is incident from the boundary, that is, the angle with the positive direction of the x-axis is 0, and the incident wind speed is defined as the wind speed at a height of 10m above the ground, with a magnitude of 10m / s, and is incident from the position x=-0m.

[0193] Figure 5The simulated beach profile, tide level, and incidence conditions are shown; Figure 6 6a-6d in the figure respectively show the changes of sediment transport rate, wind speed, vegetation height and beach profile at the initial and final moments of the simulation without considering the effect of vegetation on wind speed attenuation of the selected representative profile; Figure 7 7a-7d show the changes in vegetation height, wind speed, sediment transport rate and beach profile at the initial and final moments of the simulation of the selected representative profiles, respectively, taking into account the effect of vegetation on wind speed attenuation.

[0194] From the second embodiment, it can be seen that the simulation calculation method provided by the present application reflects a series of ecological and geomorphological processes better, such as the dynamic response process of vegetation growth and wind and sand landform erosion and accumulation in two-dimensional plane terrain, such as vegetation growth promoting wind and sand accumulation and slowing down vegetation growth.

[0195] It will be understood by those skilled in the art that, unless otherwise defined, all terms (including technical and scientific terms) used herein have the same meaning as those generally understood by those skilled in the art to which this application belongs. It should also be understood that terms such as those defined in common dictionaries should be understood to have meanings consistent with the meanings in the context of the prior art, and will not be interpreted with idealized or overly formal meanings unless defined as herein.

[0196] The meaning of "and / or" described in this application means that the situations where each exists alone or both exist at the same time are included.

[0197] The term “connection” as used in this application may mean a direct connection between components or an indirect connection between components via other components.

[0198] Based on the above ideal embodiments of the present invention, the relevant staff can make various changes and modifications without departing from the technical concept of the present invention through the above description. The technical scope of the present invention is not limited to the contents of the specification, and its technical scope must be determined according to the scope of the claims.

Claims

1. A simulation calculation method for the coupled evolution of two-dimensional coastal sand landform and vegetation growth, characterized by: The specific steps include: Step S1: Obtain key information, including coastal topography data, representative coastal sediment particle size, vegetation growth cycle, tide level, and boundary wind speed; Step S2: The wind speed at each location is obtained through the boundary wind speed and the coastal terrain data. The wind speed at each location includes the established three-dimensional coordinate system, in which the offshore direction is set as the x-axis direction, the shoreward direction is set as the positive direction, the along-shore direction is set as the y-axis direction, the downstream direction is set as the positive direction, and the ground elevation direction in the coastal terrain data is set as the z-axis direction, and the vertically upward direction is set as the positive direction; Step S3: simulating the evolution of backshore vegetation based on the vegetation growth cycle, generalizing the growth elements of the vegetation growth cycle into two parameters, namely, vegetation growth height and vegetation coverage; Step S4: Based on the tide level and representative sediment particle size information, the beach surface humidity of the shoreline is calculated, and the beach surface moisture content is further calculated; Step S5: Based on the vegetation growth height and vegetation coverage obtained in step S3, the beach surface moisture content obtained in step S4, and the wind speed at each location, the equilibrium sediment transport rate at each location on the backshore beach is calculated; Step S6: Calculate the actual sediment transport rate at each location on the backshore beach surface; Step S7: Based on the actual sediment transport rate obtained in step S6, calculate the actual sediment transport rate in the x direction and the y direction; Step S8: Calculate the change of the elevation of the coastal backshore within a time step, and then perform iterative loop calculation from step S1 until the simulation time specified by the user is completed.

2. The method for simulating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 1 is characterized by: In step S2, it is assumed that the wind at the boundary is at an angle The wind speed components in the x and y directions at a distance z from the ground and the total wind speed are: In formula (1), is the angle between the incident wind and the positive direction of the x-axis, is the surface shear stress disturbance value in the x direction, is the surface shear stress disturbance value in the y direction. Specifically, In formulas (2) and (3), is the partial derivative of ground elevation in the x direction, is the partial derivative of the ground elevation in the y direction, α and β are both functions of L / z′, L is the characteristic length of the ground elevation change, z′ is the effective roughness length, and α and β are set to 3 and 0.2 on the windward slope of the coast, and α and β are set to 2 and 1 on the leeward slope of the coast.

3. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landform and vegetation growth according to claim 2 is characterized by: In step S3, the control equations of vegetation growth height and vegetation coverage are converted into linear equations by finite element method for subsequent solution; Specifically, the control equations for vegetation growth height and vegetation coverage are: In formula (4), h v is the vegetation growth height, H v,max is the maximum growth height of vegetation, C v is the vegetation coverage, C v,max is the maximum vegetation coverage, T g For the growth cycle; Among them, the vegetation coverage rate C v At the same time, affected by the accumulation of wind and sand, the influence formula is: In formula (5), is the vegetation coverage rate at the previous moment, H v,max is the maximum growth height of vegetation, z is the ground elevation at that moment, and z pre is the ground elevation at the previous moment; Multiplying formula (4) with the weight function Φ and integrating it over the calculation area, we get In formula (6), is the vegetation growth height at the previous moment, and the weight function Φ is: make Deriving the unit equation Based on the unit matrix multiplication form of the linear equation system formula (8), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100.

4. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landform and vegetation growth according to claim 3 is characterized by: In step S4, calculating the water content of the shoreline beach surface includes the following steps: Step S41, obtaining the height of the capillary water action on the beach, the calculation formula is: In formula (9), h c is the height of the capillary water action on the beach; γ is the surface tension of water, which is 0.0728N / m; θ is the angle between the pore water and the pore wall, which is 0; ρ w is the water density; r m is the average pore radius, taking the value d m / 5;d m is the average radius of the representative particle size of sediment; Step S42, assuming that the water content at the shoreline is 25%, then: w 20mm =25%(h c -h),h=[0,h c ] (10) In formula (10), w 20mm is the water content of the beach surface layer at 20 mm, h c is the height of capillary water action on the beach, and h is the elevation of the groundwater level on the beach, which matches the tide level.

5. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 4 is characterized by: In step S5, based on the vegetation growth height and vegetation coverage obtained in step S3, the beach surface moisture content obtained in step S4, and the wind speed at each location, the equilibrium sediment transport rate at each location on the backshore beach is calculated. The calculation formula of equilibrium sediment transport rate is: In formula (11), Q sat is the equilibrium sediment transport rate; α p is a constant used to convert the height z above the ground b The wind speed is converted into the near-bottom friction velocity u * ; C is the parameter representing the particle size distribution width of sediment; ρ a is the air density, in kg / m 3 ; g is the acceleration due to gravity, in m / s 2 ; d n is the representative particle size of sediment, in m; D n is the reference particle size, in m; u zb is the height above ground z b The wind speed, u tb is the height above ground z b Critical starting wind speed; α p According to the Prandtl-Von Karman law, the calculation formula is: In formula (12), is the Karman constant, z b is the height above the ground at that moment, z′ b is the effective roughness length; u tb The calculation formula is: In formula (13), is the critical friction velocity near the bottom; u zb The calculation formula is: in zb =in * / α P (14) In formula (14), u * is the friction velocity near the bottom.

6. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 5 is characterized by: In formula (13), the critical friction velocity near the bottom is The impacts include: In formula (15), A is the empirical coefficient; g is the acceleration due to gravity, in m / s 2 ρ P is the density of sediment particles, in kg / m 3 ρ a is the air density, in kg / m 3 ;d n is the representative particle size of sediment, in m; Critical friction velocity near the bottom The influence of the beach surface water content is also included, namely: in * th,moist =α w in * (16) In formula (16), is the critical friction velocity near the bottom under the set beach water content, u * is the critical friction velocity near the bottom when the beach water content is low, a w is the conversion coefficient, α w =1+0.1(d 50,ref / d 50 ) 20mm , d 50,ref is the reference particle size, w 20mm is the water content of the beach surface layer at 20 mm, d 50 is the median particle size of sediment particles.

7. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 5 is characterized by: In formula (14), the vegetation height above the ground z b Wind speed u zb The impact is the vegetation growth height or vegetation coverage. If affected by the height of vegetation growth, Effect of vegetation growth height on the friction velocity u near the bottom * The impact is: In formula (17), ρ v is the vegetation factor, ρ v =(h v / H v,max ) 2 ,h v is the vegetation growth height, H v,max is the maximum growth height of vegetation; u *,0 is the friction velocity near the bottom without the influence of vegetation; Γ is the parameter representing the resistance of plant morphology, with a value of 4; If affected by vegetation coverage, The effect of vegetation coverage on wind speed at 10m from the surface is: you 10,veg =u 10 (1-kC ν ) (18) In formula (18), u 10,veg is the wind speed at 10 m above the surface; k is a parameter that characterizes the effect of vegetation coverage on wind speed at 10 m above the surface, and takes different values ​​for different plant species. For typical small upright or flat herbaceous dune plants, k = 1.8, and for small round stemless plants, k = 4.6; C v is the vegetation coverage rate.

8. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 6 or 7, characterized in that: In step S6, the actual sediment transport rate control equation at each location on the backshore beach is: In formula (19), C c is the sediment concentration per unit area, in kg / m 2 ξ is the conversion coefficient between wind speed and sediment movement speed, which is 0.8; T is the time required for erosion or accumulation to complete, which is 1s; C u is the saturated sediment concentration, in kg / m 2 ; p is the reduction coefficient of the source term when the sand source is insufficient; Among them, C u =Q sat / u z ; The reduction factor p is S ε is the amount of erodible sediment, in kg / m 2 ; Solve formula (19) and combine formula (19) with Multiply and integrate over the computational area, and we get make The time term is derived by backward difference to derive the unit equation In formula (20), τ s is the stabilization parameter, u is the velocity vector, is the gradient of the interpolation basis function; In formula (22), Δt is the time step, Δl is the grid length, and the average value is taken here; u x is the wind speed component in the x direction; u y is the wind speed component in the y direction; Based on the unit matrix multiplication form of the linear equation system formula (21), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100.

9. The method for simulating and calculating the coupled evolution of two-dimensional coastal sandy landforms and vegetation growth according to claim 8 is characterized by: Based on the actual sediment transport rate obtained in step S6, the actual sediment transport rate in the x direction and the y direction is obtained, that is, In formula (23), Q b,x is the actual sediment transport rate in the x direction, Q b,y is the actual sediment transport rate in the y direction; The current coastal elevation change is calculated using the mass conservation method, that is, In formula (24), ρ s is the sediment density, λ is the porosity, and z is the ground elevation; Similarly, the mass conservation equation for simulating landform evolution is obtained as follows: Based on the unit matrix multiplication form of the linear equation system formula (25), the overall equation is assembled and the generalized minimum residual method is used to solve the overall equation, and the maximum number of iterations is set to 100.

Images