Wave transmission and drag coefficients through dense cylinder arrays Implications for designing structures for mangrove restoration

At retreating mangrove coastlines, bamboo structures are built to create new habitat for mangrove colonization. Existing structures have experienced mixed rates of success due to the lack of a scientific basis in their design. Optimizing future structure designs requires investigating the effect of the bamboo poles on waves. We conse- quently conducted laboratory experiments to measure wave transformation, hydrodynamic forces, and flow velocities inside cylinder arrays, mimicking bamboo poles, with varying cylinder configurations and orientations. The experiments provided relationships for wave transmission, wave reflection, and the drag coefficients for configurations with volumetric porosities between n = 0.64 (cid:0) 0.9. Configurations with a small lateral spacing (causing higher blockage) and a relatively longer streamwise spacing (causing less sheltering) exhibit larger forces and dissipation per element. Such arrangements enable optimizing wave dissipation at locations where the wave direction has low variability over the year. Placing the poles horizontally instead of vertically increases the forces and wave dissipation per element in relatively deeper water. Based on the experiments, we developed a conceptual analytical model that predicts wave reflection and dissipation through cylinder arrays, including blockage and sheltering. The model can reproduce the influence of cylinder arrangement on wave trans- formation, and it suggests that accurate predictions of sheltering and wave reflection are important to find optimal designs. Overall, these results provide useful insights on how to model and optimize the design of structures for mangrove restoration.


Introduction
Mangrove ecosystems have gained considerable interest as a coastal protection tool during the last decades. Their protective value results from their ability to reduce waves Bao, 2011) and to promote sediment accumulation . They can thus constitute an efficient coastal defence under climate change scenarios (Temmerman et al., 2013;Cuc et al., 2015). Despite their importance for flood risk reduction, 30% of the world's mangrove forests have been lost over the last 50 years (FAO, 2007).
Mangrove removal can expose the remaining vegetation to storm waves, favoring erosion and hampering the natural restoration of the ecosystem (Winterwerp et al., 2013). Coastline retreat has been mitigated by building bamboo and brushwood structures at mangrove coastlines in South East Asia and South America (Winterwerp et al., 2005;Schmitt et al., 2013;Winterwerp et al., 2013;Cuong et al., 2015;van Wesenbeeck et al., 2015;Nguyen and Parnell, 2017). The primary function of the structures is to attenuate incoming waves. When sediment is transported landwards by the filling tide, the lower wave stirring behind the structures enhances sediment settlement, favoring coastline expansion and mangrove colonization.
Several structure designs have been implemented. Most structures consist of two or more rows of bamboo poles perpendicular to the direction of wave propagation, often combined with horizontal poles or a brushwood filling, as shown in Fig. 1. Their width ranges between 0.7 and 1.5 m in the wave direction, and their volumetric porosity varies between n ≈ 0.5 − 0.9. Despite having shown some promising results, e. g. Winterwerp et al. (2013), existing implementations do not always successfully create new mangrove habitat.
The variable performance of the structures may relate to the lack of a scientific basis in their design. Sediment settlement at the coastline, necessary to protect the mangrove habitat from erosion, requires low wave action. However, we lack models to predict wave transmission through the bamboo structures. As a consequence, existing designs are not adapted to the local wave conditions and sediment properties.
Guidelines and models developed for traditional coastal infrastructure are not directly applicable for the bamboo structures. Breakwaters also attenuate waves but since they are usually built with rubble mound or concrete, they have different hydraulic properties. When used for coastal progradation, breakwaters are designed to attenuate waves in order to slow down wave-driven alongshore currents, which is a separate physical process involving different design considerations (Villani et al., 2021). Moreover, since currents are the agent transporting sediment towards the bamboo structures, blocking them could be counterproductive for coastal accretion.
The aim of this paper is thus to develop design tools for the bamboo structures, and to identify more cost-effective designs. We consequently collected laboratory measurements of wave transformation through cylinder arrays, which mimicked the bamboo structures. The measurements compared the relative performance of different cylinder arrangements, and enabled investigating wave dissipation processes inside the structures and how to parameterize them. We also developed an empirical model based on the experimental data, which can be applied to predict the wave height behind the structures, and to identify the most efficient designs for wave dissipation.
The outline of the paper is as follows. The theoretical background of wave-structure interaction is presented in Section 2. The set-up of the experiments, the data analysis, and the development of the analytical model are explained in Section 3. The results of the experiments and the model are presented in Section 4. The limitations and implications of this work are included in Section 5, and its main conclusions are outlined in Section 6.

Wave dissipation
When a wave propagates through a bamboo structure, part of its energy is reflected seawards, as shown in Fig. 1. Wave reflection increases the flow velocity in front of the structure, and enhances scour at the toe. The remaining wave energy continues travelling into the structure and exerts hydrodynamic forces on the poles. These forces consist of several parts; (1) skin friction forces, (2) form drag forces due to flow separation behind the poles, and (3) inertia forces, associated with the acceleration of the wave-driven flow. The friction forces are often neglected, since they are much smaller than the form drag components. For a vertical cylindrical element the total in-line force F x is often parameterized using the Morison equation (Morison et al., 1950): where h is the still water depth, η is surface elevation, ρ is the water density, c D is the drag coefficient, d is the pole diameter, u is the local horizontal flow velocity, c M is the inertia coefficient, ∂u ∂t is the horizontal flow acceleration, and z is the vertical coordinate. The work done by the hydrodynamic forces over a wave cycle, ϵ v , causes wave energy dissipation, which reduces wave transmission through the structure. For an array formed by vertical cylinders, the total work done by the horizontal in-line forces, F x would be given by Eq. (2): where N is the cylinder density per unit area. Since the inertia force and the velocity are 90 • out of phase, the wave-averaged work, ϵ v , is dominated by the drag component of Eq. (1). When an incoming wave H I encounters a structure, there is a reflected wave component H R that propagates seawards and may cause scour at the toe. Another part of the wave energy is dissipated due to drag through the structure. The smaller transmitted wave height, H T , enhances sediment deposition, creating new potential habitat for the mangroves. The structures can have different configurations, such as (c) vertical bamboo poles, (d) vertical and horizontal poles and, (e) vertical poles with a brushwood filling.
Assuming negligible wave reflection and that the velocity field can be described by linear wave theory gives Eq. (3) (Dalrymple et al., 1984): where c D,b is an empirical bulk drag coefficient, which includes the effect of disturbances of the velocity field by the cylinder array, k is the wave number, ω is the wave frequency, g is the acceleration of gravity and H is the wave height. A horizontal pole exposed to waves experiences form drag forces in both the horizontal and vertical direction, since the water particles move in elliptical motions. Suzuki et al. (2019) expanded the expression for horizontal structures by incorporating the work done by the vertical drag forces, resulting in Eq. (4): The additional vertical dissipation term implies that changing the element orientation from vertical to horizontal may increase wave dissipation. The vertical velocities, drag forces, and associated dissipation, are relatively larger in deeper water conditions, characterized by a large ratio of the water depth to the wave length. The vertical velocities and their effects on wave dissipation are smaller in shallow water, which produces comparable wave dissipation rates by vertical and horizontal cylinders. The effect of element orientation on wave dissipation has been included in analytical and numerical models, e.g. Suzuki et al. (2019), but it was not tested in previous laboratory studies.

