## Abstract

Flow over complex terrain causes stress on the bottom leading to drag, turbulence, and formation of a boundary layer. But despite the importance of the hydrodynamic roughness scale *z*_{0} in predicting flows and mixing, little is known about its connection to complex terrain. To address this gap, we conducted extensive field observations of flows and finescale measurements of bathymetry using fluid-lensing techniques over a shallow coral reef on Ofu, American Samoa. We developed a validated centimeter-scale nonhydrostatic hydrodynamic model of the reef, and the results for drag compare well with the observations. The total drag is caused by pressure differences creating form drag and is only a function of relative depth and spatially averaged streamwise slope, consistent with scaling for *k*–*δ*-type roughness, where *k* is the roughness height and *δ* is the boundary layer thickness. We approximate the complex reef surface as a superposition of wavy bedforms and present a simple method for predicting *z*_{0} from the spatial root-mean-square of depth and streamwise slope of the bathymetric surface and a linear coefficient *a*_{1}, similar to results from other studies on wavy bedforms. While the local velocity profiles vary widely, the horizontal average is consistent with a log-layer approximation. The model grid resolution required to accurately compute the form drag is *O*(10–50) times the dominant horizontal hydrodynamic scale, which is determined by a peak in the spectra of the streamwise slope. The approach taken in this study is likely applicable to other complex terrains and could be explored for other settings.

## 1. Introduction

Flows over complex terrain create variable stresses on the boundary and the nonlinear generation of turbulence that can occur at multiple temporal and spatial scales. For example, flow within coral reef systems is complex because they have irregular, branching morphologies with reef topography varying at scales ranging from centimeters to kilometers (Rosman and Hench 2011). In shallow nearshore regions with rough surfaces such as reefs or rocky shorelines, bottom stress is often a significant term in the momentum balance and the primary form of dissipation loss. Thus, correct parameterization of the bottom stress is essential to understanding these flows (Monismith 2007). Bottom stress is typically represented using a bottom drag coefficient *C*_{D}. The bulk (or total) bottom drag *C*_{D,B} is a combination of two components, skin friction drag *C*_{D,τ} resulting from tangential stresses and pressure form drag *C*_{D,p} resulting from nonuniform pressure distribution on the boundary (Kundu and Cohen 2008).

Stresses on the bottom lead to the nonlinear generation of turbulence, which is ubiquitous in flows with high Reynolds numbers (Re). Turbulence causes not only mixing of tracers of environmental interest but also mixing of momentum and the formation of a boundary layer. For a well-developed turbulent boundary layer, an inertial sublayer region exists where mean velocities exhibit a logarithmic profile. Within this region, the mean velocity profile is related to the generation of turbulence by shear at the bed; the boundary layer follows the classical law of the wall, and roughness is often expressed as a hydrodynamic roughness scale *z*_{0} (e.g., Reidenbach et al. 2006; Kundu and Cohen 2008). However, on surfaces where the roughness elements take up a significant fraction of the water column, such as shallow coral reefs, it is unclear if these assumptions hold in these conditions (Lentz et al. 2017).

For flow over complex surfaces, at least three distinct types of roughness configurations exist, each with their own scaling law (Wooding et al. 1973). The first, *k*-type roughness (where *k* is a roughness height), is typical of sand grain roughness studies (Nikuradse 1933), where the horizontal spacing is typically somewhat greater than the element height. In this regime, unstable eddies form behind roughness elements that are shed into the flow above the boundary, resulting in a turbulent structure scaled to the height of the roughness elements (Grant and Madsen 1982; Perry et al. 1969; Jimenez 2004). The second type, *δ*-type roughness (where *δ* is the boundary layer thickness), results in the formation of stable eddies between roughness elements, which are typically spaced closer than their height. This causes “skimming flow” along the top of the roughness elements, and its turbulent structure is nearly independent of the roughness height (Wooding et al. 1973; Perry et al. 1969; Jimenez 2004). A third class, *k*–*δ*-type roughness, is characterized by the formation of eddies in the lee of the roughness elements and the reattachment of flow behind the elements, which requires a relatively small concentration of elements (Wooding et al. 1973). This type of roughness scales with both the height of the roughness elements, and their concentration is defined as the ratio of the frontal area of each element to the average horizontal surface area (Grant and Madsen 1982). Flow over wavy surfaces such as sand bed ripples, as well as over many coral reef forms (as will be demonstrated in the results), is of the *k*–*δ* type.

In circulation models, turbulence is often parameterized in Reynolds-averaged Navier–Stokes (RANS) approaches, which require a parameterization of the bottom roughness. New field-scale technology to measure bottom bathymetry is advancing rapidly with increasing resolution, including sonar, hyperspectral, lidar, and fluid-lensing techniques. Measured variability in bathymetry is now often available at scales much smaller than the resolution of the computational grid. Typically, the mean depth value is taken to represent depth in the grid cell, and the remainder of the data is not used.

Significant work has been conducted to measure site-specific hydrodynamic roughness values on reefs (Rosman and Hench 2011; Reidenbach et al. 2006; Rogers et al. 2017, 2016, 2015; Lentz et al. 2016; Hearn 1999; Symonds et al. 1995; Lowe et al. 2005; McDonald et al. 2006; Thomas and Atkinson 1997; Falter et al. 2004), but a clear connection to the reef surface complexity is lacking. Several studies have suggested characterizing a reef surface based on fractal dimension (Hearn 2011a,b), or spectral energy of rugosity (Nunes and Pawlak 2008; Jaramillo and Pawlak 2011), but they lack a connection to measured hydrodynamic roughness. Thus, the state of current practice for computing bottom stress is to estimate a hydrodynamic roughness value based on general site parameters, based on literature from similar systems, or calibrated from field measurements. To our knowledge, there is yet no unified theory for how turbulent flow over complex terrain connects to hydrodynamic roughness or how multiple spatial terrain scales interact to create the total drag.

