Linear and nonlinear frequency-domain modelling of oscillatory flow over submerged canopies

An analytical and experimental study of flow velocities within submerged canopies of rigid cylinders under oscillatory flows is presented, providing insights into the momentum transfer mechanisms between the different flow harmonics. The experimental dataset covers an unprecedented wide range of flow amplitudes with in-canopy velocity reductions ranging between 0.2 and 0.8 of the free stream velocity (from inertia- to drag-dominated in-canopy flow). Results from the analytical model with nonlinear drag compare favourably to the experimental data. Having application of theories for free surface waves over canopies in mind, the effects of linearization of the drag are analysed by comparing sinusoidal and nonlinear model predictions. Finally, a unified prediction formula for in-canopy velocities for sinusoidal, velocity-skewed, and velocity-asymmetric free stream velocities is presented. The formula depends on two non-dimensional parameters related to inertia and drag forces, and the unified formula allows for easy assessment of the maximum in-canopy velocity.


Introduction
The implementation of nature-based (green) coastal protection has received considerable attention in recent years.Recent reviews of the topic include Jordan and Fröhle (2022), Schoonees et al. (2019) and Morris et al. (2018), and several practical guidelines for implementing nature-based solutions in flood risk management have been developed (Bridges et al., 2022;CISL, 2022;The European Commission, 2021;The World Bank, 2017).However, nature-based solutions are still characterised by significant uncertainties in their cost and potential impact, leading to a continued demand for research in order to quantify and limit these uncertainties.
In relation to coastal engineering, the vast majority of research has focused on wave attenuation by submerged and emergent vegetation canopies, based on either analytical and numerical modelling (e.g.Cao et al., 2015;Méndez et al., 1999;Van Rooijen et al., 2020) or on direct measurements in the field or in laboratories (e.g.Anderson & Smith, 2014;He et al., 2019;Jadhav et al., 2013;Möller et al., 2014).Much of the experimental work has involved canopies of flexible vegetation (Mullarney & Henderson, 2018).In the present paper, the vegetation is represented as a canopy of rigid cylinders in order to address fundamental aspects of the hydrodynamic processes.
For waves propagating over a vegetated bed, the in-canopy velocities are important for wave attenuation and sediment dynamics (Jacobsen et al., 2019;Jadhav et al., 2013).Velocities within the canopy are lower than the wave-generated free-stream velocity above the canopy.Research thus far has established that the smallest velocity reduction occurs in the inertia-dominated regime corresponding to a 1 /S 1, where a 1 is the free-stream orbital amplitude and S is the average spacing between stems (Lowe et al., 2005;Pujol et al., 2013;Van Rooijen et al., 2020), while higher velocity reduction occurs when a 1 /S > 1, as the flow-canopy interaction becomes increasingly drag dominated.This work extends the existing empirical evidence with experiments involving high free-stream velocities reaching well into the drag-dominated regime, and investigates a closed form prediction formulation for the in-canopy velocity reduction based on non-dimensional quantities.
Four approaches can be identified for determining in-canopy velocities under wave-driven oscillatory flow, each approach having its specific application and merits: (i) oscillatory timedomain solutions with nonlinear drag and numerical integration of the momentum equation (Lowe et al., 2005;Zeller et al., 2015); (ii) analytical solutions for free surface wave propagation in submerged and emergent vegetation (Asano et al., 1992;Dubi & Tørum, 1994;Jacobsen, 2016;Kobayashi et al., 1993;Méndez et al., 1999), for which a linearization of the drag term is used through the Lorentz' principal of identical work (e.g.Sollitt & Cross, 1972); (iii) simulation of the flow field using complex depth-resolving models (Chen & Zou, 2019;Van Rooijen et al., 2020); and (iv) laboratory or field measurements (Lowe et al., 2005;Pujol et al., 2013;Van Rooijen et al., 2020).With respect to approach (ii), it is important to note that the validity of linearizing the drag has not been properly addressed in the literature, even though the linearization leads to analytical expressions for radiation stress tensors (Jacobsen, 2016;Mendez et al., 1998), Lagrangian Stokes drift (Jacobsen, 2016;Jacobsen & McFall, 2022), and vertical distributions of shear stress (Jacobsen & McFall, 2022), all of which can be readily implemented in large-scale hydrodynamic models using the procedure outlined in Jacobsen and McFall (2022) for coupling between spectral wave models and large-scale hydrodynamic models based on the shallow water equations.Models for the mean wave-induced current in canopies have been proposed by Luhar et al. (2010) and Van Rooijen et al. (2020), for which the in-canopy velocity reduction is of implicit importance.
The present work examines the oscillatory flow velocity within a rigid canopy, with particular focus on the formulation for the drag force in the governing equations.This is achieved through a combination of analytical modelling and laboratory experiments in a large-scale oscillatory flow tunnel.Two frequency-domain models, one with nonlinear drag and one with a sinusoidal velocity response akin to Airy wave theory in which the drag is linearized, are described in Section 2 and the experiments are presented in Section 3. The analytical models are validated against the experimental data in Section 4.1 and subsequently the models are compared for a wider range of hydrodynamic and canopy conditions than covered by the experiments.A practical engineering prediction model for incanopy velocity reduction is proposed in Section 4.4.Section 5 discusses implications of the results for modelling nonlinear surface waves in the presence of canopies and for sediment transport.The conclusions are presented in Section 6.