Drag coefficients
Predicting wave dissipation by the bamboo structures relies on the knowledge of the drag coefficient. (Mendez and Losada, 2004;Suzuki et al., 2012). However, most literature has investigated the drag values for a single cylinder (Keulegan and Carpenter, 1958;Graham, 1980;Obasaju et al., 1988), or for sparse cylinder arrays with volumetric porosities above n = 0.78 (Sarpkaya, 1979;Arunachalam et al., 1981;Heideman and Sarpkaya, 1985;Augustin et al., 2009;Hu et al., 2014;Ozeren et al., 2014;Chen et al., 2018;Etminan et al., 2019;Phan et al., 2019). The drag and inertia coefficients derived in those studies are often expressed as a function of the KC number, defined as the ratio of the wave excursion ξ to the cylinder diameter d (Keulegan and Carpenter, 1958). KC represents the relative importance of the drag and inertia force components, with 0 < KC ≪ 20 − 30 corresponding to inertia-dominated cases, and KC > 20 − 30 associated with dragdominated conditions (Sumer and Fredsoe, 1997). For KC > 100 the Sheltering under waves depends on the length of the wave excursion ξ compared to the spacing between cylinders s x and their diameter d. (e) If ξ/s x < 1, the wake of one element does not reach neighbouring cylinders. (f) If 1 < ξ/s x < 5 − 7 sheltering depends on ξ/s x (Heideman and Sarpkaya, 1985;Suzuki and Arikawa, 2010). (g) If ξ/s x > 5 − 7 sheltering depends on s x /d (Heideman and Sarpkaya, 1985). drag coefficients converge with the values of steady flow, and the KCdependency disappears (Keulegan and Carpenter, 1958;Graham, 1980;Obasaju et al., 1988).
The drag coefficients for arrays are often calculated by fitting Eq. (1) with velocities either measured upstream from the array, or estimated assuming a harmonic flow (Sarpkaya, 1979;Arunachalam et al., 1981;Heideman and Sarpkaya, 1985). Neglecting the hydrodynamic changes inside the arrays has resulted in considerable variability in the drag values found in literature, with values ranging between c D,b = [0, 16] for emergent rigid cylinders (Sarpkaya, 1979;Chakrabarti, 1981;Heideman and Sarpkaya, 1985).
The variability in c D,b for emergent array has been mostly attributed to the processes of sheltering and blockage. Sheltering takes place when downstream rows are exposed to the wake of upstream elements, and they thus experience lower drag forces (Heideman and Sarpkaya, 1985;Bokaian and Geoola, 2008;Liu et al., 2008). Sheltering in wave driven flows depends on whether the wave excursion is long enough to reach the next neighbouring cylinder (Heideman and Sarpkaya, 1985;Suzuki and Arikawa, 2010), as illustrated in Fig. 2 (e-g). When the ratio between flow excursion, ξ, and streamwise separation, s x , is between 1 < ξ/s x < 5 − 7, sheltering is a function of the excursion length compared to the spacing (Heideman and Sarpkaya, 1985;Suzuki and Arikawa, 2010). For larger ξ/s x ratios sheltering relates to the dimensionless spacing between cylinders s x /d (Heideman and Sarpkaya, 1985), analogously to uniform flow. Other studies found that in relatively sparse emergent vegetation, sheltering effects can often be neglected, and that the drag forces are well described by blockage. Blockage refers to flow acceleration through a cross-section of the vegetation (Etminan et al., 2017;van Rooijen et al., 2018;Etminan et al., 2019), which causes higher drag forces on the cylinders.
Previous studies often focused on one of the two processes, which dominated the drag forces in their application. However both sheltering and blockage may influence the drag forces for the bamboo structures. The influence of the distance between elements on sheltering and blockage is illustrated in Fig. 2 (a-d), with a smaller lateral spacing s y increasing blockage effects, and a smaller streamwise separation s x increasing sheltering effects.
In view of the processes of sheltering and blockage, we hypothesize that structures with small lateral distance s y (increasing blockage) and a relatively longer streamwise separation s x (decreasing sheltering) could maximize the forces, and thus the energy dissipation per element. This type of geometric arrangement could consequently reduce the material costs of a structure. However, excessively low lateral spacings could increase wave reflection and scour, hindering structure stability (Winterwerp et al., 2013). An optimum structure should thus maximize wave dissipation while minimizing reflection.