The aim of this paper is to address the fundamental question of how complex terrain is related to the hydrodynamic roughness, including what horizontal and vertical length scales are important. To address this question, we conducted an extensive field campaign to measure flow on a coral reef and developed a nonhydrostatic hydrodynamic model of the reef based on detailed bathymetric measurements using fluid-lensing techniques (section 2). The field observations are used to compute the spatially varying drag coefficients (section 3). Based on the validated model, we present a simple model to predict total drag on complex terrain using statistics of the surface (section 4). We conclude with an overview discussion (section 5) and summary (section 6).

## 2. Methods

### a. Hydrodynamic measurements and data analysis

Ofu Island, American Samoa, lies in the southwest Pacific Ocean (14.2°S, 169.7°W) (Fig. 1a). Because of the tidal range and shallow reef crest, the onshore wave-driven flow is forced by breaking waves and modulated by the tides (Koweek et al. 2015). The focus of this study is on pool 400 located on the south shore of Ofu Island within the National Park of American Samoa (Fig. 1b).

The field experiment consisted of an array of velocity, pressure, and temperature sensors deployed from 10 to 28 March 2017 designed to characterize the waves and hydrodynamics of pool 400 (Fig. 1b; appendix). Measurements of velocity data are taken in east, north, and up coordinates and transformed to a local coordinate system with *x* and *y* in the alongshore (AS) and cross-shore (CS) directions and *z* upward from mean sea level (MSL) (Table 1). The instantaneous measurements are time averaged (15 min) to give Eulerian velocity **u**_{E}(*u*_{E}, *υ*_{E}, *w*_{E}), free surface deviation from MSL *η*, and wave statistics. Because of very small waves in pool 400 (Stokes drift ≪ **u**_{E}), we assume the mean Lagrangian velocity **u** = **u**_{E}. The Lagrangian depth-averaged mean velocity **U**(*U*, *V*, *W*) is calculated by combining data at a given location assuming **u**_{E} = 0 at the bottom and averaging over the depth. Additional details are in the appendix.

### b. Bathymetric UAV measurements

An airborne survey was conducted over pool 400 in July 2016 using unmanned aerial vehicles (UAV) and fluid-lensing computational methods (Chirayath and Earle 2016; Suosaari et al. 2016; Chirayath 2016) to resolve centimeter-scale bathymetry of the reef (Fig. S1 in the supplemental material). A UAV electric quadcopter platform was custom built to host a nadir-pointing high-frame-rate video camera, relay synchronized position data and survey a region with sequential flights, each up to 20 min in duration. Videos frames were sorted into 120-frame bins and processed using the experimental fluid-lensing algorithm (Chirayath and Earle 2016) to remove refractive distortions caused by ambient surface waves. The corrected images and UAV position data were used as input frames for structure from motion to produce 2D, centimeter-scale orthophotos and a dense 3D bathymetry model. Calibration targets were distributed at varying water depths for georeferencing and bathymetry validation. Finally, terrestrial and millimeter-scale underwater gigapixel photogrammetry was performed to calibrate and verify 2D fluid-lensing reconstructions from airborne data, perform georectification, and validate derived 3D bathymetry. Where high-frequency Gaussian noise was present in the raw bathymetry as a result of resolution limitations, depth grids were low-pass filtered to 10- to 30-cm grid scale.

### c. 3D nonhydrostatic model

The Stanford Unstructured Nonhydrostatic Terrain-following Adaptive Navier–Stokes Simulator (SUNTANS) modeling framework is a three-dimensional nonhydrostatic RANS model under the Boussinesq approximation (Fringer et al. 2006). The model grid is implemented in a rectangular horizontal grid with *z*-level vertical coordinate with domain length *L*, width *W*, average depth 〈*D*〉, and model grid spacing Δ*x*_{G}, Δ*y*_{G}, Δ*z*_{G}. To compute the turbulent viscosity *ν*_{T}, a 3D *k*–*ϵ* turbulence model is employed within the generic length scale framework (Warner et al. 2005) with molecular viscosity *ν* = 1 × 10^{−6} m^{2} s^{−1}. Time stepping is third-order Adams–Bashforth, and nonlinear momentum advection is central differencing. Initial conditions are zero velocity with uniform density. To achieve an equilibrium boundary layer, boundary conditions are periodic in the streamwise (*x*) direction and free slip on the lateral boundaries and free surface. The bottom boundary condition is a quadratic drag law with roughness height *z*_{0,ur}, imposed on the bottom grid cell (Fringer et al. 2006). For reef simulations, we assume *z*_{0,ur} = 0.033 cm, a value typically used for sand. The model is forced with a uniform pressure gradient, which is adjusted in time to achieve a desired domain-averaged velocity *U*_{B} (Nelson and Fringer 2017). Model simulations are run to steady state, typically at least five flow-through periods. The bathymetry for the reef simulations is linearly detrended in *x* and *y*, to minimize edge effects, and is smoothed along the periodic edges using , where the tapering length scale *b =* 2 m. Simulations were conducted on a 64-core AMD workstation.

## 3. Field measurement results

### a. Hydrodynamic conditions

