Ocean lee waves occur on length scales that are smaller than the grid scale of global circulation models (GCMs). Therefore, such models must parameterize the drag associated with launching lee waves. This paper compares the lee wave drag predicted by existing parameterizations with the drag measured in high-resolution nonhydrostatic numerical simulations of a lee wave over periodic sinusoidal bathymetry. The simulations afford a time-varying glimpse at the nonlinear and nonhydrostatic oceanic lee wave spinup process and identify a characteristic time scale to reach steady state. The maximum instantaneous lee wave drag observed during the spinup period is found to be well predicted by linear lee wave theory for all hill heights. In steady state, the simulations demonstrate the applicability of parameterizing the drag based on applying linear theory to the lowest overtopping streamline of the flow (LOTS), as is currently employed in GCMs. However, because existing parameterizations are based only on the height of the LOTS, they implicitly assume hydrostatic flow. For hills tall enough to trap water in their valleys, the simulations identify a set of nonhydrostatic processes that can result in a reduction of the lee wave drag from that given by hydrostatic parameterizations. The simulations suggest implementing a time-dependent nonhydrostatic version of the LOTS-based parameterization of lee wave drag and demonstrate the remarkable applicability of linear lee wave theory to oceanic lee waves.
In the ocean, because lee waves occur on length scales smaller than the resolution of global circulation models (GCMs), the lee wave drag must be parameterized. Existing parameterizations are based on steady-state solutions for the flow above “linear-height” hills, wherein the height of the hill is much smaller than the wavelength of the wave (Bell 1975b; Gill 1982). Most common among the unresolved bathymetric features of the deep ocean are the abyssal hills which are statistically homogeneous on a regional scale and thus submit well to a spectral model (Goff and Jordan 1988). This allows for calculation of the lee wave drag with the spectral linear theory of Bell (1975a,b). However, because a significant portion of the abyssal hills have nonlinear heights (Nikurashin et al. 2014), spectral linear theory is not formally valid in these regions.
At the bottom of the ocean, where the background horizontal velocity U and buoyancy frequency N are nearly constant with height, lee wave dynamics can be described in terms of the lee wave Froude number J = Nh0/U and the nonhydrostatic parameter ϵ = Uk/N, where h0 and k are the hill height and wavenumber (Mayer and Fringer 2017). Using the lee wave Froude number J, nonlinear lee waves can be separated into two regimes: subcritical, J = Nh0/U < Jc, and supercritical, J = Nh0/U ≥ Jc, where Jc = O(1) is the critical Froude number. The supercritical regime is characterized by dramatic deviations of the form drag from the predictions of linear theory. Above isolated bathymetry, the flow on the downslope side of the isolated hill displays a standing hydraulic jump-like feature called a “downslope windstorm,” with an associated pressure anomaly resulting in as much as an eightfold amplification of the steady-state drag above the linear prediction (Peltier and Clark 1979; Pierrehumbert 1987). However, for supercritical lee waves above periodic bathymetry, such as the abyssal hills, the upslope of the next hill inhibits the formation of the downslope windstorm, and the flow instead displays a process called “blocking,” wherein the lowest parcels of water become trapped in the valleys (e.g., Welch et al. 2001; Nikurashin and Ferrari 2010; Winters 2016). In the presence of blocking, some water passes over the hills and generates lee wave energy. In steady state, this overtopping fluid travels along a streamline that is different from the bathymetry. The resulting wave is therefore identical to a subcritical lee wave generated by the lowest overtopping streamline, which we will refer to as the LOTS. High-resolution nonlinear numerical simulations of atmospheric lee waves demonstrate that the supercritical lee wave drag scales with the trough-to-crest height of the LOTS, which is observed to be O(U/N) (Stein 1992; Welch et al. 2001; Eckermann et al. 2010). Thus, atmospheric wave drag parameterizations posit that the wave drag over supercritical height hills is given by substituting an effective Froude number Jeff ≈ 1 for J in linear theory (Palmer et al. 1986; Pierrehumbert 1987; Lott and Miller 1997; Garner 2005). We will refer to this approach as saturation theory.
As a package, the spectral model of the abyssal hills combines with saturation theory to permit computationally cheap calculations of the lee wave drag over arbitrary height, subgrid-scale abyssal hills in ocean models. This is the approach employed in recent studies of GCMs (Nikurashin and Ferrari 2011; Scott et al. 2011; Naveira Garabato et al. 2013; Melet et al. 2015; Trossman et al. 2013, 2015, 2016; Wright et al. 2014; Yang et al. 2018, 2019). These studies conclude that globally, lee wave drag constitutes between 0.2 and 0.75 TW of work on the ocean, meaning lee waves could be of principle importance for balancing the O(1) TW of wind work at the ocean surface (Ferrari and Wunsch 2009). However, direct observations of lee waves in supercritical regions such as the Drake Passage (Cusack et al. 2017) and Kerguelen Plateau (Waterman et al. 2013) display dissipation rates one order of magnitude smaller than predicted by combining saturation theory with an energy-dissipation parameterization (St. Laurent et al. 2002). To account for this discrepancy, recent papers have implicated the interaction of geostrophic and tidal currents (Shakespeare 2020) and the absorption of wave energy aloft via conservation of wave action in vertically sheared flow (Kunze and Lien 2019). In their introduction, Kunze and Lien (2019) offer a handful of additional hypotheses, including measurement and instrument error, under sampling, wave–wave interaction, a narrow radiating bandwidth, and, finally, that lee wave saturation is incompletely parameterized.
Saturation theory, as presented above, is based on lee waves of atmospheric scale, where the forcing is hydrostatic. Assuming hydrostatic forcing implies that the wave drag is only a function of hill height, or nondimensionally, of J, and not also a function of hill length (nondimensionally ϵ) (Pierrehumbert 1987). However, in the J–ϵ regime appropriate to lee waves over abyssal hills in the ocean, nonhydrostatic effects are important, even for subcritical height bathymetry (Mayer and Fringer 2019, manuscript submitted to Ocean Modell.). Nikurashin and Ferrari (2010) and Nikurashin et al. (2014) are an exception to the existing lee wave literature in that their simulations employ high-resolution nonhydrostatic models and capture a J–ϵ regime appropriate for abyssal hill lee waves. However, both studies derive their prediction for Jeff without varying the horizontal component of the generating bathymetry. Thus, although their simulations incorporate nonhydrostatic processes, using their result for Jeff to parameterize the lee wave drag above arbitrary bathymetry is still implicitly hydrostatic because it assumes that the drag is independent of hill length.
Additionally, existing saturation based parameterizations assume both that the lee wave has reached a steady state and that the momentum radiating through the wave is deposited within the first wavelength above the bathymetry. However, for a linear lee wave in regions with negligible rotational effects, the drag is an intrinsically unsteady feature that propagates vertically with the lee wave front, such that the drag only acts on the bottom currents during the earliest period of lee wave generation (Pedlosky 2003). This vertically propagating lee wave drag should also exist in the supercritical regime because supercritical hills still launch subcritical lee waves. Although rotation complicates this by permitting lee wave–inertial oscillation interaction (Nikurashin and Ferrari 2010), the assumption of steady-state lee waves depositing their momentum within the first wavelength above the bathymetry is likely not justified everywhere in the ocean.
This paper offers a process study of an idealized lee wave to illuminate the consequence of ignoring nonhydrostatic and unsteady effects when parameterizing oceanic lee wave drag. The results are based on high-resolution, nonhydrostatic simulations of lee waves above sinusoidal hills spanning the J–ϵ range of oceanic flow over abyssal hills. The paper is organized as follows: in section 2 we review the linear theory. In section 3 we describe the numerical model and discuss methods for measuring the lee wave drag. In section 4 we analyze the temporal evolution of a subset of subcritical and supercritical simulations. In section 5 we evaluate the steady-state drag and LOTS measured in all of the simulations. And in section 6 we conclude with recommendations for nonhydrostatic and time-dependent corrections to saturation theory.
2. Linear theory and parameter space
As in Mayer and Fringer (2017), we assume a nonrotating ocean with a free-slip bottom boundary condition and an infinite depth. Although rotation is important at the scale of abyssal hills [O(1–10) km] in much of the ocean, where its primary contribution is to enforce a lower bound on the radiating bandwidth of lee waves according to f < Uk < N, the following theory and simulations are greatly simplified by omitting it. We discuss the ramifications of ignoring rotation in the conclusions. With the above assumptions, the lee wave is uniquely characterized by J = Nh0/U and ϵ = Uk/N. Mayer and Fringer (2017) demonstrate that J = Nh0/U is the lee wave Froude number formed by the ratio of perturbation advection within the lee wave to the group velocity of the lee wave. The second dimensionless number, ϵ = Uk/N, can be cast as a ratio of the wavelength of the lee wave in the direction of the wave vector, λlee = 2πU/N, to the horizontal wavelength of the bathymetry, Lhill = 2π/k. In terms of ϵ, the vertical wavenumber of the lee wave is given by
With increasing ϵ (increasing nonhydrostasy), the lee wave wavenumber vector points increasingly downstream. When ϵ > 1, the wavenumber of the hill is larger than the lee wave wavenumber, implying an imaginary vertical wavenumber and thus a decaying, or evanescent, lee wave.
Abyssal hills are periodic and quasi-anisotropic in nature, with horizontal aspect ratios of between 1:3 and 1:8 (Goff and Arbic 2010; Fig. 1), such that they appear as ridge-like corrugations on the ocean floor. As an idealization, we consider the simple sinusoidal bathymetry
The factor of 1/2 is used so that h0 is the trough to crest height of the hill. This convention follows Lilly and Klemp (1979) and Welch et al. (2001), but differs from Nikurashin and Ferrari (2010), where the authors omit the factor of 1/2. Hence, the height scale in Nikurashin and Ferrari (2010) relates to h0 as hNF = h0/2. Using Eq. (2), Mayer and Fringer (2019, manuscript submitted to Ocean Modell.) show that the dimensionless steady-state nonhydrostatic linear form drag per wavelength is given by
Equation (3) indicates that the nondimensional drag on the background current associated with generating lee waves grows in proportion to J2, but decreases with decreasing hill length (increasing ϵ) until the evanescent boundary Lhill = λlee, beyond which the steady-state drag vanishes.
and represents the largest magnitude drag that can result from a lee wave over a hill with crest to trough height h0. From this expression, a scale for the computed drag in what follows is given by
Jc = O(1) is an empirically determined critical Froude number above which the drag saturates, and Jeff is the Froude number of the resulting effective bathymetry. It is generally assumed that Jeff = Jc. As a canonical example, in their nonhydrostatic simulations of lee waves above ϵ = 0.31 sinusoidal hills, Nikurashin and Ferrari (2010) find that saturation occurs when NhNF/U = 0.7. Recalling that hNF = h0/2, this implies a critical Froude number JNF = 1.4.
All existing parameterizations assume a steady process in which the lee wave extracts momentum continuously from the currents within one lee wave wavelength (≈500 m) of the bottom. However, irrotational linear theory (e.g., Pedlosky 2003) suggests that the lee wave drag only acts on the background current in the vicinity of the wave front, which propagates with vertical group velocity
This picture is complicated by rotation, where the more hydrostatic hills, for which Uk → f, result in an infinite vertical wavenumber, a vanishing vertical group velocity, and thus an evanescent response (Klymak 2018; Kunze and Lien 2019). Nevertheless, the general observation that the lee wave drag is felt only within the vicinity of the vertically propagating wave front remains valid. Presumably, the same temporal effect is present in nonlinear lee waves because, as posited by saturation theory, there is still a subcritical lee wave generated by the LOTS, and thus still a wave front where the lee wave drag acts.
The theory presented in this section suggests that lee waves above nonlinear and nonhydrostatic bathymetry should exhibit nonlinear, nonhydrostatic, and time-dependent behavior. However, at present, saturation theory offers a hydrostatic and steady-state parameterization for the lee wave drag. To quantify significance of omitting nonhydrostatic and time-dependent effects, we turn to numerical simulations of the idealized oceanic lee wave.
3. Numerical model
a. Model setup
Our simulations employ the nonhydrostatic Navier Stokes solver SUNTANS (Fringer et al. 2006). The setup, sketched in Fig. 1, is identical to that in Mayer and Fringer (2019, manuscript submitted to Ocean Modell.), with a uniform buoyancy frequency N = 0.002 s−1 and constant volume-averaged horizontal velocity U = 0.2 m s−1 such that the lee wave wavelength is λlee = 2πU/N = 628 m. Because we are primarily interested in lee waves in the deep ocean, where the free surface plays a negligible role, the domain in all simulations is 7 km deep and employs a sponge layer throughout the upper 5 km. The vertical grid spacing is a constant Δz = 5 m beneath the sponge layer (the bottom 2 km of the domain), and stretches linearly within the sponge layer until Δz = 300 m at the surface. This is identical to the vertical grid used in Nikurashin and Ferrari (2010). The simulations begin at rest and are spun up with a forcing scheme that ensures a constant volume-averaged streamwise velocity U (Nelson and Fringer 2017). As described in Mayer and Fringer (2019, manuscript submitted to Ocean Modell.), the kinematic viscosity is a constant ν = 0.01 m2 s−1, and no explicit turbulence model is employed. Shakespeare and Hogg (2017) suggest that this value for ν is small enough to produce an effectively inviscid linear lee wave simulation (see their Fig. 3.c). Tests of our model setup with smaller viscosity (not shown) display little dependence on the Reynolds number, defined by Re = U/(νk).
Whereas in Mayer and Fringer (2019, manuscript submitted to Ocean Modell.), we only varied the length of the sinusoidal hill [Eq. (2)], in this paper we now vary both the height and length over the set h0 = [2 to 200 m] and Lhill = [4 km to 800 m]. This samples the parameter space J = [0.02 to 2] and ϵ = [0.16 to 0.8], which is characteristic of lee waves over abyssal hills (Goff and Arbic 2010; Nikurashin and Ferrari 2011). Anticipating that nonlinear-height lee waves might take longer to reach steady state, each simulation in this paper is run for 20 excitation periods, Tex = Lhill/U. Additionally, all simulations are nonhydrostatic with a constant horizontal resolution of Δx = 10 m. This avoids both the unphysical overprediction of the wave drag caused by using a hydrostatic model with nonhydrostatic-length bathymetry as well as the underprediction in wave drag caused by poorly resolving the bathymetry (Mayer and Fringer 2019, manuscript submitted to Ocean Modell.).
b. Computing the lee wave drag
This study employs two methods to compute the simulated lee wave drag. The first uses the simulated pressure p at the bottom [z = h(x)] and the analytical bathymetry h(x) to directly compute the form drag as
A second measurement of the lee wave drag is provided by integrating the vertical momentum flux through a horizontal plane above the bathymetry. Unless otherwise noted in what follows, we choose this plane to be a height zt = 15 m above the crest of the hill, and define
For a steady inviscid lee wave over a linear hill, Fform = −Fflux analytically. Viscous effects are negligible in the ocean, and still small for zt ≪ λlee even with the relatively large viscosity in our model (Shakespeare and Hogg 2017).
The relationship between the form drag and momentum flux is demonstrated by the governing area-integrated horizontal momentum equation, which is given by
is the area-integrated momentum per unit span. In equilibrium dM/dt = 0, implying Fflux = −Fform. However, the time-integrated momentum equation produces nonzero M when unsteadiness and nonlinearities lead to imbalances between Fflux and −Fform, as discussed in what follows.
4. Unsteady behavior of lee waves
a. Subcritical, weakly nonlinear case [J = O(0.1)]
Before interrogating the supercritical simulations, in this section we consider a simulation over a weakly nonlinear but still subcritical height hill. For this purpose, the 4-km-long and 60-m-tall hill (ϵ = 0.16, J = 0.6) is a good example. Snapshots of the nondimensional vorticity [ω* = ω/(JN)] and streamlines every half Tex are shown in Fig. 2. Note that the wave front travels upward with group velocity, [Eq. (8)]. For this hydrostatic hill, the nonhydrostatic effect is negligible since . Hence, Cz ≈ ϵU = λlee/Tex such that after one excitation period, the wave front has traveled approximately one wavelength into the water column. After two excitation periods, the wave is fully established within the domain shown, and little changes over the following half Tex.
The lowest parcel of water in this weakly nonlinear simulation struggles to flow over the hill during spinup, such that at t = Tex/2, in Fig. 2, there is a small jump in vorticity halfway up the hill. By the next snapshot, at t = Tex, this bottom water has traveled over the hill and no longer presents an obstacle. Nevertheless, the sharp gradients in vorticity lead to mixing and a reduction of the buoyancy frequency near the bottom. Figure 3 shows the horizontally averaged instantaneous buoyancy frequency ⟨N⟩ for this simulation at times corresponding to those in Fig. 2, as well as later times, where we define ⟨N⟩ as
In Fig. 3, the decay of the buoyancy frequency near the bottom is significant. By t = Tex, ⟨N⟩ at the bottom of the valley is roughly 25% smaller than the background N, and it continues to decay throughout the simulation. However, the decay remains confined to a region δ = U/N above the hill (indicated by the dashed horizontal line in Fig. 3). A similar bottom-confined decay of buoyancy frequency was observed in the lee wave simulations of Klymak (2018), where it was likewise attributed to the spinup of the system from rest. There is reason to believe that the phenomenon exists in nature. For example, observations in locations ripe for lee waves such as the Antarctic Circumpolar Current north of the Kerguelen Plateau (Waterman et al. 2013) and the Hoyt Hills region of the Gulf current (Zheng et al. 2012) indicate bottom layers with thickness approximately equal to U/N. Note that this decay of near-bottom buoyancy frequency is only an issue for periodic bathymetry. For an isolated hill, the mixing during spinup would momentarily reduce the stratification above the hill, but not that of upstream fluid that flows over the hill at later times. Hence the steady-state dynamics of a lee wave above a given hill are only sensitive to the local decay of buoyancy frequency during spinup when there are other hills upstream to do the initial mixing.
To approximate the effect of a decayed bottom buoyancy frequency on the lee wave above periodic bathymetry, we define an adjusted bottom-layer buoyancy frequency by averaging ⟨N⟩ from z = 0 to z = h0 + U/N with
where z = 0 is taken at the bottom of the valley. Inspection of Eq. (4) shows that the hydrostatic wave drag is proportional to the buoyancy frequency, Flin ∝ N. Thus a reduction in the buoyancy frequency over a layer δ above the bathymetry should result in smaller form drag, which we can approximate by substitution of Nadj for N in Eq. (3), giving
Note that this expression ignores the nonhydrostatic effect of Nadj on the vertical wavenumber [Eq. (1)]. We found that incorporating an adjusted vertical wavenumber madj produced unphysical predictions, especially above the more nonhydrostatic hills. Above these narrow hills, Nadj implies imaginary madj and evanescent waves (since Uk > Nadj), even though the disturbance generates a propagating wave in the region above the decayed buoyancy. In Fig. 4, we compare Fadj to the computed form drag [Eq. (9)] and vertical momentum flux [Eq. (10)] in the subcritical simulation shown in Fig. 2. The adjusted linear theory accurately predicts the general evolution of the computed drag over the course of the simulation. As demonstrated in the following section, the adjusted linear theory works equally well during steady state (see Figs. 10, 11). Hence, for subcritical simulations, we should expect a smaller drag than predicted by linear theory, Flin [Eq. (3)], due to the decay of buoyancy in the lowest layers. Additionally, in regions where rotation is important, the decay in buoyancy frequency narrows the N − f frequency band for propagating lee waves. For a field of abyssal hills with varying wavenumber, the decay of near-bottom buoyancy frequency thus has the consequence of further reducing the total lee wave drag (Kunze and Lien 2019).
The subcritical example also affords an opportunity to study the progression of the wave front as it propagates away from the hill. In Fig. 5, we display Fflux [Eq. (10)] as a function of height above the bottom at different points in time. Note that the wave front has a vertical extent of roughly one lee wave wavelength, λfront ≈ λlee, and propagates vertically at a rate of approximately one λlee for each Tex, as predicted by hydrostatic linear theory. Also note that after the wave front has reached the sponge layer, the momentum flux is nearly constant with height, with a slight decay with distance from the bottom as predicted by viscous effects (Shakespeare and Hogg 2017). The nondivergent momentum flux confirms the hypothesis of Pedlosky (2003) that in the subcritical regime after spinup, the lee wave drag is not felt by the lowest water, where the momentum flux is nondivergent, but aloft in the vicinity of the propagating wave front. From the perspective of a GCM, this supports the notion that lee wave drag parameterizations could account for unsteadiness by appropriately positioning the drag in the vertical rather than always assigning it to the bottom layer.
b. Supercritical cases (J > 1)
In this section we report on the temporal evolution of three supercritical simulations with hills of constant height, h0 = 160 m (J = 1.6), but different length, Lhill = [4, 2, 1] km (ϵ = [0.16, 0.32, 0.63]). This subset illuminates both similarities across supercritical lee waves with different values of ϵ and important quantitative discrepancies in the form drag, especially for the narrowest hill (ϵ = 0.63).
The temporal evolution of vorticity and streamlines over the first 2.5Tex for the three supercritical hills (Fig. 6) demonstrates a remarkable qualitative similarity between lee waves over supercritical hills of dramatically different widths. In all three simulations, the lowest parcels of water become blocked and induce separation on the downslope face within the first half excitation period. Within another excitation period, the shear between the overtopping flow and the blocked layer evolves into a train of vortices that populate the valley. Thereafter, these qualitative features of the blocked layer change little. Thus, the excitation period Tex again offers a reasonable time scale for the lee wave to reach steady state.
More quantitatively, Fig. 7 shows the form drag, Fform [Eq. (9)], and the vertical momentum flux through a horizontal plane 15 m above the hill, Fflux [Eq. (10)], for 0 ≤ t ≤ 6Tex. All three simulations exhibit large discrepancies between the two measures of drag during t < Tex, indicating significant momentum deposition during this period. After the first Tex, however, the simulations display nearly balanced form drags and momentum fluxes, implying that, just as in the subcritical example, the momentum within the control volume is only altered during the first Tex of the simulation [Eq. (11)]. For the two longer hills (ϵ = [0.16, 0.31]), Fform and Fflux stabilize at around 0.4Fsat [Eq. (5)], in accordance with saturation theory. However, the narrowest hill (ϵ = 0.63) stabilizes at a much smaller magnitude. Above such a narrow hill, the blocked layer quickly fills enough of the valley to render the entire overtopping streamline evanescent, a point to which we return below.
The large discrepancies between the maximum Fform and Fflux during the first Tex is unique to the supercritical simulations and results in a significant change to the momentum of the fluid within the control volume [Eq. (12)]. Comparison with the flow for t < 2Tex displayed in the top panels of Fig. 6 encourages the interpretation that the local momentum deposition serves to arrest the valley water and establish the blocked layer, as expected for supercritical periodic bathymetry (Welch et al. 2001). We can estimate the horizontal momentum change required to block the valley water as a product of the volume Vb contained within the blocked layer and the background velocity U, namely,
We compute this value at the time Tb, which we define as the first instance that Fflux = Fform, and compare it to the change in momentum within the control volume during the time period, computed by integrating Eq. (12) from t = 0 to t = Tb. For these three hills (ϵ = [0.63, 0.32, 0.16]), the comparison gives the ratios . This close equivalence confirms that the imbalance in the form drag and the vertical momentum flux during the spinup of these supercritical lee waves indeed serves to establish the blocked layer.
The observation that for the two longer hills (ϵ = [0.16, 0.31]), the maximum value of the form drag is well approximated by the prediction of linear theory, Flin [Eq. (3)], while the momentum flux instead peaks around the saturation value Fsat = ρ0U3N−1 [Eq. (5)], as predicted by saturation theory, suggests that a scale for the local momentum deposition is given by . To test this observation, Fig. 8 displays contour plots of the maximum computed vertical momentum flux and form drag for all simulations in this study, as well as the linear prediction [Eq. (3)]. The close correspondence between the maximum computed form drag and the prediction from linear theory is remarkably robust across most of the J–ϵ regime spanned by the simulations, especially for ϵ ≤ 0.5. Likewise, the vertical momentum flux displays ubiquitous saturation of O(Fsat). Conceptually, this supports the general conclusion that during the first excitation period the form drag reaches the linear value, but the background conditions imposed by U and N restrict the momentum flux to the saturation value Fsat. The extra form drag thus acts to decelerate the valley water and establish a blocked layer. However, because the momentum sink is restricted to the first excitation period of lee wave generation, it is dramatically different from the parameterizations for the nonpropagating drag found in Pierrehumbert (1987) and Garner (2005), for which the local momentum sink is assumed to persist into steady state.
Finally, in Fig. 9 we show the evolution of the LOTS for the simulations in Fig. 6 for time t ≤ 16Tex. We note that the separated portion of the LOTS in each simulation displays persistent undulations with length scales of approximately λlee. Furthermore, the LOTS evolves very slowly for t > 6Tex, indicating that the simulations indeed approach a quasi-steady state. Thus, in what follows, quantities with the overbar (e.g., ) are time averaged during 6Tex ≤ t ≤ 12Tex.
The time-variable LOTS in Fig. 9 can be used to explain the dramatic reduction in drag observed in the very nonhydrostatic ϵ = 0.63 hill after the establishment of the blocked layer. In all supercritical simulations, once the blocked layer is formed, the LOTS has two different sections, one that is defined by the unblocked portion of the bathymetry, and another that is separated from the bathymetry. For the simulation above the narrow ϵ = 0.63 hill, the separated section of the LOTS comprises more than half of the total LOTS. As discussed in more detail in the next section, the separated portion of the LOTS is evanescent and does not generate lee waves. However, for this narrow hill, the hill length is less than two lee wave wavelengths, as indicated by ϵ = λlee/Lhill > 0.5. Therefore the unseparated portion of the LOTS, which comprises less than half of the total LOTS, is also evanescent to the flow, and thus exerts no drag after the establishment of the blocked layer. This process, which we will term “evanescent masking,” has dramatic effects on the steady-state drag for all hills in which the unseparated portion of the LOTS is evanescent, as discussed in the next section.
5. Nonlinear, nonhydrostatic effects on time-averaged drag
In this section, we analyze the steady-state lee wave drag in our simulations , obtained by averaging Fflux [Eq. (10)] over the period 6Tex ≤ t ≤ 12Tex
and compare the results with linear and saturation theory. Unless otherwise noted, we use the vertical momentum flux Fflux computed just above the hill [Eq. (10)] rather than the direct computation of the form drag Fform [Eq. (9)]. The time-averaged drag as a function of J is shown in Fig. 10. In the top panel, the drag is nondimensionalized by Fsat = ρ0U3N−1 [Eq. (5)]. This makes apparent that the drag grows with J up to a maximum of roughly 0.5 Fsat when J = 1. We identify blocking as occurring when the maximum vertical excursion of the time-averaged LOTS is less than 95% of the hill height. With this definition, blocking begins at J = 1 for the narrowest hills, and is present in all simulations for which J > 1, confirming the critical Froude number Jc = 1. In the middle panel of Fig. 10, the drags are instead nondimensionalized by the predictions from linear theory Flin [Eq. (3)]. The effect of blocking in the supercritical simulations is again apparent. However, with this scaling it is evident that the time-averaged drag is smaller than linear theory predicts, even for J < Jc. The discrepancy in the unblocked cases is a consequence of the reduction of the stratification near the bottom during lee wave spinup, as discussed in section 4a. Indeed, normalizing the drag by Fadj [Eq. (15)] instead of Flin, as in the bottom panel of Fig. 10, removes this trend. It also introduces a new source of noise for the most nonhydrostatic hills, possibly due to the effect of the reduced buoyancy frequency on the vertical wavenumber, which we neglect in computing Fadj. Together, the three panels demonstrate the separation of the lee wave into three regimes based on hill height: the linear limit J ≪ 1, where Flin [Eq. (3)] is accurate, the subcritical regime J = O(0.1), where weak nonlinearities during spinup result in a decay of buoyancy near the bottom and the drag is better predicted by the buoyancy adjusted linear theory Fadj [Eq. (15)], and the supercritical regime J ≥ Jc, in which the lowest parcels of water are trapped in the valleys and the drag is limited by the radiative capacity of the fluid, Fsat = ρ0U3N−1 [Eq. (5)].
To demonstrate the dependence of drag on the hill length, Fig. 11 displays the time-averaged drag as a function of ϵ. In general, there is a monotonic decrease in the drag with increasing ϵ, as predicted by linear theory [Eq. (3)]. At moderate values of ϵ = [0.57, 0.48, 0.42, 0.35], however, the critical J = 1 and supercritical J = 1.2 height hills demonstrate interesting nonmonotonic behavior (highlighted by squares in Fig. 11). When J = 1, there is a maximum drag at ϵ = 0.48, in contradiction to linear theory. This nonmonotonic behavior can be explained as a consequence of harmonic resonance between the excitation frequency Uk and the responding buoyancy frequency N. When ϵ = 0.5, N = 2Uk, implying that the buoyancy frequency is a harmonic of the excitation frequency, amplifying the nonlinear response. On the other hand, when ϵ = 1/3 or 2/3, N = 3Uk and 3/2Uk, implying a dissonant buoyancy frequency and muted response.
In Fig. 12 we show the correlation between the average computed drag and the drag predicted by buoyancy adjusted linear theory Fadj [Eq. (15)]. We use Fadj rather than Flin to remove the effect of the reduced buoyancy frequency near the bottom on the predicted drag. From Fig. 12, it is clear that using linear theory based on the length and height scales of the bathymetry is valid in the linear and weakly nonlinear regimes, but fails as predicted by saturation theory when J = O(1). Figure 12 also suggests that the more nonhydrostatic runs (smaller circles) are more strongly affected by saturation.
As a measure of the strength of the correlation between simulations and parameterizations or theory, we use the coefficient of determination, defined as
is the sum of squared errors between the simulated values yi and their mean , and
is the sum of squared differences between the simulated values yi and the parameterized values fi. A parameterization that perfectly predicts the simulated value will give R2 = 1, while a negative value of R2 indicates that the observed values are better predicted by their mean than by the parameterization. The coefficients of determination for all steady-state parameterizations considered in this section are shown in Table 1. For linear theory using the adjusted buoyancy frequency (Fig. 12), the coefficient of determination is R2 = −53.1. The poor performance results from using the bathymetric height to scale the drag in supercritical lee waves.
Figure 13 shows the time-averaged LOTS for the simulations in Fig. 9. As posited by saturation theory, the LOTS is the effective bathymetry that generates lee waves, which can be used to compute the effective bathymetric power spectrum and thus predict the drag with the Fourier synthesis of Flin [Eq. (3)], as in Bell (1975b). Because this method employs all of the horizontal information in the LOTS, this operation represents a fully nonhydrostatic form of saturation theory. We perform this analysis on each run that exhibits blocking and plot the correlation between the average computed drag [Eq. (17)] and the drag predicted by FBell in the upper-left panel of Fig. 15. With R2 = 0.91, this combination of spectral linear theory with the LOTS is significantly better than the linear theory that uses the bathymetric height and length Fadj, for which R2 = −53.1 (Fig. 12). This demonstrates the strength of viewing the average LOTS as the effective wave-generating bathymetry.
In practice, the power spectrum of the LOTS above abyssal hills is not known a priori. Indeed, the usefulness of hydrostatic saturation theory is its prediction that, for supercritical bathymetry, there exists some average effective hill height of O(U/N) such that the effective wave generating Froude number is Jeff = O(1). To this end, we compute the average effective heights of each LOTS heff as indicated in Fig. 14, and define an effective Froude number, Jeff = Nheff/U, to form a simple prediction of the drag based on hydrostatic saturation theory
We perform this analysis on each run that exhibits blocking and compare the computed drag to the drag predicted by Feff in the upper-right panel of Fig. 15. Although using hydrostatic saturation theory to compute Feff gives a better prediction than linear theory (which has R2 = −53.1; see Table 1), with R2 = −0.296, Feff is still less predictive than the mean of the observations. It also uniformly over predicts the lee wave drag.
The discrepancy between the performance of spectral saturation theory FBell and hydrostatic saturation theory Feff suggest that there is a characteristic component of the supercritical LOTS that the latter neglects. We hypothesize that one important feature of the LOTS is the undulations downstream of the separation points with wavelengths roughly equal to λlee. If these undulations are treated as independent periodic bathymetry, they represent hills with ϵ ≈ 1 and are thus evanescent and should not contribute to the drag. The total effective height of the LOTS can therefore be partitioned into a propagating hprop and a nonpropagating hnonprop component, as illustrated in Fig. 14. Note that we define hnonprop as the maximum height of the undulations above the blocked layer rather than as the height of the jump-like feature at the beginning of the blocked layer because only the former of these features can be treated as independent periodic bathymetry. In this sense, the separation of heff into hprop and hnonprop represents an approximation of the power spectrum of the LOTS by two primary components. Only the propagating component contributes to the drag, suggesting that
The bottom-left panel of Fig. 15 compares the computed drag to Fprop. It performs much better than hydrostatic saturation theory with R2 = 0.733, indicating that the evanescent undulations are indeed dynamically relevant. However, Fprop is not as skilled as spectral saturation theory, and still tends to overestimate the drag.
The signal of an evanescent lee wave decays with height above the ocean floor at a rate determined by the vertical wavenumber [Eq. (1)]. Hence, a streamline sufficiently far from the bathymetry should contain contributions from only the propagating component of the lee wave. This observation leads to a final parameterization for the effective bathymetry, with an effective height haloft based on the maximum vertical displacement of a streamline originating at z = h0 + λlee/2 above the crest of the hill (red dotted lines in Fig. 13). The height haloft implies a new Froude number, Jaloft = Nhaloft/U, and a new parameterization for the drag
As shown in Fig. 15, this drag prediction agrees remarkably well with the computed drag, with a coefficient of determination of R2 = 0.956. Similarly excellent performance (not shown) results from selecting a streamline twice as far from the bottom. This parameterization has the added advantage that it is insensitive to the decay of buoyancy near the bottom, since the streamline upon which Jaloft is based is at a height above the observed decay. Furthermore, parameterizing the drag with the aloft streamline does not require the spectrum of the streamline, as in the calculation of FBell, only its maximum vertical excursion haloft. The combination of simplicity and accuracy in this parameterization suggests that a good estimate of the lee wave drag in the ocean may be obtained simply by measuring the vertical excursion of a streamline at a height zaloft = πU/N above the bathymetry.
Upon averaging across all of the simulations that display blocking (indicated by a tilde), we observe an average effective Froude number of , an average propagating Froude number of , and an average aloft Froude number of . Unlike most versions of saturation theory, in which Jeff = Jc, here all values for the effective Froude number are smaller than the observed critical Froude number Jc = 1, indicating that a supercritical lee wave has a smaller drag than a critical lee wave. Our values of , , and are also significantly smaller than the effective Froude number of Nikurashin and Ferrari (2010), JNF = 1.4 (this accounts for their definition of hill height, which is half of ours, i.e., h0 = 2hNF2010). Thus, the nonhydrostatic corrections to simple saturation theory embodied by Fprop and Faloft suggest that the lee wave drag over supercritical abyssal hills could be 50% smaller than currently thought. Furthermore, since Jaloft < Jc/2, these simulations suggest the most intense abyssal lee waves occur instead over critical and nearly critical abyssal hills, before the onset of blocking.
To summarize this section, Fig. 16 shows contour plots in J–ϵ space of the two measures of time-averaged drag ( and ) and the four supercritical parameterizations for drag discussed in this section (Feff, FBell, Fprop, and Faloft). Comparison of and show them to be identical (R2 = 0.97), indicating that just as in the linear regime of Pedlosky (2003), during the majority of the supercritical lee wave event over periodic bathymetry, the form drag and the vertical momentum flux are in balance, and there is no change in momentum near the ocean floor.
Figure 16 also highlights a curious ϵ-dependent feature that is not predicted by hydrostatic saturation theory, namely, the precipitous decline in simulated drag when J > 1 and ϵ > 0.5. We hypothesize that this drop off is due to evanescent masking, wherein the blocked layer grows horizontally to the point in which the attached portion of the LOTS becomes evanescent in length. Defining the horizontal length of this unseparated portion as Llaunch (see Fig. 14), we thus expect evanescent masking when
The boundary predicted upon measuring Llaunch in every simulation is sketched in the top-left panel of Fig. 16, and aligns remarkable well with the observed drop-off in drag. Evanescent masking thus represents a second nonhydrostatic process present in oceanic lee waves above periodic bathymetry which is also absent from current implementations of saturation theory. Like the evanescent undulations over the blocked layer, evanescent masking has the effect of reducing the lee wave drag above supercritical hills below that predicted by saturation theory. Indeed, for the most supercritical simulations in this study, evanescent masking renders even the hydrostatic hills practically evanescent.
In this paper, we used idealized simulations of lee waves over one-dimensional sinusoidal bathymetry to demonstrate that for oceanic flow over supercritical height bathymetry, the lowest layer of fluid becomes blocked and presents the overtopping flow with an effective bathymetry defined by the lowest overtopping streamline (LOTS) that is subcritical in height and polychromatic in wavelength. That the effective bathymetry reduces the apparent height of the hills is well known and currently employed in parameterizations of lee wave drag in global ocean models. However, both the development of evanescent undulations in the LOTS above the blocked layer and the process of evanescent masking, in which the unseparated portion of the LOTS becomes narrower than the lee wave wavelength, are novel results of nonhydrostatic lee wave processes in this paper. Recognizing their deleterious effect on the wave drag in regions of supercritical (J ≥ 1) bathymetry may help explain the discrepancy between predictions using existing parameterizations and observation of lee wave activity in the ocean discussed in Kunze and Lien (2019).
This study also offered an opportunity to observe the spinup of weakly nonlinear [J = O(0.1)] and supercritical (J ≥ 1) lee waves above periodic bathymetry. The excitation frequency, Tex = Lhill/U, emerged as the fundamental time scale of the lee wave for all simulations, regardless of hill height or length, informing both the time required for spinup and the vertical location of the drag on the background current. Furthermore, analysis of the horizontal momentum budget during spinup of supercritical lee waves permitted a characterization of time-dependent local drag associated with blocking.
The results in this paper identify four notable deficiencies in existing parameterizations of the lee wave drag. 1) In supercritical regions (J ≥ Jc), the simulations suggest that the steady-state lee wave drag will be smaller than predictions with existing saturation parameterizations due to the omission of nonhydrostatic effects of evanescent undulations and evanescent masking. 2) We observe that the drag associated with blocking is constrained to the first excitation period when the blocked layer forms, whereas predictions of the “nonpropagating” drag given by Garner (2005) and used in Trossman et al. (2013, 2015, 2016) assume the blocking component of the drag is a steady process that continually removes momentum from the flow. 3) On the other hand, in regions with large but still subcritical height bathymetry (J < Jc), our parameterizations predict larger lee wave drag because we identify the critical Froude number Jc and the effective Froude number Jaloft separately, whereas existing saturation theory assumes them to be identical. 4) Finally, our observations of upward-propagating lee wave wave fronts in every simulation suggest that all existing parameterizations need to reconsider confining the drag to within one lee wave wavelength above the bathymetry.
There are inevitably some limitations in applying our conclusions uniformly throughout the ocean. First, there is inherent uncertainty of approximately 30% in both the local height and length of abyssal hills in the Goff and Arbic (2010) abyssal hill bathymetry product. On a regional scale, this uncertainty is mitigated by the statistical homogeneity of the abyssal hills. Nevertheless, such significant local uncertainty of abyssal hill length and width suggests that some of the details in our process study, such as that of resonant responses for the J = 1 simulations, are probably not worth including in a drag parameterization.
Second, the predicted drags Fadj, FBell, Feff, Fprop, and Faloft presented in section 5 required running the simulations to find the effective Froude numbers and the adjusted buoyancy frequencies Nadj, meaning that these were not a priori parameterizations of the drag. With respect to the effective Froude numbers, the aloft Froude number Jaloft gave the most accurate prediction of the lee wave drag. Upon averaging across all supercritical simulations, we found the average aloft Froude number of . The variability in this average stems primarily from resonant responses rather than from nonhydrostatic effects. Hence our results support lee wave drag parameterizations replacing JNF = 1.4 with Jaloft = 0.4 above all two-dimensional abyssal hill bathymetry. Such a parameterization could also approximate the nonhydrostatic effect of evanescent masking [Eq. (25)] with a low-pass filter on the abyssal hill power spectrum to include only hills for which ϵ < 0.5. Using our results to parameterize Nadj is more difficult. As discussed in section 4, the decay of near-bottom buoyancy appears to depend on the progression of the valley water over the top of the hill during the first excitation period of the simulation. This suggests that Nadj is, at a minimum, a function of N, U, h0, and k. Additionally, it involves mixing, which for our simulations is a process driven by the elevated kinematic viscosity. In the real ocean, there are likely further complications such as effects of the bottom boundary layer and hysteresis from previous lee wave events. Hence the parameterization of Nadj requires further research.
Third, our study focused on single wavelength bathymetry, while the abyssal hills are better described by a linear superposition of many wavelengths. As such, our study cannot account for nonlinear interactions imposed by polychromatic nonlinear bathymetry. One hypothesis is that the flow will adjust by filling in all valleys with stretches of evanescent undulations as observed above, whereafter the remaining unblocked peaks will likely represent a subset of the original bathymetric wavenumbers. However, the hypothesis has not yet been tested with simulations. This is a significant limitation when attempting to employ the Goff and Arbic (2010) abyssal hill product.
Fourth, our study is of two-dimensional flow over one-dimensional bathymetry. In three-dimensional flow with two-dimensional bathymetry, fluid can travel around an obstacle rather than over it, and has been shown to reduce the saturation momentum flux by 40% (Nikurashin et al. 2014). However, abyssal hills are generally anisotropic, presenting the flow with a corrugation of ridges rather than a field of seamounts, which suggests that there may be little opportunity for flow around rather than over them.
Finally, our study ignores rotation. A primary consequence of rotation on lee waves is the lower bound on the bathymetric wavenumbers that produce propagating lee waves, based on the criterion f < Uk. As Uk → f, the vertical wavenumber approaches infinity, whereupon the vertical group velocity and momentum flux vanish. Hence, our conclusions about the drag for an upward propagating lee wave require modification for longer wavelength hills in strongly rotational regions. A second rotational process that complicates our conclusions is that of lee wave frequency-band narrowing, resulting from the near-bottom decay of N during spinup of nonlinear lee waves [Eq. (14)] such that the window of drag-producing bathymetry shrinks according to f < Uk < Nadj < N. If dramatic enough, or in regions of strong rotation, it could make all lee waves evanescent, snuffing out the lee wave drag (Kunze and Lien 2019). Finally, Nikurashin and Ferrari (2010) demonstrate that in the Southern Ocean where rotation is strong, saturation-level velocities in the lee wave field are sufficient to excite inertial oscillations and elicit positive feedback that results in wave breaking within the first lee wave wavelength above the bathymetry and a local deposition of momentum. Thus, in strongly rotational regions it may not be appropriate to deposit the momentum according to the vertical group velocity. However, even in this regime, the low-level breaking is a slowly developing process. Nikurashin and Ferrari (2010) observe the need for five or more days for inertial oscillations to grow over subcritical bathymetry even though the excitation period for their bathymetry was approximately just five hours. Therefore, there is ample time for the lee wave wave front to propagate away from the bottom, depositing momentum as it travels before the development of instabilities related to the inertial oscillations.
Despite the many simplifications, this paper demonstrates the remarkable applicability of linear lee wave theory to real ocean-scale flow. During spinup, Flin [Eq. (3)] offered a good estimate for the peak value of the form drag in all simulations, and was especially accurate when ϵ < 0.5, where evanescent masking is less pronounced. After spinup, the observed quasi-steady wave drag was best predicted with the same linear theory but with the hill height replaced with the maximum vertical extent of a streamline aloft Faloft [Eq. (24)]. The second best prediction came from inserting the power spectrum of the effective bathymetry into the spectral form of linear lee wave drag (Bell 1975b). These should not be foregone conclusions, since linear theory is strictly valid only for J ≪ 1, whereas the effective bathymetries in our simulations presented effective lee wave Froude numbers of Jaloft ≈ 0.4. Granted, achieving such accurate predictions required simulating the full flow to compute the adjusted bottom buoyancy frequency and the effective bathymetry. Still, our results demonstrate that even in the supercritical regime, linear theory can explain the lee wave drag. In this sense, they represent a qualified triumph of linear lee wave theory.
We gratefully acknowledge the support of ONR Grant N00014-16-1-2256 (scientific officers Dr. T. Paluszkiewicz and Dr. S. Harper).
Denotes content that is immediately available upon publication as open access.