Governing equations
In the present work, pure oscillatory (tunnel) spatially-averaged flow is considered, which means there is no organized, vertical velocity (w = 0 m s −1 ).It is furthermore assumed in the analytical derivation that the flow is unbounded along the positive and vertical z-axis and that the canopy stems are rigid.The theoretical domain is depicted in Fig. 1, where the spatially-averaged velocity profile (blue line) and the associated 2-layer solution to the bulk velocities (red lines) are shown.The authors arrive at governing equations of the same form as Lowe et al. (2005), Luhar et al. (2010), Zeller et al. (2015) (with minor differences in the canopy shear term).The present work applies these equations in the frequency domain to investigate the required number of harmonics for an accurate solution and to validate against experiments with large free stream velocity amplitudes.
The pure oscillatory and spatially averaged flow along the xaxis means that ∂/∂x = ∂/∂y = 0 when applied to any quantity.The continuity equation thus reads: which means w = 0 for all z in order to satisfy the no-slip condition at the bottom.Inserting w = 0 in the vertical momentum equation for the flow in porous media (e.g.Jensen et al., 2014), it follows that ∂p/∂z = 0, i.e. the pressure is hydrostatic over the vertical.Finally, the horizontal momentum equation takes the following form (see Appendix A.1): being solidity, as used by other authors), ρ is fluid density, ν and ν t are the molecular and eddy viscosities, and F D is the fluid drag.Above the canopy l = N = F D = 0 and n = 1.In the following, we use I for (1 + C m Nd 2 π/4)/n.The filter velocity is the average velocity over the entire cross section (water and stems), while the pore velocity is the average velocity over the water part of the cross section.The relationship reads u = nu p , where subscript p denotes the pore velocity.Far above the canopy (away from the canopy shear layer) Eq. ( 2) reduces to: where u 0 is the free stream velocity.Since ∂p/∂z = 0, substitution of − 1 ρ ∂p/∂x with the prescribed free stream acceleration (∂u 0 /∂t) is also valid within the canopy.The drag on the stems is defined as: where P(n) is a function that accounts for whether filter or pore velocities are applied for the evaluation of the drag (see discussion in Etminan et al., 2019).In the present work, the drag force in the model is based on pore velocity, because the measured drag coefficients are based on measured pore velocities (see Section 3), so P(n) = n −2 .C D is the drag coefficient.Our focus is not on predicting the vertical profile of horizontal velocity, but rather on the overall performance of linear versus nonlinear drag formulations in predicting the bulk incanopy velocity.For this reason, a two-layer solution is adopted, akin to the work by Lowe et al. (2005), where the diffusion term in Eq. ( 2) is replaced by a shear stress at the top of the canopy (the shear at the bottom is proportional to the drag and is absorbed into C D , see Appendix A.2): where U is the bulk in-canopy filter velocity, U = u 0 − U is the velocity difference between the in-canopy and free stream velocities, and C f is a friction factor for the shear layer on top of the canopy.U is defined as the vertical average over the canopy height h v ( Appendix A.2). Note that following Eq.(3), the "bulk" velocity above the canopy equals u 0 .Lowe et al. (2005) described the canopy shear based on the free stream velocity, while the present model describes the shear contribution in terms of the velocity deficit in order to account for the influence of the in-canopy velocity U on canopy shear, as outlined by Zeller et al. (2015).For very short and sparse canopies, the bottom friction may have important contributions to the momentum balance; however, bottom friction is still parameterized as proportional to U|U| without accounting for phase leads between in-canopy velocity and bottom shear stress, so an appropriate correction to C D covers this scenario in the present two-layer model.The bottom boundary layer must be resolved to accurately account for phase lead effects between in-canopy velocity and bottom shear stress, which is often impossible for large-scale hydrodynamic models based on the shallow-water equations.
The following non-dimensional parameters (indicated with a circumflex accent) are introduced: Û = U/u 1 , û0 = u 0 /u 1 , t = tω.The characteristic dimensional quantities are the amplitude of the first harmonic of the free stream velocity, u 1 , and the cyclic frequency ω = 2π/T, where T is the oscillatory flow period. 1 This gives the following non-dimensional form of the momentum equation for Û: which shows that the in-canopy velocity is governed by three non-dimensional parameters: (7) related to inertia ( I ), canopy interface friction ( f ), and drag on the stems ( D ). a 1 = u 1 /ω is the free stream orbital amplitude of the first harmonic.It is seen that the importance of the friction is inversely proportional to the height of the canopy.The parameter Nda 1 ∝ (d/S) • (a 1 /S) is the product of the two governing parameters proposed by Lowe et al. (2005), where S ∝ N −1/2 is the spacing between the stems.In the following derivations, the ˆis omitted; however, it is reintroduced in Section 4 to distinguish between dimensional and non-dimensional quantities.
Non-dimensional solutions to the free stream and in-canopy velocities of the following forms are considered: where c.c. is the complex conjugate and M is the number of harmonics.The formulation for U explicitly allows for phase lags between the harmonics U m .

Sinusoidal solution
Free surface wave models over submerged and emerging canopies based on linear wave theory are often derived with a linearized drag term, which results in sinusoidal flow velocities throughout the water column (Asano et al., 1992;Jacobsen, 2016;Méndez et al., 1999).Analogous to these linearized free surface models, a model for the in-canopy velocity is derived here for M = 1, and it is termed the "sinusoidal model" throughout.We assume that the friction and drag can be linearized in the following fashion: Here, L,f and L,D are the linearized friction and drag parameters.The equality is fulfilled by equalizing the energy dissipation over one wave period (Sollitt & Cross, 1972) separately for friction and drag: where the overline means period-averaging.This linearization based on identical energy dissipation is adopted from free surface wave models, where it corresponds to matching of the wave attenuation across a canopy, and it is employed here in order to ensure that the model can be used in conjunction with the aforementioned free surface models.The parameters L,f and L,D are found numerically.Inserting Eqs ( 8) and (9) in Eq. ( 6), it follows that: Here, subscript "1" means that only the first harmonic (sinusoidal) component is included.Therefore: The final expression shows that the in-canopy velocity is inphase with the free stream velocity for inertia dominated flows (imaginary part of G L goes to 0), while increasing influence of drag and friction gives rise to a phase lead of the in-canopy flow over the free stream velocity (maximum phase lead is π/2, when the real part of G L goes to zero).

