Numerical Modelling of Wave–Vegetation Interaction: Embracing a Cross-Disciplinary Approach for Bridging Ecology and Engineering for Nature-Inclusive Coastal Defence Systems

Abstract

Coastal areas are increasingly at risk due to climate change, necessitating innovative mitigation approaches. This study explores the integration of living environments, particularly aquatic vegetation, with conventional defence systems to provide socially acceptable and nature-inclusive coastal defence systems. Through examining the published literature, this study identifies two perspectives: engineering and ecological. From an engineering perspective, wave propagation models and simulation techniques for wave–vegetation interaction are identified. Ecologically relevant coastal and marine vegetation is presented, and based on its ecological features (morphology, biomechanics, buoyancy, and variability) a novel ecological categorization framework is developed. The results challenge the notion of a strict divide between ecological and engineering approaches. Analysis of existing wave–vegetation models reveals that many engineers consider the ecological features of vegetation-induced wave attenuation studies. However, computational limitations often lead to simplifications. Furthermore, complex models, while offering detailed ecological insight, are often limited to small-scale experimental domains. Conversely, simpler models, suitable for large-scale engineering problems, may lack ecological detail. This suggests a potential future approach numerical modelling that combines high-resolution models for small areas with large-scale, implicit engineering models operating at the ecosystem scale.

1. Introduction

The Paris Agreement, reached in 2015, stands as a formal acknowledgment and call for action against the ongoing environmental crisis. It pledged to strengthen the global response and combat climate change by limiting the temperature increase to well below 1.5–2 °C above the pre-industrial level [1]. A decade later, many countries continue to face challenges in meeting their emissions targets [2]. So much so that future projections are no longer necessary to observe the consequences of climate change. Ample evidence of sea-level rise, global temperature increase, extreme storm events, and glacial retreat has been extensively documented in the Intergovernmental Panel on Climate Change (IPCC) report [3].

Coastal areas are particularly susceptible to the impacts of climate change due to a combination of factors. Firstly, the geographic features of coastal cities in low-lying regions, nestled between river deltas and soft beachfront sediments, increase the likelihood of compounded events featuring multiple hazards occurring concurrently [4]. For example, Hurricane Idalia [5], which struck Florida’s west coast in 2023, induced hazardous flooding resulting from a combination of extreme storm surges along the beachfront and heavy precipitation leading to urban flooding upstream. Secondly, the frequency and intensity of extreme events are anticipated to increase in tandem with sea-level rise [6,7]. In quantitative terms, global probabilistic projections indicate that a storm with a 100-year occurrence probability today will occur annually by 2050 [8]. Finally, the highest population concentrations globally are located in coastal communities [9]. Projections estimate that by 2100, 52% of the global population will be at increased risk of flooding (assuming representative concentration pathway (RCP) 8.5 among other model state variables) [10].

Decision making for protection against climate change through investments in coastal infrastructure is frequently directed towards localized solutions, with conventional defence units also termed hard infrastructure [11]. Groynes are one type of these hard solutions used to interfere with sediment transport and mitigate coastal erosion [12]. Dikes and seawalls are other examples used for flood protection [13].

A more advanced form of hard infrastructure is the adaptive systems that were developed under the predict-and-employ paradigm. Under this paradigm, only when hazardous flooding conditions are predicted is the structure mechanically employed. Perhaps two of the most famous examples are the “The Maeslantkering” storm surge barrier in the Netherlands [14] and the “Mose” flood gates in the Venetian lagoon, Italy [15]. The latter systems promise a smaller ecological footprint and less interference with the marine habitat; coupled with an increased social acceptance [16,17].

A resilience paradigm in coastal defence, on the other hand, addresses previously ignored complexities, including cross-disciplinary synergies between ecosystem services and defence infrastructure, social-ecological systems, adaptability, and multifunctionality. On this basis the concept of embracing nature-inclusive designs emerged in guidelines and among policy makers at first [18,19,20,21], and then, in engineering practice [22,23]. In the near future, investments in nature-inclusive solutions across Europe and globally are expected to increase. For example, the European Commission’s research and innovation program [24] presents valuable resources, offering not only funding opportunities for researchers and developers but also fostering international collaborations to drive innovation inspired by nature.

In a broad sense, building with nature is defined as a coastal adaptation strategy that integrates living ecosystems into engineering structures, creating a mutually beneficial partnership that maximizes opportunities for the hosting environment on multiple fronts (environmentally, socially, and economically) [25].

Preservation of naturally occurring mangrove wetlands could be an illustration of NBS in coastal management [26,27]. Mangroves are salt-tolerant trees that thrive in tropical and subtropical areas along intertidal zones. Their complex and intertwining root systems promote accretion and serve as a green barrier against storm surges [28]. They also create an interaction biosphere with local communities (e.g., [29,30]). The role they played during the 2004 Indian Ocean tsunami is a testimony [31]. Today they make up part of the coastal management package in multiple countries worldwide including Bangladesh, The Phillipines, and the United States of America (USA) [32].

Vegetated foreshores with seagrasses are another nature-inclusive coastal defence system [33,34]. Seagrasses grow in light-abundant littoral zones and occupy subtidal or intertidal areas [35]. They form productive ecosystems and reproduce to form meadows that can occupy large areas of seabed in sandy or muddy sediments [36]. They stand out with the ecosystem services they offer and their resilience to climate change. Seagrasses function as feeding and nursery grounds for marine organisms [37]. Furthermore, they contribute to climate and coastal resilience through multiple functions, including carbon sequestration [38], stabilization of subsea sediments, and attenuation of wave energy within the water column [39].

A typical nature-inclusive design that incorporates seagrasses as a first line of defence preceding the conventional vertical concrete seawall is shown in Figure 1. This design leverages the ecosystem services provided by the vegetation to attenuate wave height, mitigate foreshore erosion, and enhance biodiversity. The seawall is the second line of defence in this hybrid design; it benefits from the upstream wave attenuation and features an optimized eco-friendly design with a reduced wall height.

1.1. Building with Nature: Execution Pathway

The pathway to implementing nature-inclusive solutions in coastal systems is a multifaceted and complex process that requires the engagement and contribution of multiple components. Figure 2 provides a structured overview of this process, with the horizontal axis representing the progression of time and maturity, and the vertical axis depicting the key components involved. A positive reinforcing cycle demonstrates ongoing two-way interactions and long-term learning between stakeholders and the final solution.

The process towards realizing nature-inclusive coastal defence systems begins with phase one: “Engagement and stakeholder involvement”. While the core principles and goals of NBS are universally applicable, integrating the local context proves pivotal for success. Therefore, this initial phase revolves around actively involving local stakeholders who are vulnerable to coastal hazards. These stakeholders primarily include residents of high-risk beachfront communities, government bodies, environmental organizations, and economic entities (e.g., fishermen, tourist resorts). The primary objective in this phase is to identify the problem (coastal flooding, damage to property, economic losses) and foster a collective commitment (e.g., through workshops, interviews) to embracing nature-inclusive defence systems. At this stage, establishing a trust relationship between the decision-making body and the local stakeholders is crucial; stakeholders must experience a sense of ownership and commitment throughout the execution pathway.

The second phase, “Development tools and techniques”, marks the beginning of the design process, building upon the insights gathered from stakeholder engagement (e.g., residents express desire for a defence system that preserves beach access). This phase typically involves a cross-disciplinary approach among coastal engineers, marine ecologists, policymakers, and designers. The objective is to integrate scientific knowledge in order to devise the most effective solution while addressing the societal needs and ambitions of the stakeholders (as identified in phase 1). For ecologists, the focus lies on assessing ecosystem services. This encompasses both negative and positive effects on fauna and flora, the food chain, species interactions, and the conservation of native species. Specific attention goes to the potential disruptions of and interference with the local ecosystem. For example, this includes preventing an imbalance in the food chain or the introduction of invasive species, or at least only tolerating invasive species with limited impact.

