## Abstract

Drag and turbulence in steady stratified flows over “abyssal hills” have been parameterized using linear theory and rates of energy cascade due to wave–wave interactions. Linear theory has no drag or energy loss due to large-scale bathymetry because waves with intrinsic frequency less than the Coriolis frequency are evanescent. Numerical work has tested the theory by high passing the topography and estimating the radiation and turbulence. Adding larger-scale bathymetry that would generate evanescent internal waves generates nonlinear and turbulent flow, driving a dissipation approximately twice that of the radiating waves for the topographic spectrum chosen. This drag is linear in the forcing velocity, in contrast to atmospheric parameterizations that have quadratic drag. Simulations containing both small- and large-scale bathymetry have more dissipation than just adding the large- and small-scale dissipations together, so the scales couple. The large-scale turbulence is localized, generally in the lee of large obstacles. Medium-scale regional models partially resolve the “nonpropagating” wavenumbers, leading to the question of whether they need the large-scale energy loss to be parameterized. Varying the resolution of the simulations indicates that if the ratio of gridcell height to width is less than the root-mean-square topographic slope, then the dissipation is overestimated in coarse models (by up to 25%); conversely, it can be underestimated by up to a factor of 2 if the ratio is greater. Most regional simulations are likely in the second regime and should have extra drag added to represent the large-scale bathymetry, and the deficit is at least as large as that parameterized for abyssal hills.

## 1. Introduction

Slowly varying stratified flow over topography occurs throughout the ocean either due to mean flows or eddies. Energy can be lost from the mean flow by bottom friction (usually small), by the creation of internal waves that radiate and eventually break (Nikurashin and Ferrari 2010), or by other nonlinear processes. By creating internal waves that have to break in the water column, mean flow over rough topography is one of the possible pathways by which the interior of the ocean is mixed, with long-term consequences for overturning circulations and distribution of tracers in the ocean. The same processes remove momentum from the mean flow in a way that can also affect large-scale pressure gradients and, hence, create the form stress on the large-scale flow that ultimately balances surface wind stress.

The overall goal is to be able to parameterize drag and mixing due to near-bottom stratified flows over rough topography that may not be resolved in large-scale models. Linear theory for steady stratified flow over topography (Bell 1975) allows us to calculate the rate at which energy is removed from the mean flow over a topography composed of a broad range of spectral components. Rate of energy lost implies a form drag over the topography of *F* = *u*_{0}*D*, where *F* is the radiated energy, *D* is the form drag, and *u*_{0} is the mean flow speed. Here, we will generally deal with *F*, but it is directly proportional to the drag. Tests with two-dimensional topography indicate that Bell’s theory is relevant to provide a drag parameterization for oceanic scales (Nikurashin and Ferrari 2010), though corrections need to be made if the topography varies in the cross-flow direction (Nikurashin et al. 2014), with substantially less generation of internal waves for large topography [*Nh*/*u*_{0} large, where *h* is the root-mean-squared topographic height, and *N* is the buoyancy frequency. This is called the “steepness parameter” by Nikurashin et al. (2014); we will call it the “inverse Froude number.”]. However, freely propagating internal waves are only generated for topographic wavenumbers *k* > *f*/*u*_{0}, where *f* is the Coriolis frequency. For larger-scale topography, the response is evanescent, with a vertical decay scale given by . For large wavelengths, this decay scale can reach hundreds of meters. In their work comparing to theory, Nikurashin et al. (2014) bandpassed the bathymetry so that *f*/*u*_{0} < *k* < *N*/*u*_{0}. The extra mixing and drag due to these “abyssal hills” is believed to be significant to accurate global numerical simulations (e.g., Trossman et al. 2016) and in calculations of overturning circulations (e.g., de Lavergne et al. 2016).

This leads to the central questions of this paper: How important is the large-wavenumber topography to the drag and turbulence on the near-bottom flow, and should these larger scales be included in any topographic drag parameterization? To apply Bell’s theory, the inverse Froude number *Nh*/*u*_{0} needs to be small. For a flow over bandpassed topography, this is close to being met for abyssal oceanic scales, where typical values used by Nikurashin et al. (2014) were *N* = 10^{−3} rad s^{−1}, *h* = 80 m, and *u*_{0} = 0.1 m s^{−1}, so *Nh*/*u*_{0} ≈ 0.8. However, the full bathymetric spectrum that they bandpassed from (and that we use below) has a root-mean-squared height of *h* = 305 m, so a characteristic *Nh*/*u*_{0} ≈ 3. This means that the large-scale part of the red topographic spectrum is nonlinear and not amenable to linear treatment (Bell 1975). Rather, upstream blocking and downstream hydraulic effects are predicted to be important (Baines 1995; Klymak et al. 2010), as well as an increased tendency for the flow to go around obstacles in the cross-flow direction (Nikurashin et al. 2014).

The importance of parameterizing drag and turbulence due to large-amplitude topography (*Nh*/*u*_{0} ≫ 1) has been long realized by the atmospheric science community (e.g., Bacmeister and Pierrehumbert 1988; Lott and Miller 1997). A number of schemes have been devised based on the fact that substantial portions of the flow can be blocked around an isolated obstacle, and hence the drag becomes quadratic with the flow speed and proportional to the cross-sectional area of the obstacle exposed to the mean flow (e.g., Scinocca and McFarlane 2000; Garner 2005). Recently, these parameterized drags have been shown to be important in numerical simulations of the ocean (Trossman et al. 2013, 2016) and to have promise when compared to observations (Trossman et al. 2015). These parameterizations are quadratic in flow speed and linear in obstacle height corrected by the blocking depth *h*_{b} = *u*_{0}/*N*:

Conversely, Klymak et al. (2010) found that drag for *Nh*/*u*_{0} ≫ 1 is linear in *u*_{0} and quadratic in height for a two-dimensional obstacle:

The key difference is that flow cannot go around the obstacle, so drag is set by the difference between the up- and downstream isopycnal deflections.

Here, we report a number of simulations based on those by Nikurashin et al. (2014). As described in section 3, simulations are made over three types of bathymetry from the same spectrum used by Nikurashin et al. (2014): one where the topography has been bandpassed as they present, one where the topography has been lowpassed at the same low wavenumber used in the bandpass, and a third where the full bathymetric spectrum is used (small-cap typography indicates the bathymetry in what follows). Because we need to resolve the large scales, these model runs are carried out on a very large domain. The results for simulations over these topographies at four different mean flow speeds are presented (section 4), showing that the large-scale lowpassed topography has approximately twice the dissipation of the small scale (bandpass), and all the bathymetry scales together (full) have somewhat more than the sum of the two simulations across the velocities investigated. The implications for sampling are briefly discussed (section 5) due to the localized nature of the turbulence generated by the large-scale topography. We also investigate whether the turbulence and drag from the large-scale topography will be represented in regional- and global-scale ocean models.

## 2. Form drag, dissipation, and our model setup

Our model follows Nikurashin et al. (2014) and aphysically applies a body force in the *y*-momentum equations equal to +*fu*_{0}, where *f* is the Coriolis parameter, and *u*_{0} is a “target” flow speed in the *x* direction. The idea of such a forcing is to simulate a steady surface pressure gradient in the *y* direction that would be in balance with the near-bottom mean flow:

The analog is to a channel flow in *x* with a wind stress *τ*_{w} in the *x* direction. This sets up an Ekman layer that converges on the southern boundary and diverges on the northern (assuming the Northern Hemisphere), setting up a cross-channel sea surface tilt so that ∂*η*_{0}/∂*y* < 0, and reaches steady state when the bottom Ekman layer transport (to the north) matches the surface Ekman transport (to the south). Here, we briefly review how that forcing manifests itself in the momentum and energy balances, partly because the relationship between turbulent dissipation and form drag is conceptually a bit subtle, and some communities prefer to think in terms of one or the other.

### a. Momentum budget

We omit the upper Ekman layer from our simulation, supply the surface pressure gradient via the body force, and do not have sides on our channel, so it is doubly re-entrant. The physics of this setup require a bit of justification to convince ourselves that the momentum and energy balances are the same as a channel flow. Further, the concept of form stress, or form drag, is one that is often invoked as important to large-scale circulations (e.g., Naveira Garabato et al. 2013), and it is helpful to review the concepts here.

Following Hughes and de Cuevas (2001) and Naveira Garabato et al. (2013), the vertical integral of the momentum balance can be expressed as

