Interaction of Waves with Idealized High-Relief Bottom Roughness Journal of Geophysical Research: Oceans

Considerable uncertainty exists about how to represent wave bottom friction in coastal systems like reefs where orbital excursions are similar to roughness element size. Here, interactions between waves and large bottom roughness were investigated using Large Eddy Simulations of oscillatory ﬂow over inﬁnite hemisphere arrays. Wave amplitude, period, and hemisphere spacing were varied to investigate the dependence of kinematics and dynamics on dimensionless parameters. The net effect of topography on the oscillatory ﬂow was assessed using a spatially and phase-averaged Navier-Stokes framework. Dynamics depended strongly on Keulegan-Carpenter number ( KC ), the ratio of wave orbital excursion to roughness element size. For 1 < KC < 10, ﬂow separation was weak, form drag was small, stress gradients were negligible, and the main effect of topography on the ﬂow was the inertial force associated with acceleration around roughness elements. For 10 < KC < 20, strong ﬂow separation occurred, and both drag and inertial forces were important. Phase-dependent dispersive stresses were the main mechanism for vertical momentum transfer between the canopy layer and the overlying water column. Friction factors based on the drag force increased with KC for 1 < KC < 20, different from previously proposed empirical curves, but approached these curves for high KC . Friction factors based on the total force decreased with increasing KC and were consistent with previously proposed curves. These results highlight the importance of distinguishing the total force on the bottom, the drag force that removes energy from the ﬂow, and the shear stress above the canopy layer, which were very different for the parameter range in this study. computed the rms of the time-varying momentum balance terms in equation averaged in the depth-average z / D 5 0 to z / D 5 0.5 function of f = k s to back out k s because they will necessarily fall on the curves that were used (e.g., Lentz et al., 2016; Mathisen & Madsen, 1996; Rogers et al., 2016). In sediment transport studies, sand-grain roughness k s is commonly parameterized by the median grain diameter D 50 (Fredsøe & Deigaard, 1992; Nielsen, 1992). In this study, we use the diameter of the hemisphere as the roughness scale k s for simplicity, which gives the ratio of wave orbital excursion to roughness scale f = k s 5 KC = 2 p . Our LES results show that for f = k s < 1, the friction factor representing form


