Quantifying Frontal-Surface Area of Woody Vegetation: A Crucial Parameter for Wave Attenuation

The last years, capacity of vegetation to reduce wave impact is receiving considerable attention. To predict wave attenuation processes within vegetation fields reliable estimates of vegetation parameters are needed. This proves to be difficult for woody vegetation as it consists of complex branch structures, characterized by varying branch densities, diameters and angles. State of the art physical and numerical models effectively use a single value for the diameter, bv and density, N of vegetation, which is unrepresentative for complex vegetation, such as trees. Trees can be better described by the projected frontal-surface area, Av. Hence, this work compares methods to quantify the Av in space for a pollard willow forest, and determines suitability of these methods for predicting wave attenuation using a spectral wave model (SWAN). We use data from manual measurements and Terrestrial Laser Scans (TLS), to estimate the vertical distribution of Av; and data from large-scale flume experiments performed on a willow forest to verify model sensitivity to Av inferences. As a baseline for comparison, tree models that describe the structure of the trees in various degrees of complexity are compiled. The most realistic tree model is used to quantify potential errors in TLS and basic manual measurements of N and bv. An initial comparison shows that the TLS data underestimates Av, which indicates that conducting manual measurements is more suitable to quantify a homogeneous forest. We found that the TLS suffers from shadowing effects (i.e., blockage of laser beams) and we recommend to apply a correction factor to improve its measurements. Furthermore, we identified the impact that the different methods to determine Av have on the estimation of wave attenuation using SWAN; in addition we verified the model results with data from large-scale flume experiments performed on the willow forest. The modeled sensitivity tests indicate large differences in wave attenuation and, consequently, a wide range (0.94–1.70) of bulk drag coefficients, C~D, for the various methods applied. This shows the variation of outcome between measuring methods and highlights the importance of stating the selected method for reliable frontal-surface area estimations, and consequently for reliable wave attenuation predictions.