where all quantities have been divided by *ρ*_{0}, so the units of all quantities are m^{2} s^{−2}. The quantities **A** and **B** are vertically integrated lateral stress divergences and nonlinearities in the flow, is the depth-integrated mass flux vector, *P* is the depth integral of the pressure field (expressed as an anomaly from hydrostatic), and *p*_{b} is the bottom pressure. *H* is the water depth, and *τ*_{w} and *τ*_{b} are the surface wind stress and bottom frictional stress, respectively; note *τ*_{b} is the turbulent skin stress on the water by the topography and usually acts opposite the wind stress. On a large scale, **A** and **B** are small, and if we integrate in *x* over a closed contour in a re-entrant channel, then the net pressure gradient term must drop out of the *x*-momentum balance. Further, integrating from top to bottom, the net north–south momentum fluxes *V* must disappear by continuity, and we find that the wind stress is balanced by the bottom friction and form stress:

In most of the ocean, the form stress is thought to dominate the skin friction (e.g., Hughes and de Cuevas 2001; Vallis 2006; Naveira Garabato et al. 2013). If we divide this linear balance into two vertical layers, one representing the upper Ekman layer and the other the rest of the water column, then the balances are just the Ekman balances:

In a steady-state channel flow, the two Ekman transports cancel out by continuity. Our simulation simulates the bottom layer only and ignores the surface Ekman layer. Because we have no north–south boundaries, the bottom Ekman layer never converges and has no more effect on the momentum balance of the simulation.

In the *y*-momentum balance, in the absence of *y*-direction wind stress, we have

and upper-layer stresses never enter into the problem. Our simulation does not admit a net pressure gradient in *y* because it is re-entrant in the north–south direction (as well as east–west), so instead we impose a body force in the *y*-momentum equation as *F* = *fu*_{0} to play the part of this pressure gradient, yielding

Any deviations in the *x* velocity from *u*_{0} are balanced by friction and form stress in the *y* direction. As we note below, the *y*-momentum budget is dominated (by two orders of magnitude) by the body force and the mean flow.

Our simulations capture the relevant near-bottom dynamics, with an identical momentum budget to a more realistically wind-forced channel, with the exception of not simulating the upper Ekman layer. We specify the equivalent of a mean pressure gradient in *y*, and hence the target flow speed *u*_{0}, rather than specifying the surface wind stress. This means that for the same flow speed *u*_{0}, the bottom stress will depend on the character of the bottom friction and the form stresses, rather than the flow speed adjusting so the bottom stress matches the wind stress. However, given that we are ultimately interested in parameterizing bottom drag in terms of near-bottom flow speeds, and that there are substantial other large-scale processes that affect that near-bottom flow speed beyond the surface wind stress, this is a useful exercise.

### b. Energy budget

The energy budget can be formed by taking the dot product of the momentum equations with the velocity vector. For the full water column, the volume integral just reduces to

or the network by the wind stress integrates to the net dissipation (*D*) in the fluid. If we instead just consider the flow beneath the upper Ekman layer, then the energy budget is

In our context, this means the body force does work on the flow to produce the dissipation.

It is further useful to consider how the form stress enters the energy equation. It cannot do so in a net sense, but if we put ourselves in a reference frame where the topography appears to move with velocity −*u*_{0}, then the topography does work on the fluid both as form drag work and skin friction that has to be dissipated in the interior of the fluid:

If we consider the difference with the full energy equation, we get the transfer from the large-scale forcing to the form drag:

So essentially, the work done by the surface pressure gradient (or body force) is balanced by work done by form and skin friction stress, which then drives interior dissipation or seabed friction. In the case where there is a flow speed *u*_{0} that we want to compare to the form drag, it is most useful to think of the work done to pull the topography through a quiescent fluid at speed *u*_{0} as the right-hand side of Eq. (13), which is the sum of the forces on the topography times the flow speed *u*_{0}. The work done by the form drag on the fluid and the work done by the skin friction drag have to cascade to turbulence.

In a stratified fluid, the form stress can be greatly increased over the skin stress due to internal wave drag (e.g., Peltier and Clark 1979; Durran 1986; Bacmeister and Pierrehumbert 1988). This internal wave drag is usually conceptualized as radiating internal waves that break elsewhere in the flow and deposit momentum, though the importance of local breaking has often been remarked upon. Here, we point out that the nonradiating component of the flow can also be nonlinear and, hence, plays a role in the energy and momentum budgets.

## 3. Model configuration