During the observation period, tidal variation of the free surface *η* is ±0.5 m, winds are light at less than 5 m s^{−1}, and wave forcing on the forereef consists of several long-period swell events with *H*_{rms} exceeding 1 m that propagate toward the reef crest and break; and thus *H*_{rms} in the lagoon is quite small (Figs. 2a–d). The effect of tidally modulated wave forcing creates a pressure gradient that drives a tidally varying flow in the pool center (Fig. 2e). The drifter tracks show the large-scale flow field from the reef crest exits through low areas in the reef crest and channels (predominantly near H0), and the general flow patterns are consistent over tidal phase on multiple days (Fig. 1b).

### b. Momentum balance

In shallow depths, the bottom stress is often approximated with (Grant and Madsen 1979; Feddersen et al. 2003; Lentz et al. 2017)

where *C*_{D} is a local nondimensional drag coefficient that may depend on the flow environment and bottom roughness and *ρ* is the fluid density.

Flow on shallow reefs is governed by the depth-integrated momentum equations for **U**(*U*, *V*) given by (e.g., Mei et al. 2005)

where time-varying depth , *h* is mean depth, is the mean free surface, **S** is the radiation stress tensor, *g* is gravitational acceleration, is the mean surface stress, and is a time average of function *f*. In the alongshore (*x*) direction, neglecting *V*∂*U*/∂*y*, substituting (1), assuming the water surface slope is linear, flow **q** = **U***D*, and by continuity the flow *q*_{x} does not change significantly along *x*_{1} to *x*_{2}, and taking the spatial average of (2) from point *x*_{1} to *x*_{2}, ,

From left to right, the terms in (3) will be referred to as unsteady (US), nonlinear advection (NL), pressure gradient (PG), radiation stress gradient (RSG), bottom stress (BT), and surface stress (ST), where the local *C*_{D} given by (1) becomes the spatially averaged drag coefficient . For field observations on shallow reefs, the momentum equation in integral form [(3)] is shown to have increased accuracy over the differential form [(2)] owing to its inclusion of the spatial depth variability between measurement points (Lentz et al. 2017).

In a turbulent channel flow where roughness elements are small compared to the depth, the velocity profile is well represented by the log-layer approximation:

where the shear velocity , the height above the bottom , *z*_{0} is a hydrodynamic roughness scale, *κ* is von Kármán’s constant (0.41), and Π is Cole’s wake strength (~0.2 for high Reynolds number). Integrating (4) over the depth and substituting (1) (Lentz et al. 2017),