INTRODUCTION
Current flood protection systems are under pressure due to climate change, in the form of more intense and frequent storms, heavier rainfalls, and accelerated sea level rise (Oppenheimer et al., 2019). Dikes will need to be made higher and stronger to keep up with changing climate, which will lead to increasing costs (Jonkman et al., 2013). Because of this, the interest in more costeffective solutions has grown in recent years for both coasts and rivers. Using vegetation in front of dikes has attracted a lot of attention for its contribution to the flood safety, as it adapts to rising sea levels, promotes sedimentation and is able to attenuate incoming waves (van Wesenbeeck et al., 2013). Besides the flood safety benefits, these so called "eco-engineering solutions" or "hybrid flood defenses" provide several ecological benefits, such as increasing biodiversity by providing nursery habitat for fish and breeding areas for birds (Rog et al., 2017).
Many studies were conducted on wave attenuation (oscillatory flow) through vegetation fields, such as kelp (e.g., Mendez and Losada, 2004), salt marshes (e.g., Vuik, 2019) and mangroves (e.g., Mazda et al., 2006;Quartel et al., 2007). Two processes that lead to wave energy damping through vegetation are wave breaking and wave-vegetation interactions (Vo-Luong and Massel, 2008). Both vegetation characteristics (such as vegetation density, geometry and height) and hydrodynamic conditions (such as wave-current interactions, wave orbital velocity, and wave period (Hu et al., 2014) influence these processes. The hydrodynamics around one rigid cylinder and an array of rigid cylinders is relatively well understood under oscillatory flow, and is commonly used to mimic the presence of vegetation fields in physical and numerical models. These models use a simplified description of tree structures containing a single diameter and density of the branches. Woody vegetation, such as mangroves in tropical regions and willow trees in temperate regions, have complex structures (i.e., varying branch densities, diameters, and angles). These structures not only make it difficult to capture the hydrodynamics, but also the vegetation parameters (such as frontal-surface area), which is needed to estimate the wave attenuation. This, along with the spatial and seasonal variations of vegetation, can lead to inadequate estimations of the vegetation parameters. Usually, the structure of the plants is simplified (e.g., rigid cylinders with a certain diameter and density) in analysis of field measurements and is assumed to be constant in space or roughly divided into coarse layers (Suzuki, 2011). These simplifications affect the wave attenuation from numerical or physical models, which in turn can be one of the reasons for the large range of C D values found in literature (e.g., Wu and Cox, 2015;He et al., 2019). To gain more insight into the wave energy dissipation through woody vegetation fields, it is necessary to have reliable estimates of the relevant vegetation parameters and accurate ways to measure them. We hypothesize that current methods of representing vegetation (illustrated in Figure 1) are overly simplified and that more thorough vegetation quantification including diameter tapering and branch angles will improve model predictions of vegetation and wave reduction.
The purpose of this study is to compare different methods for quantifying the frontal-surface area of woody vegetation and present their effect on the wave damping predictions. We use data from large-scale flume experiments with live willow trees under storm conditions, where significant wave heights, H m0 , up to 1.5 m were tested (van Wesenbeeck et al., 2022). Both manual measuring techniques and Terrestrial Laser Scanning (TLS) were conducted during these physical experiments, these are described in Section 3. The frontal-surface area estimation from the TLS data is compared to tree models (i.e., simplified representations of the tree structure). These tree models are developed by combining manual measurements and allometric relations, which are also included in Section 3. Afterwards, we use SWAN to calculate the wave damping corresponding to the different tree models and TLS output. These results are compared to the measured wave damping from the large-scale experiments in Section 4. Finally, methods for quantifying vegetation frontal-surface area, and recommendations for further research are suggested in Section 5.

Wave Attenuation by Vegetation
Many studies predicted the wave attenuation by vegetation using an increased bottom friction coefficient (e.g., Hasselmann and Collins, 1968;Quartel et al., 2007), while other studies (e.g., Dalrymple et al., 1984;Kobayashi et al., 1993) represented the vegetation as an array of cylinders, considering the wave forcing on these structures. The latter approach (i.e., cylinder approach) relates the wave attenuation to the plant geometry (i.e., the vegetation height, diameter and density), which is preferred over the former. Studies such as Dalrymple et al. (1984) used the timeaveraged energy conservation equation. Here, the wave energy dissipation by vegetation, ǫ v , is due to the work done by the waves on the vegetation, integrated over the submerged vegetation height. The time-averaged rate of wave energy dissipation per unit horizontal area becomes: The resulting force on the vegetation is given by the Morison equation (Morison et al., 1950): where F D is the drag force and F i is the inertia force. The Morison equation is solely suitable for determining the drag forces on slender objects (neglecting diffraction of waves around the object). We can still use this equation for predicting drag forces on trees with low canopy densities (i.e., no porosity effects namely effects of blockage and sheltering) (Etminan et al., 2019), and the ratio between the wave length and the diameter of the individual branches is small. Furthermore, most of the studies assume linear wave theory to be valid within vegetation fields; therefore, the total force on the vegetation is solely due to the drag forces as the inertia forces are out of phase with the velocity signal. Dalrymple et al. (1984) found an analytical solution for the wave height evolution through a vegetation field. His study focused on vegetation on a flat bottom subject to regular waves. This application has been extended by Mendez and Losada (2004), accounting for sloping bottom and irregular waves. The wave height decay through a vegetation field can be determined with the following equation (Mendez and Losada, 2004): with the wave attenuation coefficient, β, as where H rms is the root-mean-square wave height behind the vegetation field, H rms,in is the incoming root-mean-square wave height (in front of the forest), X is the distance inside the forest, C D is the bulk drag coefficient, N is the number of vegetation stems per horizontal area, b v is the average diameter of the individual stems, k is the wave number, h v is the height of the vegetation that is submerged, and h is the water depth. To date, the formulation of Mendez and Losada (2004) is widely applied and one can vary the vegetation parameters in space in the spectral version of this vegetation dissipation model. Suzuki (2011) adapted Equation (5) such that the vertical structure of the vegetation can be taken into account (i.e., C D , b v and N can vary over intervals along the vertical). This is implemented in numerical wave models such as SWAN (phaseaveraged model) (Suzuki, 2011). For other wave models, such as SWASH (phase-resolving model), the wave dissipation due to vegetation is based on the Morison equation as function of time, and requires also the same geometric information of the trees as vegetation input (i.e., N and b v ).
The bulk drag coefficient is an important, but an uncertain parameter (Kobayashi et al., 1993), as it generally accounts for the neglected processes. For example, we usually assume rigid cylinders, thus neglecting vegetation motion, which influences the relative velocity. Obviously, this becomes less valid for relatively flexible vegetation and can influence the wave attenuation by these fields. Furthermore, simplifications of vegetation geometry can also affect the wave attenuation processes. Nevertheless, several studies reveal relations between the bulk drag coefficient and dimensionless hydraulic parameters, such as the Reynolds number, Re (e.g., Hu et al., 2014) and the Keulegan-Carpenter number, KC, (e.g., Jadhav et al., 2013). Both Re and KC are predictors for the wake structures behind the vegetation (Sumer and Fredsoe, 1998). KC is also a measure for the relative importance between drag and inertia. These studies generally show that C D decreases with increasing KC and Re. In addition, other studies find better relations with non-dimensional parameters containing also vegetation characteristics such as the vegetation submergence ratio (Mendez and Losada, 2004) and the spacing between cylinders (Suzuki, 2011). These relations make it possible to chose a value for the bulk drag coefficient accordingly. In the case of woody vegetation with enough spacing between branches without motion, one can expect the C D values to be similar to that of rigid vertical cylinders, which holds a nearly constant value of ≈ 1.0 to 1.2 in the sub-critical flow regime (10 3 ≤ Re ≤ 10 5 ) (Sumer and Fredsoe, 1998;van Wesenbeeck et al., 2022).

Quantifying Vegetation
Knowing the vegetation densities (N) and average diameters (b v ) in space is usually sufficient vegetation input for wave damping models for large enough Reynolds numbers, keeping in mind that the bulk drag coefficient can be used for calibration. An example of two rather simple representations of a pollard willow tree is shown in Figure 1. Unlike the simplified description in these wave damping models, the realistic branch structure of woody vegetation is difficult to capture by means of a single diameter and density. It is difficult to get accurate approximations of these values, especially for woody vegetation because these strongly vary in space. Furthermore, the individual branches are also characterized by an angle (i.e., direction in which they grow) with respect to the incoming waves, which also complicates determining the 'projected' surface area. Alternatively, the frontal-surface area distribution of all the branches, A v (z) = b v .N, gives a better approximation of the entire tree structure. Hence, this research focuses on determining this parameter (van Wesenbeeck et al., 2022).
Obtaining vegetation parameters can be done in different ways, namely: by conducting manual measurements or using remote sensing techniques (e.g., optical techniques and lidar, which can be spaceborn, aerborne, or terrestrial); each having its own advantages and disadvantages.
Manual measurements can be categorized into destructive (i.e., for biomass measurements) and non-destructive measures (Nordh and Verwijst, 2004). The latter is a simple way of obtaining geometric information of trees and, accordingly, for calculating wave attenuation. For example, the density and diameters in a certain volume can be determined, after which a vertical and horizontal average value may be used. This is generally seen as a valid way to quantify vegetation and has been used in several studies (e.g., Jadhav et al., 2013;Möller et al., 2014;Ozeren et al., 2014). Some studies indicate that the vertical variation of vegetation parameters is important to consider (Ozeren et al., 2014), and this vertical variation is already incorporated into numerical models (Suzuki, 2011). Thus, a certain spatial resolution is needed depending on the vegetation type, and conducting hand measurements often becomes tedious.
Although solely hand measurements can be inefficient, we can still indirectly extent this data by using it as input for tree allometric relations. These are mainly formulations to obtain difficult tree parameters, such as tree height, through an relatively easy measured parameter such as the diameter at breast height. It is shown to work well for characterizing mangrove roots (Ohira et al., 2013). A branching method, which also uses tree allometry, is developed by Järvelä (2004) to calculate the area of trees. This method stems from the Strahlers ordering scheme, originally applied on river systems and later on tree structures (McMahon and Kronauer, 1976). This scheme starts at the smallest order ('finger-tip'), where two small bifurcations of the same order meet, they sum up to a higher order stream and so on, until it reaches the highest order stream (Order=n). Certain factors in between orders, also referred to as branching factors, are obtained and used to calculate the total area of a tree. These branching factors are mostly dependent on the tree species. For a detailed description on this branching method, we refer to Järvelä (2004).
Lastly, remote sensing techniques can be used to obtain vegetation parameters. Space-borne and airborne methods are effective for very large-scale measurements, for example identifying the state of the vegetation, but are less accurate for capturing vegetation parameters such as height and canopy structures (Srinivasan et al., 2014). On a smaller scale, optical techniques are commonly used, and provide detailed information on the geometry (e.g., Norris et al., 2017;Maza et al., 2019). This method becomes more difficult and less accurate when dealing with larger objects, such as trees. For tree and forest scale applications, Terrestrial laser scanning (TLS) shows to be a good alternative. A comparison between TLS results and the branching method of Jarvela showed a good agreement in terms of total area of trees, which was afterwards used for determining the flow resistance during floods (Antonarakis et al., 2009).

METHODS
Section 3.1 describes the large-scale physical experiments, which includes the set-up and measurements of the forest, and of the waves. Different methods for measuring and modeling the frontal-surface area in space result in different estimates for frontal-surface area distributions over height, Av(z). These methods are described in Section 3.2. Results of Av(z), which can follow from combinations of measurements and tree allometry relations, are called tree models in this work. These tree models are used to determine the required input for the wave Frontiers in Marine Science | www.frontiersin.org attenuation model, A v . Detailed data of the tree parameters and of the wave damping by these trees is used from large-scale physical experiments; hence, direct comparisons could be made between the measured wave attenuation and the calculated wave attenuation.

Physical Experiments
Full-scale experiments were conducted on live pollard willow trees under storm conditions. These experiments were carried out in a 300-m-long, 5-m-wide and 9.5-m-deep wave flume, where significant wave heights (H m0 ) up to 1.5 m were tested. A 40-m-wide forest was created with 32 pollard willow trees (Salix Alba) in front of a concrete levee. These trees formed 16 rows of 2, and were situated on a 85-m-long and 2.33-m-high platform. The platform represented a shallow foreshore, which permitted large wave height-water depth ratios to avoid having wave breaking inside the forest. A more detailed description of the experimental set-up is given in the paper of van Wesenbeeck et al. (2022).

Manual Measurements
The pollard willow trees (Salix Alba) were 15 years old with 3-year-old primary branches (i.e., branches that directly sprout from the knot of the tree). These primary branches were categorized into 3 classes (see Figure 2) based on the diameter at their base, D B (i.e., diameter at the location above the knot), namely: class 1 (D B > 50 mm), class 2 (20 < D B ≤ 50 mm) and class 3 (D B ≤ 20 mm).
The main tree characteristics were manually measured prior to conducting test series 1 (full canopy with leaves). These include the following measurements: • General tree characteristics: • Diameter and height of the trunk (D trunk , H trunk ) • Diameter and height of the knot (D knot , H knot ) • The number of branches per diameter class (N class1 , N class2 , N class3 ) • Diameter and height of the canopy (D canopy , H canopy ) After the experiments, we took the following additional measurements: • Characteristics of the primary branches: • The D B with corresponding total branch length was measured in order to create an allometric relation. Allometric relations are commonly used to estimate tree parameters, which are challenging to measure (such as the branch length), from an easier measurable parameter (such as the diameter at breast height). A random sample size of 340 primary branches out of the entire forest, which in total consists of 2,852 primary branches, was used. • Diameter decay relations were made for each branch class, as branches have a tapering form. This was achieved by measuring the diameter at increments of 1 meter along the main branch. A random sample size of 30 branches was used (10 branches per branch class). • Sketches of 9 branches were made. These sketches include the number, length and diameter of all the side branches; and they include the diameter of the main branch at vertical increments of 1 meter. • Additionally, strength and elasticity measurements were conducted on the primary branches; however, these are not analyzed in the present work.
• Distribution of branch diameter and density over the vertical, by measuring an individual tree (i.e., This tree was from the same location as the tested trees, and contained shoots of the same age), named "Tree S": • The general tree characteristics were measured (H trunk , D trunk , H knot , D knot , N class1 , N class2 , N class3 , D canopy , H canopy ) • The number of branches and their diameters were measured at 1 meter vertical increments. This was eventually used to compute the branch density (N) and the average diameter (b v ) within each vertical interval.

Terrestrial Laser Scanning
Besides taking manual measurements, we also used Terrestrial laser scanning (TLS). We scanned the forest prior and after each test series to monitor the state of the forest and to detect possible loss of tree biomass. These scans were taken with a FARO FOCUS 3D S 120 scanner from three fixed positions above the flume ( Figure 5).

Experimental Setup: Wave Measurements
The tests (T001-T042) were carried out with different wave parameters and vegetation configurations. In total 5 test series (TS) were conducted, namely: TS 1 on the trunks of the trees for low hydrodynamic conditions (i.e., low water levels and wave heights), TS 2 on full canopies with leaves, TS 3 on full canopies without leaves, TS 4 on reduced canopy density without leaves, and TS 5 without willows (calibration tests). A more detailed description of the wave conditions is given in the Supplementary Table S1.
The wave heights in front and behind the forest were measured with two types of wave gauges (i.e., radar wave gauge and resistance wave gauge). The wave damping due to the vegetation is expressed as the ratio between the incoming wave height and the wave height reduction (i.e., the difference in wave height between the tests with vegetation and the tests without vegetation), as applied in van Wesenbeeck et al. (2022): where Dr is the wave damping ratio by vegetation, H m0,veg and H m0,no veg are the measured wave heights with and without vegetation at the location behind the forest, and H m0,in is the measured incoming wave height in front of the forest.

Tree Model 1 -Primary Branch Model
The first method considers the primary branches (i.e., main shoots) and their tapering form. Information of the branching structure and side branches is not included, therewith the contribution of the side branches (secondary and tertiary) to the total A v is neglected. Allometric relations between D B and the total length of the branches (Section 3.1.1), the total number of branches counted per tree, together with diameter decay relations of the primary branches were used to predict diameter and the total number of branches at vertical increments of 1 m. With this, the frontal-surface area is obtained at vertical increments of one meter for each tree. Hereby it is assumed that the primary branches are oriented perpendicular to the knot, with the largest branches located in the center of the knot (see Supplementary Figure S1).

Tree Model 2 -Branching Model (By Example)
As in tree model 1 the side branches were neglected, a branching method similar to the one by Järvelä (2004) is used to account for the side branches. We made some adjustments based on the measured properties of pollard willows. This method will be illustrated with an example on one tree, named "tree S." Firstly, the regular pollarding practice on these trees impacts the tree structure, namely with relatively thick knot and trunk with respect to its branches, deviating from the "natural" willow structure. Therefore, the branches are regarded separately from the trunk, applying the method on the primary branches (excluding the trunk). We first identified the initial parameters for each tree. These are given in Table 1 for tree S. Besides these parameters, we also defined the minimum diameter, d min of 0.003 m and the height of the canopy, H canopy of 3.4 m.
Secondly, we used the detailed sketches of the branches to obtain three branching factors. These branching factors show the number of side branches that a mother branch can support (R B ), the diameter ratio (R D ) and length ratios (R L ) between subsequent branch orders. These factors are defined as follows: L m ; where m is the child branch, (m+1) is the mother branch, N b is the number of side branches, D is the diameter, and L is the length of the branch. The obtained branching factors differ per order ( Table 2), thus we maintained distinct values per order instead of working with the average values and are assumed constant for all trees.
Thirdly, the frontal-surface area was determined by applying the branching factors from Table 2 on the initial parameters from Table 1. The highest order branch is generally assigned to the trunk; however, in this work the primary branches are the highest order (3 rd order) as mentioned before. The steps to obtain the total frontal-area can be seen in Table 3. For example, the diameter factor, R D , is applied on the highest order branch (i.e., the base diameter), D 3 , until a diameter of the smallest order branch (d 1 ) is equal to the minimal diameter (d min ) of ≈ 0.003 m.
Fourthly, we applied a factor of 0.5 on the obtained frontal areas of Table 3 to account for the tapering form of the branches, as this is not accounted for in the initial method. Finally, we validated the branching method by applying it to the nine detailed sketches and on the measurements of the individual tree.
The method by Järvelä (2004) solely computes the total surface area of the tree; moreover, it lacks information of the  Lm ; where m is the child branch, (m+1) is the mother branch, N b is the number of side branches, D is the diameter, and L is the length of the branch.  vertical distribution (Järvelä, 2004;Antonarakis et al., 2009). The distribution cannot be determined from these methods as no angles between the branches, no tapering of the diameters, and no positions are considered. To determine the height variation of the surface area, we simulated the random structure of individual branches based on the allometric relations described above. The simulations allow for variations in the branch positions, while still following the branching rules. This results in slightly different branches from each simulation, mimicking "real" trees. The tapering form of the branches is included in the simulation by assuming linear diameter decay along its length. This tree model was validated as some tree parameters such as the angles of the side branches, the start positions of the primary branches on the knot, and the start positions of the side branches on their mother branch, do not have constant values and their distributions are not known. Thus, we validated this tree model against the hand measurements done on a single pollard willow (Tree S), keeping in mind that from these hand measurements we could also obtain b v,i at vertical increments of 1 m. The same step-size (1 m) was used to calculate b v,i from the structure as simulated with tree model 2. Furthermore, only diameters larger than 4 mm were considered, as this was the cut-off diameter applied in the hand measurements.
As mentioned above, the angles between branches were not measured. Still, this may influence the value of A v (z). Hence, a range of angles between the side branches were simulated from 10 • to 45 • with 5 • increments to determine the influence of the angles, and to chose a suitable value. The outcome is given in Figure 3. The angle between 1st and 2nd order branches, α 1 = 30 • , and the angles between 2nd and 3rd order branches, Finally, the frontal-surface area variation along the height, A v (z), was calculated. This is defined as the projected surface area of the simulated branches perpendicular to the incoming waves (see Figure 4). Thus, these three-dimensional simulations take the branch angles, positions and diameter decay into account.

Terrestrial Laser Scanning -Post Processing
Data from the scans of the leaf-less condition of the trees were used to obtain A v (z). This vegetation state coincides with the characteristics that we observe during extreme conditions, which are high hydraulic loads and leaf-less trees (winter season). Not to mention, scans during leaf-on conditions are less suitable to map the vegetation as most of the laser beams are directly blocked in the outer layers of the tree by the leaf density (see FIGURE 4 | Schematic of the coordinate system used for calculating A v with the incoming waves from the positive y-direction. P is the endpoint of the branch, and P' is the projected branch perpendicular to the incoming waves.
Frontiers in Marine Science | www.frontiersin.org   Figure 5 shows the entire forest without leaves and the 3 positions in the flume from which the TLS scans were created.
The post-processing of the point cloud data included sub-sampling, segmentation and filtering of the point cloud, performed in CloudCompare. Segmentation was done to remove excess information (e.g., measuring equipment, and flume walls) and sub-sampling was needed for the sake of work-ability of the point cloud. Afterwards, the point cloud data of each tree was analyzed separately and the frontal area (i.e., does not include for the angles of the branches relative to the incoming waves) was constructed using alpha − shape function, which is a built-in function in MATLAB. This function uses an alpha parameter to control how it creates the bounding region around the 3D-point cloud. An alpha shape value of 0.01 is used in this study, shown in

Simulating Wave Attenuation and Bulk Drag Coefficients
The Morison equation is solely suitable for determining the drag forces on slender objects (neglecting diffraction of waves around the object), but can be used for predicting drag forces on trees as well. In this study, the forest is characterized by low canopy densities of around 13 branches per m 2 . The average spacing, S becomes 28d, with d = 10mm as the average branch diameter. The density parameter as described in Nepf (1999), ad = d 2 / S 2 , is approximately 0.001, indicating a relatively sparse forest; hence, we can assume no porosity effects (i.e., blockage and sheltering). Thus, the drag coefficient can be considered similar to that of a single cylinder (Etminan et al., 2019). In addition, the ratio between the wave length and the diameter is small, following the characteristic of slender objects.
SWAN in 1D-mode is used in this study to simulate the waves through the forest. The SWAN model was defined in Cartesian coordinates with a spatial discretization of 1 meter. The model is forced with identical JONSWAP spectra to those measured at the deep water locations in the experiments (WHM01-03). SWAN was executed without accounting for wind growth, white capping, refraction, diffraction, quadruplets, triads and turbulence dissipation, as these processes are not relevant in this flume experiment. On the other hand, frequency shifting in frequency space, bottom friction with a roughness constant of 0.07 m 2 s −3 , wave breaking, and vegetation dissipation were activated. The wave energy dissipation through vegetation fields accounts for the vertical structure of the vegetation and follows the action balance formulation, as shown in Suzuki (2011). Following the same manner, the energy dissipation by vegetation becomes: where ǫ v is the averaged wave energy dissipation due to vegetation; C D the bulk drag coefficient, g the gravitational acceleration constant; k the mean wave number; α i the ratio between the depth at the top of the vertical layer i and the total water depth (α i ≤ 1), where i = 1 represents the lowest vertical layer and i = I the highest vertical layer; h the water depth; H s,in the significant wave height; and A veg,i the total frontal width of vegetation perpendicular to the waves for layer i per unit area. Varying this vegetation parameter A v,i , according to the different tree models, will result in different wave attenuation outcome.
The bulk drag coefficient, C D , is important for predictive wave models and is therefore often estimated in studies on vegetation (e.g., Suzuki, 2011;Maza et al., 2015). In this work, the SWAN model is tuned with the measurements such that the modeled wave damping corresponds to the measured wave damping behind the forest. To achieve this, the tuning parameter, namely the bulk drag coefficient, is calibrated using exhaustive search. We made a distinction between the drag coefficient of the canopy (C D,can ) and the trunk (C D,tr ) for the tests with high water levels. For these high water level tests, the drag coefficient of the trunk is set to a constant value of 1.2 (Wieselsberger, 1921). The SWAN set-up is made following the same approach as described in van Wesenbeeck et al. (2022). This is repeated for different vegetation input [i.e., A v (z)] from the tree models and TLS results.

RESULTS
This section shows general characteristics of the 32 pollard willow trees (Section 4.1), and the characteristics of the resulting tree models (Section 4.2). We analyzed the frontal-surface area, A v (z) of each tree model, and afterwards we calculated the corresponding wave attenuation and compared these results. In addition, the relation between Keulegan-Carpenter number (KC) and drag coefficient (CD) for each tree model is presented (Section 4.5). This provides us more insight into the errors of different tree models and their effects on predicting wave attenuation.

General Characteristics of Pollard Willows
From the manual measurements it follows that the smaller primary branches (Class 3) take up a large portion of the total number of branches in a pollard willow tree (64 %), while the larger primary branches (Class 1) only account for 3 % of the total number of shoots. Therefore, a larger data set could be obtained for the smaller branches. The fraction of each diameter class is given in Figure 7A with corresponding 95%-confidence interval. Furthermore, the gathered data of the D B is plotted in Figure 7B, which includes for the average value and the 95%confidence interval.
Other tree parameter such as the D b , N high and L high , are given in Supplementary Table S2.

Tree Model 1: Primary Branches
Solely considering the primary branches of the pollard willow trees is a relatively simple way of mapping the vegetation. Besides this, relatively less wave attenuation is expected from the smaller side-branches as these are more flexible than the primary branches. The manual measurements done at breast height (i.e., the number of branches per class), the allometric relation between D B and branch length, and the diameter-decay relations were the main input for this tree model. We assumed an average diameter at breast height for each class (see Figure 7B). The D B -length relation, diameterdecay relation, and a simplification of the tree without sidebranches (tree model 1), are shown in Figure 8. This model is considered as the lower limit of the frontal-surface area as it only accounts for the main branches (neglecting the side branches). It was initially expected to give a reasonable estimate of the total A v ; on the contrary, the comparison to the other tree model that follows hereafter shows that this is not the case.

Tree Model 2: Branching Method
The branching method uses RNG to describe the branch structure for all branches larger than ≈ 3 mm, which is used to create tree model 2. To statistically estimate the A v (z), trees were simulated in Python to generate random branches that follow the branching rules shown in Table 2. The parameters from Supplementary Table S2 are input for generating the branches. An example of a generated branch from each class is shown in Figure 9A. Figure 9B illustrates tree model 2, which is formed by these generated branches positioned on the knot. In Figure 11A the A v (z) as sampled from the generated trees is shown. As each realization yields a different A v (z), the output is a distribution for which the mean and the 95% confidence interval are determined.

TLS Correction Factor
The comparison of A v (Figure 11B) showed an overall underestimation by the TLS results. Though, the TLS showed an overestimation in the upper layers of the trees. This suggests that occlusion due to the canopy density is the main reason for the observed errors in the lower canopy layers. A correction factor was applied on one tree at the edge of the forest to improve the TLS results by accounting for the percentage of blocked laser beams.
A TLS,new,n = A TLS,n * f cor,n , where A TLS,n is the area obtained from the 3D-point cloud in layer n, f cor is a correction factor based on the number of blocked laser beams, N t is the total amount of laser points send out by the laser scanner, and N is the number of laser points blocked in the layer before reaching layer n. The correction increased the A v,tot with 2.45 m 2 . However, this correction could only be applied to trees near the laser scanner positions as occlusion by neighboring trees could be neglected. We applied this to tree P1 as this tree is also most likely not hit by laser beams other than the laser scanner nearest to it (see Figure 10).

Comparing Frontal-Surface Area
The three different tree models are first compared considering their frontal-surface area distributions, shown in Figure 11B.
All the tree models are set to a similar resolution (i.e., namely at vertical increments of 1 m) of the results to get a more clear comparison between the models. It shows that side branches amount to at least ≈ 25 % more A v,tot than only the primary branches. Furthermore, the TLS results underestimate the canopies total frontal-surface area compared to the other tree models. The underestimation by the TLS is mainly seen in the lower layers of the trees, while in the upper layers there is an overestimation.

Comparing Wave Attenuation
A strong relation is found between C D and the Keulegan-Carpenter number (KC). The KC number is usually calculated as follows: KC = U×T L , where U is the orbital velocity, T is the wave period and L is the average diameter of the object. For determining KC, we chose the maximum horizontal orbital velocity based on linear waves U 0max , and T m02 as a representative wave period of the wave spectrum. The definition of KC is non-trivial for realistic tree structures, as these structures are characterized by a distribution of branch diameters that varies over the height. The KC number in this work is calculated at vertical increments of 1 m and averaged over the height.
The results show a decrease of C D with increasing KC numbers (Figure 12). Furthermore, the differences in C D between the tree models are a consequence of different frontal-surface area distributions. Tree model 1, which is considered as the lower limit for the frontal-surface area, shows larger C D values relative to tree model 2, due to its lesser frontal-surface area. The largest C D values are observed for tree model 3. This is according to expectation, since the TLS method underpredicts the frontalsurface area of the trees (Section 4.4).

DISCUSSION AND CONCLUSION
Globally, considerable claims are made of how trees, such as mangroves and fresh water trees, reduce impact of waves during extreme storms and tsunamis (EJF, 2006). Most research on wave attenuation by vegetation focused on investigating hydrodynamics around the vegetation structure, whereas less attention is given to vegetation parameters and methods to quantify these structures (e.g., Huang et al., 2011;Wu and Cox, 2015). Vegetation, such as in salt marshes, is generally quantified in a simplified manner using a spatial average diameter (b v ) and density (N). For example, studies on mangrove forests used similar approaches, where single values for the density and diameter were used to represent the entire structure (e.g., Phan et al., 2019) or used single values per layer to represent the roots, trunk and canopy separately (e.g., Vo-Luong and Massel, 2008;Narayan et al., 2011;Ismail et al., 2012). More thorough quantification of especially complex canopies seems to be lacking, but is needed to obtain insight in processes that determine wave attenuation during extremes accurately. For woody vegetation (e.g., willows and mangroves) the frontal-surface area parameter, A v (z) is an useful way of quantifying this. Some studies have described vegetation using similar parameters concerning the tree frontal area (e.g., Maza et al., 2017Maza et al., , 2019. However, these solely considered damping through roots and trunks, neglecting the canopy. Previously, large-scale experiments were conducted on a well-defined 40-m-long willow forest (Salix Alba) for which the wave attenuation during storm conditions was determined (van Wesenbeeck et al., 2022), showing that most of the wave damping occurred at mid-water levels through the canopies of the trees. This is where most of the biomass and surface-frontal area was present; hence, improving our estimates of the entire vegetation structure and developing methods to quantify this will improve predictability of wave attenuation by vegetation for different water levels and wave heights. We used these large-scale experiments to obtain more insight in methods to extract A v (z) of the trees, using both manual measurements and terrestrial laser scanning (TLS). Manual measurements executed on all the branches at vertical increments of 1 meter, are relatively simple measurements, but tedious to conduct on every tree. The accuracy of these measurements depends on the size of the vertical increments (i.e., layer thickness). With this, the diameter decay along the vertical is taken into account to a certain extent. However, the angles of the branches are neglected and side branches that are within the layers are also neglected. Using simple manual measurements as input for allometric relations for trees, is an efficient way of including side branches. Based on this, a three-dimensional tree model is created (included positions and angles of the branches), which was used to determine the frontalsurface area distribution in the vertical, A v (z). For this method, it is important to point out that the initial tree parameters and branching factors used in this research are applicable for "pollard" willows with 3-year-old shoots (Salix Alba); these values are likely to differ for "natural" willows of the same species. The simulations from this method can nevertheless be used to obtain a more thorough representation of the trees for numerical and physical models.
The side branches are relatively flexible branches, leading to more extreme branch motion. This in turn, can lower the significance that these small branches may have on attenuating waves. On the other hand, it was found that the side branches could take up to 25% of the total frontal-surface area of these trees. Thus, including these branches in tree models may be necessary for wave attenuation predictions.
A more efficient method for assessment of A v (z) could be through TLS measurements. The TLS is a practical method to conduct measurements on a forest scale and has the potential to be combined with satellite data into useful parameters for wave damping models. However, the results of the TLS depend on the post-processing and A B FIGURE 10 | (A) Tree P1 with, (B) corresponding graph containing the relative error of the laser scanner as function of the penetration distance of the beam, r. The colors and the horizontal axis indicate the distance from the laser scanner.

A B
FIGURE 11 | (A) The average frontal-surface area (A v ) of the 32 trees in the large-scale experiments with 95% confidence interval using tree model 2, (B) frontal-surface area comparison of the tree models, considering the average tree without leaves in the wave flume.
in our case, the chosen α-value. We showed that the TLS underestimates the frontal-surface area, especially near the trunk region where branch densities are largest, which is mainly due to shadowing effects. However, this can be corrected by applying a factor that accounts for the shadowing phenomenon (blockage of the laser beams by high branch densities). These findings encourage further research on the TLS as it is a promising measuring technique for FIGURE 12 | Relation between C D and KC for tree model 1 (only the Primary branches), tree model 2 (primary branches including the side branches), and tree model 3 (TLS measurements).
mapping large heterogeneous forests in the field. For this, it is recommended to position the TLS in the direction of the incoming waves to investigate the blockage factor in more depth. Furthermore, a decrease in bulk drag coefficient with increasing Keulegan-Carpenter number (KC) was shown, which is a trend found also in other vegetation studies (e.g., Anderson, 2011;Ozeren et al., 2014;Etminan et al., 2019). The numerical model validation showed large differences in wave attenuation for different tree models and a corresponding large range in bulk drag coefficients, especially for lower KC numbers. This suggests that the discrepancies between the tree models in terms of frontal-surface area is significant for the wave attenuation. A range of C D between 0.94 and 1.70 for high KC numbers was obtained, emphasizing the importance of reliable frontal-surface area distributions.
Thus, having feasible methods to obtain reliable frontalsurface area of vegetation fields is relevant for reliable wave predictions. This will increase the potential of implementing hybrid flood defenses, existing of coastal or riparian vegetation and engineered structures, such as embankments or levees. Only a combination of more accurate measurements on vegetation, data and modeling can be used to systematically obtain more insight in what vegetation can actually contribute to flood risk reduction during average and more extreme conditions.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BH and BW designed the research project. The calculations and analysis were performed by SK and supervised by BH. JA carried out the swan simulations. SK carried out the tree simulations and wrote the manuscript with input from all authors. All authors discussed the final results, contributed to the manuscript revision, and approved of the submitted version.

FUNDING
This study was done within the research program WOODY, number 17194, which is financed by the Dutch Research Council (NWO) via the Open Technology Program. carried out on the Dutch national e-infrastructure with the support of SURF Cooperative under project EINF-666 Emerging technologies for Nature Based solutions in coastal waters, and we acknowledge support by GNU Parallel (Tange, 2021).