Engineers, on the other hand, focus on understanding the hydrodynamics and morphodynamics of the system. They study the interactions occurring within the nature-inspired defence system (e.g., seagrasses or mangroves), paying close attention to the dynamics of wave–vegetation interaction and the resulting wave energy dissipation. Additionally, they evaluate the sediment transport balance (focus on erosion/accretion). In designs where hybrid solutions (combination of nature-inspired with conventional grey infrastructure) are employed, engineers also consider the interaction between waves and structures. Accordingly, this phase concludes with the shaping of the final nature-inclusive design of the coastal defence, incorporating input from both ecologists and coastal engineers.

The third and final phase, “Implementation of nature-inclusive solutions”, transitions the nature-inclusive coastal defence system from design to physical reality. This phase consists of site preparation, habitat restoration, and the installation of natural barriers. “Natural barriers” in this context refer to actions such as the restoration of mangrove forests or revegetation of the site with seagrasses. Successful implementation extends beyond these initial steps and commands ongoing efforts. Firstly, it requires monitoring and adaptive management to assess the effectiveness of the implemented solutions. This includes (1) monitoring the evolution of the beach topography and nearshore bathymetry, (2) evaluating ecosystem services (e.g., water quality, maintaining biodiversity, and energy flow in the food web), and (3) assessing coastal protection against waves and storm surges (e.g., reduction in wave height, storm surge buffering). Secondly, it involves sustained collaboration, capacity building, and ongoing engagement with stakeholders, as depicted in the positive reinforcing cycle outlined in Figure 2. This approach promotes resilience and adaptive governance by raising awareness among local communities and fostering a sense of ownership and involvement.

1.2. Objectives and Structure

The present study focuses on addressing the development of NBS by exploring wave–vegetation interaction from cross-disciplinary engineering and ecological perspectives, particularly within the domain of the second phase, “Development tools and techniques” (described in Figure 2). The hypothesis under investigation theorizes that the research and development of NBS in fundamental academic research continues to be divided (polarized) between ecologists and engineers. There is a partial consideration of crucial biomechanical processes and properties by engineers and a neglect of engineering applications by ecologists.

The aim of this study is to test this hypothesis and offer recommendations for comprehensively integrating ecological and engineering perspectives in the development of NBS for coastal protection. To achieve this goal, the following objectives are identified:

• For the hydrodynamic modelling and engineering designs, the objective of this study is to identify the numerical models and techniques employed by engineers to address wave propagation and interactions with vegetation within the context of fundamental research on NBS.
• For the ecological assessment, the objective is to develop a conceptual framework to categorize the engineering tools (numerical models identified in the first objective) based on the extent to which they account for the ecological and biological nature of aquatic vegetation.

This paper is structured as follows: the methods used in this study are presented in Section 2; following that the numerical models used to simulate wave–vegetation interactions are presented in Section 3; the coastal and marine vegetation are presented in Section 4. A link between the numerical models and the ecological aspects is presented through a categorization framework in Section 5. The results obtained through applying this framework are presented in Section 6 and discussed in Section 7 by also devising recommendations for future research. This paper is concluded in Section 8.

2. Methods

The methodology of this study is guided by standards of a cross-disciplinary research process, integrating both engineering and ecological perspectives to address the research objectives. The process is structured into multiple steps as outlined below.

The conceptual chart of this study is depicted in Figure 3. The stated hypothesis is tested through the objectives presented earlier (Section 1 and also summarized in the figure), both from an engineering and an ecological perspective. The closed loop connecting these objectives with the hypothesis, labelled “YES/NO”, indicates the acceptance or rejection of the stated hypothesis. Additionally, on the right-hand side of the sketch, the study design is summarized using a five-step process.

The process begins by defining and describing the key elements present in this research. From an engineering perspective, this includes detailing the numerical models used to simulate wave propagation and wave–vegetation interactions. From an ecological perspective, this involves describing the ecological roles and physiology of vegetation with wave attenuation roles that are relevant to nature-inspired coastal defence systems. This mainly includes seagrasses, salt marshes, and mangroves.

2.1. Study Design

Keyword Identification

This step involves identifying relevant keywords that characterize the focus of this study. These keywords include terms such as “numerical models”, “wave–vegetation interaction”, and “wave dissipation” for searching published numerical models, while terms like “seagrass ecology”, “seagrass species”, “seagrass physiology”, “mangroves”, and “nature-inclusive solutions” are used to approach the ecological perspective. It is important to note that these keywords are not exhaustive; other terms with similar meanings or direct relevant are also utilized.

2.2. Systematic Literature Review

A systematic literature study is conducted using the identified keywords to search through electronic databases, mainly Google Scholar and Web of Science. Additionally, related documents, either through direct citation or subject relevance (e.g., identified by the algorithms operated by Publishers or by Scopus) are included in the data collection phase. All the information collected is exclusively derived from peer-reviewed journal articles.

2.3. Ecological Categorization Framework

The ecological framework for this study is developed by systematically identifying and organizing the characteristics of aquatic vegetation that are relevant to wave–vegetation interactions. These characteristics are categorized into distinct features that serve as the basis for classifying numerical models.

2.4. Categorization of Numerical Models

The ecological framework is applied to categorize the numerical models. This step involves matching model capabilities with ecological categories, highlighting how different models address various ecological aspects such as morphological traits (e.g., height and density), biomechanics, and variability.

2.5. Conclusions and Recommendations

Conclusions are drawn based on the results of the categorization and analysis. The study provides recommendations for future research, emphasizing the integration of engineering and ecological perspectives in numerical modelling.

3. Numerical Modelling

3.1. Wave Propagation Models

In its simplest form, the mathematical description of water waves can be achieved using linear wave theory, also known as Airy theory [40]. To account for the non-linear shape of water waves, characterized by the asymmetry between crest and trough, a more advanced description is provided by Stokes’s theory [41]. Using these analytical descriptions in combination with field data, the earliest form of wave propagation prediction using empirical methods dates back to 1947 [42]. For a more complex and accurate representation of wave propagation phenomena, the mathematical description becomes dependent on partial differential equations (e.g., the Navier–Stokes (NS) equations [41,43]) that can only be solved numerically. It is from here that the development of and reliance on numerical wave propagation models in the field of coastal engineering came to prominence.

Several numerical techniques can be used to describe wave propagation and transformation across various temporal and spatial scales. These techniques can be grouped into two exclusive categories: (1) phase-averaging models and (2) phase-resolving models.

3.1.1. Phase-Averaging Models

Phase-averaging models were developed following advances in the statistical theory of ocean waves and are based on the energy balance equation. In these models, the wave evolution is calculated based on the energy of each individual frequency component. The wave spectrum is determined by considering the effects of advection, refraction, and physical processes such as wind input, bottom boundary effects (bathymetry), and non-linear interactions (e.g., wave–wave interactions, wave growth) [44]. The range of physical processes included in these models has been continuously extended alongside the increasing computational power. These models are particularly advantageous for simulating large domains and are used to establish forecasting platforms. The reason for this is the computational efficiency resulting from the fact that phase information is omitted and there are no restrictions on the size of the computational grid in relation to the wave characteristics. However, the applicability of these models is limited in shallow waters, where wave transformations dominate and the statistical description underlying phase-averaging models no longer holds.

One of the best-known and most widely used phase-averaging models is SWAN (Simulating Waves Nearshore) [44]. The model resolves the wave kinematics based on the action balance equation, which describes how the action density spectrum (N) evolves over time and space [44], and is given by:

$$N(\sigma ,\theta )=\frac{E(\sigma ,\theta )}{\sigma }$$

(1)

where N is a normalized quantity that describes the wave energy distribution across frequencies $\sigma$ and wave propagation directions $\theta$, and E is the phase-averaged wave energy spectrum. The full wave action balance equation used to calculate the wave spectrum can be expressed as

$$\frac{d}{dt}N+\frac{d}{dx}{c}_{x}N+\frac{d}{dy}{c}_{y}N+\frac{d}{d\sigma }{c}_{\sigma }N+\frac{d}{d\theta }{c}_{\theta }N=\frac{S}{\sigma }$$

(2)