Methodology
In order to explore the effect of cylinder arrangement and orientation on wave transformation and the drag coefficients, structure prototypes consisting of arrays of cylinders were tested in a wave flume. We measured wave transformation, hydrodynamic forces and flow velocities inside the arrays. The experiments focused on regular or in-line configurations, since they simplify the analysis of the physical processes. Random arrangements would provide higher drag variability under comparable conditions, and besides this, regular arrangements may provide more efficient designs that maximize blockage and minimize sheltering.
Due to the properties of the bamboo structures, the experiments enabled investigating processes not addressed in other studies. We tested cylinder arrays that were denser and shorter (in the direction of wave propagation) compared to most experiments in the literature, since fully developed canopy conditions cannot be directly applied to the bamboo structures. Wave reflection was also analyzed, given the higher density of the structures compared to natural vegetation, and the potential detrimental effects of reflection on structure stability. The setup of the experiments and the data analysis are explained in the following sections.

Wave generation
The laboratory experiments were conducted in a wave and current flume at Delft University of Technology. The flume is 40 m long, 0.8 m wide and 0.8 m high. A wave generator with an active reflection compensation system was placed at one side of the flume and a wave absorber at the opposite end. We prescribed the second-order steering of the wave maker for all tests. Monochromatic waves were generated, with a water depth of h = 0.55 m, a wave height of H = 0.13 m, and periods of T = 1,1.25,1.5,1.75,2, and 3 s, respectively.

Physical model
The generated waves propagated through a frame with cylinders placed in the middle of the flume, as illustrated in Fig. 3 (a). The physical model consisted of a grid of 0.76 x 0.76 m, where aluminum cylinders could be introduced in different arrangements. The elements were held together by a top and a bottom plate, as shown in Fig. 3 (c). The cylinder diameter was d = 0.04 m for all experiments. The tested configurations are illustrated in Fig. 4. The configurations are named based on their lateral spacing (D, for dense with s y = 1.5d, and S, for sparse with s y = 3.0d), their streamwise spacing s x (also D or S), the number of rows, and the cylinder arrangement (with R for regular or in-line, and T for staggered). Most experiments were conducted with vertical cylinder arrangements, starting with one single row, and adding additional rows in successive experiments. For a smaller subset of configurations, indicated by an asterisk in Fig. 4, the frame was also placed horizontally in the flume to analyze the effect of cylinder orientation on wave transformation.

Instrument set-up
For each cylinder arrangement we measured the water surface elevation, flow velocities and the forces acting on individual cylinders. The locations of the instruments are presented in Fig. 3 (a). All the instruments were measuring throughout each experiment with a frequency of 100.
The water surface elevation was measured with capacitance-type wave gauges; two in front of the structures (WG1 and WG2 in Fig. 3), and two behind it (WG3 and WG4 in Fig. 3). The output of the wave gauges was in volts, and the surface elevation was obtained from linear regression, using separate calibration factors for each of the wave gauges. The accuracy of the gauges was 1% (Delft Hydraulics, year unknown). The separation between each pair of wave gauges was set to 0.25 times the wave length of each wave condition, for optimal wave reflection analysis (Goda and Suzuki, 1976).
An electromagnetic flow meter (EMF) was placed at a distance of 0.4 m upstream from the structure. The elevation of the EMF was changed between tests to provide velocity measurements at 3 elevations from the bed: z = 0.15 m, 0.25 m, and 0.4 m. The EMF measurements had an accuracy of approximately 1% (Delft Hydraulics, 1990).
The velocities inside the structures were measured with a Nortek Vectrino acoustic velocimeter (ADV), placed 0.04 m upstream from the center of a cylinder (see Fig. 3 (c) and Fig. 4). The elevation of the ADV was also varied between tests, and it measured the flow velocity at z = 0.15 m, 0.25 m, and 0.4 m. The ADV had an accuracy of approximately 1% (Nortek, 2020).
The hydrodynamic loads acting on a single cylinder were recorded with a SCAIME load cell mounted on the upper part of the element, measuring in volts with 0.017% accuracy (SCAIME, 2020). The load cells were calibrated using known weights, and fitting a linear relationship between weight and voltage output. The measured forces were calculated by multiplying the sensor output by the calibration factor, and by the acceleration of gravity. Force and velocity measurements were only collected for the vertical orientations, since we could not introduce the sensors inside the horizontal structures without removing multiple elements.

Wave transformation
The incoming and reflected wave components were separated with the method of Goda and Suzuki (1976), using the dispersion relation for non-linear waves of Kirby and Dalrymple (1986), consistently with the second-order wave paddle steering. The method was applied to the time series of WG1 and WG2 to calculate wave reflection in front of the structure, and to the measurements of WG3 and WG4 to calculate the wave transmitted through the structure. The analysis was conducted over time intervals during which (1) the propagating wave had already reached the wave gauges but (2) its reflected component from the end of the flume had not yet arrived at WG4.

Force coefficients
The drag and inertia coefficients were determined with a least-square fit method, using the depth-averaged velocity and acceleration in Eq. (1) to reproduce the measured forces (Sumer and Fredsoe, 1997). This approach thus required reconstructing the full velocity profile from the different experiments to estimate the depth-averaged quantities.
Although the instruments were automatically synchronized by the data logger during each test, the velocity measurements at the different elevations (z = 0.15 m, 0.25 m, and 0.4 m) were collected during separate experiments. Combining those measurements to obtain the vertical velocity profile required correcting for the relative time shift between tests, to ensure that the velocities along the vertical coordinate corresponded with the same phase of the wave. The time shift was calculated by maximizing the correlation between the time series of WG2 for each test with respect to the reference case, which was taken as the test with z = 0.15 m. An example of the velocity measurements before and after correcting for the time shift is shown in Fig. A.15 of the Appendix. A moving average was applied to the velocity time series over intervals of 0.25 s. For the velocity measurements of 1C, the mean flow component was removed using the detrend function in Matlab. The acceleration time series was computed from the time derivative of the velocities.
Prior to calculating the depth-averaged quantities, we estimated the values of the velocity and acceleration at the flume bottom z = 0 and at the mean water level z = h to include the velocity changes throughout the whole water column. We extrapolated those values from a hyperbolic cosine fit through the measurements at z = 0.15,0.25, and 0.4 m. The reconstructed velocity profiles from the EMF measurements of 1C are shown in Fig. A.16 of the Appendix. The hyperbolic cosine fit was made assuming that the vertical profile can be described by the main harmonic, since we generated regular waves in the flume. We also calculated the wave spectrum for all wave conditions, and reconstructed the vertical profile by adding the velocity of each harmonic. The maximum differences were up to 2% at z = h.
The depth-averaged quantities were then calculated by trapezoidal integration over the vertical. A moving average was applied over intervals of 0.25 s to the force time series of the reference case, z = 0.15 m. The drag and inertia coefficients were determined by fitting Eq. (1) over an interval of 4 wave periods in order to minimize spurious effects. The interval length of the moving average was varied to evaluate how it affected the fitted drag coefficients for 1C. Averaging over intervals of

