*2.1. The Direct Measuring Method for CD Derivation*

The force acting on a single stem can be expressed by Morison equation [23] as

$$F = F\_D + F\_M = \frac{1}{2} \rho \mathbb{C}\_D h\_{\upsilon} h\_{\upsilon} lI|ll| + \frac{\pi}{4} \rho \mathbb{C}\_M h\_{\upsilon} h\_{\upsilon}^2 \frac{\partial lI}{\partial t} \tag{1}$$

where *F* is the total acting force on vegetation which can be obtained by force measurement. *FD* is drag force, and *FM* is inertia force. *ρ* is the density of the fluid. *CD* and *CM* are the drag and inertia coefficients, respectively. *hv* is the height of vegetation in water, and *bv* is the diameter of circular cylinder. *U* is the depth-averaged flow velocity. When linear wave theory is applied, the *U* varies as a function of sine:

$$
\mathcal{U}I = \mathcal{U}\_{\text{iv}} \sin(\omega t) \tag{2}
$$

where *Uw* is the amplitude of horizontal wave orbital velocity. Following linear wave theory, *Uw* can be expressed as

$$dJ\_w = \frac{\pi H}{T} \frac{\cos \mathbf{h} [k(h + Z\_0)]}{\sin \mathbf{h} (kh)} \tag{3}$$

where *H* is wave height, *T* is wave period, *k* is wave number, *h* is water depth, and *Z*<sup>0</sup> is the vertical position of the considered point, which is 0 at the still wave level, and −*h* at the sea bed. Equation (3) was used to estimate *Uw* when the velocity measurement is unavailable. *CM* is often assumed to be equal to 2 for cylinders (e.g., [39]). To derive the time-varying *CD*, we can apply the following equation:

$$\mathcal{C}\_{D} = \frac{2F\_{D}}{\rho h\_{v}b\_{v}lI|\mathcal{U}|} = \frac{2(F - F\_{M})}{\rho h\_{v}b\_{v}lI|\mathcal{U}|} \tag{4}$$

where *FM* can be derived based on *<sup>∂</sup><sup>U</sup> <sup>∂</sup><sup>t</sup>* using the timeseries of velocity data, and other parameters (i.e., *ρ*, *bv*, *hv*, *CM*) in *FM* are known. Thus, time-varying *CD* can be obtained readily when in-phase force and velocity data is obtained.

Period-averaged *CD* is relevant to vegetation-induced wave dissipation. It was not computed as the temporal mean of the time-varying *CD*. Time-varying *CD* has great variability over one wave period [37]. Specifically, its value is infinite when the velocity is close to zero. However, those *CD* values are not relevant for vegetation-induced wave energy dissipation, as dissipation is highest at the velocity peaks. Thus, *CD* values at high velocity matter the most. To obtain relevant period-averaged *CD* values, the direct measurement method applies the technique of quantifying the power and work done by the acting force (*ε*) [37]. The time-varying power of *FD* and *FM* is evaluated as follows:

$$P\_D = F\_D \mathcal{U} \tag{5}$$

$$P\_M = F\_M \mathcal{U} \tag{6}$$

The work done by the total acting force (*F*) over a wave period (*T*) is

$$\mathcal{W} = \int\_0^T Flldt\tag{7}$$

If we substitute *F* with Equation (1), then we obtain

$$\mathcal{W} = \mathcal{W}\_{\rm D} + \mathcal{W}\_{\rm M} = \frac{1}{T} \int\_{0}^{T} \mathcal{F}\_{\rm D} \mathrm{Id} dt + \frac{1}{T} \int\_{0}^{T} \mathcal{F}\_{\rm M} \mathrm{Id} dt = \frac{1}{2T} \int\_{0}^{T} \rho \mathcal{C}\_{\rm D} h\_{\rm b} b\_{\rm b} \mathrm{Id}^{2} |\mathrm{Id}| dt + \frac{\pi}{4T} \int\_{0}^{T} \rho \mathcal{C}\_{\rm M} h\_{\rm b} b\_{\rm b}^{2} \frac{\partial \mathcal{U}}{\partial t} \mathrm{Id} dt \tag{8}$$

