NUMERICAL MODELING OF THE INTERACTIONS BETWEEN NONLINEAR WAVES AND ARBITRARILY FLEXIBLE VEGETATION

Coastal wetlands are among the natural features with the capability to dissipate wave energy and reduce storm damage. Inadequate representation of wave and vegetation characteristics in numerical models may reduce their capability in predicting wave processes over wetlands. Previous numerical wave models have typically applied simplifications on vegetation behavior. For instance, vegetation stems were usually assumed to be rigid or semiflexible and thus extreme stem deflections could not be captured. In this study, a time-domain nonlinear numerical model based on extended Boussinesq formulation is developed and coupled with a numerical model for vegetation blade dynamics that allows for arbitrary flexibility. Comparison with analytical and laboratory experiments show that the coupled model can adequately predict flow as well as vegetation blade dynamics without the need for any parameter tuning. The model is then used to obtain wave-induced forces on a stem and vegetation blade orientation. Model results indicate that the variation of the vegetative drag coefficient with wave frequency is non-monotonic.


INTRODUCTION
Aquatic marshes are ubiquitous at the land-water interface and are known to have significant effects on the hydrodynamics (Chen et al., 2007), sediment processes (Neumeier and Ciavola, 2004), and nutrient (Barko et al., 1991) and contaminant (Dixon and Florian, 1993) exchange.In addition, observations have documented the dissipative effect of vegetation on surface waves, suggesting that wetlands can mitigate damages to the upland.In turn, storms can damage wetlands and reduce their capability in dissipating storm energy.With increasing interest in using natural features such as wetlands as buffers between land and sea, it is critical to quantify the impacts of vegetation on storm waves and surge, and the impact of storms on wetlands.
Wave dissipation over vegetation has been studied in laboratory and field experiments.Knutson et al. (1982) showed that waves lose almost all their energy over a 30 m Spartina alterniflora marsh.In a study on a seagrass meadow in a fetch-limited environment, Bradley and Houser (2009) showed that higher frequecnies dissipate at a higher rate.Similar results were found in Jadhav et al. (2013) measurements over a salt marsh during a tropical storm.Laboratory experiments have utilized both rigid (e.g.Augustin et al., 2009) and flexible (e.g.Anderson and Smith, 2014) vegetation mimics to investigate wave attenuation over vegetation.These studies examined the effect of stem density and state of submergence on wave dissipation and found that dissipation over near-emergent vegetation mimics is significantly larger than dissipation over submerged vegetation.
Theoretical models have been developed to study wave-vegetation interaction with a focus placed on surface wave damping.The energy extracted from surface waves is commonly represented by a bulk drag coefficient (C D ) which depends on vegetation characteristics and flow parameters.Earlier theoretical studies of Dalrymple et al. (1984) and Kobayashi et al. (1993) derived formulations for the vegetative drag coefficient assuming vegetation can be replicated as rigid cylinders and waves are monochromatic and linear.Mendez and Losada (2004) extended the theory to narrow-banded random waves and semi-flexible vegetation.
Numerical models enable investigating wave-vegetation in a more realistic setting.For instance wave transformation and nonlinear evolution, horizontal and vertical variations in the velocity field, and interactions with currents can be simulated.In large-scale phase-averaged models, vegetation effects are captured by incorporating a bulk drag force based on empirical relationships (Loder et al., 2009;Suzuki et al., 2012).However, the parameterizations for nonlinear wave-wave interactions in these models can result in inaccuracies in shallow water.On the other hand, phase-resolving models can solve nonlinear wave-wave interactions.Boussinesq (Augustin et al., 2009;Blackmar et al., 2013) and Reynolds-averaged Navier-Stokes (RANS) (e.g.Ma et al., 2013;Maza et al., 2013;Zhu and Chen, 2015) numerical models have been used to examine wave attenuation over vegetation.The application of these models to the wave-vegetation interaction problem involves calibration for vegetation-induced drag force using laboratory data, and in general, some simplifications are applied on vegetation behavior.In particular, vegetation elements are assumed to be either rigid or semi-flexible to simplify the formulations used for vegetation motion.For instance, the vegetation model used in (Maza et al., 2013) assumes linear deformation of the vegetation blade, thus it is valid only for small deflections.The vegetation model incorporated in a nonhydrostatic model by (Zhu and Chen, 2015) assumes the stem deflection to be small enough such that the blade tip moves in the horizontal direction.Therefore, these models are not applicable to extreme deflection of vegetation blade where the tip can have significant vertical motion.Such a behavior is observed in seagrass under realistic wave and current forcing.
In order to solve the vegetation motion under wave forcing, vertical and horizontal force balance equations needs to be solved (e.g.Asano et al., 1992;Dubi and Torum, 1996;Utter and Denny, 1996;Abdelrhman, 2007).The forces applied on a vegetation stem are buoyancy, lift, form drag, skin friction, and hydrodynamic inertial forces.Recent studies have developed numerical models for vegetation stem dynamics and have studied segrass blade behavior under combined wave and current Zeller et al. (2014) and pure wave Luhar and Nepf (2016) action.In the model developed by Zeller et al. (2014), a vegetation stem is discretized into a series of rectangular plates attached by joints and torsion springs are placed at the joints to represent blade rigidity.Using this approach, the numerical model can solve motion of vegetation with arbitrary flexibility.It was shown that model prediction for blade motion compared well with laboratory measurements without calibrating for any parameter.The results suggest that the vegetative drag force depends strongly upon the ratio of the horizontal excursion of the blade tip to the wave orbital excursion.
In this study, we develop a numerical model that couples a time-domain nonlinear phase-resolving wave model with a numerical model for vegetation blade dynamics.The wave model is a based on the weakly nonlinear extended Boussinesq formulation of Nwogu (1993) and the vegetation model is a numerical model based on Zeller et al. (2014) formulation.The vegetation blade model solves the comprehensive force balance equation for a stem, and is an improvement over earlier vegetation models as it allows for arbitrary stem rigidity and it doesn't require calibration for the vegetative drag coefficient.The wave and vegetation models are dynamically coupled such that at each time step, the wave model solves for the free surface elevation and depth-varying velocity which is fed into the vegetation model.In turn, the vegetation model solves the instantaneous orientation of blade segments, computes the forces in each segment, and calculates the effective drag coefficient that represents forces applied by the blade on the flow.The effective drag force is then added to the momentum equation in the wave model.It is shown that the model results for vegetation motion compare well with laboratory experiments.The numerical model is then applied to examine vegetation blade motion, hydrodynamic forces applied on stem segments, and the variation of vegetative drag coefficient with wave frequency.