Analytical model development
Optimal structure designs are characterized by low reflection (K r ), high dissipation (K d ), and low material costs (so the lowest number of elements possible). We consequently present a simplified conceptual model to investigate how wave reflection and dissipation vary with different cylinder arrangements, in order to identify the most optimal designs.
Wave reflection through a single row of cylinders is calculated as a function of its lateral spacing s y compared to the cylinder diameter d. If the cylinders are so close that they touch each other, i.e. s y − d = 0, the incoming wave height is fully reflected, with K r = 1. In the extreme case where the lateral separation is infinite, the reflection coefficient is zero. We consequently model wave reflection with a function of the form: where c R is an empirical coefficient. Eq. (5) provides the reflection coefficient of the first row. The energy flux entering the first row after subtracting wave reflection is thus given by (1 − K r 2 )E i c g,i , where E i c g,i is the incoming wave energy flux seawards from the structure. The energy dissipated due to the drag forces acting on the first row of cylinders is obtained by introducing Eq. (3) in the wave energy flux balance (Dalrymple et al., 1984): where x is the horizontal coordinate in the direction of wave propagation. Expressing the balance as a function of the wave height results in Eq. (7): where . Solving the linear differential equation for the wave height results in Eq. (8) (Dalrymple et al., 1984): where K t is the wave transmission coefficient through the first row of cylinders, expressed as the ratio of the wave height just downstream of the row, H, to the incoming wave height, H o . s x represents the separation between two rows of cylinders center-to-center in the wave direction, and the damping factor α is given by The bulk drag coefficient used in Eq.
(3) is estimated as the product of the drag coefficient of a single cylinder, c D , times an empirical characteristic velocity, u c , representative of blockage and sheltering effects inside an array, and divided by the undisturbed flow velocity u obtained with linear wave theory: Fig. 4. Configurations tested in the experiments. The location of the ADV for each configuration is indicated by a grey cross, while the location of the force measurements is shown with a red dot. The configurations were tested by starting with a single row, and adding downstream rows in successive steps. All the configurations were tested with a vertical orientation, while the configurations with blue asterisks were also tested with a horizontal cylinder orientation. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) The subscript w denotes that since this bulk drag coefficient is implemented in the wave dissipation term, it relates to the local velocity to the power of three. The characteristic velocity is estimated as a function of several empirical factors: where the blockage factor f b is based on mass conservation through a cross-section of the structure (Etminan et al., 2017;van Rooijen et al., 2018;Etminan et al., 2019), resulting in Eq. (11): f KC is an empirical factor representing the transition of the drag coefficient between inertia and drag dominated conditions, e.g. as shown in Fig. 9 of Etminan et al. (2019).
The right term, f s , representing sheltering effects, is computed using Eq. (12): Here we assume that for highly turbulent environments the velocity reduction in the wake of a cylinder is proportional to 1 − c s /(s x /d), as shown by Eames et al. (2011) for uniform flow, where c s is an empirical parameter dependent on the turbulent intensity. This approach has also been successfully applied to predict sheltering effects for dense cylinder arrays in currents (Gijón Mancheño et al., 2021). Eqs. (6)-(12) estimate the wave dissipation caused by each row. The total wave reflection and dissipation rates of the structure can be calculated cumulatively row by row, by (1) firstly calculating wave reflection, (2) subtracting the reflected energy flux, and (3) calculating wave dissipation between each row and the downstream one.

Wave transmission
Wave transmission through the different configurations is shown in Fig. 5 as a function of KC. The transmission coefficient K t is defined as the ratio of the transmitted wave amplitude to the incoming wave amplitude. The transmission measurements range between K t = 0.4 − 1.
Overall, the transmission rates decrease for longer waves, associated with higher KC values, and for an increasing number of rows for each configuration. However, most wave height reduction takes place in the first rows of the structures.
The influence of element density on wave transmission is more pronounced than the effect of wave excursion. The densest structure, DD13R, produces the lowest transmission rates (Fig. 5, a), and the most porous structure, SS7R, the highest transmission rates (Fig. 5, d). Nevertheless, element arrangement plays an important role on the wave height reduction per cylinder. For instance, the configuration formed by dense rows with a relatively longer streamwise separation, DS7R, (Fig. 5, b), has half as many elements as the least porous configuration, but their wave transmission rates are similar. The results of Fig. 5 (b) and (d) correspond with structures that have the same number of elements, but wave transmission is higher through the staggered arrangement, SS13T, (Fig. 5 d)   , s x = s y = 3d in a regular arrangement. We also showed the measurements with a horizontal orientation for each full configuration (light blue dotted lines). The plots show that the transmission rates are mostly influenced by the structure configuration, rather than the wave excursion, and that most wave attenuation takes place on the first rows of the structure. The results also show higher wave height reduction for the horizontal arrangements compared to the vertical orientations, especially for smaller wave periods. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) separation (Fig. 5 b). Horizontal arrangements reduce wave transmission compared to vertical configurations, as shown in Fig. 5. The additional wave height reduction is largest for the smallest KC, with horizontal configurations having transmission coefficients 10 − 20% smaller, whereas for KC > 15 horizontal and vertical arrangements show similar wave transmission rates. Since the frames have the same frontal and volumetric porosity for both orientations (resulting in the same frontal area), and the measurements were collected for intermediate water conditions (for which the vertical velocities are still significant), the additional wave height reduction is likely due to the work done by the vertical drag forces (Suzuki et al., 2019).

