Two approaches are used for quantifying wave height attenuation in the Jicaral mangrove forest: the wave height attenuation rate ( r ) and the drag force approach, which are widely used in numerical models, including the Simulating Waves nearshore (SWAN) model (Booij et al., 1999). The first approach is used to evaluate the correlation between water depth ( h ) and H m 0 with the r values for the mudflat zone and the mangrove zone. In contrast, the drag force approach is used to obtain an adequate estimation of C ˜ D in the field, which will be useful for later numerical models based on this approach (Mendez and Losada, 2004).

Figure 4 shows H m 0 , T p , h and relative wave height ( H m 0 / h ) values for the bursts in which all the sensors are inundated and H m 0 in S1 is greater than 0.05 m. The recorded H m 0 values during the measurement period did not exceed 0.38 m ( Figure 4A ). The T p values ranged from 1.88 to 3.76 s ( Figure 4B ), and the maximum h value recorded at sensor S1 during spring tides was 2.38 m ( Figure 4C ). The relative wave heights did not exceed 0.36 ( Figure 4D ), suggesting that wave breaking is not the main wave dissipation mechanism since wave breaking occurs when wave heights exceed 60% of h (Horstman et al., 2014; Cao et al., 2016).

Figure 5 displays the time-averaged spectral densities for the bursts corresponding to the five largest tides with h values exceeding 2 m, which are highlighted in red in Figure 4 . The selected bursts include the tides for February 21 ( Figure 5A ), February 22 ( Figure 5B ), February 23 ( Figure 5C ) and March 4 at 2:00 and 16:00, corresponding to Figures 5D and E , respectively. The observed spectra are unimodal, with wave periods ranging from approximately 2.2 to 4 s, which are typical of local waves generated by winds, suggesting no presence of infra-gravity waves. The most energetic incident waves were recorded at sensor S1 on March 4, with maximum energy densities twice as high as those recorded in February.

The results indicate that the maximum energy density decreases by an average of 20% in the mudflat region between sensors S0 and S1, except for the spring tide case on February 23, where a slight increase in energy is recorded at S1. This increase may be attributed to local effects, such as wind-induced energy. The greatest reduction occurs along the first 24.8 m of the forest between sensors S1 and S2, with the maximum energy density decreasing by 40% on average. Conversely, for the last section of the mangrove plant, between S3 and S4, the energy dissipation is the lowest, with the maximum energy density decreasing by an average of 10%. This decrease in energy density reduction is expected due to a decrease in both mangrove forest density and wave height. In general, as waves propagate through mangroves (S1–S4), there is a notable decrease in energy density. The wave height is reduced by 34% on average, similar to the value estimated by Maza et al. (2021) for a 35-year-old Rhizophora sp. forest with the same width as the Jicaral forest and a water depth equal to 2 m.

Wave height attenuation rates are estimated by using Equation 1, where r denotes the wave attenuation per linear meter. Figure 6 displays the r values obtained between sensors S1 and S4 (mangrove zone, black dots) and sensors S0 and S1 (mudflat zone, red dots) for each recorded burst as a function of the measured water level at S1. The black and red lines represent the linear fits of the data for the mangrove and mudflat zones, respectively. Within the mangrove forest, wave attenuation tends to decrease as water depth increases, as indicated by the black line. A similar trend was reported by Mazda et al. (1997a) and Mazda et al. (2006) in Kandelia candel and Sonneratia sp. forests. In addition, Figure 6 shows that the r values in the mudflat zone tend to be constant independent of water depth.

Figure 7 shows r as a function of the incident wave height. As shown, as H m 0 increases, the wave attenuation also increases, as reported by Horstman et al. (2014). In contrast, in the mudflat, there is also no clear variation trend for r with respect to H m 0 , and in some cases, there is an increase in the wave height between sensors S0 and S1, which agrees with Figure 5C .

In general, Figures 6 , 7 present r values in the mangrove field that vary between 0.001 and 0.01 m^-1, whereas those related to the mudflat zone tend to be one order of magnitude lower (10^-4) than those for the mangrove forest, with a mean value of 0.0004 m^-1. This difference suggests that the decrease in wave height along the forest is mainly attributed to the drag exerted by the trees, as demonstrated by previous studies (e.g., Mazda et al., 1997a; Horstman et al., 2014).