Here, we use a similar model machinery to Nikurashin et al. (2014), wherein we assume a doubly periodic domain with constant stratification and a mean flow in the *x* direction forced over rough topography. The strength of the flow maintained by a body force is meant to simulate an externally imposed surface pressure gradient. The doubly periodic lateral domain of 409.6 km in the *x* direction and 118.4 km in the *y* direction is chosen to capture enough variance in the large-scale topography.

### a. Stratification, forcing, and spinup

The MITgcm was used for all simulations (Marshall et al. 1997), in a manner analogous to previous work at similar scales (Buijsman et al. 2014; Klymak et al. 2016). Background explicit vertical and horizontal viscosity and diffusivity are kept low (*K*_{ρ} = *ν* = 10^{−5} m^{2} s^{−1}), except in the presence of resolved density overturns, where the vertical viscosity and diffusions are increased in a manner consistent with the expected Thorpe scale (Klymak and Legg 2010). There is also numerical diffusivity and dissipation due to the second-order flux-limiting temperature advection scheme (tempAdvScheme = 77; see the MITgcm manual). Using a weak explicit diffusivity allows the internal wave field to evolve freely to the extent that the chosen resolution allows, rather than adding artificial damping (Shakespeare and Hogg 2017). For the work carried out here, the terms in the energy budget are all calculated, and the residual is identified as the dissipation. However, the spatial distribution of explicit dissipation (calculated from the explicit viscosities and local shears) is similar to the inverse energy budgets. The model is run in hydrostatic mode for these simulations.

The goal is to calculate drag and energy loss as a function of near-bottom *N* and *u*_{0}. To that end, the simulations are run with a constant initial stratification of *N* = 10^{−3} s^{−1}. An initial uniform velocity was set in the *x* direction (*u*_{0} = 0.02, 0.05, 0.1, 0.15 m s^{−1}), and momentum was maintained in the flow using a universal body force in the *y* direction: *F*_{B} = +*fu*_{0}. We chose *f* = +10^{−4} rad s^{−1}.

All simulations had vertical resolution of 10 m over 4000-m depth. Coarse 1-km runs over the 409.6 km × 118.4 km domain were spun up for 200 h. These runs had weak relaxation to the background stratification in a region that covered 25% of the domain in the *y* direction and the whole domain in the *x* direction. These results were interpolated onto fine-resolution 100 m × 100 m runs that had no buoyancy relaxation, and the runs were simulated for a further 20 h. There was a loss of stratification near the topography in the simulations (Fig. 1b), which is unavoidable, given the 410-km along-flow domain length (47-day transit time at 0.1 m s^{−1}) without resorting to unphysical forcing in the interior of the domain of interest. The time scales work well because it is the large-scale near-inertial internal waves that are slow to propagate in the vertical, whereas the smaller-scale waves set up quite rapidly. We diagnose the rate of change of energy in the energy budgets below, and the residual is small. We also note below that the vertical energy flux aloft is small, indicating that the model is in relative steady state.

### b. Bathymetry

The basic bathymetry used for the simulations is a stochastic version of the bathymetric spectrum used in Nikurashin et al. (2014), given by

where *H* = 305 m is the root-mean-square of the topographic height, *μ* = 3.5 is a parameter setting the high-wavenumber slope, and *k*_{0} = *l*_{0} = 1.8 × 10^{−4} rad m^{−1} are parameters that set the wavelength at which the spectrum of the topography starts to flatten out (Fig. 2; gray dashed spectrum).

Three variations on this topography are used. The full topography contains variance at all wavenumbers (Fig. 2; red line), bounded at the large scale by the domain size and at the small scales by the grid resolution. The bandpass topography (Fig. 2; blue line) is composed of wavenumbers *f*/*u*_{0} < |**k**| < *N*/*u*_{0} and the lowpass topography of wavenumbers |**k**| < *u*_{0}/*f*, which corresponds to a wavelength of 6 km.

The qualitative effect on the flow of the low-wavenumber topography is clear from an example cross section (Fig. 1a). There are peaks and valleys such that the full topography spans 1250 m of water depth. The scale of this topography strongly affects the inverse topographic Froude number *Nh*/*U*, which is approximately 0.8 for the bandpass topography but greater than 3 for the full topography. The goal of this paper is to determine the effect this large-scale topography has on the turbulence.

## 4. Results

### a. Overview of simulated flows

Example slices from the simulations illustrate the differences between the bathymetries (Figs. 3, 4). The bandpass bathymetry simulations yield a bottom-intensified (almost) steady internal wave field that has radiated energy through the domain (Fig. 3c). Directly near the seafloor, there is evidence of enhanced numerical dissipation due to the flux-limiting advection scheme (Fig. 4c) as evinced by the pixilation of the velocity field at these depths.