Wave reflection
Wave height reduction behind the structures is partly due to wave reflection. The wave reflection rates for the different configurations are illustrated in Fig. 6, where K r represents the ratio of the reflected wave amplitude to the incoming wave amplitude. The results of 1C (single cylinder) are representative of wave reflection from the end of the flume, and show values oscillating between K r = 0.02 − 0.07. We calculated K r before the propagating wave reached WG4 after being partly reflected at the wave absorber. However, small oscillations were generated by the wave maker in the beginning of the experiments, which explain the reflection rates observed for 1C.
The reflection rates in front of the structures vary between K r = 0.05 − 0.4. Wave reflection in front of the cylinder arrays increases with KC until KC = 15. Beyond this KC value, the reflection rates show a slight decrease for most configurations except for DS7R (Fig. 6, b). The highest reflection rates are measured for the least porous configuration, DD13R ( Fig. 6, a), and the smallest reflection rates for the most porous configuration, SS7R (Fig. 6, d).
Wave reflection also varies depending on the cylinder arrangement. Using the same frontal area but increasing the number of rows (Fig. 6 a  and 6 b respectively), results in higher wave reflection. This is partly due to the increase in the number of cylinders. However, DS5R and DD5R have the same number of elements and the same frontal area, while DS5R (with a longer streamwise spacing, and thus a longer structure width) experiences lower reflection rates. This suggests that increasing the streamwise spacing, and the structure width in the direction of wave propagation, reduces wave reflection. Staggering the elements also reduces the reflection rates, as it can be observed by comparing Fig. 6 (b) and (c).

Velocities and forces
Wave dissipation inside the structures is caused by the work done by the forces acting on the elements. We consequently investigated the magnitude of the forces and velocities measured inside different configurations in Fig. 7. The measured velocity signals are asymmetrical in all cases, with larger negative than positive velocities (Fig. 7 a and b). This is caused by asymmetric placement of the sensor, as illustrated in the upper sketches of Fig. 7. During the positive velocities the sensor measured the flow before it accelerated between the elements, whereas during the negative velocities it experienced the jet formed between the elements. Return currents could also increase the negative velocities, but the pronounced negative asymmetry is not observed in the measurements of the EMF, placed 0.4 m upstream from the frame. The negative velocities are thus indicative of how much the flow accelerates through the spacing between cylinders, whereas the positive velocities do not  Fig. 7. Forces and velocities measured for a single cylinder (1C), a single row (D) and a full structure (DS7R) for KC = 10 and KC = 21. The upper sketches show how the placement of the sensor affects the velocity measurements. During the positive part of the wave cycle (when the velocities were in the direction of propagation), the sensor measured the flow before it accelerated between the cylinders. During the negative part of the cycle, it received the jet formed between the elements. The lower sketches show the location of the force (coloured dots) and velocity measurements (crosses) for the different configurations. The flow velocity increases between multiple cylinders compared to the the case with a single cylinder (a,b). The higher velocities also increase the forces acting on the elements (c,d). include blockage effects. For KC = 10 (Fig. 7 a) the negative velocities for the single row are 2.5 times larger than for a single cylinder. The first row of DS7R has similar negative velocities to the single row (D). Further into the structure the negative velocities reduce due to wave attenuation. In the last row, the negative velocities are actually comparable to those measured for a single cylinder. The velocities for KC = 21 (Fig. 7 b) show a similar behaviour for the different configurations, but the increase in the negative velocities for the first row of the structure is smaller than for KC = 10, with velocities being a factor of 2 times larger than for a single cylinder.
The force measurements are shown in Fig. 7 (c-d). The force signal for KC = 10 (Fig. 7 c) has an almost 90 • phase difference with the velocity, indicating inertia-dominated conditions. This is further shown in Fig. 8, where the force coefficients were fitted for 1C (single cylinder) and D (single row with s y = 1.5d), and used to estimate the contribution of the drag and inertia components to the total force. For both 1C (Fig. 8  a) and D (Fig. 8 c) with KC = 10 the inertia force is almost equal to the total force. The forces for the single row of cylinders and the first row of the structure are approximately 2 times larger than for a single cylinder, as shown in Fig. 7 (c). The relationship between the increase in velocity and the increase in the forces is thus close to linear, which is consistent with inertia forces being linearly proportional to the acceleration.
The force signal for KC = 21 (Fig. 7 d) is in phase with the velocity, indicating drag-dominated conditions. This can also be seen in Fig. 8 (b) and (d), where the the drag component governs the total force. The forces for a single row and the first row of the structure with KC = 21 are approximately 2 times larger than for a single cylinder (Fig. 7 d). If the drag forces were fully driven by blockage between the cylinders, the factor of 2 in the velocities would result in a factor of 4 in the forces, whereas the ratio we measure is smaller. Using the velocities between the elements in Eq. (1) would consequently overpredict the drag forces for the present configurations and KC range. A similar behaviour is also found in the model results of Etminan et al. (2019). Etminan et al. (2019) observed that bulk drag coefficients of cylinder arrays, which include the effect of the velocity changes on the drag forces, increase from the value of a single cylinder from KC = 10, until larger drag values between KC = 20 − 60. For higher KC numbers, the drag forces are well represented by the velocities between the cylinders due to mass conservation. Our measurements fall on their intermediate KC range, where blockage increases the drag forces compared to a single cylinder, but its effect is still reduced compared to KC > 20 − 60. Parameterizing the characteristic velocities for the drag forces in this KC range consequently requires applying a reduction coefficient to the velocities from mass conservation.