where the terms on the left-hand side represent the partial derivatives of N with respect to time t, and with $c_x$, $c_y$ being the wave velocities in the x and y directions, respectively; and with $c_{\sigma }$ being the coefficient related to changes in relative frequency due to water depth, and $c_{\theta }$ being the coefficient related to changes in wave direction. Last, the total source term S represents all the physical processes that influence the wave energy density spectrum [44].

3.1.2. Phase-Resolving Models

The second category, phase-resolving models, solves the evolution of each individual wave over time. Within this category, different formulations of the governing equations are utilized, leading to the following sub-categories: Mild-slope equation models (MSE), Boussinesq equation models (BE), non-hydrostatic models (NH), and Navier–Stokes models (NS).

(i)

Mild-slope models:

Mild-slope models are based on mild-slope equations that describe linear waves within domains where the variation in bathymetry is gradual over a horizontal distance comparable to the wavelength [45]. These equations are expressed as:

$$\begin{array}{cc}\hfill \frac{\partial \eta }{\partial t}& =B\varphi -\nabla ·(A\nabla \varphi )\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill \frac{\partial \eta }{\partial t}& =-\mathbf{g}\eta \hfill \end{array}$$

(4)

where $\eta$ is the surface elevation, $\varphi$ is the velocity potential, ∇ is the horizontal gradient operator, t is time, $\mathbf{g}$ is gravitational acceleration, and A and B are coefficients related to the dispersion relationship. These equations are capable of simulating processes such as wave reflection at boundaries, as well as transformations such as shoaling and refraction due to changes in bathymetry, and diffraction during barrier penetration. Mild-slope numerical models are computationally efficient and are widely used in engineering for the design of coastal structures and wave penetration studies. For instance, MildWAVE [46] is a mild-slope model that supports versatile wave generation techniques [47] and is suitable for both coastal and offshore applications [48,49]. Nevertheless, it is important to recognize the general limitations of mild-slope models in regions characterized by steep bathymetric features and extreme transformations (e.g., depth-induced wave breaking cannot be captured).

(ii)

Boussinesq equation models:

Boussinesq equation models can address non-linear interactions and wave breaking phenomena that are typically beyond the scope of mild-slope models. The Boussinesq equations are based on a set of partial differential equations describing the propagation of surface gravity waves in shallow water [50]. The basic formulation consists of two coupled equations, one describing the evolution of the surface elevation and the other the depth-averaged horizontal velocity, as shown below, respectively:

$$\begin{array}{cc}\hfill \frac{\partial \eta }{\partial t}+\frac{\partial }{\partial x}\left((d+\eta )\mathbf{U}\right)& =0\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \frac{\partial \mathbf{U}}{\partial t}+\mathbf{U}\frac{\partial \mathbf{U}}{\partial x}+\mathbf{g}\frac{\partial \eta }{\partial x}& =\frac{d}{2}\frac{{\partial }^3\left(d\mathbf{U}\right)}{\partial x^2\partial t}-\frac{d^2}{6}\frac{{\partial }^3\mathbf{U}}{\partial x^2\partial t}-\mathbf{g}\frac{\mathbf{U}\left|\mathbf{U}\right|}{C^2(d+\eta )}+{\nu }_T\frac{{\partial }^2\mathbf{U}}{\partial x^2}\hfill \end{array}$$

(6)

where $\eta$ is the surface elevation, d is the water depth, $\mathbf{U}$ is the velocity vector, $\mathbf{g}$ is the gravitational acceleration, C is the Chezy resistance, and ${\nu }_T$ is the eddy viscosity coefficient. Recent developments have extended the applicability of these models to deeper waters and improved the non-linear effects by introducing several vertical layers and higher-order polynomial approximations [51]. However, the adoption and use of Boussinesq equation models in industry has been limited by their susceptibility to numerical instabilities, attributed to the presence of higher-order derivatives in both time and space.

(iii)

Non-hydrostatic models:

Non-hydrostatic models address the limitations of the previously presented models (mild-slope and Boussinesq) by solving the Navier–Stokes equations under two basic assumptions: a free surface and a constant density. Under these conditions, the equations are integrated across layers, yielding depth-averaged quantities for single-layer scenarios and layer-specific quantities for multiple-layer cases [52,53]. The governing shallow water equations in a two-dimensional medium are:

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}& =0\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill \frac{\partial u}{\partial t}+\frac{\partial uu}{\partial x}+\frac{\partial wu}{\partial z}& =-\mathbf{g}\frac{\partial \eta }{\partial x}-\frac{\partial q}{\partial x}-\frac{1}{\rho }{\mathbf{F}}_x\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill \frac{\partial w}{\partial t}+\frac{\partial uw}{\partial x}+\frac{\partial ww}{\partial z}& =-\frac{\partial q}{\partial z}-\frac{1}{\rho }{\mathbf{F}}_z\hfill \end{array}$$

(9)

where u and w are the horizontal and vertical velocity components, respectively, $\mathbf{g}$ is the gravitational acceleration, $\rho$ is the water density, q is the non-hydrostatic pressure, and ${\mathbf{F}}_x$ and ${\mathbf{F}}_z$ are external forces acting on the fluid in the horizontal and vertical directions, respectively.

These models have a reasonable computational cost while effectively capturing wave transformation processes, including wave breaking, and have fewer constraints on bathymetry and domain features such as water depth. The combination of capturing wave non-linearities, computational efficiency, and numerical stability, positions non-hydrostatic models as a favoured tool among engineers in both industrial applications and fundamental research (e.g., SWASH (Simulating Waves till Shore) [54,55]).

(iv)

Navier–Stokes models:

Navier–Stokes models are the most complex, solving the full Navier–Stokes equations with high resolution both horizontally and vertically. Assuming an incompressible fluid and a simplified mathematical description for the Navier–Stokes equations, the fluid dynamics can be described by the continuity (Equation (10)) and momentum (Equation (11)) equations, written as:

$$\begin{array}{cc}\hfill \nabla ·\mathbf{U}& =0;\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill \rho \frac{\partial \mathbf{U}}{\partial t}& =-\nabla p+\rho \mathbf{g}+\mu {\nabla }^2\mathbf{U}\hfill \end{array}$$

(11)

where $\mathbf{U}$ is the velocity vector, t is time, $\rho$ is density, p is pressure, $\mathbf{g}$ is the gravitational acceleration, and $\mu$ is the fluid viscosity.

Solving these equations allows detailed representation of wave transformation processes, turbulence, and wave breaking. However, their complexity and accuracy comes with a higher computational cost, which restricts their applicability to limited temporal and spatial domains. Common examples of NS models include the mesh-based volume of fluid (VOF) method [56,57] and the meshless smoothed-particle hydrodynamics (SPH) method [58,59]. These models offer advanced capabilities and have no limitations in terms of water depth or bathymetric features, but require careful consideration of expensive computational resources.

Each of the numerical models presented has unique features and limitations. For example, phase-averaged models offer computational efficiency but are not applicable in shallow coastal regions with wave breaking. Conversely, Navier–Stokes models are able to accurately capture violent flows, but in a very limited temporal-spatial context (limited to a few waves in a small domain). Therefore, a common strategy to overcome these limitations is coupling (connecting) different models. Through coupling, computationally efficient models are configured to resolve larger domains, and then, communicate with more intricate models at the numerical boundaries; such a setup enables high-resolution simulations of the physical processes. For example, the North Sea phase-averaged forecast model provides boundary conditions to a higher-resolution harbour model within a non-hydrostatic framework [60]. Coupling can also enhance functionality, such as coupling Mildwave with Nemoh [61] to study floating bodies [49]. Other examples include coupling wave propagation models with external structural libraries to simulate fluid–elastic structure interactions [62]. A typical cross-section of the eastern Mediterranean, shown in Figure 4, illustrates the conditions under which the models can be applied and how they can be combined. The image was constructed using a satellite base map from Google Maps and elevation data from the Marine Geoscience Data System (MGDS). Numerical wave propagation models are applied based on the water depth regime. Phase-averaging models resolve the wave propagation in the deep-water conditions over the large domain and are coupled with phase-resolving models for the intermediate and shallow waters. In this illustration, non-hydrostatic models can handle the relatively large phase-resolving area depicted, while Navier–Stokes models can simulate violent wave breaking along a small coastal strip.