*WD* and *WM* is the work done by *FD* and *FM* over a full period, respectively. As *U* is a sine function (Equation (2)), the work done by *FM* over a full wave period (i.e., second term on the right) is zero. Thus, the work done by *F* is equal to the work done by *FD*:

$$\mathcal{W} = \int\_0^T Fllt = \frac{1}{2} \int\_0^T \rho \mathbf{C}\_D h\_v b\_v l I^2 |\mathcal{U}| dt \tag{9}$$

Finally, the period-averaged *CD* can be derived based on the above equation:

$$\mathbb{C}\_{D} = \frac{2\int\_{0}^{T} F\_{\text{D}} l I dt}{\int\_{0}^{T} \rho h\_{\text{v}} b\_{\text{v}} l I^{2} |\mathcal{U}| dt} = \frac{2\int\_{0}^{T} F l I dt}{\int\_{0}^{T} \rho h\_{\text{v}} b\_{\text{v}} l I^{2} |\mathcal{U}| dt} \tag{10}$$

The in-phase time series data of total force (*F*) and velocity (*U*) can be used directly in Equation (10) to drive period-averaged *CD* values. As *W* is proportional to *U*<sup>3</sup> (Equation (9)), the integration of *FU* over a period is largely contributed to by the moments with relatively high velocity, and to a very limited extent, by the moments with low velocity. Thus, deriving period-averaged *CD* via the technique of quantifying *ε* can automatically assign large weight to the moments with high velocities in a wave period, resulting in most relevant *CD* values for wave dissipation analysis. To check the validity of the direct measuring method, we used the derived period-averaged *CD* values to reproduce the total acting forcing (*Frep*) using Equation (1), and compare it with the actual measurement. Additionally, another reproduced total force *Frep'* is included by assuming *CD* = 1. It is used as a reference for the *Frep*.

#### *2.2. Synchronized Force–Velocity Measuring System*

In-phase force–velocity data are critical to the direct measurement method. To obtain in-phase data, a synchronized force–velocity measuring system was developed, which was composed by a force sensor and an acoustic doppler velocimeter (ADV) (Figure 1a,d). Four measuring systems were deployed in the wave flume at Fluid Mechanics Laboratory at Sun Yat-sen University, Zhuhai Campus (Figure 1a). The wave flume is 20 m long, 0.8 m wide, and 0.6 m deep. A series of capacitance-type wave gauges were installed to monitor wave height changes in the wave flume. The mimicked vegetation canopy was 8 m long, and it was constructed by PVC pipes. The pipes were 0.2 m tall and their diameter was 0.02 m. The mimicked vegetation canopy was built following a stagger pattern with a density of 139 stems/m2. It was built on top of a false bottom in order to elevate the canopy so that the force sensors can be mounted underneath it. The mimicked vegetation canopy was submerged in water (water depth = 0.25 m), and it was subjected to various wave conditions. The tested wave height varied from 0.03 to 0.09 m, and the tested wave period varied from 0.6 to 1.2 s. A space-averaged *CD* can be obtained by taking the mean *CD* values of the four measuring spots in each test.

The force sensors selected were model M140 built by Utilcell, Spain (Figure 1b). As the size of this sensors is small, they can be installed at multiple locations in one flume test. The reading of the sensor is in "gram", which can be easily translated into "Newton" by multiplying the acceleration of gravity. The minimum division of the sensor is 3 × <sup>10</sup>−<sup>3</sup> N, and the maximum measuring load is 30 N. In the current experiment, the measuring frequency of the force sensors is set as 20 Hz. These sensors were chosen also because they are robust and can be easily waterproofed by sealing the cable connection point with glue. Furthermore, the sensor is small (15 cm × 4 cm × 2.5 cm), and can be easily fitted in the flume (Figure 1b). To prevent the sensor being affected by any force acting on itself, an aluminum case was put around the sensor.