Drag coefficients 4.4.1. Drag coefficients from forces
The bulk drag coefficients based on the undisturbed velocities from the EMF of 1C, are shown in Fig. 9. The drag coefficients of both the single cylinder and the single sparse row (in Fig. 9 c and d) correspond well with drag values for a single cylinder from the literature, which decrease from c D ≈ 2 for KC = 10, to c D ≈ 1.7 for KC = 21 (Keulegan and Carpenter, 1958;Graham, 1980;Obasaju et al., 1988). The single denser row (in Fig. 9 a and b), has drag higher values, approximately 2.5 times larger than for a single cylinder. These larger drag coefficients are likely due to blockage effects, as the flow contracts through the small openings between the cylinders.
The drag coefficients of the most porous configuration (SS7R) and the staggered structure (SS13T) are similar to the drag coefficient of a single cylinder, as shown in Fig. 9 (c) and (d) respectively. The drag coefficients are higher for the structure formed by rows with a small lateral spacing and a relatively longer streamwise spacing (DS7R, shown in Fig. 9 b), with c D,b = 2 − 4 at the first row. The drag coefficient decreases at the middle and last rows of DS7R, since the undisturbed velocities do not include wave attenuation through the structure. The least porous configuration (SS7R, shown in Fig. 9 a) experiences smaller drag coefficients than DS7R. This can be partially explained by the higher reflection rates of the least porous structure. However, DS7R and DD9R  have similar reflection rates for KC = 21, while the drag coefficient of DS7R is twice as large. This suggests that sheltering of downstream rows reduces the forces acting on the cylinders. Sheltering effects thus decrease the work done per element, explaining the higher wave reduction efficiency of DS7R.

Predicting wave transmission
The expressions for wave dissipation (Eqs. (3) and (4)) are based on the assumption of undisturbed flow, which does not hold for some of the configurations tested in the present study. However, we wondered whether these expressions could still provide reasonable predictions if blockage and sheltering effects inside the structure are accounted for. Etminan et al. (2019) and van Rooijen et al. (2018) observed that the bulk drag coefficient for arrays mimicking natural vegetation was well represented by mass conservation through a cross-section of the array, as shown in Eq. (13): where c D is the drag coefficient of a single cylinder, A is the total area of the cross-section of the flume, and A c is the total available flow area between the cylinders. The wave transmission predictions obtained by using the bulk drag from Eq. (13) are shown in Fig. 10 (a). We predicted wave transmission with Eq. (3) for vertical orientations and Eq. (4) for horizontal elements. Measurements and predictions are compared for the configurations tested with both vertical and horizontal orientations. The measured transmission rates were corrected to exclude the effect of wave reflection, which is not calculated by Eqs. (3) and (4). Using the drag coefficient derived from mass conservation underpredicts the measured transmission (Fig. 10 a), indicating that the work done by the cylinders is overpredicted.
We also estimated wave transmission using the bulk drag values derived from the force measurements, illustrated in Fig. 9. These predictions are shown in Fig. 10 (b). For each configuration we used the  13)). (b) Results obtained using the c D,b values from the forces measured at the first row of the structure (from Fig. 9). (c) Results obtained using the constrained velocities in the wave energy dissipation rate. (d) Results obtained using the characteristic drag velocities u c from Eq. (17) in the wave energy dissipation rate. Using mass conservation to predict wave dissipation (a and c) overestimated the wave height reduction. The best agreement with the measurements was obtained deriving the characteristic velocity from the drag forces, and introducing it in the wave dissipation with a power of 3 (d).
bulk drag coefficients fitted to the first row of the structure. The bulk drag measurements of the middle and back rows were calculated as a function of the undisturbed velocities, and they consequently include the effect of wave attenuation through the structure. Implementing them in Eqs. (3) and (4) could thus result in an underprediction of the wave dissipation. Using c D,b values derived from the forces to predict wave transmission provides a better agreement with the measurements, but it leads to an overestimation of the measured wave transmission, as it can be observed in Fig. 10 (b). Chen et al. (2018) also overpredicted wave transmission measurements when they used the same approach. This is discussed further in Section 4.5.1.

Relating drag coefficient fitted from forces and dissipation
The bulk drag coefficients derived in the present study (Fig. 9) were fitted to the forces and they are consequently related to the undisturbed velocity to the power of two: However, when the bulk drag coefficient is used as an empirical factor to reproduce wave transmission measurements, it relates to the undisturbed velocity to the power of three: where the subscript w denotes that the empirical drag coefficient is related to the wave dissipation rate. Considering the previous relationships, by using a bulk drag coefficient derived from force measurements to estimate wave dissipation we might underestimate the effect of the velocity changes inside the structure. In order to account for the power of the velocity, we related the measured forces to the drag coefficient of a single cylinder, c D , and replaced the undisturbed velocity u by the characteristic drag velocity u c , which includes the effect of the structure on the flow. The relationship between u c and c D,b is given by Eq. (16): Solving for u c results in: Eq. (17) expresses the bulk drag as a factor that multiplies the undisturbed velocity. Assuming that the magnitude of the characteristic velocity is the same for wave dissipation, but taking into account that the dissipation is proportional to the velocity to the third power, results in: This formulation expresses the characteristic velocity inside a canopy as an empirical drag coefficient that can be included in the dissipation term ϵ v . Introducing Eqs. (18) in Eq. (3) would result in: And for Eq. (4) it results in: The wave transmission predictions from Eqs. (19) and (20) are shown in Fig. 10 (c) and (d). In Fig. 10 (c) we estimated the characteristic drag velocity u c from mass conservation in a cross-section of the array, with u c /u = A/A c , where A is the total cross-section of the flume, and A c the available flow area between the cylinders. In Fig. 10 (d) we estimated the characteristic drag velocity u c from the the bulk drag measurements derived from the forces, using Eq. (17).
Using the velocities due to mass conservation to the power of three also underpredicts the wave transmission measurements (Fig. 10, c). The best agreement between predictions and measurements is obtained when using the empirical characteristic drag velocity to the power of three (Fig. 10, d). These results suggest that the bulk drag coefficients derived from wave transmission measurements and those derived from force measurements are related, but they are not directly exchangeable. Using bulk drag coefficients from forces to predict wave dissipation requires expressing c D,b as a characteristic drag velocity (as done in Eq. (17)), and introducing it in the dissipation rate to the power of 3. The results also show that the expression of Suzuki et al. (2019) as a function of the characteristic velocity provides a good agreement with the transmission rates observed for horizontal arrangements. This agreement between observations and predictions supports that the additional dissipation observed for horizontal structures is caused by the work done by the vertical drag forces, and that the drag coefficient does not experience large changes when varying the cylinder orientation.