Nonlinear solution
To obtain a nonlinear solution, the following Fourier expansions are introduced: The products | U| U = |u 0 − U|(u 0 − U) and |U|U can be given in Fourier terms by evaluating the products of the expansions in Eqs ( 8), (13a) and (13b).These products suggest that there will be an exchange of momentum between the frequencies; the momentum exchange is mathematically identical to the generation of super-and subharmonics in free surface wave theory (Madsen & Fuhrman, 2006).Inserting Eqs ( 8), (13a) and (13b) into Eq.( 6), the following momentum equation for the m th harmonic is found: where * means the complex conjugate and U m = u m − U m .
The indices are such that m = q + l (super-harmonic interaction for all q > 0 and l > 0) and m = |k − p| (sub-harmonic interaction for all k > 0 and p > 0).The function G reads: Introducing the decomposition U m = U r,m + iU i,m , where superscripts r and i refer to the real and imaginary components, and treating U r,m and U i,m as separate variables, Eq. ( 14) can be written as a 2M × 2M matrix system with a nonhomogeneous right-hand side.The matrix system is nonlinear, because the coefficients b q and B q for q = 0, 1, 2, . . .depend on all components in U.The matrix system is solved by successively inserting the solution from previous iterations until converged below a 2-norm of 10 −10 .M = 21 is used throughout the present analysis, and Eqs (13a) and ( 13b) are evaluated with 2M components to incorporate all sub-harmonic contributions.
The momentum equation (Eq.( 6)) is, as already discussed, practically identical to the two-layer models proposed by Lowe et al. (2005) and Zeller et al. (2015).However, the main difference is the solution procedure: where their models were solved in the time-domain, the present sinusoidal and nonlinear models are solved in the frequency domain, which has the advantage that the in-canopy momentum distribution (e.g. . . ) is obtained directly without additional post-processing and the super-and subharmonic interactions are quantified.This advantage will be utilized in Section 4.3, where momentum distribution, momentum exchange, and importance of super-and sub-harmonic terms are quantified for sinusoidal and velocity-skewed free stream velocities.

Experimental set-up
The experiments were performed in the Aberdeen Oscillatory Flow Tunnel (AOFT) (Fig. 2).The test section is 10 m long, 0.3 m wide and 0.75 m tall.A 7 m long false bed was installed at an elevation of 0.25 m above the tunnel floor within the test section with ramps of slope 1:4 at either end, reducing the effective flow depth within the test section to 0.5 m.Stiff, circular PVC rods, with h v = 130 mm and d = 8.3 mm, were fixed to the false bed in a regular geometric pattern over its full 672 O. E. Neshamar et al. Journal of Hydraulic Research Vol. 61, No. 5 (2023) Figure 2 Layout of the AOFT with the rigid cylinder array installed on the false bottom.The force is measured on the test cylinder as described in Neshamar (2022) Figure 3 Pictures of (a) the sparse canopy and (b) the dense canopy length.Two canopy densities were tested: a sparse canopy with N = 579 stems m −2 , shown in Fig. 3a, and a dense canopy with N = 1736 stems m −2 , shown in Fig. 3b.The corresponding porosities are n = 0.969 and n = 0.906, respectively.These array densities were selected to resemble submerged vegetation canopies, for which porosities of 0.9-0.99 were reported by Nepf (2012).One rod in the centre of the canopy ( ∼ 3.5 m from the leading edge of the canopy) was mounted on a force transducer (ATI Nano17 IP68 6-axis load cell) built into the false bed to measure forces from which the hydrodynamic force coefficients could be determined.The load cell measures forces and moments along all three axes and can be operated in three different calibration settings, with sensing ranges of 12-50 N for force and 120-500 N mm for moment; the corresponding measurement resolutions are 0.003-0.013N and 0.016-0.063N mm.
For each canopy density, experiments were conducted for five free stream velocity amplitudes with constant oscillatory flow period T = 6 s.The piston input signal was sinusoidal but the measured free-stream velocities showed slight departure from sinusoidal with higher harmonics present in the flow.The velocities were measured using a Dantec FibreFlow 2component laser Doppler anemometer (LDA), consisting of