3.2. Wave–Vegetation Interaction Models

The wave propagation models previously presented are used to address coastal engineering problems. To adapt them to nature-inspired coastal defence systems, it is necessary to incorporate the interaction between waves and vegetation. This section examines the numerical techniques used to incorporate the effect of vegetation into wave propagation models.

To include the effect of vegetation into the wave propagation models, it is crucial to examine the interaction within vegetation systems at a microscopic scale and understand how energy dissipates. Considering an individual vegetation stem situated along the wave propagation line, it creates a drag resistance, consequently extracting energy from the mean flow. This energy extraction is directly linked to the work performed on the fluid by the vegetation [63]. Under the assumption of non-breaking waves, the velocity profile of the fluid can be calculated (e.g., analytically with Stokes theory [41] or numerically), allowing the force vector acting at the fluid–vegetation interface to be estimated.

In more turbulent regimes, such as breaking waves or those with a very high Reynolds (Re) number, the same energy dissipation through turbulence applies. However, calculating fluid kinematics and energy transfer through analytical methods becomes increasingly challenging in such environments (e.g., conditions situated within breaking limit on Le Méhauté diagram [64]). Therefore, the use of numerical models to simulate the wave–vegetation interaction becomes necessary.

3.2.1. Friction Approach

An established relationship is present between wave damping and bottom friction in intermediate to shallow water conditions [65]. By benefiting from this correlation, early attempts recorded in the literature resorted to tuning the bottom friction parameter in numerical models to account for the vegetation-induced energy dissipation [66]. This factor can be tuned to enhance the boundary shear stress, impacting the velocity distribution and generating flow resistance [67,68]. In adopting this approach, one of the challenges is the calibration of the bottom friction parameter. Empirical relationships that guide the determination of the friction parameter are documented in the literature; they often include fitting parameters related to hydraulic conditions and type of sediment [69,70]. Nowacki et al. [71] demonstrate the applicability of this method by modelling the wave attenuation in a shallow water lagoon using the bottom friction equation implemented in SWAN [44], formulated as follows:

$$S_{ds,b}=-C_b\frac{{\sigma }^2}{{\mathbf{g}}^2{sinh}^{2}kh}E$$

(12)

where $S_{ds,b}$ is the bottom friction term, added to the total source term S as presented in Equation (2). Here, $C_b=C_f\mathbf{g}{u}_{rms}$, with the friction coefficient $C_f$ being the parameter used to account for the effect of vegetation.

Operating under the same generalized concept, the porous medium approach is another technique to simulate wave attenuation due to vegetation [72,73]. Here, a designated volume (area with vegetation) within the fluid domain is represented as a continuum with interconnected pores that impede the mean flow [74]. By adjusting the obstruction parameter (porosity value), the energy dissipation in the system can be simulated. Empirical correlations that are based on the vegetation characteristics, such as length, leaf width, and density, are available in the literature [72]. Nevertheless, these parameters need to be tuned using experimental or field data.

3.2.2. Cylinder Approach

Dalrymple et al. [75] document the earliest attempt to incorporate vegetation-induced energy loss through adopting an analytical damping for localized energy dissipation. This formulation was restricted to monochromatic linear waves and assumed a cluster of cylinders to represent the vegetation. The attenuation can then be fitted using the geometric properties of the cylinders and a drag coefficient (C_D) to estimate the dissipation due to the resulting force [75]. This force can be calculated using the Morison equation [76], expressed as

$$F=\frac{1}{2}\rho C_{D}D\left|u\right|u+\rho C_{M}A\frac{du}{dt}$$

(13)

where $\rho$ is the water density, D is the diameter, u is the horizontal velocity, and A is the area. The force F is considered to be the sum of a drag and an inertia component. The drag component, the first term in Equation (13), is the resistance force from friction and separation around the structure. It is calculated using a drag coefficient $C_D$. The inertia component, the second term in Equation (13), depends on the acceleration and inertia of the water displaced around the structure. It is calculated using an inertia coefficient $C_M$. The choice of the drag and inertia coefficients $C_D$ and $C_M$ can be guided through the Reynolds number $Re$ and the Keulegan–Carpenter number $KC$.

This technique, known as the cylinder approach, was widely adapted for its computational efficiency and versatility. By changing the geometric properties and drag parameters it could account for different species of seagrasses or mangroves [77,78]. The cylinder approach was extended to account for vegetation’s swaying motion by assuming small-amplitude vibrations in Asano et al. [79]. Afterwards, Méndez et al. [80] used potential flow and eigenfunction expansion to improve the formulation, also taking into account the hydrodynamic effects in the vicinity of the vegetation and validating the formulation for irregular waves. The cylinder approach continued to mature, and in 2004, Méndez and Losada [81] proposed an empirical formulation that is valid over sloping bathymetries and in areas of wave breaking.

The empirical formulations based on the cylinder approach presented above are all developed from experimental datasets in controlled wave flume environments. Up-scaling and application to larger domains becomes possible when integrating them into numerical models. Suzuki et al. [82] presented a cylinder approach-based vegetation dissipation model within the phase-averaging SWAN [44]. In this model, a spectral term representing the Méndez and Losada [81] formulation was added to the energy conservation term. With a horizontal and vertical layer-specific dissipation, this SWAN model could also account for spatial variations in vegetation characteristics. Several researchers [83,84] followed the same description to investigate the wave attenuation services of nature-inspired solutions.

A similar approach, through adding an energy dissipation term into the momentum equation, can be adopted for phase-resolving models. In this case, the force acting on the cylinder is calculated in the time domain using the Morison equation [76]. Additional terms that take into account the inertia effects are also included. For example, this implementation was presented for SWASH [85] and mild-slope models [86].

3.2.3. Coupled Approach

When comparing the friction and cylinder approaches, the latter is more representative of the physical processes in the wave–vegetation interaction problem. Firstly, it explicitly takes into account the geometrical characteristics of the vegetation, including diameter, length, and density. Secondly, the use of cylindrical mimics to represent vegetation has become common practice in experimental setups (e.g., [87,88,89]). Thirdly, there is a growing interest in incorporating the swaying motion through flexible structural solvers, since the rigid assumption tends to overestimate wave damping [90]. For these reasons, numerical platforms that represent the vegetation using a cylinder (alternative geometric shapes, such as rectangular blades, can also be utilized) and include a communication interface between waves (fluid solver) and vegetation (structural solver), adopting the coupled approach, have come to prominence (e.g., [91,92]). A complete form of the coupled approach explicitly resolves all dimensions of the wave–vegetation interaction: (1) swaying motion, (2) flow characteristics, and (3) energy transfer mechanisms. However, the numerical treatment of all three dimensions requires large computational resources, leading to cases where the coupled approach is selectively applied to address only one (or more) aspects of the interaction.

In coupled models, a common approach to modelling the swaying movement of vegetation is to utilize a structural model based on a cantilever beam configuration [93,94]. The bending characteristics of the vegetation, referred to as flexural rigidity, are defined through the section geometrical and material properties. In addressing a single dimension of the wave–vegetation interaction (i.e., (1) swaying motion), researchers frequently employ one-way coupled models [63,95,96], wherein the forcing data are communicated from the flow solver to the structural solver for the computation of deflections. This type of one-way coupling does not provide feedback from the vegetation on the fluid. In the majority of cases, the force balance on the vegetation is calculated using the Morison equation [76] in order to account for the drag and lift forces. Additional terms such as skin friction, buoyancy, gravity, and added mass are also added to the system.