For each force sensor, a mimicked vegetation stem (PVC pipe) was firmly screwed to it, so that the force acting on a vegetation stem can be detected (Figure 1b). The PVC pipes that were attached to the force sensor were not different from other pipes in the mimicked vegetation canopy. After the sensor was attached to a PVC pipe, it measured the total acting force on the pipe. For the same force, the reading is constant, regardless of the location of the acting force. This is desirable for our current experiment, in which wave-induced forces acted over the full length of the pipe. Additionally, the sensors can detect the force in both the following and the opposing direction as the wave propagation, which is ideal to measure the force generated by oscillatory flows.

**Figure 1.** (**a**) Flume experiment set-up. The numbers 1–4 indicate the locations of the in-phase force–velocity measurement; (**b**) A force sensor connected to a PVC pipe; (**c**) Mimicked vegetation canopy constructed by PVC pipes; (**d**) Instrument deployments in the flume without mimicked vegetation. "WG" stands for wave gauge, "ADV" stands for acoustic doppler velocimeter, and "FS" stands for force sensor. The dashed line indicates the ADV and force sensor are placed at the same cross-section to ensure in-phase measurements.

As a first test of the force sensors, we put known weights on one of the sensors attached with PVC pipes. This sensor was held horizontally in this test. The obtained readings were subsequently compared to the known weights. Different weights were put at three different positions on the pipe, i.e., bottom, middle, and tip. The generated forces were in both positive and negative directions. The comparison between the readings and the weights is listed in Table 1. It is clear that for the force measurement in both directions, the relative errors of the sensors are within 1%, compared to the known weights. This first test shows that the force sensors can provide precise measurements.