Adequate characterization of mangrove tree structure is essential for modeling wave-vegetation interactions and identifying mangrove tree parameters, i.e., characteristic length, that better describe wave attenuation. Hence, this study represents the vertical and spatial variation in mangrove tree characteristics by A h m , calculated from the variables N   and d r o o t . Figure 8 displays A h m values for different elevations above the ground ( h m ), where the blue line corresponds to the A h m vertical variation in Zone 1, the red line is for Zone 2 and the green line is for Zone 3. For Zone 2, the characteristics reported in Figure 2B for both SU2 and SU3 were averaged, because based on visual inspection in the field, both sample units represented areas equally distributed along the transect. Zone 1 presented the largest A h m values, which were reduced by 60% and 73% in Zone 2 and Zone 3, respectively, for h m =10 cm. Instead, as h m increases, the A h m differences between zones decrease, revealing the influence of secondary roots on low h m values.

In addition to the data collected in Jicaral, Figure 8 shows the A h m values of previous studies that have measured or parameterized mangrove geometry for Rhizophora sp., either in the field or laboratory. A detailed characterization of a mangrove forest at two sites, Kantang and Palian, which were monitored in Thailand and reported by Horstman et al. (2014), is included. Specifically, the data correspond to the areas where Rhizophora sp. trees are predominant. For Kantang, the calculated A h m values at a height of 10 cm are similar to those obtained for Zone 2 and Zone 3 at a height of 10 cm above the ground. At 5 cm, the presence of pneumatophores from Avicennia/Sonneratia trees causes A h m to increase considerably and become similar to those found in Zone 1 at a height of 10 cm, revealing the high density of roots in the Jicaral mangrove at the leading edge. Additionally, the vertical distributions of A h m , reported by Maza et al. (2019) and Kelty et al. (2022), are included in Figure 8 . In both studies, laboratory tests were conducted to quantify wave attenuation in a mangrove forest. The forests were constructed using mangrove tree representations, which were defined based on Ohira et al. (2013) parametrization for Rhizophora sp. trees. Due to the low tree density (0.0625 trees/m²) considered in the experiments of Maza et al. (2019), low A h m values are obtained compared to the data collected in Jicaral. On the other hand, Kelty et al. (2022) based their parametrization on mangrove measurements made in the Pacific of Costa Rica (Loría-Naranjo et al., 2014) where Jicaral is located, resulting in A h m values similar to those reported for Zones 1 and 2 for h m values greater than 0.5 m. Nevertheless, A h m values for h m = 0.5 m or less are lower than the A h m values measured in the field. This difference may be due to the use of the Ohira et al. (2013) parameterization, which considers only primary roots. The secondary and subsequent roots, which are primarily present at low h m values, are not considered. These roots may contribute up to 80% of the frontal area, as evidenced in the study by Yoshikai et al. (2021). Therefore, the role of stilt roots is essential for the wave damping capacity of Rhizophora trees.

Most of the models developed to obtain wave height attenuation by vegetation rely on the adequate estimation of C ˜ D . It is necessary to develop a C ˜ D formulation that considers the variation in mangrove forest characteristics and wave field conditions. Thus, field conditions are simulated using the SWAN model to obtain calibrated C ˜ D values. Figure 9 shows an example of six simulated cases, including the variation in the mangrove forest in the field and considering different wave conditions. The SWAN results are represented by continuous lines, whereas the field measurements are represented by dots. The error bars in Figure 9 indicate the variation in the SWAN results if the standard deviation of d r o o t ( Supplementary Table 1 ) is considered, following the methodology presented in section 2.2.2. In Zone 1, a maximum H m 0 difference of 8.4% and a minimum value of 3.5% with respect to the mean value are found when comparing the SWAN results using the d r o o t mean value with those using the upper and lower limits. However, in Zone 3, the maximum difference is 17.2% ( Figure 9D ), and a minimum value of 7.38% is found ( Figure 9A ). The cases represented in Figures 9D , E exhibit higher variability using the lower and upper limits of d r o o t along the transect. These cases correspond to the highest incident wave heights. This result indicates that a variation in the frontal area generates a greater change in energy dissipation in cases with high drag forces, i.e., high orbital velocities, than in cases where the wave height is lower and, consequently, the orbital velocities are also lower. The good agreement along the mangrove transect using a single C ˜ D value reveals that an adequate mangrove forest characterization that considers the spatial distribution of tree traits, combined with the varying velocity field affecting the forest, leads to a constant C ˜ D value across the forest, which agrees with the findings of Maza et al. (2017).