Expanding on the previously outlined approach, two-way coupled models have been developed to examine more than one dimension of the wave–vegetation interaction. Typically, these models combine flow solvers for the wave propagation with structural solvers for the deflection. For example, Yin et al. [97] utilized the XBeach [98] non-hydrostatic model to convey both the depth-averaged velocity in one direction and the dynamic force and swaying velocity in the opposite direction. The influence of these parameters on the wave dissipation is then incorporated into the momentum equation through a dissipation term. This is achieved using an extended version of the cylinder approach [75], whereby the fluid velocity is replaced by the relative velocity between the fluid and the vegetation. A similar strategy has been also adopted in SWASH [99]. For Navier–Stokes models Maza et al. [100] presented a similar form of this coupling.

The previously presented coupled models address the swaying motion and the flow characteristics in the system. However, the energy transfer mechanisms and dissipation in the form of vortex shedding are not considered, and the wave dissipation is implicitly captured using a modified form of the cylinder approach. Solving the full Navier–Stokes equations coupled with a structural solver can resolve the fluid flow and fluid–structure interaction explicitly. This approach allows for the capture of the force transfer without resorting to the assumptions of the cylinder approach, which depend on fitting parameters such as the drag coefficient. Furthermore, it directly computes the loss of energy in the system without the necessity of a modification in the governing momentum or energy balance equations. This explicit capturing of energy comes at a computational cost. El Rahi et al. [91] utilized the meshless smoothed-particle hydrodynamics (SPH) code DualSPHysics, coupled with a finite element analysis (FEA) structural solver, to address vegetation dynamics. The method explicitly resolves all the dimensions of the wave–vegetation interaction and has no dependency on any fitting parameters of empirical formulations. However, the effectiveness of this method may vary depending on specific conditions as it employs an embedded approach to represent vegetation within a larger dummy boundary envelope. It is particularly suited to scenarios where buoyancy-restoring forces have minimal impact and the vegetation has a rectangular cross-section (limited to one geometric shape) with a thickness significantly smaller than the wavelength of the oscillatory water waves. In another coupled SPH model, Paquier et al. [101] addressed the wave-flexible vegetation problem using a two-way coupling between a fluid solver in combination with the slender rod theory. However, constrained by the computational cost, they adjusted the dimensions of the vegetation within the numerical model while modifying the section properties. Other examples of coupled models that have been developed to simulate all the turbulence around flexible elements in fluid environments are presented in Tschisgale and Fröhlich [102], Girfoglio et al. [103], and El Rahi et al. [104].

4. Coastal and Marine Vegetation

Vegetation in coastal areas can be grouped according to species, characteristics, and habitat. This classification yields two groups that are relevant to the subject of wave dissipation: (1) emergent rigid vegetation, such as mangroves; and (2) flexible vegetation, such as seagrasses. Figure 5 illustrates the distribution of both groups of vegetation types. Flexible vegetation, such as seagrasses, occupies the submerged region and extends into the emergent zone. Rigid vegetation, such as mangroves, predominantly occupies the emergent zone.

4.1. Mangroves: Emergent Rigid Vegetation

Mangroves are terrestrial plants that have adapted to saline waters (halophytic) and occur in intertidal zones of tropical and subtropical regions [105]. Their leaves require exposure to air, which confines their habitat to shallow intertidal waters, allowing them to emerge during all tidal cycles [106]. Due to their dependence on fresh water, mangroves typically inhabit wetlands with abundant freshwater sources. This is a feature of estuaries and river deltas, where inland fresh water discharges into the sea. These areas are particularly favourable for mangrove growth as the freshwater inflows bring a rich abundance of nutrients. However, certain species, belonging to the genera Avicennia and Ceriops, exhibit slower growth rates when exposed to freshwater and are notable exceptions to this pattern [107]. Furthermore, permeable soils with porous characteristics are essential for proper drainage and nutrient supply. Field studies conducted in various locations, including Northern Australia [108] and Japan [109], have demonstrated the significant impact of soil characteristics on mangrove growth (e.g., Avicennia marina exhibits better growth in soils with high pH levels [109]). Typically, the soil is composed of well-sorted silts mixed with other organic materials. The largest mangrove habitats are located in South and Southeast Asia, spanning countries along the shorelines of the Indian Ocean such as Indonesia, as well as regions along the great barrier reef in Australia. Following that, the Atlantic shorelines of Brazil in South America and Mexico in Central and North America host significant populations [110].

Mangroves are an integral part of the ecosystem, acting as a nursery for a vast array of species, including oysters and sponges [111]. They also provide a nursery for juvenile fish species. Sasekumar et al. [112] documented more than 100 species of fish and around 10 species of prawns spending at least part of their life cycle in the mangrove forests of Selangor, Malaysia. Additionally, mangrove ecosystems provide nesting grounds and shelter for birds and terrestrial animals. This includes, but is not limited to, bird species hosted by the Queensland mangroves in Australia [113], and crocodiles (Crocodylus porosus) in Sri Lanka [114]. Beneath the sediment bed, the roots of mangroves provide a habitat for many benthic organisms (e.g., gastropods [115]). These organisms significantly contribute to the ecosystem by performing multiple functions. For instance, they connect mangrove organic matter to the broader food web by breaking down and processing fallen mangrove leaves [116]. Mangroves also deliver economic and ecosystem services to local communities and countries. Economically, they are harvested for their wood [117] and serve as hotspots for eco-tourism [118]. Furthermore, their coastal protection services, which involve attenuating wave energy and preserving shorelines, have been extensively documented in the literature [32,119,120].

Mangroves are characterized by a woody trunk which supports a crown of branches and leaves [111]. They are anchored into the sediment bed by an intertwining root network that emerges sparsely from the soil for the uptake of air [111]. There are more than 70 documented species of mangroves, of which two are predominantly present in mangrove habitats and provide coastal protection services [121]. The two most common mangrove species are Rhizophora mangle, commonly known as red mangroves, and Avicennia germinans, known as black mangroves [122].

Red mangroves, Rhizophora mangle, are among the most prevalent species that grow in muddy areas with access to fresh or highly diluted waters [123]. They typically thrive in low-lying regions of swamps where brackish waters (i.e., with reduced salinity) are present for at least part of the year. Their occurrence is globally distributed and includes but is not limited to the shorelines of Mauritania in South Africa, the Atlantic coast of Florida, and the Gulf of Mexico in northern America [121]. They are characterized by a tall emerging root system, known as prop roots [123]. These prop roots have a dual function: above the surface exchanging air and below the surface absorbing nutrients [123]. This species is particularly effective in providing coastal protection due to its tall and rigid structure, which can reach up to 30 m in the tropics (as observed in Puerto Rico [124]).

Black mangroves, Avicennia germinans, exhibit a preference for drier habitats that are not fully submerged in water and can tolerate hypersaline conditions [125]. They are recognizable by their numerous rigid emerging roots, which resemble columns and are known as pneumatophores [126]. The black mangrove structure is rigid and reaches heights of up to 15 m, and it also provides coastal protection services [125].

4.2. Salt Marshes and Seagrasses: Flexible Vegetation

The group of flexible vegetation can be divided into two main types of plants. The first type, salt marshes, includes marine halophytes which are found in intertidal zones. The second type consists of seagrasses, which are flowering marine plants that thrive in shallow water zones.

The term salt marsh refers to salt-tolerant terrestrial plants that thrive in the upper intertidal zone [127]. These plants, including flexible, non-woody species such as grasses and herbs, require exposure to air [127]. They typically grow in areas with muddy substrates and soft sediment beds, forming in coastlines and estuaries [128]. The development of salt marshes is characteristic of sheltered conditions as they can recede quickly under the force of energetic erosive waves. They are distributed across all latitudes; however, they predominantly flourish in temperate zones and are mostly replaced by mangrove ecosystems in tropical zones [129]. Mcowen et al. [130] estimated the global distribution of salt marshes, and identified the largest areas in North and Central America followed by Australia. As halophytes, their distribution extends from mean sea level to the high water level, allowing them exposure to air during the tidal cycle [131].

Within salt marsh ecosystems, plant species diversity is considerable and is primarily influenced by competition at higher elevations and physiological tolerances at lower elevations. For example, in the Tagus estuary of Portugal, Simas et al. [132] observed a correlation between species distribution and salinity. Spartina is the dominant genus in low-salinity areas, while Halimione portulacoides and Arthronemum frauticosum dominate in higher-salinity areas.