For the velocity measurement, we used four ADVs. They were deployed at the same cross-section as the force sensors to obtain roughly synchronized signals. The cases wave0312 (3 cm wave height and 1.2 s wave period) and wave0709 (7 cm wave height and 0.9 s wave period) were selected to conduct velocity profile measurements in order to obtain information on velocity structure and depth-averaged in-canopy velocity. The profile was obtained by manually adjusting the vertical measuring location of the ADVs in repeated sequences. For most cases, the velocity measurement was taken at half of the water depth, which was a representative value of the depth-averaged in-canopy velocity. The accuracy of this treatment is acceptable when the submergence ratio (i.e., *h*/*hv* = 1.25) is small [24,37], and it can significantly reduce the labor involved. Two of the ADVs were made by Nortek (*Vectrino* see http://www.nortekusa.com/usa/products/acoustic-doppler-velocimeters/vectrino-1), and the other two were made by SonTek (*MicroADV*, see https://www.sontek.com/argonaut-adv). These four ADVs are common instruments in fluid mechanics labs. Their basis measurement technology is coherent Doppler processing. They measure 3D water velocity of a small cylinder (i.e., within 1 cm3) that is a few centimeters away from measuring probes in the water. They can measure at frequencies as high as 64 Hz, which are desirable for the direct measuring method. In our experiment, the ADV data acquisition followed their respective user manuals. The measuring frequency was set as 40 Hz in order to accommodate the measuring frequency of the force sensor (i.e., 20 Hz). The obtained data are filtered through a low-pass filter to remove high frequency spikes following a similar method described in Strom and Papanicolaou [40].

**Direction Known Weights (g) 1st Reading <sup>a</sup> (g) 2nd Reading <sup>b</sup> (g) 3rd Reading <sup>c</sup> (g) Mean Reading (g) Absolute Error (g) Relative Error** + 5 g 5.10 5.00 5.00 5.03 0.03 0.60% 10 g 10.10 10.10 10.00 10.07 0.07 0.70% 20 g 20.10 20.00 20.00 20.03 0.03 0.15% 50 g 49.90 49.90 49.90 49.90 0.10 0.02% − 5 g −5.00 −5.00 −5.10 −5.03 −0.03 −0.60% 10 g −10.00 −10.10 −10.10 −10.07 −0.07 −0.70% 20 g −20.00 −20.00 −20.00 −20.00 0 0% 50 g −50.0 −49.9 −49.9 −49.93 0.07 0.14%

**Table 1.** Comparison between known weights and force sensor reading.

a,b,c The 1st, 2nd, and 3rd time readings were taking when the weights were put at the bottom, middle, and tip of the testing pipe, respectively.

## *2.3. Automatic Alignment Algorithm*

Although the force and velocity measurements were deployed at the same cross-section of the wave flume, it did not ensure perfectly in-phase data. In fact, small time lags commonly existed between the obtained original force and velocity time series. These time lags were induced by small misalignments between force and velocity measurement, and/or by intrinsic delays of these electronic devices. To reduce the time-lag between velocity and force measurement, the original force and velocity time series should be realigned.

According to the Morison equation [23], velocity (*U*) and drag force (*FD*) should be in phase, which can be used to evaluate the time-lag between velocity and force measurement. A flow chart for the data realignment is shown in Figure 2. The inputs are the timeseries of force (*F*) and velocity (*U*). As we can assume that *CM* =2[39], the inertia force can be calculated based on the velocity. Then, the drag force (*FD*) can be computed by subtracting inertia force (*FM*) from the total force (*F*). Subsequently, we can determine the phase shift (Δt) between the velocity and drag force peaks. Lastly, this phase shift (Δt) will be recorded and used to adjust the velocity timeseries, aiming to obtain more in-phase velocity and force data. The obtained new velocity and force data will be used as input in the same loop. This loop continues 30 times, and we chose the minimum phase shift (Δt) and the resultant velocity and force timeseries as outputs for *CD* derivation. The automatic alignment algorithm is provided in the Appendix A as a MATLAB script. To verify the 30 loop count criterion, a sensitivity analysis is conducted by changing the loop count to 10, 20, 30, and 50. The resulting phase shifts with those loop counts are subsequently compared.

**Figure 2.** Flow chart to realign velocity and force data signals. The algorithm is provided in the Appendix A as a MATLAB script.

#### **3. Results**

#### *3.1. Velocity Profiles in the Vegetation Canopy*

Figure 3 shows the velocity profiles were measured at location 3 in case wave0312 (3 cm wave height and 1.2 s wave period) and wave0709 (7 cm wave height and 0.9 s wave period). For the rest of the tested cases, velocity was measured at the half water depth as a proxy of depth-averaged in-canopy velocity. Figure 3a shows that in the case wave0312, the velocity profiles are rather uniform in the vertical direction. The depth-averaged in-canopy velocity amplitude is 0.052 m/s, whereas *Uw* measured at the half water depth is 0.049 m/s. The difference between these two is small. The velocity profiles in wave0709 have greater vertical gradient, i.e., higher velocity at the top and lower velocity near the bottom. Overall, the difference between *Uw* measured at the half water depth (0.106 m/s) and the amplitude of depth-averaged in-canopy velocity (0.115 m/s) is small. Therefore, it is acceptable to use *Uw* measured at the half water depth as a representative value of the depth-averaged in-canopy velocity.

**Figure 3.** (**a**) Measured velocity profile for case wave0312 with 3 cm wave height and 1.2 s wave period, *Umin* is the highest wave orbital velocity in negative direction (oppose to wave propagation), *Umax* is the highest wave orbital velocity in positive direction (same as wave propagation), and *Uw* is the amplitude of wave orbital velocity; (**b**) Measured velocity profile for case wave0709 with 7 cm wave height and 0.9 s wave period.

#### *3.2. Wave Height and Wave Orbital Velocity in the Mimicked Vegetation Canopy*

Reductions of wave height (*H*) and magnitude of wave orbital velocity (*Uw*) through mimicked vegetation canopy can be observed in Figure 4. The wave height reduces continuously from the canopy front to the end. The final wave height reduction rate was 55% (Figure 4a). The shown wave orbital velocity is obtained by ADV measurement at location 1–3. The ADV measurement at location 4 failed during the experiment. The shown *Uw* is obtained by using Equation (3) based on an average wave height between x = 6–8 m in Figure 4a. With the reduced wave height, the magnitude of wave orbital velocity also reduces from 0.155 m/s to 0.095 m/s from the beginning to the end of the canopy. The reduced wave orbital velocity leads to variations in acting force on vegetation stem as well as in vegetation drag coefficient (*CD*), which are shown in the following sections.

**Figure 4.** (**a**) Spatial variations of wave height (*H*) through the vegetation canopy (i.e., x = 0–8 m) indicated as green bars; (**b**) Spatial variations of the magnitude of wave orbital velocity (*Uw*) through the vegetation canopy. The first 3 *Uw* data points (x = 0–6 m) are obtained from ADV measurement, whereas the last data point at x = 7 m is obtained by using Equation (3) based on an average wave height between x = 6 m and 8 m in panel (**a**). The shown test case is wave0712 with 7 cm wave height and 1.2 s wave period.

### *3.3. Data Alignment and Time-Varying CD*

The total acting force (*F*) and velocity measured at four locations in the wave flume are shown in Figure 5. Waves reach these four locations at different moments. The velocity data at the last measurement point were missing due to the failure of the ADV measurement. It is clear that the acting force (*F*) and velocity reduce as waves pass through mimicked vegetation canopy. It is also apparent that there are time lags in synchronized force–velocity measurements at all locations. However, it should be noted that these seeming time lags (i.e., time shifts between *F* and *U*) are not the real time lags (i.e., time shifts between *FD* and *U*), which lead to errors in *CD*.

**Figure 5.** (**a**–**d**) Raw velocity and total force data measured at four locations (1–4) in the order of the wave propagation in the mimicked vegetation canopy; (**e**) The shaded area in panel (**c**) is blown up in panel (**e**) for detailed analysis, where inertia force and drag force were derived based on non-synchronized data. The time shift (Δt) between the *FD* and *U* is 0.26 s. The shown test case is wave0712 with 7 cm wave height and 1.2 s wave period. The velocity data in panel (**d**) is not included, due to the ADV measurement failure at location 4.

Following the previous study [37], we only tracked the first 2–3 full wave periods after the start-up, but before waves reached the back end of the flume, to avoid possible influence of wave reflection. Figure 5e shows the force and velocity data of the chosen first two wave periods after the start-up. Based on the original data and Equation (1), both *FD* and *FM* can be estimated as shown in Figure 5e. However, it should be noted that as the total force and velocity data are not in phase. The derived *FD* and *FM* are the first estimation. The derived peaks of the *FD* are even higher than the total force, which is not possible. This result highlights the necessity of obtaining in-phase data. Judging from the peaks between *U* and estimated *FD* peaks, the time lag between those two signals is 0.26 s.

In order to eliminate the time lags between these signals, the realignment algorithm was applied to obtain synchronized velocity and force data (Figure 6a–d). It is clear that after the realignment procedure, the time lag between *U* and *FD* is largely reduced (Figure 6a,b), but it cannot be completely eliminated, as small shifts in signal peaks and troughs still exist. Based on the phase shifts of the two peaks and two troughs (as indicated by the red double arrow lines), the time shift of the realigned

*U* and *FD* is significantly reduced to 0.003 s, which is only 1% of the original time shift before the realignment. Note that this optimized time shift is the mean shift of the tested two wave periods (i.e., shifts of two peaks and two troughs). Apart from the reduced time shift, the magnitude of *FD* is also reduced to be lower than the total force (*F*), which is in line with the original Morison equation. It is noted that the peak drag force before alignment is about twice as large as the peak drag force after the alignment. Thus, if the original *F* timeseries were used, the derived drag coefficient would be considerably overestimated.

**Figure 6.** (**a**–**d**) Time-varying *U*, *FD*, *FM*, and *F* data over two wave periods. The gray lines are before realignment, and the red ones are after realignment. The vertical double-arrowed lines in panel (a) and (**b**) indicate the synchronization status before and after the realignment; (**e**) Time-varying *CD* derived based on realigned *U* and *FD* data. The shown test case is wave0712 with 7 cm wave height and 1.2 s wave period.

Based on the synchronized velocity and force data and Equation (4), the time-varying *CD* can be derived (Figure 6e). The time-varying *CD* values vary periodically with the changing velocity signals. The values are small when the velocity is large, but they reach to infinity when the velocity is close to zero. The reason for the unrealistically large *CD* values is because the *FD* are divided by very small velocity values in Equation (4). These unrealistically large *CD* values are not useful in modelling wave dissipation by vegetation, since they are associated with period with very low velocity when the energy dissipation is very limited.

To verify the loop count criterion applied in the realignment algorithm, a sensitivity analysis of the loop is conducted (Figure 7). It is clear that with the increase of loops, the phase shifts can be sufficiently reduced. When 30 loop count is applied, the phase shifts can be reduced to minimum, i.e., 0.003 s for the case 1 and almost zero for the case 2. If the number of loops is increased to 50, the time shifts cannot be further reduced. Thus, the results of sensitivity analysis indicate that applying 30 loop count in the realignment algorithm is a reasonable criterion.

**Figure 7.** Sensitivity analysis of the number of loops used in the realignment algorithm. Case 1 is the case wave0712 with 7 cm wave height and 1.2 s wave period, and case 2 is the case wave0512 with 5 cm wave height and 1.2 s wave period.

#### *3.4. Deriving Period-Averaged CD*

Based on the realigned force–velocity data, we can also derive period-averaged *CD* by quantifying the power and work done, *FD* and *FM*. It is clear that the time-varying power of *FD* is always positive, and its magnitude varies in phase with the velocity magnitude (Figure 8). For ideal sinusoidal velocity signals, *PD* at the velocity peaks should be equal to the troughs. However, in our test case (and in real field conditions), the wave orbital velocity is asymmetrical: higher in the positive direction and lower in the negative direction. Thus, *PD* is larger near the wave peaks and smaller near the wave trough. The difference between peaks and troughs are much more apparent in *PD* compared to the difference in velocity. It is because that *PD* is to the third power of velocity. Small asymmetry in velocity will be greatly magnified in *PD*. The variation of *PM* is different from that of the *PD*. It is clear that *PM* varies between positive and negative values. The zero-crossings in *PM* occur when velocity is zero or when velocity is at its maximum in both directions, i.e., when *FM* is zero in Equation (5). Note that there are small fluctuations at the peaks of the *PD*. These fluctuations may be induced by the small phase shifts between the *FD* and *U* timeseries.

**Figure 8.** Time-varying *U*, *PD*, and *PM* over two wave periods. The shown test case is wave0712, with 7 cm wave height and 1.2 s wave period. The averaged work done by drag force (*WD*) and inertia force (*WM*) over the shown two periods is 2.0 <sup>×</sup> <sup>10</sup>−<sup>3</sup> *<sup>J</sup>* and <sup>−</sup>2.0 <sup>×</sup> <sup>10</sup>−<sup>4</sup> *<sup>J</sup>*, respectively. The derived period-averaged *CD* is 3.25.

By integrating *PD* and *PW* over one wave period, we can obtain the work done by drag force (*WD*) and inertia force (*WM*) as shown in Equation (8). The averaged *WD* and *WM* over the two periods in Figure <sup>8</sup> is 2.0 × <sup>10</sup>−<sup>3</sup> *<sup>J</sup>* and −2.0 × <sup>10</sup>−<sup>4</sup> *<sup>J</sup>*, respectively. In case of ideal sinusoidal velocity signals, *WM* should be exactly zero. Due to the asymmetric wave velocity, the *WM* is not zero, but the magnitude of *WM* is fairly small, i.e., one tenth of the *WD*. Since the magnitude of *WM* is considerably small compared to *WD*, the assumption that *WM* can be ignored in the deriving period-averaged *CD* (in Equation (8)) is still valid, and the period-averaged *CD* in the shown case is derived as 3.25.

Following the same method, period-averaged *CD* at three functional measuring locations of all the tested cases are listed in Table 2. It shows that relatively large deviations in *CD* values exist among different measuring locations. The *CD* values generally increase from locations from 1 to 3. Previous studies have shown that *CD* values increase with the reduced velocity (i.e., Reynolds number) [24,37]. The obtained increase of *CD* values may be related to the reduction of wave orbital velocity from the front to the end of the vegetation canopy, as shown in Figure 4b Thus, the spatial variation in *CD* values is in-line with previous studies. Additionally, it is noted that the cases with larger wave height and wave period (i.e., higher wave orbital velocity) generally lead to smaller spatially averaged *CD*, which is also in agreement with previous studies [24,37].


**Table 2.** Period-averaged *CD* in all the tested cases.