The flow over the full and lowpass bathymetry shows the impact of including the large-scale bathymetry (Figs. 3a,b). First, there are approximately 100-km-scale regions of barotropic acceleration and deceleration due to the inhomogeneous nature of the bottom form drag. Further, there are thick regions of acceleration and deceleration near the topography of the same order as the strength of the forcing. These regions scale in thickness roughly as Δ = *πu*_{0}/*N* ≈ 300 m for the flows here (see below, where we change *u*_{0}) and represent the thickness of the active layer of the flow near the bathymetry (e.g., Klymak et al. 2010). These regions are the same order as the bathymetric scale, so the flow is substantially nonlinear, with blocking, steering, and hydraulic responses all expected phenomena.

Also of note is the existence of radiating internal waves in lowpass solutions, despite there being no topographic variance for wavenumbers *k* > *f*/*u*_{0} (Fig. 4b). These waves have relatively high amplitudes and are due to the local flow over the bathymetry being faster than the mean flow, so that *f*/*u* > *f*/*u*_{0}; note that the clearest examples of wave packets are seen in regions of enhanced near-bottom velocities (Fig. 3b).

The flows over the large-scale topography have variance out to low wavenumbers, as confirmed by looking at lateral temperature spectra above the topography (Fig. 5). The bandpass run, as expected, has variance that is largely confined in the freely propagating regime. This energy drops, particularly for large scales, with distance from the topography (Fig. 5b). For the runs with large-scale topography, there is substantial variance at large scales (Fig. 5a; red and purple lines). By 2300-m depth, this variance also drops, and it drops selectively at the smaller scales. While there is undoubtedly some nonlinearity in this response, it is consistent with Bell’s solutions that predict an exponential decay of the linear response with height above the topography proportional to the horizontal wavenumber *k*.

### b. Laterally averaged energy budgets

We form an energy budget of the simulated flows that we integrate laterally (and then vertically) to determine the important terms and take a residual to get the dissipation in the model. We linearize the potential energy term, which given the relatively small vertical displacements is an acceptable approximation, and use a Boussinesq approximation, so that

where **u** is the velocity vector, *g* is the acceleration due to gravity, is the square of the background buoyancy frequency, is the density anomaly relative to the background density profile *ρ*_{b}(*z*), and *ρ*_{0} = 1000 kg m^{−3} is the reference density.

The background density profile *ρ*_{b}(*z*) is calculated from the modeled density field following Tseng and Ferziger (2001). Akin to the sorting procedure proposed by Winters et al. (1995), the cumulative distribution of water area is calculated as a function of density and matched to the cumulative distribution of water area as a function of depth. Interpolating from one distribution to the other gives the depth of each density in the distribution, and mapping back to depth gives a profile of density that is sorted by depth.^{1} The background density gradient changes (slightly) in the hour of simulation that we form our energy budget over, allowing us to calculate the change in background potential energy in the model (*E*_{B}) between model snapshots as

where *V* is the total volume of water, and *A*(*z*) is the area of model domain at each depth that is in the water.

The horizontally averaged energy budget is then

where *w* is the vertical velocity, is the pressure anomaly compared to a reference pressure profile , is the energy input via the body force, and *D* is the residual due to dissipation and changes in the background potential energy *E*_{b}. The first two terms on the right-hand side each integrate to zero in the vertical, but serve to redistribute energy vertically in the water column. The divergence of the nonlinear vertical advection of energy is nonzero and largely in opposition to the wave pressure work divergence. The net effect is that the dissipation is inferred to occur slightly above where the energy is put into the system by the body force.

#### 1) Energy changes with topography

The vertical integral of the energy budgets shows that the simulations are not perfectly in steady state, with the rate of change term (Fig. 6; red line) varying between 5% and 30% of the body force term (Fig. 6; purple line). We consider this an uncertainty in the energy budget—if allowed to run to full steady state, it is possible that the dissipation would increase, or that the rate of conversion (as represented by the body force) would go down. Below, we use this to assign an uncertainty to the dissipation estimate, putting the “dissipation” between the body force and the residual.