Wave Model
Governing equations: The wave model is a based on the weakly nonlinear extended Boussinesq equations formulated by Nwogu (1993).These equations improve the dispersive properties of classical Boussinesq equations (Peregrine, 1967) for long waves, hence are applicable to deeper waters.The system of extended Boussinesq equations for a homogeneous fluid is comprised of the following continuity and momentum equations, where η is the free surface elevation, u α = (u α , v α ) is the velocity at a reference level z α = −0.531h,h is the still water depth, g is the gravitational acceleration, and ∇ = (∂/∂x, ∂/∂y).In this study, we incorporate a vegetation-induced drag force in the momentum equation which is denoted by f v .
Numerical scheme: The above system of equations is solved numerically using a predictor-corrector method.The approach is similar to that outlined in Wei and Kirby (1995).First, the velocity vector is written in terms of its (u, v) components, then the velocity components and and the free-surface elevation are written as follow: where f (x, y, t) is a source function for generating surface waves, and follows the formulation developed by Wei et al. (1999).Variables U and V are given by, where b 1 = (z α /h) 2 /2, and b 2 = z α /h.E, F, G, F 1 and G 1 are functions of η, u, and v and are given by, where The first and second order spatial derivatives in Eqs (9-12) are discretized using fourth and second order finite difference schemes, respectively.x, y, and t are descretized as x = i∆x, y = j∆y, t = n∆t.The time marching algorithm is a predictor-corrector approach.The predictor step is the third-order explicit Adams-Bashforth scheme.U n+1 is obtained in the predictor scheme as: η and V are written in the same manner but differ in that η n+1 will not include the terms in the above equation that contain F 1 , and V n+1 will include G 1 instead of F 1 .The corrector step is the fourth-order implicit Adams-Moulton scheme.In the corrector step, U n+1 is given by, Similar to the predictor step, η n+1 and V n+1 follow the above equation but η n+1 will not include the last two terms on the left hand side and G 1 replaces F 1 in the equation for V n+1 .To ensure convergence, the corrector step is run iteratively and an error is defined as, where ζ is η, u, or v.The asterisk refers to an intermediate result of the corrector step, and the iteration continues until the corrector step yields an error smaller that a user defined valued for all the variables.
We use e = 10 −6 here.Comprehensive details of the approach are provided in Wei and Kirby (1995), and Wei et al. (1995) in which the mode is extended to fully nonlinear waves.Boundary conditions: Appropriate boundary conditions are required to solve the governing equations.We assume that the domain is rectangular, has constant depth, and is enclosed by vertical impermeable walls.The flow at the solid vertical boundaries (x = 0, x = L x , y = 0, y = L y ) satisfies the impermeability condition: where n is the unit vector normal to the boundary.In addition the the above equation, the conservation of mass in the domain needs to be satisfied.It can be shown that this condition can be satisfied by applying the following conditions at the solid vertical boundaries: In order to minimize reflection from the vertical boundaries, sponge layers are placed adjacent to the boundaries at both x and y directions.The formulation of the sponge layers follows that of Wei and Kirby (1995), and their thickness is chosen to be three times the surface wavelength.