which is valid for *D*/*z*_{0} > 10. Given a velocity measurement bounded by two depth measurements and a horizontal profile of depth, (3) and (5) can be solved to determine an average *z*_{0} for the profile (Lentz et al. 2017).

The mean water surface offset *η*_{0} between the offshore and each site in the lagoon is based on solving (3) for low forcing conditions (see appendix) and shows a steep increase from the shallow outlet at K0 and a more gradual increase within the pool from site G0 to A0 (Fig. 3a). The time-averaged flow magnitude |**q**| generally increases along the pool to a maximum at H0 where much of the flow exits the pool (Figs. 1b and 3b).

The results for the spatially averaged 〈*C*_{D}〉 from the momentum balance vary along the transect between *O*(0.01) and *O*(0.1) (Fig. 3c). Comparing to the results using a logarithmic fit to the velocity profile in (4) and the Reynold’s stress method ( appendix, Eq. A2), 〈*C*_{D}〉 are similar magnitude but fall outside the 95% confidence limits (Fig. 3c). The likely explanation for this is that 〈*C*_{D}〉 is averaged over *O*(100) m, while the log fit and Reynold’s stress methods are local measurements. The spatially averaged roughness scale *z*_{0} from the momentum balance varies from *O*(1) to *O*(10) cm along the transect, and results from the local log fit are similar (Fig. 3d). Time-average velocity profiles normalized by the depth-averaged velocity show a logarithmic profile within the midwater column (Fig. S2). This confirms the applicability of (4) for determining a local roughness value. At the two deeper sites (B0, D0) the profile is logarithmic over the measured depth, but at the two shallower sites (F0,H0) the profile deviates from a logarithmic profile below 0.3 and above 0.6 , likely caused by the large roughness elements relative to depth. It is possible that some of the deviations are due to inaccurate reference height , which on reefs is sometimes accounted for with an additional offset variable in (4) (Rosman and Hench 2011).

The magnitude of the terms in (3) shows the leading terms are between pressure gradient (PG) and bottom stress (BT), while the NL and RSG terms may become of secondary importance in some locations. (Fig. 3e). Thus, in this study site, (3) can be simplified to

For this work, we ignore the effects of waves on creating a rougher boundary layer and influencing *C*_{D}. For this site, with mean flows *U* of *O*(0.1) m s^{−1}, the relatively small wave heights in the pool (Fig. 2) create wave excursion distances *A*_{b} of *O*(0.2) m, and wave velocities of *O*(0.1) m s^{−1}. Considering the physical scale of roughness elements *h*_{b} is typically 1 m, the relative excursion *A*_{b}/*h*_{b} < 1, and ratio of wave velocities to mean currents is *O*(1) or less. For flows in this regime, the effect of waves may create a rougher boundary layer height (Lentz et al. 2017). However, considering the high uncertainty for some of the field sites in the results for *z*_{0} (Fig. 3d), we consider this a second-order effect.

## 4. Theory and model results

### a. Theoretical approach

To better understand the detailed flows in this environment and extend beyond the limited field observation sites, we model the complex flows on the reef using the detailed bathymetric measurements. We apply the SUNTANS modeling framework in periodic boundary conditions forced by a uniform pressure gradient acceleration *F*_{PG} in an unstratified fluid in a nonrotating frame. The momentum and continuity model equations are

and

where **u** = (*u*, *υ*, *w*) is the velocity vector, *p* is the dynamic pressure composed of a nonhydrostatic and hydrostatic component (*p* = *p*_{nh} + *ρ*_{0}*gη*), and background density *ρ*_{0} = 1000 kg m^{−3}.

Upon time (*t*) averaging, domain volume (∀) averaging, invoking Gauss’s theorem, and assuming steady state, the *x* component of (7) becomes the following:

where *f*_{τ,b} is the time-average bottom shear stress and *f*_{p,b} is the stress from mean pressure acting on the bottom (indicated by the subscript *b*). Note this is an analogous statement to (6), that is, that a mean pressure forcing balances stresses on the bottom. However, (9) separates the stresses on the bottom into shear stress and pressure stress, which creates form drag. The *F*_{τ,b} is the acceleration due to the mean velocity parallel to the bottom:

where the bottom following coordinate system (**n**, **s**) are the unit normal and unit parallel bottom surface vectors, respectively. The *F*_{p,b} is due to the bottom pressures acting normal to the bottom and computed from (9). Normalizing (9) by density times the mean depth over twice the kinetic energy gives drag coefficients due to the total (bulk) flow, shear stress, and form drag in the streamwise direction:

### b. Model validation and results

While RANS models have shown effectiveness in computing hydrodynamic flows and drag coefficients at a range of scales, previous studies note some limitations in computing the details of turbulence and eddy formation especially for bluff bodies where the details of flow separation are important (Rodi 1997; Lübcke et al. 2001). To validate the hydrodynamic model, a series of three idealized simulations on a cylinder, cube, and sinusoidal bottom are conducted with known experimental results (Table 2) in addition to validation on the Ofu reef presented next (Table 3).

The simulations are conducted at Re similar to those expected on the field site, where mean velocities are on average *O*(0.1) m s^{−1} (Fig. 2). The horizontal width *w* of the dominant bottom features are *O*(1) m (Figs. 4 and 5). Thus, at the site, the average Re_{w} = (*U*_{B}*w*)/*ν* = 1 × 10^{5} but can vary from 6 × 10^{4} to 1 × 10^{6}.

To compare to results from the literature, *C*_{D,B} for the cylinder and cube simulations are computed by , with total drag force and projected frontal area *A*. For the cylinder simulation Re_{w} = 1 × 10^{5}, and *C*_{D,B} is within 5% of the published value (Kundu and Cohen 2008) (Table 2). For higher Re_{w} = 5 × 10^{6}, where we expect transition to a turbulent boundary layer and decrease in the net drag, the model result overpredicts *C*_{D,B} by 42% (Roshko 1961) (Table 2). For the cube simulation, the model result for *C*_{D,B} is within 5% of the published value and not a function of Re_{w} (Bearman 1971) (Table 2). For the sinusoidal bathymetry Re_{D} = (*U*_{B}*D*)/*ν* = 4.2 × 10^{5} and using (11) to compute 〈*C*_{D,B}〉, which is within 18% of the experimental result (Gong et al. 1996) and within 36% of the LES simulation result (Salvetti et al. 2001) (Table 2).

Thus, the model accurately computes the drag coefficients for a range of rounded and sharp surfaces at Re typical of the study site. For forms with rounded features such as a cylinder at Re_{w} greater than 5 × 10^{5}, the model may overpredict the drag coefficient. However, the typical Re_{w} on this reef is 1 × 10^{5}, below this threshold. In addition, many of the forms on the reef have sharp edges and the drag coefficient is not expected to vary with Re.

The reef surface is modeled using a similar approach as the idealized domains (Table 3). Model results at site D0 with rounded isolated roughness elements show strong horizontal variability in the flow field near the bottom with separated wakes in both the horizontal and vertical (Figs. 4a,d). The aerial image shows the roughness elements at this site are Porites corals surrounded by sand (Fig. 4b). The pressure field shows corresponding high pressures in the front of the roughness elements (some from the free surface and the remainder from nonhydrostatic pressure), with negative pressures in the lee-separated wake regions, leading to a net form drag (Figs. 4c,f). Eddy diffusivity is highest in the interior and is highly turbulent *O*(10) cm^{2} s^{−1} (Fig. 4e), similar to field measurements on other reefs (Reidenbach et al. 2006). The separated wakes are qualitatively similar in flow direction and extent of separation to detailed field observations on other reefs (Hench and Rosman 2013).

Model results at site F0, which has more irregular bathymetry, are similar, but the wakes behind the roughness elements and resulting pressure field are larger and dominated by the large-scale features (Figs. 5a,c,d,e,f). The aerial image shows this site has smaller irregular coral shapes with some sand patches (Fig. 5b).

The model results for 〈*C*_{D,B}〉 generally show good agreement with the field observations, although with a high degree of scatter (Fig. 6). The high degree of scatter is likely due primarily to the spatial scales of the different results. The model simulations are conducted at 20–32-m spatial scale owing to model time and bathymetric limitations. The field results are averaged over *O*(100) m for the results from the momentum equation, but are local measurements for the logarithmic fit and Reynolds stress methods. Thus, while there is scatter in these results, the validation confirms that this model reasonably computes 〈*C*_{D,B}〉 on these reef surfaces and can be applied to additional model scenarios and reef locations.

### c. Drag

Given shallow steady unidirectional flow over a periodic wavy bathymetric surface, the relevant variables are a resisting drag force over area (stress) *f*_{D}, depth-averaged flow *U*_{B}, average depth *D*, characteristic bottom height *h*_{b}, characteristic bedform wavelength in the direction parallel and normal to flow (*λ*_{x}, *λ*_{y}), fluid density *ρ*, and fluid viscosity *ν*. Per Buckingham pi theorem a common way to nondimensionalize this problem is the following (Grant and Madsen 1982; Wooding et al. 1973):

At high Reynolds numbers (Re_{h} = *h*_{b}*U*_{B}/*ν* > 10^{6}), the flow is in the rough turbulent regime and the *C*_{D,B} should be independent of Re_{h}, while in the turbulent transition regime, *C*_{D,B} should vary like (Schlichting 1979). For high relative depth (*D*/*h*_{b}) the flow profile is a well-developed log layer and (4) applies, while at low relative depth the roughness elements take up a significant part of the water column and the flow approaches a porous media flow (Rosman and Hench 2011; Lentz et al. 2017). At small bedform steepness in the direction of flow (*h*_{b}/*λ*_{x} < 0.1) there is no wake separation behind the bedforms and the drag is primarily due to skin friction, while at high bedform steepness wake separation is pronounced and the drag is primarily due to pressure form drag (Grant and Madsen 1982; Schlichting 1979). The last parameter *h*_{b}/*λ*_{y} is a spanwise steepness factor and governs how much of the flow is blocked; *C*_{D,B} is at a maximum for high steepness factors (2D bedforms normal to flow) as a result of high pressure form drag, and *C*_{D,B} decreases to skin friction for low steepness factors approaching 0 (2D bedforms aligned with flow).

We approximate the complex reef surface as multiple wavy bedforms superimposed together. The spatially averaged drag is then given by the spatial average of (12):

where 〈*D*〉 is the volume-averaged depth and Re_{h} = 〈*h*_{b}〉*U*_{B}/*ν*, and we have assumed *z*′/〈*D*〉 ≪ 1 for the second term to be uncorrelated. The terms in (13) can be approximated assuming equivalent parameters to a wavy sinusoidal surface using the root-mean-square of the surface and the bottom slope:

where is the deviation of the bottom surface, *D* is detrended (no mean or mean slope), and the spatial root-mean-square operator is applied over the model domain. For this study, (14) and (15) are applied to the bathymetry at the same grid resolution as the model (Δ*x*_{G}). The model grid depths are taken as the average of all subgrid-scale points from the high-resolution bathymetry (Δ*x*_{B}).

The model is used to explore the variability of the terms in (13). The Reynolds number Re_{h} is varied over the range of observed values from 10^{4} to 10^{5} by changing *U*_{B} and, as expected, has very little effect on *C*_{D,B} (Fig. 7). Since the flow is in a turbulent regime, the irregular forms of the reef create separated boundary layers (Figs. 4 and 5), which are not dependent on Re_{h}.

The relative depth 〈*D*〉/〈*h*_{b}〉 is varied by changing 〈*D*〉 at sites D0 and F0 over the range of observed values. Both sites show decreasing 〈*C*_{D,B}〉 with increasing 〈*D*〉/〈*h*_{b}〉, consistent with the theory in (5) (Fig. 8). This trend is primarily a result of 〈*C*_{D,p}〉 changing, while 〈*C*_{D,τ}〉 remains very small over these simulations. Equation (5) matches the model results well for high relative depth (site D0) while somewhat underpredicting the drag at smaller relative depth (site F0).

The variation in relative slope 〈*h*_{b}/*λ*_{x}〉 is achieved by simulating many different model sites, while holding Re_{h} and 〈*D*〉/〈*h*_{b}〉 constant (the same as site D0 base case) (Fig. 9). The 〈*C*_{D,B}〉 increases linearly for 〈*h*_{b}/*λ*_{x}〉 > 0.05 because of an increase in *C*_{D,p}. For idealized sinusoidal forms, flow separation and high form drag occur at *h*_{b}/*λ*_{x} > 0.1 (Salvetti et al. 2001), but the reef model results show a dominance in form drag for 〈*h*_{b}/*λ*_{x}〉 > 0.05, likely since the irregular reef bedforms create more flow separation at a lower average slope. At very low slopes, 〈*C*_{D,p}〉 trends toward zero and 〈*C*_{D,τ}〉 becomes the dominant term in (13). These simulations are of a highly smoothed spatially downsampled reef (Δ*x*_{G} > 3 m) with the model grid resolution the same as the bathymetry. Thus for this reef, form drag due to pressure 〈*C*_{D,p}〉 is the dominant term at all sites.

In Fig. 9, there is some variability in 〈*h*_{b}/*λ*_{y}〉 that cannot be independently varied. But, because most reef surfaces at this scale are symmetric, 〈*h*_{b}/*λ*_{x}〉 and 〈*h*_{b}/*λ*_{y}〉 are typically of similar order (Fig. 10a). Thus, 〈*h*_{b}/*λ*_{y}〉 is not expected to be a significant parameter in predicting 〈*C*_{D,B}〉, a simplification also used in studies on sediment bedforms (Wooding et al. 1973; Grant and Madsen 1982). However, for forms that are not symmetric it could become an important factor.

Finally, in all the simulations we have assumed unidirectional flow, and the factor 〈*h*_{b}/*λ*_{x}〉 is not dependent on the sign of the flow. However, it is well known that, because of flow separation, a bluff body will have less form drag if it has an abrupt front and gently sloping lee rather than a gentle front and abrupt lee (Kundu and Cohen 2008). We vary the flow direction on a square domain (32 m × 32 m) at site D0. The 〈*C*_{D,B}〉 varies slightly with flow direction by up to 7% (Fig. 10b). As a result, we consider flow sign a minor effect for these reef surfaces, which are primarily symmetric.

### d. Hydrodynamic length scales

For a well-developed turbulent flow, the characteristic hydrodynamic length scale for *k*–*δ*-type roughness (Wooding et al. 1973) should be proportional to the following (Grant and Madsen 1982; Nielsen 1992):

where *h*_{b} is the roughness height and *h*_{b}/*λ*_{x} is the roughness spacing or the roughness slope. A reasonable parameterization for (17) for periodic bedforms is given by *z*_{0} = *a*_{1}*h*_{b}(*h*_{b}/*λ*_{x}), and previous experimental results found *a*_{1} equal to 0.27 for sand ripples (Nielsen 1992), 0.36 for triangular ripples (Jonsson and Carlsen 1976), 0.50 for the atmospheric boundary layer (Lettau 1969), and 0.83 (Swart 1977) and 0.92 (Grant and Madsen 1982) for oscillatory flow over ripples. For complex terrain, we propose to approximate the roughness height by 〈*h*_{b}〉 and the roughness concentration by 〈*h*_{b}/*λ*_{x}〉, and (17) becomes

where is the contribution from the unresolved hydrodynamic roughness, at scales smaller than the measurement resolution. Substituting spatial averages into (5),

Thus, (18) and (19) contain the important nondimensional parameters in (16) concluded from the model results. While the current study uses downsampled bathymetry with a grid spacing equal to the model grid (Δ*x*_{G} = Δ*x*_{B}) to estimate *a*_{1}, in practice, for estimating hydrodynamic roughness this method should apply (18) to bathymetric data at the highest available resolution, which should be higher than the model resolution. For irregular surfaces, a vertical step in the high-resolution bathymetry would give the singularity 〈*h*_{b}/*λ*_{x}〉 = ∞ and the unphysical result 〈*z*_{0}〉 = ∞. In this case, the fundamental assumption of a wavy surface is violated, and some spatial smoothing may be required to apply (18).

An analysis of all model runs on the Ofu reef with a best fit of 41 model simulations to (18) and (19) gives *a*_{1} = 0.38 with a coefficient of determination *R*^{2} of 0.79 (Fig. 11). The results in Fig. 11 are based on results within 6 × 10^{3} < Re_{h} < 7 × 10^{4}, 1.1 < 〈*D*〉/〈*h*_{b}〉 < 6.3, and 0.7 < 〈*h*_{b}/*λ*_{x}〉 < 0.29. However, (18) is likely valid as long as the flow is turbulent; additionally the upper limit of 〈*D*〉/〈*h*_{b}〉 approaching a pure boundary layer, and the lower limit of 〈*h*_{b}/*λ*_{x}〉 approaching a flat bottom are likely not limiting. The results for 〈*z*_{0}〉 compare well to other studies on ideal sinusoidal forms (Salvetti et al. 2001; Zilker and Hanratty 1979; Buckles et al. 1984; Henn and Sykes 1999) and triangular ripples (*a*_{1} = 0.36) (Jonsson and Carlsen 1976) but are somewhat higher than 〈*z*_{0}〉 for sand ripples (*a*_{1} = 0.27) (Nielsen 1992).

The vertical profiles of horizontal velocity at sites D0 and F0 show high horizontal variability in *u* but relatively similar profiles in horizontally averaged 〈*u*〉_{H} when normalized by *U*_{B} above the mean depth (Fig. 12a), where is upward from the mean depth. Values of 〈*u*〉_{H} below the mean depth tend toward zero within the reef canopy and are highly variable depending on the reef bathymetry. Scaling all reef model simulations by and scaling all horizontally averaged velocity 〈*u*〉_{H} by *U*_{B} collapses the data into self-similar profiles (Fig. 12b). The profiles are logarithmic within the middle of the water column, as predicted by theory in (4). Thus, while the local velocity is highly variable, a spatial horizontal average 〈*u*〉_{H} is well approximated by a log-layer parameterization above the mean depth.

A final question of interest is what horizontal spatial scales are most important on the reef for creating drag. A series of simulations at sites D0 and F0 are conducted by varying the bathymetry grid resolution Δ*x*_{B} from 10 cm to 24 m, where the 24-m resolution is flat. Model grid resolution Δ*x*_{G} is held constant at 10 cm, and a cubic spline method is used to interpolate from the coarser bathymetry grid to the finer model grid. The total drag 〈*C*_{D,B}〉 increases and then levels off with increasing grid resolution with 〈*C*_{D,p}〉 the primary component of this increase (Figs. 13a,b). At the highest modeled resolution, 〈*C*_{D,B}〉 approaches the field observed values, but as noted previously, there are inherent differences in the spatial averaging of these different methods. For the modeled bathymetry, the important result is that 〈*C*_{D,B}〉 is approaching an asymptotic value, where enhanced grid resolution is contributing little to the drag.

We use Fourier analysis of the bathymetry in the *x* direction with wavelength *λ* and the spectral energy density (SED) *S*(*λ*)averaged in the *y* direction applied at different bathymetric grid resolutions Δ*x*_{B}. SED of the depth shows highest spectral density for the large-scale features and decreasing density for the small-scale features (Figs. 13c,d). SED of the streamwise bottom slope shows increasing magnitude with smaller grid scales and a peak wavelength *λ*_{x} of approximately 4.5 and 1.1 m at sites D0 and F0, respectively (Figs. 13e,f), the same scale at which 〈*C*_{D,p}〉 begins to increase (Figs. 13a,b). Qualitative inspection of the coral forms in the plan view shows the spacing is of similar order (Figs. 4 and 5). Thus, appears to be a good quantitative metric to identify the dominant horizontal length scales creating hydrodynamic drag. The grid spacing required to adequately resolve the wake effects and thus accurately model 〈*C*_{D,p}〉, appears to be approximately 11 (F0) to 45 (D0) times smaller than the spectral peak wavelength of (Fig. 13).

The mean dynamic pressure field at site D0 shows decreasing differences in for larger grid scales (Fig. 14), which leads to the decrease in 〈*C*_{D,p}〉 (Fig. 13a). However, while there are some enhanced details in the near bed at 10-cm resolution, the primary pressure features and magnitude remain similar between 10- and 30-cm resolution (Figs. 14a–c). The physical reason for this is that the large-scale features dominate the flow field, and smaller features are contained within the wakes and boundary layer created by the dominant scale features, a characteristic feature of *k*–*δ*-type roughness (Figs. 4 and 5). Thus, once the dominant hydrodynamic features are resolved with *O*(10–50) grid cells, enhanced resolution has little effect on the drag.

## 5. Discussion

Based on field observations from a shallow reef in American Samoa, we compute spatially averaged drag coefficients at several sites, which are averaged over *O*(100) m. Based on the momentum balance, we conclude that within the pool the primary terms in the depth-averaged momentum equation are between pressure gradient and bottom stress, a result which is well known in these conditions (Monismith 2007). The pressure gradient develops as a result of breaking waves on the forereef, which are tidally modulated, and drive the unidirectional alongshore flow in the pool. This relatively simple flow in the absence of waves provides the ideal setting to better understand the role of bottom bathymetry of the complex reef in creating hydrodynamic roughness.

Using high-resolution UAV fluid-lensing techniques, we mapped the bottom bathymetry of the pool and used this data to model the reef using the SUNTANS nonhydrostatic hydrodynamic model. The results from the model for bottom drag coefficient 〈*C*_{D,B}〉 are in excellent agreement with idealized test cases and generally in good agreement with the field observations, with some scatter likely due to variations in spatial averages from the different methods. The choice of using a RANS model allows for modeling field-scale sections of the reef with moderate computational cost but drastically simplifies the turbulent properties of the flow. Work to extend this analysis using a large-eddy simulation (LES) code could address some of these issues and potentially improve the model results but would require very large computational resources. However, for the purposes of this study, primarily to predict total drag, the selected model appears to be adequate.

The modeled hydrodynamic flows on the reef illustrate the complexity of flows over these surfaces. Wakes form in both the horizontal and vertical, which create imbalances in the pressures, and generally reattach behind the larger reef forms, confirming that many reefs can be characterized as having *k*–*δ*-type roughness. These pressure differences lead to form drag 〈*C*_{D,p}〉, which is the primary cause of total drag in the model results. Drag due to shear stress 〈*C*_{D,τ}〉 is a minor term except for nearly flat model geometries, a condition that is very rare in reef environments.

For parameterizing the bathymetry and dimensional analysis, we approximate the reef surface as a superposition of wavy sinusoidal bedforms with potentially four nondimensional parameters governing the flow. These include a roughness Reynolds number , relative depth , slope in the streamwise direction 〈*h*_{b}/*λ*_{x}〉, and slope in the spanwise direction 〈*h*_{b}/*λ*_{y}〉. We approximate these terms by spatially averaging over the model domain, and based on variation of parameters in the model, we conclude that the bulk drag is a function of the relative depth and slope in the streamwise direction only. The direction of flow may be important for certain forms that are not symmetric, the effects of which are not considered in this method.

Roughness on wavy bedforms such as reefs are characterized by *k*–*δ*-type roughness, which scales with a roughness height and a roughness concentration. Drawing on previous work on sediment bedforms and idealized studies, we propose a parameterization for the bottom roughness scale 〈*z*_{0}〉 [(18)], which is a function of the roughness height based on the spatial rms of the surface 〈*h*_{b}〉, the roughness concentration based on the spatial rms slope 〈*h*_{b}/*λ*_{x}〉, the coefficient *a*_{1,} and the unresolved roughness . Combining this equation with an expression for the spatially averaged drag coefficient 〈*C*_{D,B}〉 [(19)], a best fit to all model runs, we obtain an *R*^{2} = 0.79 and coefficient *a*_{1} = 0.38. The results for *z*_{0} compare well to previous studies on idealized sinusoidal and triangular bedforms (*a*_{1} = 0.36).

The result in (18) is likely applicable for a range turbulent flows over complex terrains with *k*–*δ*-type roughness where relative depth and average slope 〈*h*_{b}/*λ*_{x}〉 < 0.29. For very shallow depths , the horizontally averaged log-layer approximation likely breaks down and (18) and (19) are no longer valid (Lentz et al. 2017). For very deep cases as → ∞, the flow becomes a pure boundary layer and (18) is valid, but the depth-averaged approach in (19) is likely not applicable as other forces (e.g., surface, buoyancy, rotation) may become important above the boundary layer.

As 〈*h*_{b}/*λ*_{x}〉 increases, flow will not reattach between elements, and a larger-scale boundary layer will form dependent only on 〈*h*_{b}〉, characteristic of -type roughness. Thus, the transition from the present study limits 〈*h*_{b}/*λ*_{x}〉 = 0.29 to this new regime will likely occur around 〈*h*_{b}/*λ*_{x}〉 = 0.7 (idealized packed spheres), and further studies into this transition would be of great interest. For very small 〈*h*_{b}/*λ*_{x}〉, flow may not separate behind roughness elements, leading to a smaller than predicted by (18). This transition is pronounced for idealized sinusoidal bedforms for 〈*h*_{b}/*λ*_{x}〉 < 0.1 but is not observed in the reef model results as irregular forms enhance flow separation. An alternate interpretation for the variation in roughness regimes described above is to instead divide regimes into *k*- and *d*-type roughness, which account for roughness density (or solidity) (Jimenez 2004) and should be generally equivalent to varying 〈*h*_{b}/*λ*_{x}〉.

For geometry that is highly resolved, represents the contribution from shear drag and should tend toward zero with increasing resolution. For coarsely resolved geometry, represents a parameterization of the form drag, which must be estimated.

While the local velocity profile on the reef is highly variable, the spatially averaged flow is logarithmic above the mean depth and justifies using a log-layer approximation for modeling larger-scale flows. Additionally, while only the slope in the streamwise direction is considered in (18), the spatial averaging does account for some variations in the spanwise direction, since the spatial average is applied in both horizontal directions over some areas. Waveforms that are parallel to flow would yield 〈*h*_{b}/*λ*_{x}〉 = 0 and thus *z*_{0} = 0, while flow perpendicular to the same bedforms would give nonzero 〈*h*_{b}/*λ*_{x}〉 and *z*_{0}.

Finally, we investigate which horizontal spatial scales are relevant to creating total drag on the reef. Based on model simulations where we hold the model resolution constant but change the bathymetry resolution using an average of the high-resolution subgrid-scale bathymetry, we find that the drag approaches an asymptote for small scales, implying that a moderate horizontal length scale dominates the drag. This also implies that the very finescale features do not contribute much to the overall pressure field and thus the total drag, the result of the large-scale features creating wakes that isolate the small features from the main flow. This dominant horizontal spatial scale predicted by energy spectra of the slope showed a maximum around the same horizontal length scale of the large roughness features of the reef. A model resolution of *O*(10–50) grid cells per this dominant spatial scale is required to resolve the flow separation and accurately model the form drag on the reef.

## 6. Summary

We conducted extensive field observations on a coral reef and developed a nonhydrostatic hydrodynamic model for flows over the complex terrain. We present a simplified method to estimate the hydrodynamic roughness *z*_{0} on a complex surface, with only the rms of the depth and streamwise slope as well as a coefficient *a*_{1} required as inputs. Moving forward, the method presented here is likely applicable to other wavy bedforms with *k*–*δ*-type roughness. Further work could explore how the coefficient *a*_{1} may vary depending on other flow conditions and geometries and the transition to *k*-type roughness for high slopes. Additionally, surface waves are common in shallow flows; thus, understanding how oscillatory flows interact with the mean hydrodynamic drag over complex surfaces is a logical next step. Finally, while this method well predicts the average stress to model the large-scale flow, it likely does not well approximate the local bed shear stress, which is often used to model sediment transport and should be investigated further.

## Acknowledgments

This work was supported by a grant from the National Science Foundation (OCE-1536502; “Collaborative research: Wave driven flow through a shallow, fringing reef”) to SGM, JJA, and CBW. The data from this study are available by request from the corresponding author. Suggestions from two anonymous reviewers, and discussions with Geno Pawlak significantly improved this manuscript. We wish to acknowledge the field team: Annie Adelson, Emma Reid, Benjamin Hefner, Ron Instrella, and Will Roderick. Assistance with modeling was provided by Oliver Fringer, Yun Zhang, Eric Mayer, Kurt Nelson, and Wenhao Chen. This work was conducted under permits from the U.S. Department of the Interior National Park Service, National Park of American Samoa, and the American Samoa Department of Marine and Wildlife Resources.

### APPENDIX

#### Instrumentation, Sampling, and Initial Processing

Pressure, velocity, drifter, conductivity, and temperature measurements were made on the reef, and a weather station was located 1.8 km to the southwest (Table A1). Monitoring stations with larger instruments were attached to polyethylene plates and secured to the reef in areas of dead coral or sand. Monitoring stations with only small sensors were secured directly to dead corals or rubble.

The measured pressure *p* is filtered to obtain *η*_{hyd}, the free surface deviation due to hydrodynamic processes; see the supplemental material. The net offset *η*_{0}(*x*) at each site between *h* and MSL is determined from (3):

where , and Δ*η*_{0} is the intercept in the least squares sense for the smallest 10% of the sum of the RHS; that is, as forcing approaches zero, water slope should be flat (Monismith et al. 2013; Lentz et al. 2016). The terms in parentheses in (A1) are described in section 3b. The mean surface stress is approximated by a quadratic drag law, , where *ρ*_{a} is air density, *C*_{Da} is the wind drag coefficient, and **U**_{10} is wind velocity (Smith 1988). At the forereef sites, *η*_{0} = 0, and *η*_{0} at the first reef site onshore of the surfzone (K0) is estimated from the wave energy flux using the method of Vetter et al. 2010. Wave analysis is conducted using standard spectral methods (Dean and Dalrymple 1991); see the supplemental material.

An alternate method of computing the mean bottom stress is computed from the turbulent Reynolds stress, which is assumed constant within the inertial sublayer (Reidenbach et al. 2006), , using the measured turbulent velocities **u**′ from the ADVs and removing wave effects (Benilov and Filyushkin 1970). Combining with (1) gives

## REFERENCES

*Water Wave Mechanics for Engineers and Scientists*. Advanced Series on Ocean Engineering, Vol. 2, World Scientific, 353 pp.

*Encyclopedia of Modern Coral Reefs: Structure, Form and Process*, D. Hopley, Ed., Springer, 563–573.

**34**,

*Fluid Mechanics.*4th ed. Academic Press, 872 pp.

*Theory and Applications of Ocean Surface Waves*. Advanced Series on Ocean Engineering, Vol. 23, World Scientific, 1071 pp.

*Coastal Bottom Boundary Layers and Sediment Transport.*World Scientific, 323 pp.

*Boundary-Layer Theory.*7th ed. McGraw-Hill, 817 pp.

## Footnotes

Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JPO-D-18-0013.s1.

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).