The vertical integral of the energy budget shows the clear difference between the simulations over the different topography types. The bandpass simulation has the least amount of energy loss from the mean flow of 11.7 ± 0.3 mW m^{−2} (Fig. 6c), concentrated near the seafloor in a thin layer approximately 300 m thick. This corresponds very well with the 10 mW m^{−2} of energy conversion found by Nikurashin et al. (2014) for similar parameters. However, when just the lowpass bathymetry is used, there is substantially more dissipation (21.0 ± 4.0 mW m^{−2}) than the bandpass simulation. This is because the flow is blocked or accelerated in regions due to the large-scale topography.

The full topographic simulation has more dissipation than the other two simulations combined (39.5 ± 0.8 mW m^{−2}; Fig. 6a). Dissipation extends higher into the water column than the bandpass simulation, partially because the topography extends higher, but also because turbulent features are on the order of *πU*/*N*_{0} ≈ 300 m high (Fig. 4a). Overall, the high dissipation covers about 800 m of depth. (Note that we have not presented dissipation vs height off bottom, which is not possible to do with the residual budgets we are making here.)

The background potential energy changes significantly in these simulations, with a total between 15% and 40% of the energy budget residual. These simulations are not direct numerical simulations, so the exact ratio of dissipation to irreversible buoyancy flux should not be taken very seriously. However, it does indicate that vertical mixing can be substantial in these flows, as indicated by the erosion of the near-bottom stratification (Fig. 1b).

Finally, we compute the form drag work as , which should be equal to the work done by the body force in steady state. Within the moderate errors of our calculations, this number agrees with the body-force calculations (Fig. 6). The moderate discrepancies can again be considered a measure of the uncertainties in our budgets.

#### 2) Energy changes with mean flow speed

The same model setup was used to simulate flows with mean flows of *u*_{0} = 0.02, 0.05, 0.10, and 0.15 m s^{−1}. Note that the three types of topography were not changed in these runs, so the bandpass topography, bandpassed between 600-m and 6-km scales, would not exactly match the lower and upper bounds of the permissible internal wave band set by *f*/*u*_{0} < *k* < *N*/*u*_{0}. With this in mind, it is clear that stronger forcing leads to stronger dissipations with all three topographies (Fig. 7). As with *u*_{0} = 0.10 m s^{−1}, the bandpass topography has significantly less dissipation than the lowpass and full topography, by about a factor of 2 and 3, respectively. These differences drop as higher mean flow speeds are simulated because the bandpass simulations have a steeper power law with *u*_{0} than the other two topographies.

The power law of the dissipation versus *u*_{0} in the bandpass simulations of 1.95 is very close to the theoretically expected value of 2 (Fig. 7; light blue). The power laws for the lowpass and full simulations are less than this value, with the lowpass power law being 1.7 and the full power law slightly steeper at 1.76 (Fig. 7; purple and red, respectively). For the lowpass case, a likely reason for the lower power law is that as *u*_{0} increases, the amount of blocking (and, hence, the strength of downstream hydraulic jumps) decreases as *Nh*/*u*_{0} decreases. Hence, the flow becomes more linear, and the nonlinearity that drives the nonpropagating drag decreases. The regime in these runs is *Nh*/*u*_{0} ≈ 25, 10, 5, and 3, so we are between the classical “linear” wave drag regime and the strongly nonlinear wave drag regime.

The dissipation in the full simulation exceeds the sum of the lowpass and bandpass dissipations, so the different wavenumbers in the topography obviously interact. The simplest explanation of this is that the large-scale elements in the full simulation locally accelerate or decelerate the flow over the small-scale elements. While the corresponding mean is close to *u*_{0}, the mean of , so the dissipation due to the small-scale topography is greater than in the bandpass case. We can see evidence for this in the snapshots (Fig. 3a), where there appears to be more internal wave activity emanating from regions where the flow has been accelerated versus regions where it has been decelerated (recall that the velocity anomaly is plotted, so deep blue colors are getting close to zero velocity, not negative velocity). Most remarkably, the dissipation is much closer to quadratic with target flow speed *u*_{0}, consistent with a linear drag law, in contrast to quadratic drag laws in atmospheric parameterizations of nonpropagating wave drag.

## 5. Summary and discussion

Simulations of mean flow over topography that include large scales exhibit significant energy removal from the mean flow (and hence, form drag) over and above that exerted by steady flow over the smaller scales. This is despite the fact that flow over the large scales cannot emit propagating internal waves because the topographic wavenumbers *k* < *f*/*u*_{0}. In steady state, this nonpropagating part of the wave field would not extract energy from the mean flow under linear internal wave dynamics because the disturbances would be evanescent, and other linear processes like topographic steering are similarly inviscid. However, at the large scales, *Nh*/*u*_{0} > 1, and the flow is significantly nonlinear and dissipative, with flow blocking upstream of topography and breaking motions downstream. This means that at the parameter ranges considered here, the dissipation due to the full-spectrum topography flows is more than 3 times as much as just the flow over high-wavenumber topography, a result that holds across a broad range of velocities over the same topography.

### a. Scaling energy loss from mean flow

It would be desirable to predict the dissipation due to the large-scale topography. Because of the nonlinearity and three-dimensionality, this does not appear simple to do, and certainly there is no linear theory to appeal to like that used for the high-wavenumber flows. As noted above, attempts in the atmospheric science literature appeal to bluff–body drag to get a quadratic drag law proportional to the height and width of individual topographic features (e.g., Lott and Miller 1997; Garner 2005). However, our simulations are not consistent with that scaling, yielding a linear drag law instead.

A rough estimate for large *Nh*/*u*_{0} flows can be derived from the isolated topography case (Klymak et al. 2010), where the form drag integrated along an obstacle can be approximated by