Vegetation model
The vegetative drag force is computed using a numerical model that solves the forces applied on a vegetation stem by the flow and obtains the instantaneous orientation of a flexible blade.The model is based on the formulation developed by Zeller et al. (2014) where a stem is modeled as a series of rectangular plates which are connected at joints by torsion springs.The Young's modulus of elasticity can vary along the blade, and no limitation is applied on stem motion.Thus, the model allows for arbitrary flexibility of the vegetation stem, and can be used to simulate the motion of extremely flexible seagrass blades.Figure (1) shows the blade orientation as computed by the numerical model for a typical case.In this example, a 15 cm stem is divided into 10 segments which are attached by joints that include a torsion spring and a rotary dashpot.The stem is subject to a wave with velocity 0.173 m/s at the free surface and uniform current velocity of 0.132 m/s.The water depth is 1.0 m, and the width and thickness of the blade are 1.00 cm and 0.25 mm, respectively.Vegetation-induced damping is represented by a drag force per unit volume, f v , in the momentum equation (Eq.2), and is calculated as, where C D is the drag coefficient that represents the total energy dissipation in the domain, and A v is the frontal area of the vegetation per unit volume.The total vegetation-induced dissipation in the flow can be written as (e.g.Dalrymple et al., 1984), Where ϵ v is dissipation rate per unit volume, the overbar denotes the time-averaged quantity, T is the wave period, L v is the vegetation stem lenbgth, ρ w is water density, w is the width of the rectangular plate that mimics vegetation blade segments, N v is the vegetattion density (number of stems per unit area), and b is the spacing between plants.Note that although Eq. ( 19) is written for a rectangular plate, it can readily be applied to cylindrical elements by replacing w with cylinder diameter (see Dalrymple et al. (1984)).
The most accurate way to evaluate ϵ v is to compute turbulence production using a turbulence resolving approach, i.e. direct numerical simulation (DNS).However, with the present computational power, this approach is suitable only for very small domains.Therefore, vegetation-induced turbulence is parameterized in large domains.In the present blade dynamics model, the turbulence is parameterized using a comprehensive set of forces and moments applied on a vegetation stem.These forces include inertial forces, lift, form drag, skin friction, buoyancy, and moments applied by torsion springs which determine the rigidity of the stem (Zeller et al., 2014).Therefore, the total dissipation can be parameterized as, where F D , F I , F f , F L , and F B are the drag, inertial, skin friction, lift, and buoyancy force, respectively, u rel is the relative velocity of the blade with respect to the fluid, and the dimension along the stem is denoted by s.As discussed by Zeller et al. (2014), the dominant forces are the drag and skin friction forces.Following Abdelrhman (2007), the drag force from blade segment i is formulated as, where C D i is the drag coefficient, θ i is the angle of attack of the fluid, u i is the fluid velocity at the segment center, and θ i is the angle between the segment and the vertical direction.The angle of attack is given by, where i and k are unit vectors i then horizontal and vertical directions, respectively.The segment drag coefficient uses a parameterization for flow encountering a rectangular plate and is given by, where θ i is the angle between segment i and the vertical direction.The skin friction force is given by, where C f is the friction coefficient and is estimated to be 0.02 (Abdelrhman, 2007).Details of the formulations used for the inertial and lift forces are provided in Zeller et al. (2014).As mentioned, the blade segments are connected via joints which are equipped with torsion springs and rotary dashpots.
The torque applied by torsion springs on blade segments determines blade rigidity and depends on the Young's modulus of elasticity as well as orientation and geometry of blade segments: where β i = nE i I i /L and θ 0 i is the initial angle between segment i and the vertical direction.In addition to torsion springs, rotary dashpots are placed at joints to damp spurious oscillations.The torque applied by the rotary dashpot at joint i is given by, where γ is the damping coefficient.To calculate the effective drag coefficient (C D ) in Eq. ( 18), the forces and moments applied by the flow on each blade segment need to be calculated, and the forces and moments, as seen in Eqs.(21-25) depend upon fluid velocity at the segment centers as we as orientation and velocity of the segments.The instantaneous horizontal depth-varying velocity at segment centers is computed using the wave model as (Wei et al., 1995), The model accounts for the effect of vertical velocity but its impact on drag coefficient is found to be negligible in the present parameter space.As detailed by Zeller et al. (2014), a system of ordinary differential equations is derived based on the vertical and horizontal force balance and moment balance for each blade segment.The ODE system is then solved numerically and the instantaneous orientation of the segments is obtained.Once forces are computed, the total dissipation rate is calculated using equation ( 20).The effective drag coefficient which represents the total energy dissipation, is then obtained by equating Eqs ( 19) and ( 20).The vegetation and the wave model are dynamically coupled such that the drag coefficient computed using the vegetation model will be passed to the wave model at each time step, the wave model computes the updated flow field, and the new flow filed is passed to the vegetation model at the next time step.The velocity of the segment i, ẋi , is computed using the θ.Therefore, the effective drag coefficient is obtained at each time step using instantaneous wave field and blade orientation.