Introduction
Complex topography is a common feature in many natural systems. Flow over complex topography is of great interest in applications such as flow over sand dunes (Charru et al., 2013;Parson et al., 2004), river flows over complex bed topography (Anderson, 2013;Lane et al., 2002;Nicholas, 2001), and wind over terrain (Kim et al., 2000) and urban areas (Britter & Hanna, 2003;Coceal & Belcher, 2005). Predicting flow over complex topography has been the goal of many research efforts, as it is critical to understand processes such as geomorphic evolution (Livingstone et al., 2007), mixing and transport of mass and momentum in coastal regions with complex topography ( € Ozg€ okmen et al., 2004), pollutant dispersion in the atmosphere (Sl adek et al., 2007), and heat transport in urban areas (Kr€ uger et al., 2011).
Unlike systems described above where the flow is effectively steady, flow in the coastal ocean is often dominated by unsteady water motions associated with surface waves. As waves propagate, they lose energy via interactions with the bottom topography. Understanding interactions between water waves and bottom topography is important for predicting wave transformations and quantifying turbulent mixing and transport ( € Ozg€ okmen et al., 2004), which in turn affect the larger-scale circulation patterns. Understanding wavetopography interactions is also important for estimating transport of dissolved and particulate materials, such as nutrients and sediments (Boehlert, 1988;Carroll & Dixon, 2002).
In coastal wave and circulation models, the horizontal spatial resolution, which is typically O(km), and O(100 m) for high-resolution simulations, is much larger than scales of the topography (e.g., Lambrechts et al., 2008;Lowe et al., 2009;Tarya et al., 2010). Interactions of the fine-scale topography with flows cannot be resolved in these models, and these subgrid processes need to be parameterized. In ocean models, bottom friction is typically represented using a quadratic law with a wave friction factor (for waves) or a drag coefficient (for current). In systems like terrestrial forests and urban canopies, effects of spatial heterogeneity are dealt with by averaging the governing equation in space. Spatial averaging leads to drag, averaged vertical mixing, and dispersion terms associated with spatial heterogeneity, and these terms are then parameterized. While a spatial averaging approach has become standard for representing steady flows in terrestrial systems (e.g., Finnigan, 2000), a formal spatial averaging approach has not yet been applied to unsteady flows associated with waves, which are important in the ocean.
The interaction of surface waves with the seafloor is typically parameterized using a quadratic drag relation with a wave friction factor computed from a roughness length-scale that is often tuned to match simulation results with measurements. While the wave friction factor was originally proposed to represent friction over smooth bottoms and sand grain roughness (e.g., Dean & Dalrymple, 1991), it has become the standard way to parameterize effects of unresolved bottom topography on waves over large areas of the seafloor (e.g., a model grid cell, Cavaleri et al., 2007). It has also been used to represent frictional effects on waves in systems with high relief, complex topography, such as coral reefs (e.g., Lentz et al., 2016;Lowe et al., 2009;Monismith et al., 2015). A large number of laboratory and field studies have shown that for turbulent environmental flows over sand-grain roughness, the wave friction is controlled primarily by the Keulegan-Carpenter number (KC5U b T=D, where U b is the amplitude of the near-bed wave orbital velocity, T is the wave period, and D is the grain diameter) representing the ratio of the wave orbital excursion to the roughness element size (Dean & Dalrymple, 1991). However, the majority of this work has focused on situations in which the wave orbital excursion is large compared with the roughness element size (KC > 10). In this parameter range, a turbulent wave boundary layer forms above the roughness elements, the near-bed eddy size is directly related to the geometry of the roughness elements, and friction factors collapse onto a single curve and decrease monotonically with increasing KC. A number of theoretical and empirical formulas have been proposed for the dependence of friction factor on KC (Dean & Dalrymple, 1991;Grant & Madsen, 1982;Nielsen, 1992); however, these curves diverge for low KC (Lentz et al., 2016;Rogers et al., 2016). The applicability of these models in low KC situations is also questionable, because the basic assumption that roughness elements are small compared with the wave boundary layer thickness is violated. There is considerable uncertainty about the behavior of the friction factor, and processes that control it, in this part of the parameter space; however, this KC range is often of interest in the coastal ocean, particularly in systems like coral reefs. Some more recent work has investigated physical processes that control wave friction at lower KC relevant to coastal systems like reefs (KC $ 1210) using laboratory experiments and field measurements with idealized geometries (Lowe et al., 2005(Lowe et al., , 2007. Laboratory experiments illustrate that interactions between obstacles and unsteady flows can alter both the amplitude and phase of oscillatory velocities within the canopy layer (Lowe et al., 2005;Luhar et al., 2010;Rosman et al., 2013). An important transition occurs around KC $ 1. For KC ) 1, the primary force exerted on water is the drag force, whereas for lower KC, the inertial force is largest. Although the inertial force is in quadrature with the flow speed, and so does not remove energy from waves, it affects the magnitude of the flow velocity in the layer containing solid obstacles and can thereby affect drag on the flow. In these previous studies, inertial and drag coefficients have been assumed to be constant; however, measurements of forces on cylinders in oscillatory flows suggest that both the inertia and drag coefficients can vary with flow parameters when KC is of order unity (Dean & Dalrymple, 1991;Keulegan & Carpenter, 1958). Additionally, in the parameter range where wave orbital excursion is similar to obstacle size (KC $ 1), the flow field is likely to be extremely heterogeneous. For example, in field measurements on a coral reef, Hench and Rosman (2013) observed that spatial flow patterns around a single coral colony depend strongly on the frequency in the velocity spectrum; there was a strong wake signature in longer period flow variations (KC > 1) and no wake signature in shorter period variations (KC $ 1). The implications of three-dimensional spatial flow structure for the dynamics of the larger-scale oscillatory flow have not been formally investigated.
Flow over complex topography often results in vortex shedding, flow recirculation, and wakes, and is composed of a broad spectrum of turbulence length-scales. Large eddy simulations (LES) can resolve large-scale turbulent motions and have been used to study flows over aquatic vegetation (Chakrabarti et al., 2016;Stoesser et al., 2009), sand ripples (Zedler & Street, 2001), and urban geometries (Xie & Castro, 2009). LES can also accurately predict turbulent wakes behind solid obstacles, as well as interactions between wakes and obstacles. In contrast RANS (Reynolds Averaged Navier-Stokes) models fail to capture the organized large-scale unsteadiness and coherent structures (Chang & Scotti, 2004). Complex interactions of flow with the topography play dominant roles in redistributing mass and momentum in the flow and are believed to be the driving mechanism for scalar transport over regions with complex topography. High-resolution LES models can resolve these processes and provide a detailed description of the instantaneous unsteady turbulent flow, pressure fields, and provide insight into the underlying dynamics.
Herein, we investigate the interaction of waves with high-relief roughness. The focus is the part of the parameter space where wave orbital excursions are similar to the size of individual bottom features (KC $ 1-10), which is important in coastal systems like coral reefs but currently poorly understood. We conducted LES simulations of oscillatory flow over idealized topography composed of an infinite array of hemispheres. Hemispheres were chosen for the roughness elements due to their three-dimensionality and their similarity to both sand-grain roughness used in many previous studies and shapes of individual coral colonies on reefs. Geometric and flow parameters were varied to investigate how the kinematics and dynamics of the flow field vary with important dimensionless parameters. The implications of processes observed at the individual roughness element scale for the large-scale oscillatory flow are analyzed using a volume and phase-averaged momentum budget. We explore the dependence of drag and inertial forces, as well as turbulent and dispersive stresses, on important dimensionless parameters. This work provides insight into fundamental processes that need to be considered when parameterizing interactions between waves and seabeds with large, high-relief roughness.

Background: Double-Averaging Framework for Flow Over Complex Topography
In order to investigate the effects of wave-topography interactions on large-scale hydrodynamics and wave dissipation, we employ a spatial averaging approach. Spatial averaging of the flow properties is applied over volumes with horizontal length scales greater than the size of a single solid obstacle but fine enough in the vertical directions that the vertical structure of the dynamics can be resolved. Previous work in systems such as forest canopies (Finnigan, 2000), streambeds (Nikora et al., 2007), and aquatic vegetation canopies (Nepf, 2011), has applied a spatial-averaged framework to steady flows. Here, we extend equations developed in those studies to oscillatory flows.
The instantaneous streamwise velocity is decomposed as u5 u c 1 u w 1u 0 , with u c the current, u w the wave component, and u 0 the turbulent fluctuations. In this study, there is no current, so the decomposition can be simplified as u5 u w 1u 0 . For monochromatic waves, the velocity component associated with waves can be computed from the phase-average over N waves, defined as u w ðtÞ5 where T is the wave period, x is the wave frequency in rad/s, and xt is the wave phase.
An intrinsic spatial average (Raupach & Shaw, 1982) is used, where the average is taken over the volume occupied by fluid between altitudes z 1 and z 2 , with z 1 5z2Dz=2 and z 2 5z1Dz=2, and Dz is the height of the averaging volume.
L x and L y are the dimensions of the averaging volume in the horizontal directions, and the spatial average is indicated by angle brackets. V f is the volume of fluid within the averaging volume. The velocity field can therefore be decomposed as Double prime u 00 w represents the deviation of the phase-averaged velocity from its spatial average (u 00 w 5u w 2hu w i). The spanwise velocity (v), vertical velocity (w), and pressure (p) are decomposed and spatially averaged in a similar way.
Double-averaged forms can be derived for conservation equations of mass, momentum, and turbulence properties such as turbulent kinetic energy. The spatially averaged momentum equation is with new terms arising from spatial averaging u i is the fluid velocity in the i direction (with i 5 1,2,3 corresponding to the streamwise (x), spanwise (y), and vertical (z) direction, respectively), / is the solid fraction, defined as the volume of solid within the averaging volume divided by the total averaging volume. q is the density of the fluid and l is the dynamic viscosity of the fluid. n is the unit surface normal vector. s ij is the sum of Reynolds stress, dispersive stress, and viscous stress. The Reynolds stress represents momentum transport due to turbulent fluctuations. The dispersive stress is the spatial analog of the Reynolds stress, and it represents momentum transport due to spatially averaged effects of persistent spatial variations in the flow. The pressure force (f Pi ) and viscous drag terms (f Vi ) come from integration of the pressure gradient and viscous stress terms, respectively, over the surface of the solid obstacles. In the coastal ocean, the pressure force term is typically significantly larger than the viscous drag term, therefore, the viscous drag term can be neglected.
In previous studies that considered only steady flow, the pressure force f Pi was taken to equal the form drag and parameterized using a quadratic drag law (e.g., Finnigan, 2000). In oscillatory flow, f Pi has two components, the form drag and an inertial force that arises because the fluid is accelerating as it moves over the solid surface. Following the Morison equation (Morison et al., 1950), the total in-line force exerted on a solid body in oscillatory flow can be expressed: where i is a unit vector in the flow direction. The first term on the right hand side is the inertial force component and the second term is the drag force component. The inertial force is in phase with the local flow acceleration du/dt and the drag force is in phase with the flow velocity. V is the submerged volume of the body and A is the frontal area perpendicular to the flow direction. C M is the inertia coefficient and C D is the drag coefficient. The force applied on the fluid by the solid surface is equal and opposite to the force applied on the solid surface by the fluid. Therefore, by dividing equation (8) by the fluid mass in the averaging volume, the streamwise pressure force term f P;x can be therefore obtained The inertial force can be further decomposed as where the first term is due to the unsteady pressure gradient that drives the flow, and the second term is the added mass force that occurs due to resistance to acceleration of flow around the surface of the hemisphere. The added mass coefficient can be computed C a 5C M 21 following equation (10).

Model Setup
Numerical modeling experiments were designed to understand oscillatory flow over a bottom covered with idealized high-relief roughness elements. The high-relief bottom was specified as an infinite array of regularly spaced hemispheres of the same height, aligned through their centers, on an otherwise flat horizontal bottom ( Figure 1). Simulation domain horizontal dimensions were set equal to the distance between centers of adjacent hemispheres S. The height H of the domain was 2 m for all simulations. Periodic boundary conditions were applied in both horizontal directions, so interactions between flow structures generated from adjacent hemispheres could be captured. The top boundary was sufficiently far from the hemisphere that a free-stream region existed. A smooth wall boundary condition (noslip and no-penetration) was applied at the bottom boundary, including the surface of the hemisphere, and a free-slip boundary condition was applied at the top boundary.
For oscillatory flow over high-relief roughness elements, the length scale of the largest turbulent eddies is set by the length scale of individual roughness elements (Barr et al., 2004). The minimum domain size needed to accurately simulate oscillatory flow over the hemisphere array therefore scales with hemisphere diameter. Smaller domains were used because using large domains would entail a considerable increase in computational expense. The influence of the domain size was tested by comparing the flow statistics between simulation results from a large domain and small domain. We carried out a simulation for KC 5 12 with the domain size doubled in each horizontal direction (four hemispheres in the domain). Velocity fields agreed with results with the smaller domain within 8%, and drag and inertial forces agreed within 10%.
A block structured grid was generated using utilities in the OpenFOAM package (Weller et al., 1998). The blockMesh utility was used to generate the grid, which was then adjusted via the blockMeshBodyFit tool to fit the mesh to the surface of the hemisphere. A finer grid was applied near the surface of the hemisphere to better resolve turbulence ( Figure 2 and Table 1). N total is the total number of grid points and D r min is the size of the first grid cell (smallest) right above the surface of the hemisphere.
A fixed time step was picked for each case based on the CFL condition (CFL < 1), which guarantees the distance that any information travels during one time step within the mesh is smaller than the mesh grid size. All simulations were initially run for 10-15 wave cycles to ensure flow reached quasi-steady state. After the initial spin-up, each simulation was continued for another 10 wave cycles and these data were used for analysis.   Table 1.

Large Eddy Simulation Model
Simulations were conducted using OpenFOAM (version 3.0). Oscillatory motion was driven by an oscillating horizontal pressure gradient that was introduced in the momentum equation. The pressure gradient term is similar to that used to generate the oscillatory motion in a U-tube. Sufficiently high above the hemispheres and the wave boundary layer, the momentum equation reduces to A sinusoidal free-stream velocity was used with u 1 5U 1 sin ðxtÞ, where U 1 is the amplitude of the oscillating free-stream velocity. The pressure gradient driving the flow can be derived as 2dp 1 =dx5qU 1 xcos ðxtÞ. The governing equations can be written u i is the resolved fluid velocity, p d is the dynamic pressure, and s ij is the stress term including both turbulent stress and viscous stress. d ij is the Kronecker delta. Note p d is the nonhydrostatic, dynamic pressure part, the total pressure is p5p d 1ðx2x 0 Þ dp1 dx . LES turbulence closure was used to calculate the turbulent stress. We initially carried out a series of runs using a two equation k2x turbulence closure model (Wilcox, 1991). The RANS model failed to predict the flow instabilities around flow reversal, which are important to the generation of turbulence when KC is small. Furthermore, the RANS model overpredicted the drag force because the turbulent wake was poorly developed. These results are consistent with previous work that has shown RANS models may fail to generate satisfactory predictions of flow around bluff bodies under certain conditions when flow separation and vortex formation are not well represented (L€ ubcke et al., 2001;Rodi, 1997). We therefore used LES for all cases presented herein. In LES models, the turbulent stress term represents the momentum fluxes by smallscale processes occurring at length scales that cannot be adequately resolved on the computational mesh. The wall-adapting local eddy-viscosity (WALE) model was used as the subgrid closure. The stress term is where S ij is the strain rate tensor of the resolved scale, m is the kinematic viscosity of the fluid, and the eddy viscosity m t is given by Ducros et al. (1998) The constant C w 50:325 and V c is the volume of the grid cell.

Simulations Conducted
A series of simulations (22 total) were carried out with different parameter combinations to explore how the flow kinematics and dynamics varied with important dimensionless parameters (Table 2). With the geometric parameters D (hemisphere diameter), S (hemisphere spacing) and H (water depth), and flow parameters Table 1 Summary of Computational Grid Information Dr min (m) 1 2 0.5 7.36 3 10 6 6.9 3 10 23 2 1 0.5 5.68 3 10 6 2.6 3 10 23 3 0.75 0.5 4.40 3 10 6 1.7 3 10 23 T (wave period), U 1 (freestream velocity amplitude) and m (kinematic viscosity), the important dimensionless groups derived based on the Buckingham p theorem are H/D plays an important role when a steady current exists in the system because Froude number (Fr5U= ffiffiffiffiffiffi gH p , and g is the gravitational acceleration) is one of the most important parameter in steady shallow water flows (Grubi sić et al., 1995). For cases with waves only, the wave boundary layer and roughness element size control the length scales of turbulent eddies and the water-depth plays a negligible role in the dynamics as long as H is large enough such that a freestream region exists. The effect of water depth is therefore not discussed in this study. For all simulations, the diameter of the hemisphere D was 0.5 m. Keulegan-Carpenter number (KC), defined as the ratio of wave excursion to the size of the object, was varied by changing either the wave period T or the amplitude of the freestream velocity U 1 . The Stokes number (b) in this study is only a function of wave period T because D was the same for all cases. b can be interpreted as the ratio of the characteristic length scale of bottom topography to the laminar wave boundary layer thickness (Stokes boundary layer thickness d s 5 ffiffiffiffiffiffiffiffiffiffi mT=p p ). A Reynolds number can be formed from the product of KC and b, but it is not an independent dimensionless parameter. Three different spacings S were used to vary S/D and explore the effects of interactions between adjacent hemispheres.

Flow Kinematics
Two contrasting cases are described here in detail to illustrate distinct flow features that occur at different KC, one with KC 5 2 (U 1 50:05 m/s, T 5 20 s, case 6) and the other with KC 5 12 (U 1 50:3 m/s, T 5 20 s, case 9) (Figure 3). In case 6, because the orbital excursion (f5U 1 =x) is small relative to the hemisphere diameter, there is no flow separation and flow follows the surface of the hemisphere closely. The dynamic pressure distribution is smooth without fluctuations, and is nearly antisymmetric at the flow reversal and symmetric at the flow peak. When no flow separation occurs, streamlines attach to the surface of the hemisphere, and flow is expected to show symmetry under flow peak and flow reversal. In contrast, in case 9, strong flow separation occurs. Small-scale fluctuations in the pressure field coincide with turbulent velocity  (Figures 3i and 3j). A region of recirculation behind the hemisphere occurs when the flow separation is strong (Figure 3p). By comparing results from the two cases, it is clear that characteristics of the flow field vary strongly with KC and the degree of flow separation plays an important role in the dynamics of oscillatory flow over a three-dimensional obstacle array.
The development of the flow field in one wave cycle can be seen more clearly in the vorticity fields ( Figure 4). Flow separation in the low KC case (case 6, KC 5 2) is very weak, and vortices follow the hemisphere for the majority of the wave cycle (Figures 4a-4h). During the decelerating phase (Figures 4e-4h), the growth of vorticity around the hemisphere can be clearly identified, with negative vorticity in the upper half and positive vorticity in the lower half. Because KC in this case is small with f ( D, vortices are not able to grow strong enough to evolve into a fully developed wake, and only some small eddies are shed from the hemisphere. These small eddies at the downstream side of the hemisphere arise mainly from flow instabilities caused by the centrifugal effect when an accelerating flow moves over a curved surface, similar to G€ ortler instabilities (Saric, 1994) observed in oscillatory flow past a rippled bed (Tseng & Ferziger, 2004). These small eddies are generated primarily around the flow reversal when an unstable velocity profile develops ( Figure 5). As the flow peak is approached, flow reattaches to the surface of the hemisphere. During the accelerating phase (Figures 4a-4d), small eddies that were generated during flow reversal are advected back to the hemisphere and interact with the hemisphere.
Flow separation in the high KC case (case 9, KC 5 12) is much more intense than in case 6. At flow reversal (Figure 4i), strong interactions between eddies and the upstream (left) side of the hemisphere occur. These small-scale eddies are the remains of the turbulent wake from the previous cycle. The eddies are advected back and impinge on the hemisphere, affecting the pressure on its surface (Figure 4j). The flow reattaches to the surface of the hemisphere at the upstream side when these eddies have advected past the hemisphere (Figure 4 k). Because the wave excursion is greater than the size of the domain, the eddies can be transported further downstream and interact with the adjacent hemisphere. During the decelerating phase (Figures 4m-4p), the flow does not change direction so the intensity of the flow separation continues to grow, as indicated by the size of the wake zone and number of small eddies in the wake. Just before flow reversal when the excursion is close to its maximum (Figures 4o and 4p), strong interactions of the wake from the adjacent hemisphere with the local hemisphere can be observed as eddies encounter the hemisphere, and get stretched, squeezed, and rotated.
Spatially averaged velocity profiles have an inflection point near z/ D 5 0.5 at the top of the hemisphere (Figures 6a and 6c). Laboratory and field measurements have also shown that a shear layer exists at the top of submerged vegetation and coral canopies (Hench & Rosman, 2013;Monismith, 2007;Nepf, 2011). The spatially averaged velocity below z/D 5 0.5 is affected by the presence of the hemisphere, and the velocity is smaller than it would be over a flat bottom. At flow reversal (xt 5 0), the velocity below the top of the hemisphere changes direction before the flow in the freestream because of the drag exerted on the flow by the hemisphere. The phase lead is more pronounced in the high KC case (case 9) due to strong flow separation (Figure 6c).

Turbulent dissipation rate () was approximated by the dissipation rate of the unresolved scales in LES
In case 6, there is almost no flow separation and turbulence is weak; therefore, remains small, on the order of 10 28 m 2 =s 3 (Figure 6b). is largest at flow reversal because of flow instabilities. For all wave phases, the averaged turbulent dissipation rate h i shows a peak near z/D 5 0.2, and an almost linear decrease to zero at z/D 5 0.5. The dissipation rate above the top of the hemispheres is almost zero, indicating that for low KC there is negligible transport of turbulent energy into the overlying water column.
For case 9, the averaged turbulent dissipation rate is on the order of 10 25 m 2 =s 3 (Figure 6d), which is about three orders of magnitude greater than case 6 due to strong flow separation. h i is relatively constant below the top of the hemisphere and drops sharply to zero at a distance z/D of 0 to 0.1 above the top of the hemispheres. h i begins increasing after the flow reversal, and reaches a maximum after the flow peak, when the separation is also the largest (Figure 6d). The high dissipation zone is closest to the bottom at this wave phase, and extends highest in the water column after flow reversal (xt 5 p=4), suggesting transport of turbulence into the overlying water column.
To examine the time-variation of the energy dissipation in more detail, we examined the depth-integrated dissipation rate The normalized depth-integrated velocity U avg =U 1 5 Ð H 0 huidz=U 1 H (Figure 7a) does not vary significantly among cases. However, the normalized depth-integrated dissipation rate D=U 3 1 (Figure 7b) shows strong dependence on KC. Flow separation occurs more intensely in case 5 (T 5 10 s, U 1 5 0.53 m/s, KC 5 10.6) and case 9 (T 5 20 s, U 1 5 0.3 m/s, KC 5 12), which have similar KC but different b. In these cases, the peak dissipation rate occurs around xt50:6p, after the flow peak but while the flow separation region is still growing. For case 8 with an intermediate KC (T 5 20 s, U 1 5 0.15 m/s, KC 5 6), the phase lag between the freestream velocity and the depth integrated dissipation rate is more significant, with the peak around xt50:75p. For case 7 (T 5 20 s, U 1 5 0.1 m/s, KC 5 4), the normalized D peaks near the flow reversal because the flow separation is weak and generation of turbulence is dominated by flow instabilities that occur around flow reversal.  (Jeong & Hussain, 1995). G€ ortler vortices are indicated by the black ellipse.

Drag and Inertial Forces
The total in-line force acting on the hemisphere in the streamwise direction, which is the sum of drag and inertial forces, was computed from the simulated pressure and velocity fields as where i is the unit vector in the streamwise direction, n is the unit vector normal to the surface, and V hem 5 1 12 pD 3 is the volume of the hemisphere. The first term on the right hand side of equation (22) is the sum of the drag force and added mass component of the inertial force. The drag force is the force component that is in phase with the freestream flow velocity and 908 out of phase with the flow acceleration. The added mass force component represents the resistance to flow acceleration around the hemisphere and is in phase with the flow acceleration, but 908 out of phase with the freestream flow velocity. The last term on the right hand side of equation (22) is equivalent to the Froude-Kryolv force, caused by the unsteady pressure gradient that drives the flow. The inertial force is the sum of the added mass force and the Froude-Krylov force.
In planar oscillatory flows, the time-averaged drag force exerted on an object can be parameterized using a quadratic law that relates the root mean square (rms) drag force to a characteristic velocity scale such as U 1;rms and hemisphere frontal area A hem 5 1 8 pD 2 via a constant drag coefficient C D F D;rms 5 1 2 qC D A hem U 2 1;rms The root mean square of the free stream velocity U 1;rms is defined as for a sinusoidal wave. Similarly, the inertial force is parameterized using a constant C M that relates the rms inertial force to the rms acceleration. F I;rms 5qC M V hem xU 1;rms 5qC a V hem xU 1;rms 1qV hem xU 1;rms The inertial force is the sum of the added mass force and the Froude-Krylov force. The total in-line force acting on each hemisphere can therefore be parameterized as The force coefficients in equation (26) were determined using Fourier averaging (Keulegan & Carpenter, 1958). The drag coefficient was calculated from The added mass coefficient C a was computed as C M 21.
For case 6 with small KC (KC 5 2), the inertial force dominates, and the total pressure force following equation (22) is in phase with du 1 =dt (Figures 8a-8d). At the flow peak (Figure 8c), the pressure distribution is almost symmetric, which leads to a net pressure force around zero. During flow reversal (Figure 8a), the net pressure force is largest. In contrast, for case 9 (KC 5 12) with high KC and significant flow separation, the drag force is significant. The pressure field is not symmetric at the flow peak; therefore, there is a net pressure force at this phase. Around flow reversal, the pressure distribution is more symmetric than for case 6, so the pressure force is smaller during this part of the wave cycle.
The inertial force (F I ) varies with the local acceleration of fluid (du/dt) and is simply a cosine function of xt, with amplitude proportional to the inertia coefficient C M (Figure 9a). When the free stream velocity is small and flow separation is weak, the normalized drag force (F D ) is small (Figure 9b). The normalized drag force has a flat peak under the wave crest/trough, and it is about zero near the flow reversal. For case 8 (U 1 5 0.15 m/s, KC 5 6), the inertial force still dominates the dynamics. The normalized drag force starts to show two peaks in a half wave cycle, which are similar in magnitude, one around xt 5 0:3p and the other around xt 5 0:7p. For case 9 (U 1 5 0.3 m/s, KC 5 12), the flow separation becomes intense, and drag and inertial forces are of equal importance. The normalized drag force shows a large peak under the wave peak/trough, meaning the drag force is largest when the free stream velocity is largest. There is another peak around xt 5 0:16p, when eddies shed from the previous wave cycle are advected back and interact with the hemisphere ( Figure  4f). When these small-scale eddies impinge on the hemisphere, they create pressure fluctuations on the surface of the hemisphere. This results in an increase in the drag force. For case 10 (U 1 5 0.5 m/s, KC 5 20), the normalized drag force shows a similar pattern as case 9.
When KC increases, the drag coefficient increases because the intensity of the flow separation increases (Figure 10a). KC is varied by either changing the freestream velocity U 1 while keeping the wave period T constant (solid curves in Figure 10a), or changing the wave period T while keeping the freestream velocity U 1 the same (black dashed curve in Figure 10a). The slope of the C D versus KC curve decreases around KC 5 6, and C D seems to approach an asymptotic constant value of about 0.24. The drag coefficient in our study strongly depends on the intensity of flow separation. When the wave excursion is small (e.g., case 6), flow changes direction before the vorticity field grows strong enough to form strong flow separation, flow  follows the surface of the hemisphere and therefore the drag coefficient is small. When the wave excursion is large enough, flow separation is more developed and the drag coefficient increases with KC due to increased flow separation. The same trend occurs for different wave periods. With the same KC, the drag coefficient decreases with increasing Stokes number b ((D 2 =mT) or increasing wave period T). This is because with a longer wave period, the vorticity field around the hemisphere has more time to grow, which results in stronger flow separation.
Somewhat surprisingly, drag coefficient values are relatively insensitive to the spacing S/D between the hemispheres (Figure 10c). When KC is small, there is almost no variation in C D with S/D, and there is no significant interaction between flow structures formed on adjacent hemispheres. As KC increases (KC % 12), interactions of eddies shed from adjacent hemispheres with the local hemisphere become important, and these interactions are more intense in closely spaced arrays. However, these interactions only seem to have a minor effect on C D , even for the highest KC (KC 5 20) and closest spacing (S=D51:5) we simulated. Hemispheres were never within the separation bubble of neighboring hemispheres and little flow sheltering was observed. For S/D used in this study, the total projected frontal area per unit bottom area (k f ) is below or around 0.15, the value beyond which flow sheltering typically becomes important (Jim enez, 2004). Flow sheltering is therefore expected to play a minor role in the dynamics studied here.  The added mass coefficient C a decreased with increasing KC but did not vary with Stokes number ( Figure  10b). This pattern was the same for all three spacings simulated (Figure 10d), C a values were smaller for more closely spaced arrays, because S limits the volume of fluid affected by the presence of a single hemisphere.

Spatially Averaged Momentum Balance
An equation describing the effects of the obstacle array on the spatially averaged velocity field can be obtained by subtracting the far-field momentum balance equation (11) from the spatially averaged momentum balance equation (4).
@ðh u w i2u 1 Þ @t 5 @s xz @z 2f res ; where h u w i is the streamwise component of the double-averaged velocity, s xz is the shear stress, and f res is the resistance force per unit fluid mass exerted on the flow by topography. In the double-averaged framework, the stress term s xz is the sum of the spatially averaged turbulent stress s turb xz and the dispersive stress s disp xz . In LES the turbulent stress is the sum of a resolved component and a subgrid-scale (modeled) component. The spatially averaged turbulent stress is therefore calculated as where the first term is the resolved part of the Reynolds stress, and the second term is the subgrid stress from unresolved scales. The dispersive stress term is s disp xz 52h u 00 w w 00 w i. The resistance force f res can be obtained by dividing the in-line force exerted on each hemisphere by the mass of the fluid in the spatial averaging volume: The Froude-Krylov force associated with the oscillating pressure gradient that drives the flow does not appear in these equations because the oscillating pressure gradient is removed when the far-field momentum balance is subtracted. The resistance force is therefore the sum of the drag force and the added mass component of the inertial force For cases with small KC, the inertial force dominates the dynamics (Figures 11e-11h). The added mass force balances the acceleration term and the drag force is one to two orders of magnitude smaller. The stress term @s xz =@z is 2 orders of magnitude smaller than f a , and is therefore negligible.
For high KC cases, both inertial and drag forces play important roles in the dynamics (Figures 12 and 13). At flow reversal (xt50), the added mass force term f a is at its maximum. The drag force term f d is not zero and changes its sign around z/D 5 0.2, because the velocity in the canopy layer changes direction before the freestream velocity does. At xt5p=4, the drag force f d increases and the added mass force f a decreases, and they are of the same magnitude. At the flow peak (xt 5 p=2), the added mass force f a reduces to zero and the drag force f d is at its maximum. During flow deceleration, f d decreases and f a increases. In this case, stress terms, which represent vertical transport of momentum, although they are still much smaller than the drag and inertial force terms, are significant. The dispersive stress term is about 10 times greater than the Reynolds stress term. Flow around the hemisphere is highly three-dimensional. The vertical velocity is positive at the upstream side of hemisphere where the horizontal velocity is large and negative at the downstream side in the wake where the horizontal velocity is small. The dispersive stress term, which is the correlation of u 00 w and w 00 w , is therefore substantial. The dispersive stress plays an important role in vertical momentum transport especially between z/D 5 0.5 and 0.6, where pressure force terms (both the drag and inertial force terms) are zero.
To investigate how the relative importance of momentum budget terms varies with dimensionless parameters, we computed the rms of the time-varying momentum balance terms in equation (29) averaged in the canopy layer. The depth-average from z/D 5 0 to z/D 5 0.5 was calculated first, and then the rms was computed over all wave phases to obtain the magnitude of each term. Figure 14a shows the normalized momentum budget terms for all cases with S/D 5 2. The relative acceleration term is largest for both low and high KC, and a minimum around KC 5 10. At low KC, flow separation is weak (or no flow separation). In this KC range, the momentum boundary layer thickness at the upstream side decreases as KC increases, and the relative acceleration term therefore decreases as KC increases. At large KC, the increased strength of flow separation leads to the increase of the relative acceleration term with increasing KC. The drag force f d progressively becomes more important as KC increases and flow separation becomes stronger. Momentum budget terms were also compared across runs with several different wave periods, or equivalently different Stokes number (Figure 14a). The deviation of normalized drag force values among runs with the same KC, particularly evident for KC 5 4 and 6, occurs due to a dependence of drag force on Stokes number. This dependence is also seen in Figure 10a, where C D values for runs with the same KC are larger when Stokes numbers are smaller (longer period waves). The dispersive stress term initially increases with KC then drops as KC increases further. The dispersive stress term also shows almost no dependence on Stokes number b. The turbulent Reynolds stress is nearly zero for low KC cases and increases sharply for KC ! 12 when strong flow separation develops.
The main patterns are similar for different spacings S/D. The turbulent stress is more significant when S/D is smaller because wakes occupy a larger fraction of the fluid volume, resulting in a more significant spatially averaged vertical shear (Figure 14b). The drag force term is more significant for closely spaced arrays while the inertial force is more significant for sparser arrays, and this difference increases as KC increases. As S/D decreases, wake regions occupy a larger fraction of the total fluid volume in the canopy layer, and the phase difference between the spatially averaged canopy layer velocity and the freestream velocity increases. Since the inertial force is in phase with the freestream acceleration, the inertial force becomes relatively smaller and the drag force becomes relatively larger as the phase shift increases.

Momentum Balance
The pattern of how the normalized momentum budget terms vary with KC ( Figure 14) can be explored by scaling terms. We use notation introduced by Britter and Hanna (2003) where the obstacle frontal area per unit plan bottom area is Figure 11. Normalized terms in the spatially averaged momentum budget for low KC case 6 (KC 5 2) at four different wave phases (columns): xt 5 0 (flow reversal), xt 5 p=4, xt 5 p=2 (flow peak), and xt 5 3p=4. Plots are (a-d) acceleration term and (e-h) drag and added mass force terms. 33) and the solid volume fraction is A hem 5pD 2 =8 is the frontal area of an individual hemisphere, A T 5S 2 is the corresponding area of the bottom, / is the solid volume fraction, V hem 5pD 3 =12 is the volume of a single hemisphere, and V T 5S 2 D=2 is the volume of the layer containing the hemisphere. For the idealized topography used in this study, k f 5pD 2 =8S 2 and /5pD 2 =6S 2 for the layer between the bed and the top of the hemisphere.
The drag force (f d ) and added mass force (f a ) terms can be written and The rms of each force term is jf j5 The drag force term increases in importance as KC increases for two reasons. First, the rms of the dimensionless drag force term jf d j=U 1 x increases linearly with KC due to the quadratic form of the drag relationship. Second, the drag coefficient C D increases with KC because as the orbital excursion becomes larger relatively to the hemisphere, the flow separation becomes progressively stronger. If C a did not vary, then jf a j=U 1 x would not vary with KC ( Figure 14). The decrease in this term with increasing KC occurs because C a decreases with KC as the volume of the fluid affected by the presence of the hemisphere decreases with increasing KC.
The relative importance of these two terms can be estimated by a defined as At low KC, the drag coefficient C D is smaller than the added mass coefficient C a , therefore a ( 1. The added mass force dominates the dynamics, which suggests the following simplified momentum balance @ðh u w i2u 1 Þ @t 52 C a / 12/ @u 1 @t The added mass coefficient C a decreases as KC increases, therefore, for low KC the velocity in the canopy layer is in phase with the freestream velocity, and its magnitude is smallest for small KC.  At large KC (5-20), the drag force and the added mass force are both important. As KC increases, the intensity of the flow separation grows, the drag force term increases in importance relative to the added mass force term, and the flow in the wake is expected to deviate from the freestream velocity significantly. Because f d is 908 out of phase with the freestream acceleration, the drag force decreases the magnitude of the velocity in the canopy layer relative to the freestream velocity and also slightly shifts its phase. For the KC range in this study, the change in magnitude of the freestream velocity was small (less than 6%) and the phase shift was also small (less than 38). For larger KC, the within canopy velocity would decrease and the phase shift would increase with increasing KC until the shear stress terms began to dominate.

Friction Factor
The force per unit plan area (s eff ) exerted by the bottom under waves is typically represented using a wave friction factor, defined by Jonsson (1966): where s eff is the maximum force per unit bottom area in a wave cycle. There is some discrepancy in the literature as to how s eff and f w should be defined and interpreted. Traditionally, s eff is defined as the bottom shear stress of the wave boundary layer above the roughness elements (e.g., Grant & Madsen, 1979Jonsson, 1966;Smyth & Hay, 2002). However, f w is often computed from the dissipation of wave energy along a flume or the work done in driving an oscillatory flow. In wave dissipation studies, s eff represents the component of the total force exerted by the bottom that does work on the oscillatory flow (e.g., Bagnold & Taylor, 1947;Carsten et al., 1969;Simons et al., 1988). The inertial force is in quadrature with the water velocity and does not do work; therefore, in these studies, s eff is implicitly defined as the drag force per unit plan bottom area. In another class of studies, s eff is measured from the total force per unit area acting on a shear plate, and therefore includes both drag and inertial forces (e.g., Kamphuis, 1975;Mirfenderesk & Young, 2003). For wave orbital excursions (f) large compared with the roughness element scale (k s ), the inertial force is small compared with the drag force and the shear stress immediately above the roughness elements balances the drag on those elements; therefore, the three definitions of s eff are equivalent. This was not the case in the present study.
In our simulations, shear stresses were very small compared with forces associated with the pressure field around the hemisphere, and the inertial force was similar to, or larger than the drag force. The friction factor f w was calculated from our simulation results in two different ways for comparison with previous studies. First, we computed the friction factor associated with the drag force alone. This friction factor was computed from the drag coefficient C D using This friction factor parameterizes only the component of the total force that removes energy from waves. Second, we computed the friction factor associated with the total force per unit area exerted on the bed, which is the sum of drag force and inertial force. The rms of the total force is used to calculate the friction factor and computed as where a5 3 8 ffiffi 2 p p CDKC Ca 5 3CD 4 ffiffi 2 p Ca f D and f5U 1 =x is the wave orbital excursion. Friction factors derived from analysis of our LES cases are compared with laboratory data from previous studies, as well as empirical and theoretical curves proposed in the literature (Figure 15). We only included data from studies that measured k s and f w independently (Grant & Madsen, 1982;Mirfenderesk & Young, 2003;Sleath, 1986) and omit those that used expressions for f w as a function of f=k s to back out k s because they will necessarily fall on the curves that were used (e.g., Lentz et al., 2016;Mathisen & Madsen, 1996;Rogers et al., 2016). In sediment transport studies, sand-grain roughness k s is commonly parameterized by the median grain diameter D 50 (Fredsøe & Deigaard, 1992;Nielsen, 1992). In this study, we use the diameter of the hemisphere as the roughness scale k s for simplicity, which gives the ratio of wave orbital excursion to roughness scale f=k s 5KC=2p. Our LES results show that for f=k s < 1, the friction factor representing form drag (f w;d ) increases with f=k s . For wave orbital excursions that are small compared with roughness element size (f=k s < 0:5), flow separation is weak, roughness elements do not significantly affect the drag on the flow, and the friction factor is close to the value for a smooth wall. This is consistent with Van Rijn (2007) who found that bed forms with length scale much larger than the wave orbital excursion do not contribute significantly to drag. As f=k s increases, flow separation grows stronger, the drag on the flow increases, and the friction factor increases. The pattern is markedly different for the friction factor representing the total force (f w;t ). f w;t decreases monotonically as f=k s increases, and values from our LES study agree closely with the laboratory results of Mirfenderesk and Young (2003), and the empirical formula proposed by Nielsen (1992).
The contrast between our friction factor estimates highlights that, when f=k s is of order unity or smaller, it is extremely important to differentiate between the total force on the bed, the force that does work on the flow, and the shear stress above the roughness elements. Depending on the application, different friction factor definitions may be appropriate. For example, if one is interested in the initiation of sediment motion, one would consider the total force per unit area, parameterized by f w;t . If one is interested in the dissipation of wave energy, one should consider only the force component that is not in quadrature with the orbital velocity, parameterized by f w;d . While the two estimates are similar when f=k s ) 1, they are orders of magnitude different when f=k s is of order unity.
Most models that are commonly used to estimate f w (e.g., Grant & Madsen, 1979;Jonsson, 1966;Nielsen, 1992) were developed for f=k s ) 1, when the wave boundary layer thickness is large compared with roughness element height, shear stress at the top of the canopy layer is similar to form drag per unit area exerted by the bed, and the inertial force is small compared with form drag. For fully developed turbulent flow, the thickness of the turbulent wave boundary layer can be estimated d % ju Ã T=p, where u Ã 5 ffiffiffiffiffiffiffiffiffi f w =2 p U 1 is the friction velocity, and j is the von Karman constant. The turbulent wave boundary layer thickness is therefore larger than the roughness element height when f w > 1 2 j 22 ðf=k s Þ 22 , where k s represents the roughness element height (nonshaded region in Figure 15). None of the assumptions outlined above are true for f=k s of order unity or smaller. In this parameter range, models for the friction factor must explicitly account for form drag and inertial forces (e.g., Lowe et al., 2005Lowe et al., , 2007Monismith et al., 2015).
Because the details of flow separation, and the resulting forces, vary depending on the geometry of roughness elements, it is expected that friction factors for f=k s < 1 depend strongly on the three-dimensional structure of the bottom.
This study investigated the forces imposed on the flow by an array of identical, evenly distributed hemispheres, and results may differ for other geometries. Schlichting (1936) carried out a series of laboratory experiments to study the drag over roughness elements of different shapes. The roughness influences of a variety of objects, such as ribs (Mirfenderesk & Young, 2003), sand grains (Grant & Madsen, 1979;Mirfenderesk & Young, 2003), and ripples (Sleath, 1986) have also been examined in other studies. These studies reveal that forces exerted on oscillatory flows vary with the height, shape, density, and distribution of the roughness elements, which are all important to the characteristics of the surface roughness (van Rij et al., 2002;Van Rijn, 2007).
Results presented in this paper rely on identification of a single representative roughness length-scale. Natural topographies are typically multiscale and composed of irregular shapes, making identification of a single roughness length-scale difficult. The very high friction factors inferred from measurement of wave dissipation on reefs (Lentz et al., 2016;Monismith et al., 2015;Rogers et al., 2016) can likely be attributed partially to the distributed drag layer above the true bottom, and partially to the presence of multiple topographic length scales. This remains an outstanding problem to be further explored. Figure 15. Friction factor f w versus f=k s (i.e., the ratio of wave orbital excursion to the roughness length scale). Solid line represents the theoretical solution by Grant and Madsen (1982). Dashed curve is the empirical formula by Nielsen (1992). The parameter space is divided by the curve f w 50:5j 22 ðf=k s Þ 22 , below which the roughness height k s is larger than the wave orbital excursion. The inertia force dominates the dynamics when f=k s < 1, and the drag force dominates the dynamics when f=k s ) 1.