Local and spatial averaged velocities
Figure 4 shows the x -y positions of the in-canopy measurement positions; the positions form a row of approximately evenly-spaced points spanning the periodic domain along y.
The effective number of velocity profiles is doubled because of a "mirroring" effect: since the array geometry is symmetric, and the free-stream velocity is approximately sinusoidal, phase-averaged velocities around each cylinder will be opposite in magnitude along x but otherwise identical during the first and second halves of the flow cycle.Consequently, measurements of the "local velocity" u l at a given "real" measurement coordinate can represent an additional measurement at the "mirror" position by inverting the measured velocity time-series and phase-shifting by 180 • .For example, u l recorded at a position immediately in front of a cylinder in the first flow half-cycle is equal to −u l recorded at a point immediately behind the cylinder in the second half-cycle.Since phase-averaged vertical velocities within the array are negligible, in-canopy flow is effectively divided into thin horizontal "layers" with negligible mass transfer between layers.From continuity, the spatially-averaged filter velocity, u e (z, t) (subscript e refers to experimental value), is therefore found by averaging velocities along y at x = −S x /2, where S x is the cylinder centre-to-centre spacing along x.The spatially-averaged filter velocity was obtained by interpolating the velocity measurements (real and mirrored) onto an evenly-spaced grid spanning the periodic domain along y, and averaging the result.This procedure results in an approximation of u e , which assumes the overall flow to be fully periodic and symmetric.Due to the mirroring, the calculated u e is symmetric in time and contains no even-numbered harmonics.Neshamar (2022) describes how this method was validated experimentally, by comparison with measurements made along two different y-transects (at x = −S x /2 and x = −S x /4).
The depth-and spatially-averaged in-canopy filter velocity, U e , is found through vertical integration: whereby the bulk pore velocity follows: U e,p = U e /n.