MODEL RESULTS
A typical model output is illustrated in Figure ( 2). Figure 2 (a), (b), and (c) respectively show the computational domain around the vegetation blade, the horizontal velocity profile at the blade location as calculated using equation ( 26), and a zoomed in view of blade that shows deflection under surface waves.In this example, wave amplitude, wave period, and water depth are, a = 0.10 m, T = 4.00 s, and h = 1.00 m, respectively.The model is run in the 1D mode with spatial resolution dx = 0.25 m, and temporal resolution dt = 0.04 s.The Ursell number is U r = 0.36.This case includes a single vegetation stem with length L v = 0.40 m, which is divided to n v = 10 segments.The thickness and the width of the blade are b = 0.25 mm and w = 1.00 cm, respectively.Water density is ρ w = 1000 kg/m 3 , and the vegetation density is assumed to be ρ v = 920 kg/m 3 .The Young's modulus of elasticity of the blade is E = 0.29 GPa, and the moment of inertia, I, in the segment connected to the seabed is assumed to be 30 times the moment of inertia in other segments to represent increase stiffness at the root.This increase in I is arbitrary and is not essential for a successful simulation.The damping coefficients of rotary dashpots at all the joints are the same and are γ = 10 −3 N.m/rad.Since the present numerical model is composed of a wave and a blade dynamics models, the wave and vegetation models as well as the performance of the coupled model should be validated.The Nwogu (1993) formulation has been widely validated previously.In particular, Wei and Kirby (1995) show that their numerical model performs well in predicting wave height and shape in several wave propagation problems that involve solitary and random waves.Since our wave model uses the same formulation and numerical scheme as Wei and Kirby (1995), we do not repeat all the validation tests for the wave model and only test the model performance for one-dimensional solitary wave propagation.Figure (3) compares the free-surface elevation at t = 50 s and x = 171 m as computed by the model for a solitary wave with the analytical solution.C is the analytical solitary wave speed.In figures (3)a and b, the nonlinearity parameter, δ = H/h, is 0.10 and 0.20, respectively.The water surface elevation computed from the wave model shows a slight reduction in waveheight indicating a small numerical diffusion.Furthremore, the phase speed of the numerical solution is slightly smlaller than that of the anlytical solution.The deviation from the anlytical solution in both waveheight and phase speed increases with δ.These properties of the numerical solution are consistent with the results of Wei and Kirby (1995).The blade dynamic model was validated with laboratory experiments by Zeller et al. (2014).The study showed that there is a satisfactory agreement between the numerical blade model prediction for the location of the blade tip with laboratory measurements.The study used velocity induced by combinations of currents and waves in a laboratory flume as the model input.Our approach to validating the coupled model is similar, but differ in that we use the wave model to generate the velocity input to the vegetation model.To quantify the relative blade rigidity, a dimensionless rigidity parameter is defined as, where U max is the maximum wave orbital velocity in the water column over a wave period.This parameter normalizes the blade rigidity with the drag force.Zeller et al. (2014) discuss that the time of maximum blade deflection corresponds to the point of maximum turbulence production, and thus maximum wave dissipation.Therefore, the coupled model should accurately predict the maximum deflection to accurately compute the vegetative drag force.
To validate the coupled wave-vegetation model, we compare model prediction for blade defletion for different values of λ r with experimental results of Zeller et al. (2014).We note that the laboratory experiments were conducted in a flume that generates waves and currents simultaneously, and that currents in some experiments were of the same order of magnitude as the maximum wave velocity.However, the internal wavemaker in the present wave model can only generate waves.Therefore, we increase the wave amplitude such that the maximum orbital velocity is equal to the maximum combined wave and current velocities from experiments.Since the waves are in shallow water region and wave orbital velocity is relatively uniform in depth, this approach is not expected to cause a large error in estimating the maximum horizontal velocity in the laboratory experiments.Similar to Zeller et al. (2014), we use the elevation of the blade tip at the point of maximum deflection as normalized by the stem length, z * t = z t /L v , to evaluate blade deflection.A comparison between the present model results for maximum blade deflection and laboratory measurements is shown in Figure (4).The results show a satisfactory agreement.As expected, an increase in blade rigidity decreases blade deflection.The numerical model is initialized using the laboratory data on wave and current velocities and wave period.Wave velocities, current velocities, and wave periods in the experiments used for validation vary between 0.212−0.232m/s, 0.0274−0.274m/s, and 3.91 − 5.11 s, respectively.The vegetation parameters are the same as those in Figure (2) but differ in that γ = 10 −5 .We note that some discrepancy between the present numerical model and laboratory results is inevitable.One source of discrepancy can be due to bottom boundary layer effects which is not captured in the numerical model but affects the velocity profile in the laboratory experiments.Therefore, the forces applied by the flow on the blade may differ in blade segments in the vicinity of the bed.As discussed above, vegetation-induced wave dissipation has a strong dependency on the deflection of vegetation stems.If the deflection is large, the vegetation frontal area decreases and the drag coefficient decreases as a result.The present numerical model can be used to obtain the instantaneous position of blade segments under wave action.Figure 5  Storm damages can reduce wave dissipation effect of wetlands.Since the present numerical model can give forces applied on each vegetation stem segment, it can potentially be used to predict damages induced by storms.In particular, the model can identify the segments at which the stem may break off, or if the stem can be uprooted under a particulate wave forcing.The present model can be employed to assess the impact of important wave and vegetation parameters on vegetative drag coefficient.Figure (7) shows the variation of C D with wave period is non-monotonic.As seen, in the present parameter set, an increase in wave period results in an initial decrease in the drag coefficient while C D starts to increase at larger wave periods.The wave and vegetation parameters used here are the same as those in the example discussed in Figure (2) except that L v = 0.15 m, and a = 0.024 m.