Salt marshes are valued as nurseries for fish and other marine organisms [133]. For instance, in a single salt marsh estuary in Georgia, USA, Gorden Rogers et al. [134] identified multiple fish species including the Atlantic croaker (Micropogonias undulatus), the silver perch (Bidyanus bidyanus), and the southern flounder (Paralichthys lethostigma). Furthermore, along the French Atlantic coast, the nursery function of salt marshes was identified for the European sea bass (Dicentrarchus labrax) [135]. Salt marshes are also active in filtering heavy (toxic) metals from estuarine circulation [136]. Furthermore, they are recognized for their role in acting as a protective barrier between coastal communities and extreme weather conditions, such as waves and storm surges (e.g., [137,138]).

Seagrasses are marine flowering plants that have a quasi-global presence in the form of seagrass meadows in coastal waters worldwide [36]. Seagrasses are constituted by a group of approximately 60 different species grouped into two families, Patamogetonaceae and Hydrocharitaceae, dominated by three main genera: Halophila, Zostera, and Posidonia [36]. Typically, populations grow into monospecific meadows constituted of a single species. For example, the Mediterranean coasts of Spain, France, Italy, Greece, and Tunisia host Posidonia oceanica meadows [139], while Zostera noltii is prevalent in North Atlantic and Baltic cities [140]. In tropical waters, meadows are reportedly more diverse [141].

While most seagrass species are able to tolerate a wide range of salinity levels, only a select few can withstand air exposure [36]. Species capable of withstanding air exposure, such as Zostera noltii meadows on the west, east, and northwest African coastlines, can establish intertidal meadows that extend from shallow littoral waters to the high tidal marks along the coast [142]. With regard to the depth at which they can survive, all seagrass species are constrained to the boundaries of the light extinction depth, where sunlight penetrates sufficiently for photosynthesis [143]. The water depth limit correlates with the physiological adaptations of each species and the water quality of the hosting environment. For example, Posidonia oceanica and Halophila depiciens can grow in depths of up to 40 m, whereas Cymodocea serrulata and Syringodium isoetifolium are typically found in depths of around 10 m [143]. With regard to their global distribution, the southern temperate oceans are dominated by species belonging to the genus Zostera, the tropical Atlantic by Thalassia, and the clear Mediterranean waters by Posidonia [144]. Other species, such as Ruppia maritima, exhibit a global distribution across both the Northern and Southern Hemispheres [144].

Seagrass meadows support the hosted coastal ecosystem through multiple ecosystem services [145]. Although few marine species feed directly on seagrasses due to their low nutritional value and their being difficult to digest, green turtles (Chelonia mydas) are known to graze on seagrasses [146]. Futhermore, seagrass meadows serve as nursery and feeding grounds for multiple juvenile fish species, such as herring (Clupea harengus) and pollack (Pollachius pollachius) [37], by supporting a rich epiphytic fauna and flora on their leaves [147]. Additionally, in the sediment of the seagrass bed a large diversity of benthic organisms can serve as a food source. Furthermore, seagrasses develop roots and rhizomes in the sediment, stabilizing the sediment bed and preventing erosion [148].

5. Ecological Framework

The simplification of the wave–vegetation interaction process is necessary to represent the coastal vegetation within numerical wave propagation models. Numerous numerical techniques (described in Section 3.2), like the implicit bottom friction and explicit CFD direct modelling of energy transfer, are available for this purpose. However, these numerical techniques often prioritize interaction dynamics and neglect the ecological and physiological features of the vegetation. This section explores vegetation characteristics which are relevant to wave attenuation. These features are organized within the “ecological framework”, which is illustrated by a conceptual diagram in Figure 6.

5.1. Plant Morphology

Plant morphology is the study of the geometric characteristics and external structure of plants [149]. One key morphological trait that is relevant to wave–vegetation interaction is length. For seagrasses and saltmarsh vegetation, leaf length is measured from the sediment bed to the tip of the leaf. On average, the typical length of a leaf is approximately 20 cm. For example, the mean leaf length for some common North American species is around 17 cm for Halodule wrightii, 41 cm for Thalassia testudinum, and 23 cm for Zostera marina [150]. For the common Mediterranean Posidonia oceanica species, the leaf length averages about 25 cm [151]. Mangroves on the other hand, typically measure on the order of meters (height is the common term for mangroves rather than length). Several species, including Avicennia marina, Ceriops tagal, and Bruguiera gymnorhiza, have been documented to reach heights ranging from 2.5 to 20 m, depending strongly on environmental conditions [152,153]. In rigid mangrove vegetation, tree height and wave-interacting frontal area remain relatively constant. However, in flexible seagrasses and salt marsh vegetation, canopy height can vary. Canopy height is defined as the height in the neutral position. When bending due to gravity, the canopy height is often less than the extended length.

Local populations of seagrasses and saltmarsh plants form meadows, while mangroves form forests. The population characteristics are described by (shoot) density, which refers to the number of plants in a given unit area. Typically, (shoot) density is lower for mangroves. For instance, in the coastal areas of Indonesia, the maximum density is approximately 1200 trees per ha, and the minimum is around 100 trees per ha, for Rhizophora apiculata and Bruguiera cylindrica, respectively [154].

Seagrasses, on the other hand, often exhibit much higher densities that strongly depend on the species. Pioneer species colonizing bare sediment have relatively low densities, whereas climax vegetation with stable, long-term presence has much higher densities and is more effective in wave attenuation. For instance, along the Spanish Mediterranean coast, Posidonia oceanica has been reported with a shoot density ranging from 150 to 1500 shoots per m² [155]. Along the Atlantic coast of Canada, occurrences of Zostera marina have been observed with shoot densities of up to 1000 shoots per m² [156].

A strong correlation between plant morphology, in particular shoot density, and wave attenuation has been established in the literature [157,158]. For example, Stratigaki et al. [88] demonstrated that the wave height reduction increases almost two fold when the shoot density is doubled. In separate research, focusing on pollard willow trees, Kalloe et al. [159] established a correlation between density and wave attenuation through the development of a novel density estimation technique based on frontal surface area.

5.2. Biomechanics

In flexible vegetation, the plants exhibit a swaying motion under the action of the waves [137]. This deflection and elastic reconfiguration observed in swaying plants is known as biomechanics [160]. This dynamic problem can be solved by treating the plant as a cantilever structure with a simplified geometry. In this case, the root system acts as the fixed end, while the tip of the plant as the free oscillating end (e.g., [91,99,101,161]). The bending resistance is then described by the modulus of elasticity (E) and the cross-sectional area’s moment of inertia (I).

In order to resolve the dynamic response within structural solvers, assumptions are made regarding uniform section properties and internal compositions. However, given the complex and intertwined root and shoot systems in plants, the assumption of a uniform flexible structure is a significant simplification. In addition, plants can be composed of different cell types with diverse chemical compositions. Furthermore, the epiphytic overgrowth of seagrass plants by epiflora (and to a lesser extent epifauna) alters the hydrodynamic interactions with the flow (increased drag) and the bending stiffness (i.e., modification of section properties) of the vegetation [162]. Therefore, biomechanical properties are always idealized and simplified in engineering applications.

An illustration of the required simplifications to study biomechanics from an engineering perspective can be presented for Posidonia oceanica. The modulus of elasticity (E) of natural Posidonia oceanica plants is reported to be in the range of 4.7 GPa [163]. However, within the engineering studies of Folklard [163] and Stratigaki et al. [88], the mimics with a simplified geometry had a modulus of elasticity (E) of 0.51 GPa and 0.903 GPa, respectively. Other researchers selected an appropriate value of E that ensures a swaying motion, but without a direct link to a species. For example Luhar and Nepf [95] used an HDPE blade with E equals 0.93 GPa and Reis et al. [87] used a rubber cylinder with E equals 0.82 × 10⁻³ GPa.

5.3. Buoyancy