The C ˜ D values that best represent the wave height attenuation in the mangrove forest are shown in Figure 10 with their respective K C numbers with green dots. Furthermore, Maza et al. (2019) and Kelty et al. (2022) data for C ˜ D vs. K C are included, purple and yellow dots respectively. These data are used because they tested Rhizophora sp. representations under random wave conditions, similar to the conditions found in the field. For consistency, the d r o o t values reported by these two studies are also employed as the characteristic length scale in the K C formula. By considering a common relationship that allows us to obtain C ˜ D as a function of K C , that is, by knowing the incident wave conditions and the characteristic diameter of the forest, the data in this study, combined with those from Maza et al. (2019) and Kelty et al. (2022), are fitted to an equation in the form of C ˜ D = ( a K C ) b + c , resulting in Equation 9:

The proposed C ˜ D   equation based on K C as a function of d r o o t for Rhizophora sp. yields an acceptable result with a correlation coefficient ( ρ ) equal to 0.7, a root mean square error ( R M S E ) of 0.44 and a coefficient of determination ( R 2 ) of 0.5 ( Figure 10 ).

In addition, C ˜ D values from van Wesenbeeck et al. (2022) and Cao et al. (2016), obtained in the laboratory and field, respectively, are also included in Figure 10 to compare their results with those reported in this study. The values of C ˜ D are shown because they originate from woody vegetation but are not considered in Equation 9 since their geometric characteristics greatly differ from those of Rhizophora sp., leading to different wave-vegetation interaction patterns and, therefore, different drag force values. In the study by van Wesenbeeck et al. (2022), a unique large-scale laboratory experiment was conducted with woody vegetation, where the tested vegetation was pollard willow forest (Salix alba). They tested willows with full leaves and canopies, willows without leaves and full canopies, and willows without leaves and reduced canopies. Here, the cases in which the forest has all its leaves and full canopy are considered. In addition, the C ˜ D − K C relationship obtained by Cao et al. (2016) is also included as a continuous light blue line in Figure 10 . The curve corresponds to Avicennia marina, which has a root system different than that of Rhizophora sp., and is mainly characterized by pneumatophores. Figure 10 reveals that the C ˜ D values found for Rhizophora sp. are similar to those produced by pollard willow despite the differences in their respective structures. Similarities in C ˜ D values are obtained predominantly for large K C numbers when trees are exposed to greater wave hydrodynamics. However, further investigations are needed in these cases, when mangrove forests are exposed to high K C numbers, i.e., extreme events. In contrast, Cao et al. (2016) reported low C ˜ D values compared to those in this study since, unlike Rhizophora sp., Avicennia marina presents pneumatophores as roots instead of as stilt roots. When the water depth exceeds the pneumatophore height (approximately 12 cm), the maximum orbital velocities are not directly affected by the pneumatophores, leading to less drag force and consequently lower C ˜ D values.

Equation 9 provides an estimate of C ˜ D for regions predominantly covered by Rhizophora sp. forest and affected mainly by shortwave conditions since the recorded field conditions do not include long waves ( Figure 5 ). On the other hand, the influence of current tidal is explored by differentiating the C ˜ D values according to ebb or flood tidal cases ( Supplementary Figure 1 ). The influence of tidal currents in C ˜ D − K C formulation is not evident; this may be attributed to low current velocities. Calleja et al. (2022) reported current tidal values below 0.025 m/s for the inner part of the Gulf and zones near to the coast when evaluating the suitability of aquaculture farms in the Gulf of Nicoya. The reported current velocities are lower than those from studies that have reported an influence of current in wave attenuation when interacting with waves and vegetation (e.g., Yin et al., 2020; Hu et al., 2021).

The proposed Equation 9 is obtained by combining laboratory and field data, resulting in an integrated C ˜ D formulation with d r o o t as the characteristic length scale. To the best of our knowledge, this is the first attempt to combine laboratory and field data to develop a unique empirical equation C ˜ D for Rhizophora sp. forest.