CONCLUSION
In this paper, we outline the development of a numerical time-domain wave-vegetation interaction model.The wave model is a nonlinear phase-resolving model based on the extended Boussinesq formulation of Nwogu (1993).It is shown that the wave model accurately predicts solitary wave propagation over a long distance.The vegetation model is based on the blade dynamics fomulation of Zeller et al. (2014).The vegetation model receives flow information, particularly the depth-varying vertical and horizontal wave orbital velocity and solves the instantaneous position of blade.The vegetation model divides the blade in several segments that are connected by joints equipped with torsion springs that determine the blade rigidity and rotary dashpots.The numerical blade dynamics model does not impose any limitation on the rigidity of the segments, hence it can simulate blade motion of arbitrary flexibility under wave action.The numerical blade dynamics model accounts for a comprehensive set of forces applied on vegetation including the drag, inertial, lift, buoyancy, skin friction, and gravity forces.The blade model computes the instantaneous drag coefficient which will be used to obtain the effective drag force in the momentum equation of the wave model.The wave and vegetation models are dynamically coupled ensuring a two-way feedback.It is shown that the results of the coupled model for blade motion compares well with laboratory experiments.
The model is employed to examine the instantaneous position of the blade segments, the forces applied by waves on the segments, and the variation of the vegetative drag coefficient with wave period.Although comparison with experiments were not conducted due to lack of data on wave forces in stem segments, model results are qualitatively consistent with expectation.For instance, the forces applied on the lowermost blade segment were found to be significantly larger than forces applied on the free end of the blade.Furthermore, it is shown that the vegetative drag force does not vary monotonically with wave period.In the example discussed here, the drag coefficient is largest at the highest and the lowest ends of the range of the wave period.
Although the present model includes a single stem, extension to multiple stems is straightforward.and is a subject of future work.Future work will also include extension of the wave model to a fully nonlinear formulation.Furthermore, the performance of the model in predicting wave dissipation will be compared with laboratory experiments.