Summary and Conclusions
Our LES simulations of oscillatory flow over hemisphere arrays illustrate that the dynamics of wave-driven flow over high-relief topography are strongly dependent on Keulegan-Carpenter number (KC5U 1 T=D). When KC is small (1 < KC < 10), there is no strong flow separation because the wave excursion is small compared with the roughness element size. Flow instabilities (G€ ortler vortices) on the curved hemisphere surface, which occur near the flow reversal, contribute significantly to the generation of turbulence. In this parameter range, the inertial force dominates flow dynamics, form drag is small, and both the Reynolds stress and the dispersive stress are negligible. The inertial force becomes less important and form drag becomes more important as KC increases. For 10 < KC < 20, strong flow separation occurs, and turbulence is generated primarily in the wake. In this parameter range, both inertial and drag forces have important roles in the momentum budget. The wave dispersive stress, newly derived here, is the main mechanism for vertical momentum transfer between the canopy layer and the overlying flow, although the stress gradient term is still small compared with the drag and inertial force terms.
If the spacing between obstacles (S) is sufficiently small (Keulegan-Carpenter number based on the spacing KC S 5U 1 T=S > 2p), turbulent structures generated in wakes interact with adjacent solid obstacles. When these structures impinge on the hemisphere, they cause fluctuations in the pressure on its surface, and hence fluctuations in the resulting force. However, for the parameter range investigated in this study, flow sheltering effects are negligible and the dynamics are dominated by interactions of flow with each obstacle individually. The drag coefficient C D was insensitive to the spacing between the hemispheres. In contrast, the added mass coefficient C a showed strong dependence on the spacing and decreased as S/D decreases, because the volume of fluid surrounding each hemisphere is limited by S/D.
Friction factors computed from the drag force alone (f w;d ) were 1-2 orders of magnitude smaller than friction factors computed from the total force (f w;t ), which includes the drag force and the inertial force. Values of f w;d , which are appropriate when quantifying wave dissipation, increased as the ratio of orbital excursion to hemisphere diameter (f=k s ) increased, due to the increasing strength of flow separation. Our values of f w;d approached previously proposed empirical and theoretical curves at high f=k s , but deviated from those curves for f=k s Շ1. Values of f w;t , which are appropriate when quantifying the total force per unit area exerted on the bottom, decreased as f=k s increased, due to the decreasing magnitude of the inertial force. Our f w;t values for closely spaced hemispheres (S/D 5 1.5, 2) matched previously proposed empirical curves reasonably well for all f=k s . These results highlight the importance of distinguishing between the total force on the bottom, the drag force that removes energy from the flow, and the shear stress above the canopy layer when f=k s is of order unity or smaller. Models developed for f=k s ) 1 that assume a fully developed turbulent wave boundary layer above the roughness elements are not appropriate for small f=k s ; in this parameter range, models for f w must explicitly account for form drag and inertial forces.