In submerged flexible plants, the buoyancy force plays a pivotal role in the force balance equation and acts as a stabilizing force in a neutral posture, and as a restoring force during the swaying cycle. The positive buoyancy in seagrasses is attributed to air-filled cavities, referred to as lacunae [164]. The magnitude of the buoyancy force depends on the density gradient between the fluid environment and the vegetation. For example, in the case of Posidonia oceanica, which has a reported density in the range of 900 kgm⁻³ [163], the positive buoyancy-restoring force is relatively small. In contrast, Zostera marina is reported to have a density, at the low end, in the range of 700 kgm⁻³ [165].

5.4. Variability

At the local population scale among plants, variations in properties can be observed [166]. Shoot densities are affected by seasonality, with a peak during the warmer months and a minimum during the colder winter months. For example, the shoot density of Zostera noltii in the Netherlands is reported to increase during spring due to warmer temperatures and increased light availability [167]. Following that, a decline in biomass is recorded due to unfavourable environmental conditions and grazing. This pattern is also recorded within seagrass species in temperate regions, for example, for Thalassia testudinum in the South Florida region [168]. Regarding mangroves, seasonal variation, quantified through wood production, is also documented for multiple mangroves including Avicennia germinans and Rhizophora mangle [169]. Furthermore, other factors such as salinity also have an influence on the growth of mangroves [170].

At the single-plant scale, seasonal changes also impact plant biomechanical behaviour, with higher flexural rigidity observed during warmer months compared to colder winter months. In a study conducted on the salt marsh species Spartina anglica, exposure to a warmer water temperature induced an increase in flexural rigidity [171]. Additionally, the lifecycle of plants plays a role, with plants tending to be stiffer during the maturing phase and becoming more flexible during ageing [166]. To add an additional layer of complexity, all of that natural variability can be potentially impacted by environmental changes due to climate change [172].

6. Results: Categorization of Numerical Models

Numerical models addressing wave–vegetation interaction are identified through a literature review of peer-reviewed journal articles. Subsequently, these models are grouped into five categories based on the ecological framework. Each category and its properties are presented in the classification matrix shown in

Table 1. Classification matrix of numerical models according to the ecological framework.

	Morphology	Biomechanics	Buoyancy	Variability
Category A	No	No	No	No
Category B	Yes	No	No	No
Category C	Yes	Yes	No	No
Category D	Yes	Yes	Yes	No
Category E	Yes	Yes	Yes	Yes

.

The following standards are employed in identifying the applicability of the individual ecological features. An ecological feature is considered included if at least one characteristic is present. For morphology, this means vegetation length or shoot density. For biomechanics, it is flexural rigidity. For buoyancy, it is the buoyant force. For variability, it includes spatial variation in shoot density, seasonal variation, or variation in flexural rigidity. Furthermore, buoyancy is considered included only if it is represented as a restoring force, which is part of the force balance equation during the swaying cycle in flexible vegetation.

Numerical models classified under “category A” disregard all ecological features of the vegetation. Instead, these models focus on wave dissipation, achieved by subtracting energy from the system based on either tuning a friction parameter against experimental data or using empirical formulations to guide the selection. For instance, Blackmar et al. [173] used a phase-resolving numerical model validated against experimental data to capture wave attenuation. The experiments were conducted using scaled mimics comparable to Schoenoplectus pungens arranged in multiple configurations and shoot densities. In this study, the numerical model implicitly simulated the wave attenuation without any consideration of the features of the vegetation. Attenuation was achieved by calibrating a bottom-friction factor against experimentally recorded surface elevations. Similarly, Augustin et al. [68] employed an iterative process to determine the friction factor for waves interacting with mimics comparable to Spartina alterniflora. In a separate study, validated against field observations, an enhanced bottom-friction approach within a phase-averaging model was employed to simulate wave attenuation in coastal saltmarshes in Chesapeake Bay [174].

In “category B”, numerical models address morphology by incorporating the geometric features of vegetation. These models employ the cylinder approach to simulate drag forces using the Morison formulation [76], which requires input parameters such as the length and diameter of the vegetation. Subsequently, wave attenuation is modelled as a dissipative sink of momentum or energy. For example, Li and Yan [175] calculated wave dissipation in a Navier–Stokes model by introducing an additional source term to account for vegetation-induced turbulence. This term, derived from the cylinder approach, incorporates morphological variables such as width, thickness, length, and shoot density. A similar method for drag calculation and momentum dissipation is described by Marsooli and Wu [176]. A number of studies have described the implementation of a similar cylinder approach in various wave propagation models, including mild-slope models [177], non-hydrostatic models [85], and phase-averaging models [82].

By treating vegetation as flexible elements swaying under the action of waves, numerical models classified under “category C” address both morphology and biomechanics. These models employ a coupled approach that combines a fluid solver with a structural solver. There are two types of such models.

The first type computes drag forces using the cylinder approach, and then, communicates the hydrodynamic loading to the structural solver. For example, Marjoribanks et al. [178] coupled a Navier–Stokes model with a structural solver. In this approach, a simplification of the fluid forces was assumed, only including drag forces to compute deflections. The deflected shape was then communicated back to the fluid solver, and a mass blockage was assigned to the corresponding computational grid cells to dissipate energy. Mendez and Losada [80] tackled the swaying motion similarly, considering only the horizontal hydrodynamic forces.

The second type of coupled models abandons the cylinder approach for computing force balance and directly uses Navier–Stokes equations. However, these models still assume simplifications, not in force transfer but in the geometric properties of the vegetation. For example, El Rahi et al. [91] used an embedded method to represent thin vegetation within a larger dummy boundary, while Paquier et al. [101] increased the thickness of the vegetation and decreased the flexural rigidity to maintain flexibility. Both techniques were aimed at maintaining achievable computational runtime at the expense of simplified hydrodynamic forcing calculations. These calculations do not take into account buoyancy and hydrostatic forces.

Under specific hydrodynamic conditions during the swaying cycle, where the restoring force due to buoyancy is more dominant than the restoring force due to stiffness, a complete representation of the forcing environment becomes necessary. Numerical models under “category D” address this by including morphology, biomechanics, and buoyancy. At the scale of an individual vegetation stem, Luhar and Nepf [95,179] presented a complete description of the forces using a coupled approach based on an extended form of the Morison equation, which includes drag, inertia, added mass, skin friction, and buoyancy. This approach, validated against experimental data, was also adopted in Zeller et al. [63] and Zhu et al. [96]. At the scale of a vegetation field represented by an array of flexible mimics, van Veelen et al. [93] coupled the previously described formulation with a non-hydrostatic model to account for wave damping through energy dissipation. At a regional scale, this has also been implemented and validated within phase-resolving models [180].

Models classified under “category E” extend the numerical treatment of vegetation by incorporating variability among the other ecological features. Marjoribanks and Paul [181] emphasized the significance of this variability by proposing a numerical approach that could accommodate variable thickness and stiffness. The model was developed with the specific aim of simulating the behaviour of salt marsh vegetation, which is characterized by a high degree of thickness and stiffness in the vicinity of the sediment bed, and a reduction in these properties as it extends towards the stem. This was achieved by configuring a structural model that discretizes individual stems into discrete elements connected by springs. In order to address the biomechanical behaviour of the system, the local flow conditions were used to calculate the hydrodynamic forces, including drag, inertia, and buoyancy. The model was validated against experimental data of waves interacting with flexible mimics representing two species: Spartina alterniflora and Spartina anglica.

In a separate study, Wu et al. [182] addressed the variability in the vertical biomass distribution of vegetation, such as Schoenoplectus pungens. They simulated wave attenuation over an array of non-uniform vegetation using two techniques. In a phase-averaging model, a friction approach was employed to support variation in vertical density by discretizing the domain into multiple vertical layers and applying a variable stem density at each layer. In a phase-resolving model, the cylinder approach was extended to include a factor representing variation in biomass in the vertical layers.

7. Discussion and Future Research Directions