Figure 1 :
Figure 1: Numerical model of vegetation stem.Dots show joints that connect segments.

Figure 4 :
Figure 4: Comparison of numerical models results for dimensionless vegetation blade tip against rigidity parameter.
(a)  and (b) show the horizontal and vertical positions of the centers of the segments i s = 1, i s = 5, and i s = 10, where i s = 1 refers to the segment at the free end and i s = 10 refers to the segment connected to the seabed.The wave amplitude in this simulation is a = 0.034 m, L v = 0.15, h = 0.40 m, and the rest of the parameters are the same as those used in Figure (2).

Figure 5 :
Figure 5: (a) Instantaneous vertical (a) and horizontal (b) location of blade segment center.i s = 10 refers to the lowermost segment, a = 0.034 m, n s = 10, L v = 0.15 m, h = 0.40 m, and the rest of the parameters are the same as those in figure 2.
Figure (6) (a) and (b) respectively show the horizontal and vertical computed by the model in segments i s = 1, 5, and10.Since the surface wave is in shallow water range, the horizontal forces are significantly larger than the vertical forces.The vegetation and wave parameters in are the same as those in Figure (5), and the forces are computed at segment centers.

Figure 6 :
Figure 6: Instantaneous (a) horizontal and (b) vertical forces in i s = 1, 5, and10.Forces are computed at segment centers, and the parameters are the same as those in Figure (5).

Figure 7 :
Figure 7: Variation of vegetation drag coefficient with wave period.Wave and vegetation parameters are the same as those in Figure(2) except that a = 0.024m, L v = 0.15 m, and E = 0.30 GPa.