Force coefficients
Hydrodynamic force on a cylinder is determined from the Morison equation, where the spatially-averaged pore velocity 674 O. E. Neshamar et al. Journal of Hydraulic Research Vol. 61, No. 5 (2023) (u e,p = u e /n) is used as the characteristic velocity: where F x (z) is the in-line force per unit length at elevation z.
Note that force varies with z due to variations in u e,p ( Fig. 10); the coefficients C D and C m are assumed constant with z.
The added mass coefficient C m can, in principle, be estimated directly from force measurements.However, since the associated force is much smaller than the drag force under the present conditions, there is large uncertainty when estimating C m from the measurements.Instead, the established value of C m ≈ 1 is used in the present work (e.g.Sumer & Fredsøe, 1999).
Since the load cell measures moments with a higher resolution compared to the forces, the drag coefficient C D is calculated from the moment measurements.A force in the x-direction produces a moment about the y-axis, and the total moment applied to the load cell M y can be calculated by multiplying the inline force (Eq.( 17)) by the "lever-arm" and integrating vertically, i.e.: where M y,D and M y,m are the components of moment due to drag and added mass, respectively, and z offset = 10 mm is the vertical distance between the front face of the load cell and the tunnel floor.To determine C D , M y,m is first calculated from Eq. (18c), with C m = 1, obtaining the flow acceleration at each z using a central difference method, and performing the numerical integration along z by interpolating linearly between measurement points.The result is subtracted from the total moment measured by the load cell, M y,LC , to produce a "measured" drag-related moment M y,D,LC = M y,LC − M y,m .Using an initial value of C D = 1, the actual value of C D is then determined by matching the "predicted" M y,D time-series (Eq.( 18b)) to the "measured" M y,D,LC time-series using a least-squares method (see e.g.Sumer & Fredsøe, 1999).Figure 5 shows an example comparison between M y measured by the load cell and the calculated M y obtained using the method described above, which for this case (S4) results in a fitted drag coefficient of C D = 1.37.The figure shows that the two curves match very well in both phase and shape, indicating that the given values of C m and C D are accurate.The same method was used to calculate C D for all flow cases, and each case produced a similarly close match between measured and calculated moments as shown for the example case in Fig. 5.
Figure 6 shows the C D values obtained from the moment measurements for all flow cases, presented as a function of Reynolds number, defined as R e = U p,max d/ν, where U p,max is  2019) hypothesized that C D for a cylinder within an array is governed by two effects, namely "blockage" and 'sheltering'.The blockage effect describes the locally increased flow velocity in the constricted space around the cylinder, which results in increased drag on the cylinder.Conversely, sheltering occurs when a cylinder is "hidden" in the wake of an upstream cylinder, and drag is reduced due to a lower incident velocity.Of the two mechanisms, their simulations showed the blockage effect to be dominant, particularly for high canopy densities.To account for the blockage effect, they proposed that C D for a cylinder in an array can be quantified using the constricted velocity, u c , defined as the velocity averaged over the constricted cross-sections around the cylinder.For the array layout shown in Fig. 4, u Figure 7 shows a comparison between the present data and Eq. ( 20).While there is some scatter, the results from the present isolated-cylinder and canopy experiments fall within approximately 10% of predicted values from Eq. ( 20 It should be noted that most wave attenuation studies define an empirical "bulk drag coefficient" C D,b , obtained from the measured wave attenuation using the method first described by Dalrymple et al. (1984).Henry et al. (2015) show examples of different C D,b formulations.In these studies, C D,b tends to decrease significantly with increasing R e , because the calculation does not account for velocity reduction within the canopy (relative to the free-stream velocity), which increases with increasing R e .In contrast, C D in the present work is based directly on the incident in-canopy velocity, and, as Fig. 6 shows, C D does not appear to vary significantly with R e .

Model validation
A comparison between the experimental results and the two analytical models is shown in Fig. 8 for all 10 cases.The measured free stream velocity is used to force the analytical models, where the first three free stream harmonics are included in the nonlinear model and only the first harmonic in the sinusoidal model (Table 1).Based on the laboratory measurements, model calculations use C m = 1, C D = 1.3 for the sparse canopy and C D = 2.3 for the dense canopy, while C f is tuned for each case to achieve the best fit between the experimental data and the nonlinear model.The resulting non-dimensional quantities I , D and f are given in Table 1.It can be observed that the nonlinear model accurately captures both the phase, amplitude and overall temporal variation of the in-canopy velocity.The sinusoidal model is limited to a single harmonic (M = 1) and cannot reproduce the higher-order in-canopy velocity components, but it does accurately capture the phase and magnitude of the incanopy velocity.The accuracy of the models is evaluated by calculating the normalized deviation to the experiments:  2005) differs to the one proposed in the present work, which explains the small differences between δ NL and δ LKM .Identical results are obtained when applying the resistance formulation from Section 2.1 in both nonlinear models.
Figure 8 shows a comparison between measured velocities U and predictions obtained from sinusoidal and nonlinear models.Results from Lowe et al.'s (2005) time-domain model (henceforth LKM, 2005) are also shown in Fig. 8: there is practically no difference between the present nonlinear frequency domain model and LKM (2005), except for a small phase shift and difference in amplitude, which are attributed to the difference in canopy resistance formulation between the two nonlinear models.
The experimental data and nonlinear model predictions show a rapid flow acceleration followed by slow deceleration after the maximum velocity (Fig. 8).This is explained by the fact that the free stream velocity has its maximum at t = 1.5 s while the in-canopy velocity peaks much earlier, which means the external pressure gradient ∂p/∂x = −ρ∂u 0 /∂t remains negative for some time after the maximum in-canopy velocity is reached, providing a forward driving force during this time.
Figure 9 shows the frequency composition of flow velocities for six representative cases (S1, S3, S5, D1, D3, and D5).The 676 O. E. Neshamar et al. Journal of Hydraulic Research Vol. 61, No. 5 (2023) Figure 8 Comparison between the spatially-averaged in-canopy velocity, U, from experiments and predicted values from sinusoidal, nonlinear, and LKM (2005) models.The free stream velocity, u 0 , is shown with axis on the right of each panel.Left: Results for sparse canopy.Right: Results for dense canopy figure shows the magnitude and phase of the 1st, 3rd and 5th harmonic of free-stream and in-canopy velocities for each case, together with model predictions obtained using the sinusoidal and nonlinear models, all normalized by the magnitude of the 1st harmonic of the free-stream velocity.Predictions obtained using LKM (2005) are also shown, and they are practically identical to the results from the nonlinear frequency-domain model.All three models predict the magnitude and phase of the 1st harmonic with relatively good accuracy.For each case, the figure shows that the free-stream velocity contains a 3rd harmonic with magnitude of between 1-5% of the first harmonic, while the 5th harmonic has a magnitude of < 1% of the 1st harmonic.The results demonstrate the nonlinearity introduced by the canopy, since for the first four cases shown, the 3rd harmonic of the in-canopy velocity is increased relative to the free-stream.This nonlinearity is most pronounced for case S5, where the 1st harmonic reduces by |U 1 |/|u 1 | ≈ 0.6 and the 3rd harmonic increases by |U 3 |/|u 3 | ≈ 2.5.For the other two cases (D3 and D5), the 1st harmonic is significantly reduced (|U 1 |/|u 1 | ≤ 0.33) while the 3rd harmonic remains near-constant.The canopy also introduces a 5th harmonic component with magnitude of between 0.5-2% of the 1st harmonic.In all cases, the nonlinear model, as well as LKM (2005), are shown to predict the magnitude and phase of the 3rd harmonic with good accuracy.The 5th harmonic is also approximated, albeit with less accuracy due to its small magnitude.
A comparison between the predicted (sinusoidal and nonlinear models) and measured velocity profiles is shown in Fig. 10.The two-layer models are found to capture the bulk behaviour accurately.The boundary layer formation on top of Figure 10 Measured and predicted vertical velocity profiles.Top panels: Root-mean-square velocity, u rms .Bottom panels: Phase difference between in-canopy and free stream velocities based on the first harmonic only the canopy is not captured, as this would require a higher vertical resolution (e.g.Chen & Zou, 2019;Zeller et al., 2015).The velocity in the canopy leads over the free stream velocity (bottom panels) and the sinusoidal and nonlinear models both predict almost identical phase leads (φ 1 = arg U 1 − arg u 1 ) with good match to the measurements.The experimental measurements did not resolve the wave boundary layer at the bottom of the flow tunnel, though it is hypothesized that the phase lead will further increase towards the bed akin to observations from laminar and turbulent wave boundary layers over smooth and rough beds (Jensen et al., 1989).It is furthermore hypothesized that the additional phase lead will be smaller than that for undisturbed wave boundary layers, because the canopy has already introduced a phase lead with respect to the external pressure gradient.These hypotheses should be addressed experimentally and theoretically, as they might prove important for a description of in-canopy sediment transport.The accuracy of the nonlinear model in terms of predicting the phase lead, momentum distribution, and in-canopy velocity reduction factors is incorporated in the results presented in Sections 4.2-4.4.