Here, a reasonable value for an obstacle height is *h*_{m} = 350–500 m. Turning this into an energy density requires assuming a spacing between obstacle peaks, which is approximately Δ*x* ≈ 100 km (Fig. 1a). This rough calculation yields a dissipation of 17–33 mW m^{−2} for *u*_{0} = 0.1 m s^{−1}, which brackets the 23 mW m^{−2} simulated for the lowpass topography simulation. The power law of this parameterization scales as for large *Nh*/*u*_{0}, but as *Nh*/*u*_{0} gets smaller, the correction terms start to dominate, flattening the power law, as observed as *u*_{0} was increased (Fig. 7). Of course, the assumptions that go into Eq. (18) are not valid at lower *Nh*/*u*_{0}, so the analogy breaks down. Future effort should be aimed at parameterizing the large-scale drag and/or dissipation using a priori calculations.

Our scenario is a closed system, and while flow can locally go “around” the topography, it cannot do so indefinitely because there is always more downstream topography. So, in this sense, it is more akin to periodic bathymetry than the isolated bathymetry used to derive the drag laws used in the atmospheric literature. Which situation is more appropriate in the ocean or atmosphere probably depends on the exact location. However, we are not aware of simulations in either environment that actually test the drag laws on three-dimensional topography, as we have attempted to do in this paper.

### b. Lateral inhomogeneity

The reason we had to use such a large modeling domain is that the dissipation caused by the large-scale topography is spatially inhomogeneous (Fig. 8). Sensitivity tests on the 1-km coarse models indicated that dissipation results converged as the domain approached 100 km × 400 km, so that was what was used for the simulations above. The reason for that can be readily seen (Fig. 8a), where the regions of strong dissipation are on the scale of 50–100 km in the along-flow direction and about 25 km in the cross-flow direction. Subsets of these regions have a strong potential to be biased.

The distribution of the dissipation has important implications for oceanographic sampling. If an experiment were deployed similar to DIMES (e.g., St. Laurent et al. 2012) where 34 vertical profiles could be accomplished, then the mean dissipation rate could be substantially biased, as shown by a Monte Carlo sampling of the sample means (Fig. 8b; blue line). The median of this distribution is 0.6 times the mean of the dissipation in the simulation, so the expected bias is relatively large and biased low, but there is a significant fraction of the Monte Carlo means that are many times the actual mean. The problem is worse if the sampling is deterministic or biased toward sampling hot spots, and of course a few individual moorings could be placed anywhere in this region and get answers that are either far too low or far too high (e.g., Waterman et al. 2013; Brearley et al. 2013).

### c. Resolution dependence: Do coarser models have the dissipation?

The argument for including a parameterization for drag and mixing due to abyssal hills is that numerical models cannot simulate these scales, and hence, the extra drag and mixing should be added. In this paper, we argue that the drag and dissipation due to the large-scale (stochastic) topography is significantly larger. However, it is entirely possible that coarse models already simulate the turbulence and drag due to this large-scale part of the topographic spectrum because they are (partially) resolving them.

We ran the same simulations over a range of lateral resolutions (1.0, 2.5, and 4.0 km) and vertical resolutions (between 10 and 307 m). We used the lowpass bathymetry, so there is no topographic roughness at scales smaller than 6 km. We used the same nominal forcing noted above with flow speeds of *u*_{0} = 0.1 m s^{−1}. Lateral and vertical model resolution strongly affect the resulting dissipation, and interestingly, the two effects counter each other (Fig. 9). For a very high lateral resolution by global model standards (1 km), the agreement with the 100-m resolution runs drops as vertical resolution is coarsened, with the energy loss plateauing relatively quickly at about 0.6 of the loss in the high-resolution run. To compare with the energy budgets above, this is 15 kW m^{−1}. If the model has a well-parameterized abyssal hill dissipation (i.e., bandpass) of 12 kW m^{−1}, then the total dissipation inferred in this coarse model is 27 kW m^{−1}, compared to 40 kW m^{−1}, or 67%. This is not terrible disagreement, and the use of the abyssal hill parameterization will definitely help get the correct dissipation and drag.

For coarser simulations, Δ*x* = 4 km, the energy loss is exaggerated at fine vertical resolutions (Fig. 9; blue line) and then suddenly drops to much less than the fine-resolution runs for Δ*z* > 200 m. This occurs as the root-mean-square of the topographic slope (vertical colored lines in both plots). Therefore, the guidance for numerical modelers running at regional-scale resolutions is somewhat ambiguous. If models resolve the root-mean-square of the topographic slope (approximately 2*π*/*k*_{0}, where *k*_{0} is the topographic bandwidth parameter described above), then they will slightly overpredict mean flow stratified drag. If the vertical resolution is too coarse, then they will not resolve this drag, and they will need to increase their drag significantly more than that just due to the parameterizations for abyssal hills.