Model results
Applying the empirical model presented in Section 3.2.3 requires defining the empirical coefficients for the drag forces, f KC , and c s , and the empirical coefficient for wave reflection, c R .
f KC was calculated using a linear fit through the laboratory measurements, resulting in Eq. (11): We fitted c s such that we could reproduce the sheltering effects observed in the present experiments. The bulk drag coefficients of downstream rows include the effect of wave attenuation, and using them would overestimate sheltering effects. However, due to flow reversal under waves the elements of the first row also experience sheltering during half of the wave cycle. c s is thus obtained by calculating the ratio of the c D,b value measured at the first row of the full configurations (SS7R, DS7R and DD13R) to c D,b values of the single rows (S and D), resulting in c s = 0.796.
The factor c R was obtained by fitting Eq. (5) to the reflection measurements of configurations S and D, resulting in c R = 41.81. We subtracted the energy reflected from the end of the flume from the measurements, since it is not accounted for in Eq. (5).
The model can reproduce the trends observed in the reflection and transmission measurements, as shown in Fig. 11. The maximum differences between modelled and measured wave heights are 0.019 m for the transmitted components, and 0.020 m for the reflected components. The deviations for different wave periods are probably linked to neglecting the influence of the wave length on sheltering and wave reflection.
Including sheltering and reflection is important for design optimization, since both processes influence how wave dissipation varies with cylinder density. For instance, the original formulation of Dalrymple et al. (1984) gives lower transmission rates for higher cylinder densities. The inclusion of blockage in their formulation would enhance this trend further, since higher densities would also lead to higher dissipation per element. Sheltering would have the opposite effect, reducing wave dissipation per element if the streamwise separation s x is small enough. Wave reflection would also decrease the wave energy available for dissipation. However, the relative influence of the previous processes is not known. We thus compared the effect of wave reflection and sheltering on wave transmission, and illustrated the results for s x = 1.5d and s x = 10d in Fig. 12. We analyzed structures with a total width (in the wave direction) of b = 0.76 m, a cylinder diameter of d = 0.04 m, and varying lateral spacing s y , with s y /d = 1.1 − 10. The wave conditions were set to H = 0.13 m, T = 3 s, h = 0.55 m, which corresponded with KC = 21. Fig. 12 (a) shows that both sheltering and wave reflection have a small effect on wave transmission for s x = 10d. Wave reflection becomes non-negligible for densities larger than N = 20 elements per m 2 (Fig. 12  c), but it also has a small effect on wave dissipation (Fig. 12 e). The effects of sheltering are more pronounced for s x = 1.5d, with the dissipation being reduced almost by half for N = 100 elements per m 2 (Fig. 12 f). The influence of wave reflection on wave dissipation is larger for s x = 1.5d than for s x = 10d, but still smaller than sheltering. This comparison suggests that accurate descriptions of sheltering may be important to predict wave transmission. Including wave reflection has a relatively smaller effect on the wave transmission predictions, but precisely assessing its magnitude is necessary to ensure the stability of the designs.
We also assessed the relative performance of different configurations in Fig. 13, with downstream spacings varying between s x = 1.2 − 10d. The remaining model parameters were set equal to those of Fig. 12. Fig. 13 shows that the same amount of wave transmission can be reached over a fixed structure length with different cylinder densities. For instance K t = 0.5 can be obtained with N v = 45 elements/m 2 (for s x = 10d), and with N v = 308 elements/m 2 (for s x = 1.2d). However, the reflection rate is lower for N v = 45 elements/m 2 , with K r = 0.13, compared to N v = 308 elements/m 2 , with K r = 0.44. Using less elements in sparsely placed rows consequently increases wave dissipation and reduces wave reflection compared to a denser and more homogeneous structure. The same trends were obtained with cylinder diameters of d = 0.02 and d = 0.08.  5), and without including sheltering. The yellow lines are obtained accounting for reflection and including sheltering in the characteristic velocity (Eq. 10). Sheltering effects are very small for the configuration with a large streamwise separation (a,c,e), whereas it reduces the wave dissipation coefficient up to 50% for the configuration with smallest spacing (f). Wave reflection has a relatively smaller effect on wave transmission, but high reflection rates could hinder structure performance in the field. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Limitations of the experimental data
Our work provided insights on the factors affecting the drag coefficients inside dense cylinder arrays within a limited range of conditions, with KC = 10 − 21 and two values of cylinder spacing. Obtaining measurements for a wider range of KC values and spacings is recommended to develop more generic parameterizations. Moreover, additional physical processes could modify the bulk drag coefficients in the field compared to the values of this study.
The combination of waves and local currents could be one factor influencing the bulk drag coefficients. Coexistent currents generally decrease the drag coefficient, with a more pronounced drag reduction for a higher magnitude of the current compared to the orbital velocities (Sumer and Fredsoe, 1997). This effect is attributed to stronger currents sweeping turbulence away from the cylinders, suppressing turbulence enhancement by waves.
Despite the drag coefficient reduction, currents can also enhance wave dissipation by increasing the flow velocities and the total work (Eq. 3). Hu et al. (2014) observed that the generation of wave-driven return currents in pure wave flows increased flow asymmetry, and thus wave dissipation. Their study showed that relatively small currents counteracted the wave-driven return flows, reducing wave dissipation, while large currents increased the total work. Wave-driven currents had a negligible effect on both drag coefficients and wave dissipation for the conditions tested in the present work. Wave-current interaction effects are also expected to be small in Demak, where the structures are placed in shallow waters where wave orbital velocities are one order of magnitude larger than the mean flow. However, this factor could differ at other sites. Element roughness, due to irregularities from the bamboo or barnacle growth, could also influence the drag coefficients. Roughness generally increases the drag coefficient to higher values. However, it can also cause a drag reduction for Re ≈ 10 4 , as shown in Fig. 4.20 of Sumer and Fredsoe (1997). The net effect of roughness on the drag coefficient will thus depend on the local flow and material properties. Changes in diameter due to degradation of the bamboo could also gradually decrease the drag forces on the poles. Our work suggested that c D,w did not change significantly for horizontal elements, but wave dissipation was higher for horizontal arrays than for vertical arrays. The increase in wave dissipation was attributed to the work done by the vertical velocities in relatively deeper water. This additional dissipation term for horizontal elements could also be relevant for modelling aquatic vegetation. Neglecting the vertical drag for horizontal roots, such as those of red mangroves, or for horizontal branches, would lead to having to fit higher values of c D,w to compensate the lack of one dissipation term. However, this process will only be significant for relatively short waves compared to the water depth.