The estimation of wave height attenuation based on drag-force models requires the characterization of mangrove forests to determine the frontal area exposed to flow action. In the laboratory, great efforts have been made to represent and to characterize in detail artificial mangrove trees using techniques such as laser scanning and photogrammetry (e.g., Maza et al., 2019; Kelty et al., 2022). In the field, detailed trunk and root diameter measurements (e.g., Horstman et al., 2014; Zhou et al., 2022) are often performed at the forest scale to capture mangrove traits due to limitations in employing techniques used in laboratory. When measuring trunk and root diameters in the field, significant effort is required to obtain similar vertical resolution of the tree frontal area as in the laboratory. Nevertheless, the results obtained for A h m in the laboratory and the field are comparable ( Figure 8 ). On the other hand, in the laboratory, the spatial heterogeneity of the frontal area along mangrove forests is not considered; in contrast to field studies, it is possible to capture the spatial variability in mangrove forest traits. The detailed mangrove forest characterization, as performed in Jicaral, allows for a reliable representation of field conditions while accounting for mangrove spatial distribution, which is essential for the adequate estimation of wave energy dissipation through the forest.

To estimate C ˜ D as a function of vegetation properties, it is important to determine an appropriate characteristic length scale that represents mangrove geometry and how it influences flow conditions. For Rhizophora sp. under random wave conditions, Maza et al. (2019) and Kelty et al. (2022) established the equivalent mean diameter and D B H as characteristic length scales, respectively. The characteristic length scale is closely related to Rhizophora sp. tree morphology, specifically to the mangrove tree elements that interact with waves. The parametrization of mangrove trees conducted by Maza et al. (2019) and Kelty et al. (2022) clearly revealed a main trunk (as illustrated in Figure 11A ) and primary roots that greatly contributed to the forest capacity to attenuate waves. For example, primary roots contribute between 75 and 85% of the total frontal area at depths of less than 75% of the highest root height. In the present study, Rhizophora mangle trees grew prostrate, making it difficult to distinguish between a main trunk and aerial roots (Jiménez and Soto, 1985; Duke et al., 1998). Consequently, roots are also the main contributors to forest attenuation capacity, and the main trunk and its associated D B H are not easily identified (as illustrated in Figure 11B ). Therefore, D B H is not the appropriate characteristic length scale for this type of tree; instead, the characteristics of the roots seem to be more representative of the mangrove geometry. A similar concept was described by Maza et al. (2017), where the characteristic length scale was depth dependent, leading to the consideration of roots or trunk diameter depending on the water level and its relationship with the vertical distribution of the trees. Therefore, this study highlights d r o o t as the appropriate characteristic length scale for obtaining a C ˜ D formulation that can predict wave height attenuation in mangrove forests. For this purpose, the Maza et al. (2019) and Kelty et al. (2022) root diameters are used to obtain Equation 9 rather than the equivalent diameter and D B H , respectively.

Wave attenuation rates have been used in numerous studies worldwide to evaluate the wave attenuation capacity of distinct mangrove species. When comparing r values, it is important to contrast the particular site conditions between studies (e.g., wave conditions, tidal range, forest characteristics, mangrove species). Supplementary Table 3 present the site characteristics, vegetation species and wave conditions for field studies conducted worldwide (e.g., Mazda et al., 1997a; Brinkman, 2006; Zhou et al., 2022). As can be evidenced, this study present similar wave conditions ( Figure 4 ) to those reported previously with most of the wave conditions not exceeding 20 cm and wave periods ranging between 2 and 6 s. The rates of wave height attenuation obtained along the forest in Jicaral are within the ranges reported by Quartel et al. (2007), with r varying between 0.002 and 0.011 m^-1; Mazda et al. (2006) reported values of approximately 0.0012-0.006 m^-1; and Horstman et al. (2014) reported values of 0.0061 and 0.012 m^-1 for Kantang and Palian, respectively. Studies conducted in the Colombian Caribbean region (Vanegas et al., 2019; Sánchez-Núñez et al., 2020) reported r values higher than those obtained in this study; these values vary between 0.06 and 0.29 m^-1 in a Rhizophora mangle forest, the same species that is present in Jicaral. This difference in r may be attributed to local effects. Vanegas et al. (2019) reported wave attenuation by a mangrove patch of 3.5 m in length in a shallow zone located 1 m from the shore, where dissipation mechanisms, such as breaking and bottom friction, together with the vegetation effect, may lead to high r values. Sánchez-Núñez et al. (2020) reported high slopes at the leading edge of a forest, which produce reflection and wave breaking, influencing r . However, the slope measured at the leading edge of the forest in Jicaral was very small, and it did not induce significant wave transformation processes. Despite the variability in   r values among studies, these values tend to decrease as water depth increases ( Figure 6 ) and increase as the incident wave height increases ( Figure 7 ), as shown in previous studies (e.g., Mazda et al., 2006; Horstman et al., 2014; Zhou et al., 2022).