Comparison between sinusoidal and nonlinear drag models
In this section, we compare sinusoidal and nonlinear model predictions of in-canopy flow for a wide range of conditions: n = [0.9,0.999], , and C f = {0, 0.02}.The comparison is done for 4700 conditions in total.In all cases, the free-stream velocity is purely sinusoidal (u m = 0 for m = 2, 3, . ..).Two assessment measures are used in the comparison: (i) magnitude of the first harmonic, and (ii) the phase between u 1 and U 1 .
The magnitude of the first harmonic Û1 was found to be within ±5% of each other for all cases (not shown).This is of considerable practical importance, because it suggests that the linearization applied in free surface models is sufficiently accurate (e.g.Jacobsen, 2016;Méndez et al., 1999).It could be that the vertical velocity component in the case of real waves influences the accuracy of the sinusoidal model; however, in the absence of a valid nonlinear wave theory for canopies, it seems to be a good engineering approximation to apply a wave theory based on a linearized drag force.
The phase lead in the oscillatory flow case has limited direct physical implications; however, for real waves the vertical distribution of the phase lead between the horizontal and vertical velocities is important for residual wave-induced stresses (Deigaard, 1993;Deigaard & Fredsøe, 1989;Guannel & Özkan-Haller, 2014;Jacobsen, 2016;Jacobsen & McFall, 2022;Philips, 1980).Since there is no nonlinear free surface wave theory for canopies, the oscillatory solution is used to assess the implications of using a linearized wave theory to calculate period-averaged wave stresses.The two analytical models, with C f = 0.0, are compared in Fig. 11 (top panel) and it is seen that there are some differences between the two.However, the differences are small (less than 3 • ) and are seen only at high D for which the phase lead is already considerable (greater than 45 • ).
Figure 11 (bottom panel) shows the phase leads predicted by the two models with C f = 0.02.It is seen that the models predict lower phase leads for C f = 0.02 compared to C f = 0.00.Furthermore, for large values of D , the results for C f = 0.02 show decreasing phase lead with increasing D , which is attributed to the increasing importance of f for large in-canopy velocity.For D > 2, the predicted phase leads branch off and the different branches are found to coincide with constant values of NC D , where the lowest branches correspond to small values of NC D and the highest branches to large values of NC D .The branching off therefore relates to the relative importance of friction versus in-canopy drag.The experimental phase leads are included in Fig. 11 and are seen to align better with the model predictions for C f = 0.02 than for C f = 0.0.The importance of C f on the in-canopy velocity reduction is addressed below.For the remainder of the work, only the nonlinear model is applied.(Svendsen & Jonsson, 1982), which leads to skewed free stream velocity signals.
The free stream velocity signal has a secondary bump in the trough of the velocity signal for |u 2 |/|u 1 | > 0.25, which is non-physical for a regular wave signal.
Results obtained for C f = 0.0 in Fig. 12 (all panels) collapse well when plotted against D , demonstrating that D is an appropriate descriptive parameter for the higher harmonic content of the in-canopy velocity.The results for C f = 0.02 and |u 2 |/|u 1 | = 0 exhibit the same branching as seen for the phase lead, whereby branches correspond to different values of NC D .The experimental data are also shown (top panel) and there is reasonable correspondence with the model results.However, the experimental ratio for |U 3 |/|U 1 | is larger than the model ratio.This is explained by the presence of small amounts of higher harmonic content in the measured free stream velocity (Table 1), e.g.|u 3 |/|u 1 | 5% for case D1, which is the "outlier" in Fig. 12.There is an overall good match between model and experiments for |U 5 |/|U 1 |.It is furthermore seen that |U 3 |/|U 1 | becomes larger when friction is introduced (C f = 0.02) and D is not the sole descriptive parameter for the momentum exchange between  13, where the secondary bump can be seen in the bottom panel.The free stream velocity signals in the middle and bottom panels are velocity-skewed, though the in-canopy velocity signal becomes increasingly asymmetric (sawtooth-shaped) for increasing D and |u 2 |/|u 1 |.The influence of the time variation in U on sediment transport is discussed in Section 5.3.
The importance of the sub-harmonic contributions is depicted in Fig. 14, where the sum of self-induced and superharmonic contributions are shown as solid markers and the sub-harmonic contributions are shown as open markers.For |u 2 |/|u 1 | = 0 (top panel), it is seen that the sub-harmonic contributions are small for m = 1 for all values of D , while for m = 3 the sub-harmonic contributions are 10% or more relative to the self-induced and super-harmonic contributions for  is practically unaltered for m = 3 and that the sub-harmonic contributions are small for m = 2.For m = 1 the sub-harmonic contribution exceeds 10% for D > 1.The importance of subharmonic momentum on m = 1 is relevant to consider, if a solution to nonlinear, free surface waves in canopies is sought: Any nonlinear wave will have significant velocity-skewness, so the relative importance of drag contributions will have more resemblance to the bottom panel than to the top panel and sub-harmonic contributions may not be negligible.
The existence of higher in-canopy harmonics (even for u m = 0, m = 2, 3, . ..) means that there must be a shear transfer of higher harmonic velocities into the free stream, since there is otherwise no momentum balance at the canopy interface.This higher harmonic shear component could be a source for eddy viscosity generation and shear layer instabilities in the wave boundary layer on top of the canopy.Lowe et al. (2005) presented a solution for, and qualitative discussion of, the in-canopy velocity reduction.However, they considered two independent nondimensional quantities (d/S and a 1 /S) instead of their product: D ∝ (d/S) • (a 1 /S).All the results from the nonlinear computations above are presented as the ratio max Û/ max û0 in Fig. 15 (top panel), where it is important to note that the use of Û and û0 includes all higherorder velocity components.Overall, there is limited scatter in max Û/ max û0 at given D , except for small D (inertiadominated flow).The scatter for small D follows directly from the momentum equation, which reads as follows for vanishing The functional fit in the bottom panel is given by Eq. ( 23) drag:

In-canopy velocity reduction
This gives a minimum inertia-related velocity reduction by −1 I < 1, which is incorporated in the bottom panel, where the ratio max Û I / max û0 is shown to be a function of D only.Note that the figure contains data for both C f = 0.00 and C f = 0.02, so the friction on top of the canopy has no noticeable effect on the velocity reduction, which is contrary to the findings for the phase lead and momentum distribution on individual harmonics.The bottom panel shows that when normalized by the minimum inertia-related velocity reduction, D is the sole nondimensional parameter that determines the velocity reduction.Figure 15 also contains the experimental data and there is a good match from inertia-to drag-dominated regimes.
Based on these results, an analytical expression is presented that can be used for an engineering estimate of the in-canopy velocities for regular, oscillatory flows: max Û max û0 This function is valid for horizontal, oscillatory free stream flows with a first and second harmonic content, where the two harmonics are in phase, and |u 2 |/|u 1 | ≤ 0.25.It is furthermore valid for 0.9 ≤ n and 0.00 ≤ C f ≤ 0.02.It remains to be shown whether this expression can be used for the prediction of in-canopy velocities for real waves with a considerable vertical velocity component.The inertia-related velocity reduction ( −1 I ) on the left-hand side is reasonably estimated by setting C m = 1 for large Keulegan-Carpenter numbers (Sumer & Fredsøe, 1999).It is also found that: following the same functional fit.This is a relaxation (use of " ") of the exact identity for purely sinusoidal motion for which max Û = √ 2 Ûrms and max û0 = √ 2û 0,rms .

Development of higher-order free surface wave theory
Analytical models of dissipative, free surface waves in canopies are found in the literature and they are all restricted to first order with a linearized drag formulation (Asano et al., 1992;Dubi & Tørum, 1994;Jacobsen, 2016;Kobayashi et al., 1993;Méndez et al., 1999), i.e. with a sinusoidal time variation of all variables.The present analysis has shown that accurate results should be expected from linear wave theories with respect to the in-canopy velocity, so the estimated radiation stress tensors (Jacobsen, 2016;Mendez et al., 1998), Lagrangian Stokes drift (Jacobsen, 2016), and vertical shear stress distribution (Jacobsen & McFall, 2022) should be reasonably approximated and can be applied in large-scale, practical engineering models, which rely on the second-order terms from Airy wave theory.However, for waves in non-vegetated, shallow water with high nonlinearity, the use of linear wave theory breaks down for the horizontal momentum flux, which was illustrated quantitatively by Dean and Bender (2006), so development of a higher-order wave theory for in-canopy wave propagation is Figure 16 Velocity-asymmetric free stream velocity.Left column: maximum free stream velocity of 0.075 m s −1 and varying r = {0, 0.375, 0.75}.Middle column: same as left with maximum free stream velocity of 0.70 m s −1 .Right column: fit with Eq. ( 23) for ratio of maximum velocities and ratio of standard deviations, respectively needed for special applications.The present analysis suggests that there is limited gain in developing a second-order free surface wave theory, because the theory would still rely on a linearized drag (no third-order contributions proportional to e.g.B 1 U 2 , B 2 U 1 ).Secondly, the present analysis indicates that a higher-order wave theory ought not to be based on a perturbation approach, but rather on a direct solution of the closure coefficients, as done with stream function wave theory (Fenton, 1988;Rienecker & Fenton, 1981), because the sub-harmonic contributions may not be negligible to the first harmonic and thus invalidate a perturbation approach for large values of D .

Velocity-asymmetric free stream velocity
The bulk of this work investigated velocity-skewed free stream velocities, so here the transformation of a velocity-asymmetric free stream velocity into a canopy is briefly presented.The free stream velocity by Abreu et al. (2010) is applied: where 0 ≤ r ≤ 0.75 gives the degree of velocity-asymmetry for ϕ = 0, and u w = (max u 0 − min u 0 )/2.u 0 is decomposed as per Eq. ( 8) and the harmonic amplitudes are given as input to the nonlinear frequency-domain model.Examples with u w = {0.075,0.70} m s −1 , N = 1736 stems m −2 , d = 0.0083 m, C m = 1, C D = 2.3, T = 6.0 s, and r = {0, 0.375, 0.75} are shown in Fig. 16.It is seen that there is a significant difference between the positive and negative part of U and qualitatively it is concluded that U becomes increasingly velocity-skewed for increasing r (increasing free stream velocity asymmetry).This holds both for D 1 and D 10. Here, " " is applied, since |u 1 | varies with r for a fixed u w .
U was computed for a range of u w ∈ [0.01, 2.5] m s −1 and the two forms of the velocity reduction Ûrms I /û 0,rms and max Û I / max û0 are depicted against Eq. ( 23) in Fig. 16 (right panel).It can be seen that only the form based on the rms-velocities collapses with the functional fit.