A growing body of evidence suggests that wave attenuation is significantly influenced by the ecological features of the vegetation. For instance, research indicates that stiff salt marsh vegetation is more effective at attenuating waves compared to flexible seagrasses, which tend to become more streamlined by reconfiguring under wave action [183]. Furthermore, shoot density and height, both subject to seasonal variability, have been correlated with the potential for wave attenuation [184].

Numerical models across all categories acknowledge the influence of ecological factors on wave attenuation. However, due to limitations in computational resources, certain model assumptions are necessary.

Recall that implicit methods implemented within “category A” include the most simplifications, yet they are computationally efficient and well suited for applications on a regional scale and for extended periods (e.g., wave attenuation over salt marshes, North Norfolk, England [185]). Recall also that explicit methods that incorporate multiple ecological features (e.g., categories C, D, and E), are more computationally demanding and are applied within smaller temporal and spatial domains, at scales ranging from a single plant [91] to an array of vegetation stems under the action of regular waves [93]. Additionally, consider the distinction between numerical models that are limited to vegetation dynamics and those that also resolve the wave attenuation. Models that focus on vegetation dynamics are employed to investigate reconfiguration patterns and understand the correlations between hydrodynamics and drag. Nevertheless, the numerical interaction between the fluid and the vegetation is one-directional, meaning that the resulting wave attenuation is not captured (e.g., [63,96,179]). Comparison of numerical simulations against experiments is based on two markers: the horizontal force component and the swaying distance. The numerical models can replicate the experiments with varying degrees of accuracy. For example, El Rahi et al. [91] reported a root mean square error (RMSE) ranging from 0.006 N to 0.020 N for the horizontal force component. Additionally, the swaying distance was numerically captured with an error of 0.001 m. Similarly, Luhar and Nepf [95] presented numerical results that fall within the accuracy margin of the experimental measurement setup, which is 0.001 N (considering a 10% accuracy inherent to the load cell).

With regard to morphology, there is a lack of numerical models that accurately represent the complex shapes of vegetation, including intricate intertwining shoots, branches, and leaves. To simplify these complexities, basic geometric shapes, such as cylinders and blades, are typically used. These simplifications facilitate the quantification of hydrodynamic forces and the simulation of system dynamics. Some studies have attempted to compute energy dissipation around vegetation directly; however, this approach requires an enormous number of computational nodes and has been limited to validation cases involving a few waves interacting with simplified rigid cylinders. For example, researchers have employed direct numerical modelling to simulate (i) tsunami interaction with mangroves represented by rigid cylinders [186], and (ii) energy dissipation through vortex shedding from an isolated rigid cylinder [187].

By moving beyond the simplified cylindrical representation, a comprehensive description of mangrove morphology in numerical models has been explored by our research group (Coastal Engineering Research Group (CERG), Ghent University). For example, Figure 7 illustrates flow kinematics within a fully discretized willow tree. The morphological characteristics of the willow tree are extracted from an STL file. Regular waves are generated, with the velocities within the vegetated zone represented using a colour map. Although this model is still limited by the number of waves and lacks a validation case, it provides a promising indication of how future models might be developed.

Concerning variability, particularly under “category E”, there is a scarcity in numerical models that address this ecological feature. Among the few models identified in the literature, one is limited to variability in flexural rigidity [178] and the other to variability in density along the vegetation length [182]. These models can handle morphological variability within the canopy and better describe the biomechanics. Nevertheless, future numerical models should also address seasonality in order to provide a more comprehensive understanding.

With regard to applicability to nature-inclusive coastal defence systems, some models within “category A” are used to simulate realistic conditions and address specific species (e.g., Spartina alterniflora in Chesapeake Bay [174]). Nevertheless, the majority of the models in other categories replicate experimental setups with some consideration of vegetation species but little consideration to the geographic zone, which is also crucial for deriving the hydrodynamic conditions. Note how the idealized laboratory experiments with regular waves presented in Luhar and Nepf [95] were used to validate a large number of models (e.g., [91,96,101,188]).

Furthermore, Temmerman et al. [27] emphasized the critical role of natural resilience in maintaining vegetation ecosystem services in the face of stress induced by storms. Additionally, they introduced the concept of recovery time following storm damage and investigated the impact of anthropogenic stressors, such as climate change, on vegetation functionality. These processes remain complex and currently beyond the scope of numerical models.

This study acknowledges that the use of computationally efficient models operating at large spatio-temporal scales often necessitates simplifications in order to address wave–vegetation interaction. These models incorporate implicit formulations for energy dissipation and include no or few ecological attributes. However, it is important to clarify that these simple models do not necessarily yield results that are divergent from those of more complex models that include a greater number of ecological attributes. Simple models are designed to operate over large domains for extended durations (i.e., irregular sea states), while complex models are typically applied to smaller domains and for shorter time frames (i.e., few regular waves). Consequently, due to the distinct spatio-temporal scales at which these models operate, definitive comparisons of the results for wave propagation over vegetation are not appropriate. Nevertheless, in ideal computational scenarios, it is reasonable to assume that more complex models yield more accurate results for two reasons. Firstly, the energy transfer at the fluid–vegetation interface is directly solved and the physics are captured. Secondly, complex models are able to represent reality by including more ecological features.

Interest in the numerical modelling of wave–vegetation interaction is expected to grow significantly in the near future. Therefore, understanding the functions and limitations of models within each category is crucial for drawing valid conclusions and establishing effective links across the various categories and scales. As more complex models emerge, it is anticipated that the number of categories will converge into two groups: (1) computationally efficient, large-domain implicit models; and (2) computationally expensive, small-domain explicit direct models.

In future research, explicit direct models should pay greater attention to seasonal variability in properties and account for the complexities in the morphological description of vegetation. This approach will enable a beneficial coupling between the two models, where the high-resolution model, running for a limited duration under design conditions, can define the parameters for the larger model operating at the ecosystem scale.

Lastly, with the increased computational power provided by high-performance computing (HPC) and the rise of computational applications on highly parallelized graphical processing units (GPUs), computational limitations will change in the near future [189,190]. For wave–vegetation interaction numerical models, this implies that the boundaries for high-resolution and low-resolution models are dynamic. Currently, high-resolution models are applied on a stem scale, but future applications could extend to the meadow scale or even the ecosystem scale. Additionally, while regular waves are currently employed in these models, realistic climate scenarios could be used in the near future to run high-resolution models for extended periods.

8. Conclusions

This study examined the implementation of nature-inclusive coastal defence systems from both engineering and ecological perspectives. From an engineering perspective, wave propagation models and wave–vegetation interaction techniques were identified. From an ecological perspective, types of coastal and marine vegetation relevant to nature-inspired coastal defence systems were presented. An ecological framework was developed based on the vegetation features. This framework includes five categories, which refer to morphology, biomechanics, buoyancy, and variability.

The results of the categorization of the numerical models challenges the hypothesis that the development of NBS in fundamental research is divided between ecologists and engineers. Evidence from published numerical models show that coastal engineers acknowledge the importance of incorporating the ecological features of vegetation while investigating wave attenuation. Furthermore, multiple wave–vegetation interaction models accommodate some ecological features, if not all of them. However, the complexity and the expensive computational costs often necessitate simplifications.

The analysis of the results also revealed that the complexity of the numerical models is inversely proportional to their ability to solve engineering problems. Complex models are generally limited to validation cases based on idealized laboratory studies, while simpler models address practical engineering challenges.

Moreover, there is little evidence indicating that future advancements will result in the development of complex numerical models that are able to cover all five ecological features and directly compute the energy dissipation, all while operating at a large spatio-temporal scale. Therefore, a mutually beneficial link between high-resolution small-scale models and large-scale implicit models could be a promising avenue for future research.

In order to advance the fundamental understanding of wave–vegetation interactions and the implementation of nature-inspired coastal defence systems, it is important to give equal importance to both simple and complex models, since their applicability depends on the conditions and the specific research question being addressed. Furthermore, it is necessary to expand on the current models by describing emergent processes including stressors such as climate change, and additional ecological features such as survivability and recovery from storm damage. Finally, this study aimed to explore wave–vegetation interaction; however, future research could be oriented to embrace additional processes such as morphodynamics and ecosystem services.