Several studies have been conducted to obtain C ˜ D values in the field; these studies are based on different approaches. Mazda et al. (1997a), Quartel et al. (2007) and Zhou et al. (2022) obtained C ˜ D values for different species using the formulation reported by Mazda et al. (1997a), which is based on an equation developed to describe the effect of bottom friction in swell waves under shallow water conditions. The formulation considers the incident and final wave heights, the distance and the mean water depth between both points but does not consider the vegetation characteristics of the studied species. On the other hand, Vanegas et al. (2019) and Sánchez-Núñez et al. (2020) considered vegetation properties by employing the effective vegetation length scale defined by Mazda et al. (1997b), which resulted in high C ˜ D values compared to those of other studies ( Supplementary Table 3 ). These C ˜ D values were obtained using the methodology presented by Mazda et al. (1997b), which was developed for unidirectional flow conditions and is based on a free surface gradient produced along the forest, which may explain the high coefficients obtained. However, none of these studies reported a predictable formulation for estimating C ˜ D by combining hydrodynamics and vegetation properties. In this study, the use of SWAN model allows us to obtain wave energy dissipation produced by vegetation using mangrove trees charactersitics, i.e., the drag force approach (Mendez and Losada, 2004), rather than considering vegetation wave dissipation as an enhancement of bottom friction. Resulting in a predictive formula which has into consideration the inherent spatial vegetation properties and hydrodynamics conditions in the field. Additionally, the calibration was conducted by adjusting the SWAN results against the observed data at multiple locations within the mangrove forest. This approach extends beyond utilizing only two reference points, one offshore and one onshore, as conducted by previous authors. Furthermore, the following approach allows us to obtain C ˜ D values associated with only the vegetation effect and by numerically resolving other dissipation processes (such as wave breaking and bottom friction), allowing to quantify their contribution to the total dissipation along the forest ( Supplementary Figure 2 ). In contrast, aforementioned studies, based on field observations and analytical solutions, included all dissipation processes measured in the field when obtaining C ˜ D values.

The mangrove tree structure in Jicaral mainly consists of roots and trunks, which are difficult to identify, followed by a canopy layer ( Figure 3 ). If the mangrove canopy is not considered, the wave height attenuation may be underestimated when h reaches the canopy height (Maza et al., 2021; Zhang et al., 2023). Based on field observations of the canopy height in the study area, the canopy is expected to influence wave attenuation when the water depth at S1 exceeds 2 m. Based on the cases selected to obtain C ˜ D ( Figure 10 ), 35% of these cases exceed the water level threshold (black dots in Figure 12 ), indicating a possible overestimation of C ˜ D , as the frontal area of the canopy was not characterized nor included in the SWAN model. When the water depth exceeds 2 m, the SWAN model uses a constant frontal area calculated from the geometric tree characteristics measured at a height of 1.5 m above the ground. In this study, cases where the flow reaches the canopy level are still within the scatter ( Figure 12 ). This behavior suggests that when the water depth reaches the canopy level, the interaction between waves and the lower section of the canopy is similar to the interaction between waves and roots. Figure 3 shows that there is no identified interface between roots and the canopy; mixed layers occur at low canopy elevations due to the inherent prostrate growth of Rhizophora mangle, in contrast to the typical three-layer morphology used in several laboratory studies (He et al., 2019; Zhang et al., 2023) and other field campaigns (e.g., Horstman et al., 2014). Therefore, the turbulent structures generated by roots and branches seems to result in similar energy dissipation rates. Nevertheless, further research should be considered to characterize the detailed influence of the canopy on the resultant wave attenuation.