### d. Concluding remarks

In the linear limit, large-wavenumber topography should not generate radiating internal waves and, hence, should not have drag or play a part in the energy budget. However, such topography is actually the dominant term in the energy budget because of nonlinearity and turbulence in an ocean-relevant regime. Coarse-scale models are challenged to get this large-scale nonpropagating drag right, though exactly whether they over- or underpredict the drag depends on the ratio of vertical to horizontal resolution. The dissipation due to the large-scale bathymetry is very spatially inhomogeneous, leading to a strong observational challenge. Finally, the extra drag scales linearly with flow speed, in contrast to being assumed quadratic in atmospheric parameterizations and ocean sensitivity tests (Trossman et al. 2013, 2015). Further studies are planned to try to predict the low-wavenumber drag over the stochastic topography better and to incorporate it into numerical models so that drag is not double counted.

The importance of getting this parameterization correct is an open question. Models that are freely running (with no data assimilation) must have bottom stresses that are adequate to balance the externally imposed surface wind stress. If the bottom drag is conceptualized as quadratic with the near-bottom flows, then having too small a drag coefficient means having near-bottom flows that are too fast. This then has implications for the dissipation of energy in the flow and, hence, mixing. Such feedbacks are found in Trossman et al. (2013, 2016), who found that adding the wave drag parameterization to the radiating part of the spectrum indeed dropped the near-bottom flow speeds. Large-scale form stress will still be similar to maps provided by Naveira Garabato et al. (2013), but the activating mean flows will be even more reduced.

## Acknowledgments

Thanks are given to Raffaele Ferrari, Brian Arbic, Ken Hughes, John Scinnoca, Norm McFarlane, and Kurt Polzin for helpful discussions on this manuscript, particularly Kurt Polzin for introducing me to the term “nonpropagating wave drag.” The input from two anonymous reviewers also greatly improved this effort. This work was funded under the Office of Naval Research Flow Encountering Abrupt Topography Defense Research Initiative (Grant N00014-15-1-2585), and NSERC Discovery Grant 327920-2006. Thanks are also given to Odessa Murray, who manages the High-Performance Computing accounts for ONR. Example model setup and processing scripts are available at http://web.uvic.ca/~jklymak/leewaves17.

## REFERENCES

*Topographic Effects in Stratified Flows*. Cambridge University Press, 498 pp.

*Atmospheric and Oceanic Fluid Dynamics: Fundamentals and Large-Scale Circulation*. Cambridge University Press, 964 pp.

## Footnotes

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

^{1}

The computational advantage is that a histogram is far easier to form than sorting and redistributing the fluid, particularly in a convoluted geometry, though the grid cells all must be the same size.