Model limitations
The model presented in Section 3.2.3 can qualitatively reproduce the influence of cylinder arrangement, but it should be further developed for its application in detailed designs. For instance, we assumed that the expression of White (1991) for the drag coefficient of a single cylinder remains applicable for very small s y values. Considering an analogy with a cylinder close to a wall in uniform flow, vortex shedding could be inhibited for very small lateral separations between cylinders (Sumer and Fredsoe, 1997). This would in turn reduce the drag coefficient compared to the values of White (1991), but this process has not been investigated for wave flows.
The wake flow model represented by Eq. (12) does not describe the changes in flow velocity as a function of s y , which would be necessary for Fig. 13. (a) Wave transmission rates, (b) wave reflection rates, and (c) wave dissipation rates for configurations with a downstream spacing between s x /d = 1.2 − 10, plotted with different shades of blue (with lighter colours indicating longer downstream separation). For each s x value the lateral spacing varies between s y /d = 1.1 − 10. The same wave transmission can be achieved with different cylinder densities, but higher densities are associated to higher wave reflection rates. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) modelling staggered and random arrangements. The model is also limited for turbulent flow, since Reynolds numbers in the field are of O ( 10 3 − 10 4 ) . For applications where viscous effects are significant, the velocity deficit in the wake will decrease compared to the results of Eq. 12 (Eames et al., 2011). Moreover, c s values may vary for different structure and wave properties. Predicting c s for any cylinder arrangement, given its geometry and the local wave conditions, requires a turbulence model that reproduces turbulence enhancement by waves. For high KC values, the model developed by Gijón Mancheño et al. (2021) for dense cylinder arrays in a current could be applied to estimate c s . For low KC values, the model should be expanded to include the effect of flow reversal on the turbulent intensity.
Lastly, the present formulation for the wave reflection factor c R neglects the influence of varying wave properties, since our wave transformation measurements were mostly influenced by structure configuration. However, this assumption should be verified for KC values outside the range tested in this study.

Implications for design optimization
Many structures in the field consist of vertical poles with a dense brushwood filling, with volumetric porosities down to n = 0.5. In view of the present work, a more porous structure with a different arrangement could be more optimal for wave dissipation. If waves approach the coastline from a relatively constant direction, placing the elements in dense rows with a relatively longer streamwise spacing could maximize the dissipation per element. If the direction of wave incidence has considerable variability over time, combining several structures with different orientations (based on the most frequent wave directions) or using staggered arrangements may be preferable. Aspects such as the construction procedure or soil properties also influence design optimization, but our results already provide some insights on how to obtain more cost-effective designs.
The effect of the structures on coastal accretion will also depend on the local sediment properties. After additional model development, implementing Eqs. 5-12 in morphodynamic models would enable assessing the impact of wave dissipation on erosion mitigation. Future work will investigate the influence of the structure location on the coastline position, in order to find designs that protect and expand the mangrove habitat.

Conclusions
The present study investigates wave transmission and hydrodynamic forces acting on groups of cylinders with volumetric porosities between n = 0.64 − 0.9. Laboratory experiments show that the forces in groups of cylinders with KC = 10 − 21 are proportional to the blockage effect taking place between the cylinders, but force and wave dissipation predictions would be overestimated by assuming that the drag term is defined by mass conservation only. This seems to suggest that there is a transition region between the behaviour of a single cylinder and to the conditions of steady flow for this KC range.
When we use bulk drag coefficients fitted to force measurements to calculate wave transmission with Dalrymple et al. (1984) we underestimate wave transmission, as also observed by Chen et al. (2018). However introducing the bulk drag coefficient to the power of 3/2, to account for the fact that the velocity goes to the power of three in the dissipation, provides better results. This suggests that bulk drag coefficients from the forces could be applied to estimate wave transmission, but the difference in the velocity scale should be accounted for.
The experiments also compare the efficiency of different structures in terms of wave dissipation. Placing the bamboo poles horizontally instead of vertically can provide additional wave dissipation for deep water conditions. Configurations with a small lateral spacing and a longer streamwise spacing increase blockage and decrease sheltering, maximizing the forces and dissipation per element. Adding more rows and placing them closer together increases reflection and sheltering, reducing the dissipation per element. Staggered arrangements experience less blockage than a configuration with the same number of elements placed in dense rows with large streamwise separation, but they could have an extra advantage at locations with strongly varying wave direction.
We also present a conceptual analytical model to predict wave reflection and dissipation through cylinder arrays, including sheltering and blockage effects. The model can qualitatively reproduce the influence of structure configuration on wave transmission and reflection, and it suggests that accurate predictions of sheltering and wave reflection are important for the optimization of future designs. Overall, these experiments provide useful insights for optimizing the designs of structures for mangrove restoration.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.