Sediment transport
Sediment transport within canopies under oscillating flow and waves has been reported in a few studies (Tinoco & Coco, 2014, 2018;Yang et al., 2016), who report that the sediment transport is linked to the in-canopy turbulent kinetic energy.Yang et al. (2016) found that the critical velocity within the canopy for incipient sediment transport decreased with decreasing porosity (increasing solidity), while the critical turbulence level was estimated to be constant.Furthermore, Tinoco and Coco (2018) measured suspended sediment concentrations and showed positive correlations between sediment concentrations, turbulence levels and solidity.Nakayama and Kuwahara (1999), Yang et al. (2016) and Tinoco and Coco (2018) proposed formulations for turbulent kinetic energy in repeated arrays that are proportional to the square of the in-canopy velocity, with the proportionality factor depending on the porosity (solidity).A combination of the proposed empirical formula for the in-canopy velocity (Eq.( 23)) and the work by Nakayama and Kuwahara (1999) can be used to predict in-canopy turbulent kinetic energy, though validation against experimental data is required for a validated, predictive formulation.A validated model would allow for a prediction of the in-canopy suspended sediment concentrations based on D and I as sole input parameters along with the characteristics of the sediment (density and median diameter).
The in-canopy redistribution of momentum to higher harmonics may be important in determining period-averaged (net) sediment flux, which is required for practical engineering applications.It is documented for non-vegetation wave boundary layers that the direction of net sediment transport under nonlinear free stream velocities depends on the sediment fall velocity and the velocity asymmetry and skewness (Fuhrman et al., 2013; of the University of Aberdeen technical staff for the experimental work, especially that of Fluids Laboratory Technician Roy Gillanders.The experimental dataset is available on https://dx.doi.org/10.5281/zenodo.4560141.

Disclosure statement
No potential conflict of interest was reported by the author(s).1)).Consequently, focus is given to the derivation of the horizontal momentum equation Eq. ( 2).The horizontal momentum equation reads as follows (Jensen et al., 2014, with  where horizontal diffusion terms are already ignored, since all derivatives along x and y vanish for spatially-averaged quantities.In the same way, the first two convective terms can be removed.The third convective term vanishes, since w = 0 given the bottom boundary condition and the continuity equation (Eq.( 1)).Finally, the diffusion term proportional to ∂w/∂x vanishes (again, since w = 0 for all z).This results in the horizontal momentum equations in Eq. (2).

A.2 Bulk two-layer model
The bulk two-layer model is found by a vertical averaging of Eq. ( 2) over the canopy height: and introduction of the bulk velocity: The first two terms in Eq. (A3) become I ∂U/∂t and 1 ρ ∂p/∂x, respectively, as per Eq. ( 5): the first term, because ∂/∂t can be taken outside the integral, and the second term, because p is independent of z.The only horizontal shear terms acting on the control volume bounding the vegetation are on the seabed and the top of the canopy, which are lumped into a single shear term τ = C f | U| U (see Section 2.1 for a discussion on this parameterization).However, since Eq. ( A3) is averaged over the canopy height it follows that: Finally, the drag term reads: Here, α U is the shape factor accounting for the slight deviations from the near-constant in-canopy velocity profile (Fig. 10); α U is in practice lumped into C D , when C D is derived from the experimental results.

Figure 1
Figure 1 Definition sketch for oscillating (tunnel) spatially averaged flow with submerged, rigid vegetation.The analytical domain is vertically unbounded.u 0 is the bulk free stream velocity, U(t) is the bulk in-canopy velocity and the blue u(t, z) is the velocity profile including viscous effects

Figure 4
Figure 4 Periodic representation of the (a) sparse and (b) dense array geometries viewed from above.The dashed line at the bottom of each figure indicates symmetry, and the dotted lines around the other three sides indicate periodic boundaries.The figures show the coordinates where vertical profiles were measured (•), as well as corresponding "mirrored" coordinates based on symmetry of geometry and flow characteristics (X) as described in Section 3.2

Figure 5
Figure5Measured and calculated moment around the y-axis M y , case S4.The free-stream velocity u 0 is also shown for reference

Figure 6
Figure 6 Drag coefficients obtained from moment measurements as a function of Reynolds number

Figure 7
Figure 7 C D,c obtained from moment measurements as a function of constricted Reynolds number, with predicted values from Eq. (20) shown for comparison Figure7shows a comparison between the present data and Eq.(20).While there is some scatter, the results from the present isolated-cylinder and canopy experiments fall within approximately 10% of predicted values from Eq. (20), demonstrating that the equation can provide a good estimation of C D in the absence of direct measurements.Nevertheless, in the present work, model calculations are based on measured C D -values as shown in Fig. 6, taking the average measured values of C D = 1.30 for the sparse canopy and C D = 2.30 for the dense canopy.It should be noted that most wave attenuation studies define an empirical "bulk drag coefficient" C D,b , obtained from the measured wave attenuation using the method first described byDalrymple et al. (1984).Henry et al. (2015) show examples of different C D,b formulations.In these studies, C D,b tends to decrease significantly with increasing R e , because the calculation does not account for velocity reduction within the canopy (relative to the free-stream velocity), which increases with increasing R e .In contrast, C D in the present work is based directly on the incident in-canopy velocity, and, as Fig.6shows, C D does not appear to vary significantly with R e .

Figure 9
Figure 9 Comparison of the harmonic behaviour of the experimental results, sinusoidal, nonlinear, and LKM (2005) models.The label "Mode" refers to the n = 1, 3, 5 harmonics in the nonlinear models and experiments, and the n = 1 harmonic in the sinusoidal model.Top rows: The odd harmonic velocity magnitudes |U 1 |, |U 3 |, and |U 5 |.Bottom rows: The odd phases relative to the first harmonic free stream phase: φ n = arg U n − arg u 1 for n = 1, 3, 5

Figure 11
Figure 11 Comparison of the arg U 1 relative to arg u 1 = 0 for the sinusoidal and nonlinear models.Top: C f = 0.00.Bottom: C f = 0.02.The range of NC D is 13-2546

Table 1
Overview of the ten experimental cases Note: Cases S1-S5 have a sparse canopy and cases D1-D5 have a dense canopy.Case names are applied in Figs8 and 10. u 1 , u 2 , and u 3 are free stream velocity amplitudes.N w and N a are the number of averaging points within and above the canopy, respectively.I , D , and f are non-dimensional resistance quantities.The normalized deviations to the experiments for the sinusoidal (δ S ), nonlinear (δ NL ), and Lowe et al. (2005) (δ LKM ) models are defined in Section 4.1.
the present work's formulation of drag and