## Abstract

Baroclinic instabilities are important processes that enhance mixing and dispersion in the ocean. The presence of sloping bathymetry and the nongeostrophic effect influence the formation and evolution of baroclinic instabilities in oceanic bottom boundary layers and in coastal waters. This study explores two nongeostrophic baroclinic instability theories adapted to the scenario with sloping bathymetry and investigates the mechanism of the instability suppression (reduction in growth rate) in the buoyant flow regime. Both the two-layer and continuously stratified models reveal that the suppression is related to a new parameter, slope-relative Burger number *S*_{r} ≡ (*M*^{2}/*f*^{2})(*α* + *α*_{p}), where *M*^{2} is the horizontal buoyancy gradient, *α* is the bathymetry slope, and *α*_{p} is the isopycnal slope. In the layer model, the instability growth rate linearly decreases with increasing *S*_{r} {the bulk form *S*_{r} = [*U*_{0}/(*H*_{0}*f*)](*α* + *α*_{p})}. In the continuously stratified model, the instability suppression intensifies with increasing *S*_{r} when the regime shifts from quasigeostrophic to nongeostrophic. The adapted theories are intrinsically applicable to deep ocean bottom boundary layers and could be conditionally applied to coastal buoyancy-driven flow. The slope-relative Burger number is related to the Richardson number by *S*_{r} = *δ*_{r}Ri^{−1}, where the slope-relative parameter is *δ*_{r} = (*α* + *α*_{p})/*α*_{p}. Since energetic fronts in coastal zones are often characterized by low Ri, that implies potentially higher values of *S*_{r}, which is why baroclinic instabilities may be suppressed in the energetic regions where they would otherwise be expected to be ubiquitous according to the quasigeostrophic theory.

## 1. Introduction

Baroclinic instabilities are ubiquitous throughout the global ocean. They release potential energy stored in horizontal density gradients, creating unsteady and evolving motions in the flow at approximately the first baroclinic deformation radius. As such, the Burger number (Bu) associated with baroclinicly unstable flow is Bu = *R*_{d}*L*^{−1} = RoRi^{1/2} ~ 1 (Eady 1949; Stone 1966, 1970). Beyond this, there are two general categories of baroclinic instabilities associated with the Rossby number (Ro) and the Richardson number (Ri) of the flow: mesoscale instabilities within the large-scale, geostrophic fronts with Ro ≪ 1 and Ri ≫ 1, and submesoscale instabilities with Ro ~ 1 and Ri ~ 1 (Boccaletti et al. 2007). Examples of the submesoscale instabilities include instabilities within the ocean mixed layers (Boccaletti et al. 2007; Fox-Kemper et al. 2008), along the fronts of mesoscale eddies (Callies et al. 2015; Brannigan et al. 2017), in the bottom boundary layers (Wenegrat et al. 2018), and, rarely, in certain coastal buoyancy-driven flows (Hetland 2017); they are potentially more nonlinear and energetic in nature than can be accurately described by quasigeostrophic (QG) theory.

The dynamics of QG baroclinic instabilities over a flat bottom has been well studied since the seminal papers by Eady (1949) and Phillips (1954), few studies have investigated the nongeostrophic scenario with a flat bottom (Stone 1966, 1970) and the QG scenario with a sloping bottom (Blumsack and Gierasch 1972; Mechoso 1980), but fewer studies have considered both the sloping-bathymetry effect and nongeostrophic effect. The submesoscale baroclinic instabilities in bottom boundary layers and coastal buoyancy-driven flows are influenced by these two effect (Wenegrat et al. 2018; Hetland 2017), and are less understood than the submesoscale instabilities in the surface mixed layers and the mesoscale instabilities, which are not influenced by the effects at the same time. One unique feature is the instability suppression—the submesoscale instabilities over sloping bathymetry exhibit weaker growth compared to the QG theories (Blumsack and Gierasch 1972; Mechoso 1980). Based on linear stability analysis, Wenegrat et al. (2018) shows that the growth of the instabilities in bottom boundary layers decreases as the regime shifts from QG to nongeostrophic. Also, although coastal buoyancy-driven flows are often associated with stronger lateral density gradients than open ocean fronts, they are seldom observed to be associated with instabilities; in particular, baroclinic instabilities are seldom observed in river plumes, even though lateral buoyancy gradients within the fronts are strong (Horner-Devine et al. 2015). External forcing agents, such as winds and tides, could suppress baroclinic instabilities through mixing processes. However, it has been demonstrated that a rotating plume without external forcing can be very stable over many rotational periods (Fong and Geyer 2002; Lentz and Helfrich 2002; Hetland 2005; Horner-Devine et al. 2006; Hetland 2017); this implies that the suppression can be due to the intrinsic inhibiting effects of the front. Baroclinic instabilities can enhance dispersion of water borne particles (Thyng and Hetland 2018), and decrease predictability in numerical simulations (Marta-Almeida et al. 2013). Better understanding the instability growth would improve our understanding on numerical predictability and be helpful for further investigating other submesoscale processes (e.g., symmetric instabilities and frontogenesis), material (e.g., nutrients and sediments) transport, mixing, and turbulence in the bottom boundary layers and coastal buoyancy-driven flows.

This paper explores nongeostrophic baroclinic instability theories adapted to the scenario with a sloping bathymetry. The two-layer model is adapted from Sakai (1989), and the continuously stratified model is an existing model adapted from Stone (1971) by Wenegrat et al. (2018). Both models are used to investigate the suppression of nongeostrophic baroclinic instabilities over sloping bathymetry, which is not revealed in the classical QG theories (Blumsack and Gierasch 1972; Mechoso 1980). In particular, this paper attempts to seek a nondimensional parameter for indicating the instability suppression and understand the underlying mechanism controlling the suppression in the buoyant flow regime.

## 2. Theory

### a. Layered model of nongeostrophic baroclinic instability

Phillips (1954) transformed the continuously stratified fluid to a two-layer system and constructed the layered model of QG baroclinic instabilities. Sakai (1989) investigated the ageostrophic instabilities using an ageostrophic version of the Phillips model. We adapt the Sakai model to the scenario with sloping bottom and surface (hereinafter referred to as the adapted Sakai model). The schematics of the Sakai model and the adapted Sakai model are shown in Fig. 1. The adapted Sakai model is a rotating two-layer channel with the sloping topography and background flow in the thermal wind balance. Considering the time scale as *f*^{−1}, the horizontal length scale as the Rossby deformation radius $Rd=\u2061(1/2)g\u2032H0/f$, and the vertical length scale as *H*_{0} (see the primitive equations and scale analysis in appendix A), the dimensionless form of the equations governing the perturbations is

subject to

where the subscript 1 denotes the variables in the upper layer, subscript 2 for the lower layer, *u* is the along-slope velocity, *υ* is the across-slope velocity, *p* is the pressure, ±*Y*_{max} are the across-slope boundaries, $H1=1\u2212(\delta +1)Rib\u22121/2y$ is the thickness of the upper layer, and $H2=1+(\delta +1)Rib\u22121/2y$ for the lower layer. Also, $Rib\u2261Frb\u22122=(g\u2032H0)/(2U02)$ is the bulk Richardson number, where $Frb=[U0/\u2061(1/2)g\u2032H0]$ is the bulk Froude number, *g*′ = (Δ*ρ*/*ρ*_{0})*g* is the reduced gravity, *H*_{0} is the dimensional average thickness, and *U*_{0} is the dimensional background flow (*U*_{0} in the upper layer and −*U*_{0} in the lower layer). The term *δ* ≡ *α*/*α*_{p} is the slope parameter, the ratio of the bottom slope *α* and the isopycnal slope *α*_{p} = 2*U*_{0}*f*/*g*′; both *α* and *α*_{p} are taken to be positive in this study. Here, we use the sign convention of *δ* as Hetland (2017) and Wenegrat et al. (2018), the opposite sign convention compared to the original one in the study of Blumsack and Gierasch (1972), such that positive *δ* represents the common case of buoyancy-driven flow over sloping bathymetry, where the isopycnal and bathymetric slopes are opposite. In this study, we will only focus on the scenarios with positive *δ*, which is referred to as the buoyant flow regime. The across-slope boundary $Ymax\u2261\Delta H/[\u2061(\delta +1)Rib\u22121/2]$ ensures that the total change of the layer thickness is Δ*H* [Δ*H* = 0.5 is set as in Sakai (1989)].

Substituting a wave-like solution $\varphi =\varphi \u02dc\u2061(y)ei\u2061(kx\u2212\sigma t)$, where $k=Rdk*$ is the dimensionless along-slope wavenumber and $\sigma =f\u22121\sigma *$ is the dimensionless wave frequency, into Eqs. (1) and (2) yields an eigenvalue problem (see appendix A). The growth rate of instabilities *σ*_{i} = Imag[*σ*] is obtained numerically through linear stability analysis. We use the software package DEDALUS for the stability analysis; DEDALUS is an open-source PDE solver written in Python that uses spectral methods (Burns et al. 2016). Figure 1 shows the normalized growth rate in the flat-bottom case (*δ* = 0) and the sloping-bottom case (*δ* = 0.5); see section 3a for more details.

### b. Continuously stratified model of nongeostrophic baroclinic instability

Eady (1949) developed a continuously stratified QG framework to study baroclinic instabilities. Blumsack and Gierasch (1972) adapted the Eady model to the scenario with a sloping bottom (hereinafter referred to as the BG model). Mechoso (1980) extended the BG model to the scenario with a sloping surface. Stone (1966, 1970, 1971) extended the Eady model to the nongeostrophic limit and constructed the continuously stratified model of nongeostrophic baroclinic instabilities. We use the modified model (Wenegrat et al. 2018), which adapts the Stone model to the scenario with sloping bottom and surface (hereinafter referred to as the adapted Stone model). The schematic of the BG, Mechoso, and adapted Stone models are shown in Fig. 2. In the adapted Stone model, the coordinates are rotated to align with the sloping bottom. The background buoyancy has a constant vertical gradient *N*^{2} and a constant horizontal gradient *M*^{2}, and the background flow is constrained by the thermal wind relation. Considering the time scale as *f*^{−1}, the horizontal length scale as *U*/*f*, and the vertical length scale as *H* (see the primitive equations and scale analysis in appendix B), the dimensionless form of the equations governing the perturbations in the rotated coordinates is

subject to

where *u* is the along-slope velocity, *υ* is the across-slope velocity, and *w* is the slope-normal velocity, *p* is the pressure, *b* is the buoyancy, *θ* is the bottom slope angle, $u\xaf=(z/cos\u2061\theta )$ is the background along-slope velocity, $b\xaf=(cos\theta \u2212\epsilon Ri\u22121\u2009sin\u2061\theta )z\u2212(\delta +1)Ri\u22121\u2009cos\u2061\theta \u2009y$ is the background buoyancy, and *z* = 0 (*z* = 1) is the sloping bottom (surface). Ri = *N*^{2}*f*^{2}*M*^{−4} is the Richardson number, and *δ* ≡ *α*/*α*_{p} = *αN*^{2}*M*^{−2} is the slope parameter, the ratio of the bottom slope *α* and the isopycnal slope *α*_{p} = *M*^{2}/*N*^{2}. The term *ε* = *f*^{2}*M*^{−2} is the nonhydrostatic parameter (Stone 1971).

Noting that the surface is also assumed to be tilted, the adapted Stone model is intrinsically suitable for baroclinic instabilities in a tilted bottom boundary layer, but seems not directly applicable to the situation with a flat surface that is not parallel with the tilted bottom (e.g., coastal buoyancy-driven flow). One way is to adapt the Stone model to the scenario with a flat surface, but this presents two challenges for theoretical approaches: first, making the uniform depth assumption like Eq. (4) is theoretically invalid; second and subsequently, assumed solutions with a wave form in the across-slope direction become theoretically invalid. Consequently, instead of adapting the Stone model to the flat-surface scenario, we will address the feasibility of the adapted Stone model (Wenegrat et al. 2018) in the flat-surface cases, for example, baroclinic instabilities in a coastal buoyancy-driven flow over a continental shelf.

Assuming a wave-like solution $\varphi =\varphi \u02dc\u2061(z)ei\u2061(kx\u2212\sigma t)$, where $k=(U/f)k*$ is the dimensionless along-slope wavenumber and $\sigma =f\u22121\sigma *$ is the dimensionless wave frequency (symmetric instabilities are excluded), then substituting it into Eqs. (3) and (4) yields the eigenvalue problem of the adapted Stone model (see appendix B). DEDALUS is also employed to calculate the growth rate of instabilities *σ*_{i} = Imag[*σ*]. Figure 2 shows the normalized growth rate based on the BG model (Blumsack and Gierasch 1972), the Mechoso model (Mechoso 1980), and the adapted Stone model (Wenegrat et al. 2018); See section 3b for more details. To keep consistent with the normalization in the QG models (the BG and Mechoso models), the dimensionless variables *σ*_{i} and *k* in the adapted Stone model are converted to the normalized variables $\sigma ^i$ and $k^$. The normalized wavenumber $k^$ is defined by normalizing the dimensional wavenumber $k*$ by the Rossby deformation radius *R*_{d} = *NH*/*f*; that is, according to the scaling relations,

The normalized growth rate $\sigma ^i$ is defined by normalizing the dimensional growth rate $\sigma i*$ by the advection time scale *R*_{d}/*U*,

With regard to the normalization in the adapted Sakai model, the wavenumber $k*$ is normalized by *R*_{d}, so $k^\u2261Rdk*$ is equivalent to the *k* in the adapted Sakai model, while the dimensionless $\sigma i\u2261\sigma i*f\u22121$ still needs to be converted to $\sigma ^i$ as Eq. (6). Hereafter, we will use the normalized variables ($\sigma ^i$ and $k^$) to describe the properties of baroclinic instabilities.

## 3. Results

In this section, we will address the suppression of nongeostrophic baroclinic instabilities over sloping bathymetry from the temporal perspective—the growth rate of instabilities—and discuss the suppression mechanisms in the layered and continuously stratified models. The suppression problem will be addressed in the dimensionless space; the dimensionless parameter space is Ri–*δ*, where Ri represents the nongeostrophic effect and *δ* represents the effect of sloping bathymetry.

### a. Suppression of instabilities in the layered model

Based on the adapted Sakai model, baroclinic instabilities start to be suppressed (meaning the growth rate of the instabilities is reduced), when the bottom slope increases (slope parameter *δ* increases from 0 to 0.5, see Fig. 1). Also, instabilities are found to be suppressed with decreasing bulk Richardson number Ri_{b}. The suppression of instabilities has opposite dependencies on Ri_{b} and *δ*. To demonstrate this, the maximum growth rate as a function of Ri_{b} and *δ* is calculated. The maximum normalized growth rate $\sigma ^i,max$ is defined as the maximum of the normalized growth rate across all wavenumbers $k^=k*Rd$ for given Ri_{b} and *δ*. The left panel of Fig. 3 shows the $\sigma ^i,max$ in the parameter space of Ri_{b} and *δ*. The maximum normalized growth rate $\sigma ^i,max$ exhibits opposite dependencies on Ri_{b} and *δ*, with the suppression as Ri_{b} decreases and *δ* increases. Baroclinic instabilities are sufficiently suppressed at low Ri_{b} and high *δ*. To understand the suppression of baroclinic instabilities, we will identify the primary parameter linked with the suppression and explore the controlling mechanism.

A new parameter, the slope-relative Burger number *S*_{r} ≡ (*M*^{2}/*f*^{2})(*α* + *α*_{p}), is considered as the coefficient controlling the suppression of instabilities in the nongeostrophic case. Compared to the horizontal slope Burger number *S*_{H} ≡ (*M*^{2}/*f*^{2})*α* (Hetland 2017), the slope-relative Burger number *S*_{r} uses a bottom slope relative to the isopycnal slope *α*_{p}. Thus, if the slopes are aligned, with the bottom parallel to the isopycnals, then *S*_{r} = 0. The slope-relative Burger number uses the layer thicknesses in the adapted Sakai model; $H1,2*=H0\u2213(\alpha +\alpha p)y*$ may be used to write a horizontal Burger number as (*M*^{2}/*f*^{2})(*H*/*L*) = (*M*^{2}/*f*^{2})(*α* + *α*_{p}). The slope-relative Burger number can be written in terms of the Richardson number Ri and the slope parameter *δ* as

where *δ*_{r} is the slope-relative parameter,

with an interpretation similar to the slope-relative Burger number *S*_{r}. For clarity, hereafter 1 + *δ* will be written as *δ*_{r}.

The bulk form of the slope-relative Burger number is

The distribution of *S*_{r} in Ri_{b}–*δ* space (shown in the center panel of Fig. 3) exhibits a similar pattern as the maximum normalized growth rate $\sigma ^i,max$. Furthermore, $\sigma ^i,max$ linearly decreases with increasing *S*_{r} (see the right panel of Fig. 3); *S*_{r} contains the opposite dependencies of $\sigma ^i,max$ on Ri_{b} and *δ*. In addition, instabilities are sufficiently suppressed for $Sr\u22730.1$. The linear relation between *S*_{r} and $\sigma ^i,max$ does not hold in the QG limit [e.g., Ri_{b} ~ *O*(100)].

A strength of the Sakai model is to interpret instabilities in the framework of wave resonance based on the physical wave coordinates (Sakai 1989). For instance, Kelvin–Helmholtz instabilities are interpreted as the resonance of interacting gravity waves. In this study, we will only focus on the interactions between Rossby waves to interpret the suppression of baroclinic instabilities. The physical wave coordinates, then, only consist of the Rossby waves; since the Rossby wave coordinates are the Fourier basis (see appendix C), the projection onto the physical wave coordinates is equivalent to conducting the Fourier transform (Pedlosky 2013). The details of the Rossby wave resonance in the adapted Sakai model can be found in appendix C. The resonance rate $R*$ in the adapted Sakai model is, then,

where $\epsilon n*=1/[2Rd2\u2061(k*2+ln*2)+1]$ is the interaction coefficient, $k*$ is the dimensional along-slope wavenumber, and $ln*=(n\pi )/(2Ymax*)$ (*n* = 1, 2, 3, …) is the dimensional across-slope wavenumber. If *δ*_{r} = 1 (the flat-bottom case), Eq. (10) is reduced to $R*=Imag\u2061(k*U01\u2212{2/[Rd2\u2061(k*2+ln*2)+1]})$, which reproduces Eq. (28) in Sakai (1989). The left panel of Fig. 4 shows the comparison between the nonzero $\sigma ^i,max$ and the normalized resonance rate $Rib1/2f\u22121R*$, where $R*$ is calculated by only considering the resonance of the lowest mode (*n* = 1) and using the $k*$ corresponding to the maximum growth rate. The resonance rate closely follows the growth rate, especially at high growth rates (presumably, the discrepancy at low growth rates is due to the only consideration of the resonance of the lowest mode). This implies that the maximum growth of instabilities is highly related to the resonance of the lowest mode of Rossby waves. Furthermore, the right panel of Fig. 4 shows that the resonance rate linearly decreases with the increasing *S*_{r}, and they are highly correlated (*r*^{2} = 0.999, *p* = 0.0); it indicates that *S*_{r} controls the growth of instabilities by influencing the Rossby wave resonance.

The underlying mechanism is that *S*_{r} influences the wave resonance by modifying the properties of the Rossby waves. According to appendix C, the normalized Doppler-shifted phase speed of the Rossby waves in both layers is

where $c^1,2=\u2213[\delta r/\u2061(k2+ln2+0.5)]$ are the normalized phase speed, *k* is the normalized along-slope wavenumber, and $ln=[(n\pi )/(2\Delta H)]\delta rRib\u22121/2$ (*n* = 1, 2, 3, …) is the normalized across-slope wavenumber. The slope-relative parameter *δ*_{r} could be considered as the normalized topographic beta, which is $\u2061(f/H0)|\u2061(dH1,2*/dy)|=f\u2061(\alpha +\alpha p)/H0$ normalized by *fα*/*H*_{0}. Noticing that $Sr\u22121\u221d(\delta r/ln2)$ (part of $c^1,2$), the phase speed $|c^1,2|$ would decrease with increasing *S*_{r}; higher *S*_{r} would indicate larger difference of the Doppler-shifted phase speed $c^D$, thereby larger difference of the Doppler-shifted frequency, and hence weaker wave resonance. Figure 5 shows that the difference between $c^1,2D$ does, in fact, increase with increasing *S*_{r} (using the lowest mode and the *k* corresponding to the maximum growth rate). As the schematics of Fig. 5 show, low *S*_{r} represents the scenario where $c^1,2D$ are close to each other so that the Rossby waves could interact efficiently to create strong wave resonance, while high *S*_{r} represents the scenario where $c^1,2D$ are different so that the wave resonance would be weaker. In sum, *S*_{r} influences the Rossby wave resonance (thus the instability growth) by altering the Doppler-shifted phase speed of the Rossby waves, thereby controlling the growth of instabilities.

### b. Suppression of instabilities in the continuously stratified model

The suppression of instabilities through the reduction of instability growth rate is also found in the adapted Stone model (Wenegrat et al. 2018). Particularly, the suppression is significant when the regime shifts from QG to nongeostrophic. To demonstrate the suppression, the QG models, that is, the BG model (Blumsack and Gierasch 1972) and the Mechoso model (Mechoso 1980), are briefly reviewed here. In the BG model, the background state consists of a density field with constant vertical and horizontal buoyancy gradients and a velocity field constrained by the thermal wind relation. By assuming that the Rossby number Ro = *U*/(*fL*) ≪ 1 and Burger number Bu = (*NH*)/(*fL*) ~ 1 (leading to the Richardson number Ri = Bu^{2}Ro^{−2} ≫ 1), the equations governing perturbations can be reduced to a single PDE about the QG potential vorticity. Then, the growth rate of the instabilities is analytically obtained by solving the associated eigenvalue problem; the normalized growth rate of the BG model is then

where $\sigma ^i=Ri1/2\sigma i*f\u22121$ is the normalized growth rate, $k^=k*Rd$ is the normalized wavenumber, and *δ* is the slope parameter. On the other hand, the Mechoso model relaxes the flat-surface assumption of the BG model and takes into account a surface that could be tilted as the bottom would be. Here, we only consider the scenario where the surface is parallel with the sloping bottom; the normalized growth rate of the Mechoso model is, then,

Figure 2 shows $\sigma ^i$ of the QG models and the adapted Stone model. The tail of the unstable modes is significantly smaller and weaker in the adapted Stone model than the QG models—the adapted Stone model exhibits a significant reduction of instabilities and predicts lower growth rates than the QG models.

Maximum growth rates predicted by the QG and adapted Stone models have different dependencies on the Richardson number Ri. The maximum normalized growth rate $\sigma ^i,max$ is defined as the maximum of the normalized growth rate across all wavenumbers, $k^$, for a given Ri and *δ*. The Ri–*δ* space is spanned within 1 ≤ Ri ≤ 5 and 0 ≤ *δ* ≤ 0.6, which follows Hetland (2017). This parameter space covers most scenarios of energetic coastal fronts (Hetland 2017) and a large variety of energetic deep ocean BBLs (Wenegrat et al. 2018). The maximum normalized growth rate $\sigma ^i,max$ for the adapted Stone model is denoted as $\sigma ^NG$, $\sigma ^BG$ for the BG model, and $\sigma ^M$ for the Mechoso model. Figure 6 shows $\sigma ^NG$, $\sigma ^M$, and $\sigma ^BG$ in the parameter space of Ri–*δ*. Here we implicitly associate QG theories to high Ri conditions, because it is unknown how well QG theories describe flow at low Ri conditions. The growth rates $\sigma ^M$ and $\sigma ^BG$ decrease with *δ*, but do not vary with Ri, since the QG normalized growth rates are independent of Ri [see Eqs. (12) and (13)]. In contrast, $\sigma ^NG$ exhibits a clear dependency on Ri, decreasing as Ri decreases.

A dimensionless number is sought for indicating the suppression of instabilities in the adapted Stone model. One requirement of the number is that it can reduce to the QG limit (Blumsack and Gierasch 1972; Mechoso 1980) and the nongeostrophic limit with flat bathymetry (Stone 1970). The following discussion reviews the dimensionless numbers $\delta r\u22121$ and (1 + Ri^{−1})^{−1/2}, which controls the instability growth in the QG limit and nongeostrophic limit with flat bathymetry, respectively. First, as demonstrated in Pedlosky (2016), the slope-relative parameter *δ*_{r} is the only dimensionless number appearing at the bottom boundary condition of the Eady problem [Pedlosky 2016, Eq. (5b)]; it involves the vertical shear of the background flow that is supplemented by the topographic production of vertical vorticity by the perturbed across-slope flow. Figure 7 (left) shows $\sigma ^M$ and $\sigma ^BG$ as functions of $\delta r\u22121$; $\delta r\u22121$ has robust relations with $\sigma ^M$ (*r*^{2} = 0.999, *p* = 0.0) and $\sigma ^BG$ (*r*^{2} = 0.993, *p* = 0.0)—$\delta r\u22121$ is the controlling number in the QG limit. Second, as demonstrated in Fox-Kemper et al. (2008), the maximum normalized nongeostrophic growth rate in the flat-bottom case is a linear function of (1 + Ri^{−1})^{−1/2},

which is obtained based on the asymptotic analysis of Stone (1970). Figure 7 (right) shows the comparison between the numerical $\sigma ^NG|\delta =0$ and (1 + Ri^{−1})^{−1/2}; (1 + Ri^{−1})^{−1/2} has a robust relation with $\sigma ^NG|\delta =0$ (*r*^{2} = 0.999, *p* = 0.0). Equation (14) is also shown in Fig. 7 (right); the offset between $\sigma ^NG|\delta =0$ and $\sigma ^NG,FlatAsym$ is presumably due to the approximation and truncation of the asymptotic analysis of Stone (1970). But both the analytical and numerical solutions suggest that (1 + Ri^{−1})^{−1/2} is the controlling number in the nongeostrophic limit with flat bathymetry.

Noticing that *S*_{r}|_{Flat} = Ri^{−1}, the controlling number (1 + Ri^{−1})^{−1/2} might be the reduction of (1 + *S*_{r})^{−1/2} in the flat-bottom case. The multiplication of $\delta r\u22121$ (indicating the QG growth) and (1 + *S*_{r})^{−1/2} (as the nongeostrophic modification) would be a physically intuitive form of the dimensionless number in the adapted Stone model, because it can simply reduce to $\delta r\u22121$ as Ri → ∞ (the QG limit; Blumsack and Gierasch 1972; Mechoso 1980) and (1 + Ri^{−1})^{−1/2} as *δ* → 0 (the nongeostrophic limit with flat bathymetry; Stone 1970). The dimensionless number $\delta r\u22121\u2061(1+Sr)\u22121/2$ in the Ri − *δ* space is shown in Fig. 8a, which has a similar distribution as the $\sigma ^NG$ (shown in Fig. 6, left). Figure 8c shows the comparison between $\delta r\u22121\u2061(1+Sr)\u22121/2$ and $\sigma ^NG$; $\sigma ^NG$ is highly related to the form indicated by a robust linear regression (*r*^{2} = 0.991 and *p* = 0.0). Given that the QG growth is controlled by $\delta r\u22121$, the modification of (1 + *S*_{r})^{−1/2} on $\delta r\u22121$ stands for the suppression of instabilities occurring when the regime shifts from QG to nongeostrophic (as shown in Fig. 6).

The fact that $\sigma ^NG$ increases with increasing $\delta r\u22121\u2061(1+Sr)\u22121/2$ indicates that the instability suppression intensifies with increasing *δ*_{r} (or *δ*) and *S*_{r}. Energetics is used to understand these dependencies. According to appendix D (following Wenegrat et al. 2018), the eddy potential energy (EPE) and eddy kinetic energy (EKE) equations (the primed variables are the dimensional perturbations) are

and

where $HBFc=\upsilon \u2032b\u2032\xaf\u2009cos\u2061\theta \u2061(M2/N2)$ is the horizontal buoyancy flux (HBF) contributed by the cross-slope motion, $HBFn=w\u2032b\u2032\xaf\u2009sin\u2061\theta \u2061(M2/N2)$ is the HBF contributed by the slope-normal motion, $VBFc=\upsilon \u2032b\u2032\xaf\u2009sin\u2061\theta $ is the vertical buoyancy flux (VBF) contributed by the cross-slope motion, $VBFn=w\u2032b\u2032\xaf\u2009cos\u2061\theta $ is the VBF contributed by the slope-normal motion, and other terms are specified in appendix D. For baroclinic instabilities, HBF_{c} and HBF_{n} are the driving energy sources [HBF_{n} ≪ HBF_{c} because (HBF_{n}/HBF_{c}) ~ *ε* tan*θ*], which convert the mean available potential energy (APE; stored in the sloping isopycnals) to EPE, as the classic nonrotated energetics would suggest. However, the two VBF terms exhibit opposite roles in terms of transferring energy. VBF_{n} represents the stratifying effect of the instabilities; it is the energy transfer, induced by the slope-normal motion, from EPE to EKE and hence is an energy source of EKE and a sink of EPE. VBF_{c} represents the EPE gained by overcoming the cross-slope gravity during the stratifying process; so it is an energy sink of EKE and a source of EPE. Based on the scaling relation in appendix B and the nonhydrostatic parameter *ε* ≪ 1, $VBFc=\upsilon \u2032b\u2032\xaf\u2009sin\u2061\theta $ can be scaled as

which is the ratio of two different types of EPE sources, that is, the EKE conversion and APE conversion to EPE; larger *δ* indicates larger EKE conversion back to EPE and hence stronger suppression. On the other hand, $VBFn=w\u2032b\u2032\xaf\u2009cos\u2061\theta $ can be scaled as the ratio of the EPE sink and all EPE sources, that is,

larger *S*_{r} indicates smaller EPE conversion to EKE and hence weaker instability growth.

It has been suggested by previous studies that the horizontal slope Burger number *S*_{H} = *δ*/Ri is relevant to the instability suppression (Hetland 2017; Wenegrat et al. 2018). Wenegrat et al. (2018) demonstrate that VBF_{c}/VBF_{n} can be scaled as *S*_{H}; larger *S*_{H} indicates larger reduction of VBF and hence weaker instabilities. The difference of scaling is that *S*_{H} scales the ratio between the VBF components rather than the components as Eqs. (17) and (18) do. Figure 8b shows *S*_{H} in the Ri − *δ* space, and Fig. 8d shows the comparison between *S*_{H} and $\sigma ^NG$. The growth rate $\sigma ^NG$ does, generally, decrease with increasing *S*_{H} but in a complex way. Although both *S*_{H} and $\delta r\u22121\u2061(1+Sr)\u22121/2$ contain the growth dependency on the VBF_{c} (EKE sink) and VBF_{n} (EKE source), the later formula seems to be more accurate in term of indicating growth rate.

## 4. Discussion

### a. Application to coastal buoyancy-driven flow

The adapted theories are directly applicable to baroclinic instabilities in a bottom boundary layer (Wenegrat et al. 2018), where both the bottom and surface can be assumed to be tilted and parallel. However, these theories may not seem as directly applicable to surface-intensified baroclinic instabilities that are usually associated with tilted bottoms but flat surfaces (e.g., instabilities in coastal fronts). In this section, we will address, by scale analysis, the applicability of the adapted Stone model to the baroclinic instabilities in coastal buoyancy-driven flow over sloping bathymetry.

The following scale analysis demonstrates that the adapted Stone model can be used as an approximation of the flat surface case so long as the slope Burger number $S\u2272O\u2061(10\u22121)$ and the horizontal slope Burger number $SH\u2272O\u2061(10\u22121)$. The scaling relations in appendix B are used in the scale analysis. In the rotated coordinates, a flat surface in the dimensional form can be written as $z*=H+\alpha \u2061(y*\u2212yc*)$ is the depth, *H* is the depth at the center of the front $yc*$, and *α* is a constant bottom slope. Considering the motions within the length scale of baroclinic instabilities that are on the order of the Rossby deformation radius *R*_{d} = *NH*/*f* (Eady 1949; Stone 1966, 1970), the dimensional flat surface can be scaled as

where *S* ≡ (*N*/*f*)*α* = *δ*Ri^{−1/2} is the slope Burger number. For the situations characterized by *δ* ~ *O*(10^{−1}) and Ri ~ *O*(1) (which are the cases in this study), the slope Burger number *S* is *O*(10^{−1}) so that the uniform fluid depth *z* = 1 will be a reasonable assumption. In other words, if *S* ~ *O*(10^{−1}), motions with lateral displacements *O*(Rd) will span a depth range of *H* ± 0.1*H*, and hence the uniform fluid depth will be a first-order approximation.

With the uniform fluid depth approximation, the rigid-lid boundary condition at the flat surface in the rotated coordinates is $w*|z*=H=\alpha \upsilon *$. The dimensionless form is, then,

where *S*_{H} ≡ *αM*^{2}*f*^{−2}= *δ*Ri^{−1} is the horizontal slope Burger number (Hetland 2017). For the situations characterized by *δ* ~ *O*(10^{−1}) and Ri ~ *O*(1) (which are the cases in this study), the horizontal slope Burger number *S*_{H} is *O*(10^{−1}) so that the no-flow boundary condition $w|z=1=0$ will be an acceptable approximation. In summary, if a front over a sloping bottom satisfied $S\u2272O\u2061(10\u22121)$ and $SH\u2272O\u2061(10\u22121)$, the assumptions of the uniform fluid depth and the no-flow boundary condition will be first-order accurate, and hence the adapted Stone model can be used to represent the surface-intensified baroclinic instabilities formed within the front.

### b. Numerical simulations of instabilities in coastal buoyancy-driven flow

A series of existing idealized numerical simulations, first analyzed by Hetland (2017), are used to examine the feasibility of the adapted Stone model. The idealized model domain is a 260 km (along-slope) × 128 km (across-slope) continental shelf with the uniform bathymetric slope *α* = 10^{−3} across all simulations. The depth increases from 5 m onshore to 133 m offshore. The model grid has 1-km uniform horizontal resolution and 30 layers in the vertical direction. The boundary conditions are periodic along-slope, open (with a sponge layer) offshore, and closed (no-slip) at the coast. The *k*–*ε* turbulence closure scheme is used to calculate the vertical mixing, and bottom friction is defined using a specified bottom roughness and a log-layer approximation. The model is unforced and run as an initial-value problem. The initial density field is a coastal buoyant front with a constant vertical stratification *N*^{2} over the whole shelf and a constant lateral buoyancy gradient *M*^{2} within the offshore distance of *W* = 50 km. The initial current field is configured in the thermal wind balance with the density field. The initial density and current fields of the base case (Ri = 2.0 and *δ* = 0.1) are shown in Fig. 9. Initial fields are varied among simulations to cover a range of situations of instability formation.

The idealized simulations were configured in the parameter space of the Richardson number Ri and the slope parameter *δ*, and all the simulations used in this study are listed in Table 1. All the simulations were run with same stratification *N*^{2} and same bottom slope *α*, but with different lateral buoyancy gradients *M*^{2} and Coriolis parameters *f* that are determined by each combination of Ri and *δ*. Note that the simulations listed in Table 1 are part of the simulations conducted in Hetland (2017); simulations with the Richardson number Ri = 1 or Ri = 10 or the slope parameter *δ* > 0.5 or the horizontal slope Burger number *S*_{H} > 0.2 were excluded. Ri = 1 is around the boundary between baroclinic instabilities and symmetric instabilities (Haine and Marshall 1998; Boccaletti et al. 2007), so the simulations with Ri = 1 are excluded to ensure only baroclinic instabilities can form. The simulations with Ri = 10 are excluded to ensure the nongeostrophic regime, and to minimize the influence of bottom friction in these simulations with long instability growth rates. The instability formations in the simulations with *δ* > 0.5 are excluded because they are generally weak and also strongly influenced by bottom friction. The simulations with *S*_{H} > 0.2 are excluded to ensure $SH\u2272O\u2061(10\u22121)$ and $S\u2272O\u2061(10\u22121)$ so that the adapted Stone can be applied to the flat-surface case. One example of the development of instabilities in the base run is shown in Fig. 10.

The growth rate of instabilities is estimated in each simulation and then compared to the QG model and the adapted Stone model. Note that the spatial scale of the instabilities increases in the offshore direction (see Fig. 10) because the deformation radius, Rd = *NH*/*f*, increases with increasing depth offshore. We only focus on the growth of the largest instabilities at ~50 km offshore with the depth *H* ~ 50 m, because they are the most energetic and dispersive components (Thyng and Hetland 2017, 2018).

EKE is used to quantify the growth rates of instabilities. Given that the domain is periodic in the along-slope direction, the velocity field associated with the instability eddies is calculated by subtracting the along-slope background velocity from the original velocity field. So, EKE is dominated by the largest eddies. Then, the EKE can be determined by integrating the kinetic energy of the eddy flow field over the whole domain. Last, the EKE is normalized by the initial domain-integrated mean kinetic energy MKE_{Initial}. Figure 11 (top) shows the normalized EKE, EKE/MKE_{Initial}, of all the simulations listed in Table 1.

The EKE in each case appears to increase exponentially from the start (Fig. 11, top), but eventually the rapid increase is retarded by friction and finite amplitude effects. To isolate our results from these frictional effects and compare our results more directly with the theories that do not consider the influence of friction, we truncate the EKE time series where it reaches half of the maximum, removing the later part that is potentially influenced by friction. The truncated time scale for each simulation is listed in Table 1, and the truncated EKE time series are shown in Fig. 11 (middle). We take the base case (Ri = 2.0, *δ* = 0.1) as an example to show the comparison between the simulated growth rate and the theoretical predictions. The truncated EKE time series of the base run is shown in Fig. 11 (bottom), and the best exponential function to fit it has a growth rate of 1.73 day^{−1} (*r*^{2} = 0.996). On the other hand, the maximum dimensional growth rate for the base case is estimated based on the adapted Stone model (Wenegrat et al. 2018), the flat-bottom Stone model (Stone 1971), and the QG model (Blumsack and Gierasch 1972); they are $\sigma NG*=1.82\u2009day\u22121$, $\sigma NG,Flat*=2.25\u2009day\u22121$, and $\sigma BG*=2.50\u2009day\u22121$, respectively. The theoretical predictions as exponential functions with the estimated growth rates are shown in Fig. 11 (bottom); as expected from the growth rate estimates, the adapted Stone prediction tracks the simulated increase in EKE closer than the flat-bottom Stone prediction and the QG prediction.

A regressed estimate of simulated EKE growth rate, $\sigma R*$, is calculated for each simulation listed in Table 1, in order to compare with the theoretical predictions. The calculated values of $\sigma R*$, $\sigma NG*$, $\sigma BG*$, and $\sigma NG,Flat*$ for each simulation are listed in Table 1. The upper panels of Fig. 12 show the comparison between the regressed growth rates and the theoretical predictions. We can see that $\sigma NG*$ values better follow $\sigma R*$ than $\sigma BG*$ and $\sigma NG,Flat*$. Two-sided *t* tests are conducted to see if the regressed and predicted growth rates are statistically equivalent. The two-tailed *p* value for the adapted Stone theory comparison is *p* = 0.42, for the QG theory *p* = 0.04, and for the flat-bottom Stone theory *p* = 0.03. We cannot reject the null hypothesis for the adapted Stone theory if a *p* value threshold of 5% is used, indicating the regressed and predicted growth rate distributions are indistinguishable. However, we can reject the null hypothesis in the tests for the QG theory and the flat-bottom Stone theory; the regressed and predicted distributions are distinct, as apparent from the offset from the 1:1 line in Figs. 12b and 12c. Moreover, the lower panels of Fig. 12 show the distribution of growth rate errors in the Ri–*δ* space. Compared to the adapted Stone model, the QG model and the flat-bottom Stone model have higher growth rate errors, particularly in the low Richardson number cases (Ri = 2.0 and 3.0); this implies that these models are not able to accurately describe the development of the submesoscale baroclinic instability eddies under energetic flow situations. However, under the conditions of $SH\u2272O\u2061(10\u22121)$ and $S\u2272O\u2061(10\u22121)$, the adapted Stone model accurately predicts the growth rate of the instabilities in the energetic flow situations. Consequently, the numerical simulations complement the scale analysis in section 4a and validate the applicability (under certain conditions) of the adapted Stone model to the flat-surface situations.

Both the nongeostrophic effect and sloping bathymetry are important for indicating the instability growth. The overestimation of growth rates in QG theory is presumably due to the missing of the nongeostrophic suppressing effect. Applying the flat-bottom Stone theory misses the suppressing effect of the sloping bathymetry. From the perspective of energetics, assuming flat-bottom (*δ* = 0) will overestimate the total VBF (by underestimating the EKE sink VBF_{c} and overestimating the EKE source VBF_{n}) and hence underestimate the suppression. From the perspective of wave interactions, assuming the flat-bottom will facilitate/overestimate the wave interaction and hence underestimate the suppression as well.

Furthermore, the adapted Stone theory might overestimate the suppression in the coastal case by assuming a sloping surface, as the comparison between the BG and Mechoso models suggests (see Fig. 2). But the nonlinear simulations indicate that the adapted Stone theory is most accurate among these existing theories, which suggests that the overestimation is not prohibitively large. On the other hand, there is little improvement when moving from the QG theory to the flat-bottom Stone theory; the underestimation of the suppression by assuming a flat-bottom seems not negligibly small.

### c. Linking S_{r} to potential vorticity

Above, the parameter *S*_{r} is linked to instability growth rate. Here we show how this parameter is related to normalized potential vorticity. The reference potential vorticity in this case is *fN*^{2}, the case with no horizontal density gradients and thereby no associated balanced flow. In the adapted Stone model, the normalized Ertel potential vorticity (Thomas et al. 2016) can be written as

For the adapted Stone model, the potential vorticity is uniform throughout the domain, but in the adapted Sakai model, the sloping interface creates a gradient in potential vorticity. In this case, we consider the change in potential vorticity, ΔPV, across an inertial radius *L*_{i} = *U*_{0}/*f*

where *β*_{T} = [(*α* + *α*_{p})*f*]/*H*_{0} is the topographic beta supplemented by the isopycnal slope. Thus, if we consider a parcel going from high to low potential vorticity, the normalized potential vorticity will be $1\u2212\Delta PVPV0\u22121=1\u2212Sr$, where PV_{0} = *f*/*H*_{0} is the reference potential vorticity starting location on the high potential vorticity side of the inertial circle. Thus the two interpretations are somewhat similar.

## 5. Conclusions

Layered and continuously stratified models of nongeostrophic baroclinic instability over sloping topography are explored in the buoyant flow regime to study the suppression mechanism. The layered model (i.e., the adapted Sakai model) reveals a new parameter, the slope-relative Burger number *S*_{r} ≡ (*M*^{2}/*f*^{2})(*α* + *α*_{p}) {the bulk form is *S*_{r} = [*U*_{0}/(*H*_{0}*f*)](*α* + *α*_{p})}, which is a dimensionless parameter that controls the suppression of instabilities in the nongeostrophic limit. The continuously stratified model (i.e., the adapted Stone model) shows that, when the regime shifts from QG to nongeostrophic, the growth of instabilities is inhibited with increasing *S*_{r}.

In the adapted Sakai model, the instability growth rate linearly decreases with increasing *S*_{r}. The physical mechanism behind *S*_{r} is explored based on the wave resonance theory (Sakai 1989). In the physical wave coordinates consisting of the Rossby waves, baroclinic instabilities are interpreted as the Rossby wave resonance; supported by the adapted Sakai model, the maximum normalized growth rate of instabilities is found to be nearly 1:1 to the resonance rate of the Rossby waves. The slope-relative Burger number *S*_{r} modifies the Doppler-shifted phase speed of the Rossby waves, alters the wave resonance, and hence influences the growth of instabilities.

In the adapted Stone model, the instability growth rate linearly decreases with decreasing $\delta r\u22121\u2061(1+Sr)\u22121/2$. This relation indicates that the suppression intensifies with increasing *δ*_{r} (or *δ*) and *S*_{r}. Supported by the energetics, larger *δ*_{r} represents larger VBF_{c}, an EKE sink, and hence stronger suppression; large *S*_{r} represents smaller VBF_{n}, an EKE source, and hence weaker instability growth. Given that $\delta r\u22121$ indicates the QG growth, (1 + *S*_{r})^{−1/2} represents the instability suppression when the regime shifts from QG to nongeostrophic.

The adapted models are intrinsically applicable to baroclinic instabilities in a bottom boundary layer but has one limitation when applying to coastal buoyancy-driven flow—the sloping-surface assumption. Supported by scale analysis, idealized numerical simulations of coastal buoyancy-driven flow are used to test the feasibility of the adapted Stone model in flat-surface situations. The comparison of the numerical results and theoretical predictions indicates that the limitation is not prohibitive if the slope Burger number $S=(N/f)\alpha \u2272O\u2061(10\u22121)$ and the horizontal slope Burger number $SH=(M2/f2)\alpha \u2272O\u2061(10\u22121)$, which complements the scale analysis. The numerical results also show that neither the BG model (Blumsack and Gierasch 1972) nor the flat Stone model (Stone 1971) can accurately describe the growth of the submesoscale baroclinic instabilities in a coastal buoyancy-driven flow; it implies that both the nongeostrophic effect and sloping bathymetry are important. The scale analysis and numerical simulations validate the feasibility of the adapted Stone model in coastal situations, and hence the inhibiting mechanism provided by the sloping-surface theories can be conditionally used to interpret the suppression of instabilities in the buoyancy-driven flow of coastal zones.

## Acknowledgments

This work was funded by NSF (Grant OCE-1851470), Texas General Land Office (Contract 18-132-000-A673), and Gulf of Mexico Research Initiative through the CSOMIO consortium (Contract R01984). Lixin Qu was supported by a graduate research fund from Texas Sea Grant and a scholarship from China Scholarship Council. We thank Jacob Wenegrat, Leif Thomas, and two anonymous reviewers for very helpful suggestions and codes when preparing this manuscript.

### APPENDIX A

#### Adapted Sakai Model

The following derivation follows Sakai (1989) but is modified to account for the presence of sloping bottom and surface. Considering a rotating two-layer channel with sloping bottom and top and currents in the thermal wind balance, the linearized equations of the perturbations are

subject to

where $u*$ is the perturbed along-slope velocity, $\upsilon *$ is the perturbed across-slope velocity, $p*$ is the perturbed pressure, $h*$ is the interface displacement, *g*′ = (Δ*ρ*/*ρ*_{0})*g* is the reduced gravity, *f* is the Coriolis parameter, $\xb1Ymax*$ are the across-slope boundaries, and the upper-and lower-layer variables are denoted by the subscripts of 1 and 2, respectively. The background flows in the upper and lower layers are set to *U*_{0} and −*U*_{0} for simplicity (Sakai 1989). Also, $H1*=H0\u2212\alpha y*\u2212\alpha py*$ is the thickness of the upper layer, and $H2*=H0+\alpha y*+\alpha py*$ for the lower layer, where *α* is the bottom slope and *α*_{p} = 2*U*_{0}*f*/*g*′ is the isopycnal slope.

Considering the time scale as 1/*f*, the horizontal length scale as the Rossby deformation radius $Rd=\u2061(1/2)g\u2032H0/f$, and the vertical length scale as *H*_{0}, the scaling relations about the variables in Eq. (A1) are

where $Rib=(g\u2032H0)/(2U02)$ is the bulk Richardson number and *δ*_{r} ≡ (*α* + *α*_{p})/*α*_{p} = *δ* + 1 is the slope-relative parameter [*δ* ≡ *α*/*α*_{p} = (*αg*′)/(2*U*_{0}*f*) is the slope parameter]. The dimensionless form of Eq. (A1) is, then,

subject to

Assuming an ansatz of the form $\varphi =\varphi \u02dc\u2061(y)ei\u2061(kx\u2212\sigma t)$, substituting the ansatz into Eqs. (A4) and (A5) yields the eigenvalue problem as followed (dropping tilde accents for clarity),

subject to

### APPENDIX B

#### Adapted Stone Model

The following derivation follows Stone (1966, 1970, 1971) but is modified to account for the presence of sloping bottom and surface. The coordinates are rotated to align with the sloping topography as Wenegrat et al. (2018), and the derivation is essentially equivalent to Wenegrat et al. (2018) but with a different coordinate orientation. Dimensionally, the equations describing a rotational fluid field with the Boussinesq approximation in the rotated coordinates are

subject to

where *θ* is the slope angle, $u*$ is the along-slope velocity, $\upsilon *$ is the across-slope velocity, and $w*$ is the slope-normal velocity, $p*$ is the pressure, and $b*=g\u2061(\rho 0\u2212\rho )\rho 0\u22121$ is the buoyancy (*ρ* and *ρ*_{0} are the seawater density and the reference, respectively). Rigid-lid boundary conditions are applied at the surface ($z*=H$) and bottom ($z*=0$).

Considering the horizontal velocity scale as *U*, the time scale as *f*^{−1}, the horizontal length scale as *U*/*f*, and the vertical length scale as *H*, the scaling relations about the variables in Eq. (B1) are, then,

subject to

where Ri = *N*^{2}*H*^{2}*U*^{−2} = *N*^{2}*f*^{2}*M*^{−4} is the Richardson number, *δ* = *αN*^{2}*M*^{−2} is the slope parameter, and *ε* = *fHU*^{−1} = *f*^{2}*M*^{−2} is the nonhydrostatic parameter.

Considering a background current flowing only in the along-slope direction, constrained by the thermal wind relation, the background state can be described as

where *M*^{2}/*f* is the thermal wind shear in the nonrotated coordinates and $b\xaf*$ will be derived in appendix D. The background state can be normalized as, then,

where *S*_{r} = *δ*_{r}Ri^{−1} is the slope-relative Burger number. Assuming *u*′, *υ*′, *w*′, *b*′, and *p*′ are the small perturbations from the background state, the equations governing the perturbed state can be linearized by neglecting the product of small terms as (dropping the primes for clarity and neglecting the viscosity and diffusion terms)

Assuming an ansatz of the form $\varphi =\varphi \u02dc\u2061(z)ei\u2061(kx+\lambda y\u2212\sigma t)$, substituting the ansatz and the background state into Eq. (B8) yields the eigenvalue problem as followed (dropping tilde accents for clarity),

subject to

### APPENDIX C

#### Rossby Wave Interactions in the Adapted Sakai Model

The following derivation follows Sakai (1989) but is modified to account for the presence of sloping bottom and top. The details about the interaction theory and associated derivation can be found in section 4 and appendix A of Sakai (1989). We will only focus on the interactions between Rossby waves. In the adapted Sakai model, the physical wave coordinates consisting of the Rossby waves in the upper and lower layers are, then,

where $k*$ is the dimensional along-slope wavenumber, $ln*=n\pi /2Ymax*$ (*n* = 1, 2, 3, …) is the dimensional across-slope wavenumber, $\sigma *$ is the dimensional wave frequency, $Rd=\u2061(1/2)g\u2032H0/f$ is the deformation radius, and *δ*_{r} = 1 + [(*αg*′)/(2*U*_{0}*f*)] is the slope-relative parameter. Assuming an ansatz of the form $\varphi =\varphi \u02dc\u2061(y*)ei\u2061(k*x*\u2212\sigma *t*)$, one eigenmode in the mathematical coordinates can be projected onto the Rossby wave coordinates as follows:

where *A*_{n} and *B*_{n} are the magnitudes in the physical wave coordinates, $d1n2\u2261\u222bE1nT\u22c5e1n\u2009dy$ (**E**_{1n} is the complex conjugate of the adjoint vector $[H1*u1n*,H1*\upsilon 1n*,\u2061(1/g\u2032)p1n*]$) and the same for *d*_{2n}. The interactions between Rossby waves can be described by

where $\epsilon n*={1/[2Rd2\u2061(k*2+ln*2)+1]}$ is the interaction coefficient (invariant form with the presence of sloping bottom and top). Consequently, the resonance rate $R*$ can be obtained by eliminating *A*_{n} and *B*_{n} in Eq. (C3):

which can be reduced to the flat-bottom case, $R*=Imag(k*U01\u2212{2/[Rd2\u2061(k*2+ln*2)+1]}}$ [Sakai 1989, Eq. (28)], by setting *δ*_{r} = 1.

### APPENDIX D

#### Energetics in the Rotated Coordinates

The following energetics closely follows Wenegrat et al. (2018) but with a different coordinate orientation. The coordinates are rotated to align with the sloping topography. The relation between the nonrotated and rotated coordinates is

where the tildes denote the nonrotated coordinates. The background buoyancy has a constant vertical gradient $\u2202b\xaf/\u2202z\u02dc=N2$ and a constant horizontal gradient $\u2202b\xaf/\u2202y\u02dc=\u2212M2$. The cross-slope and slope-normal gradients are, then,

So, the background buoyancy in the rotated coordinates could be set as $b\xaf=N2\u2061(cos\u2061\theta z\u2212sin\u2061\theta y)\u2212M2\u2061(sin\u2061\theta z+cos\u2061\theta y)$. The background along-slope velocity is in the thermal wind balance with the background buoyancy, that is, $u\xaf=(M2/f)\u2061(z/cos\u2061\theta )$.

By linearizing the dimensional governing equations, Eq. (B1), we can get the dimensional equations describing the perturbations as below (dropping asterisks for clarity)

where the bars denote the background variables, the primes denote the perturbation variables, and $d\xaf/d\xaft=(\u2202/\u2202t)+u\xaf\u22c5\u2207$ is the material derivative with the background velocity.

By multiplying *b*′/*N*^{2} to the buoyancy equation in Eq. (D3) and taking the Reynolds averaging, the eddy potential energy (EPE) equation can be obtained as

where HBF_{c} is the horizontal buoyancy flux (HBF) contributed by cross-slope motion, HBF_{n} is the HBF contributed by slope-normal motion, VBF_{c} is the vertical buoyancy flux (VBF) contributed by cross-slope motion, VBF_{n} is the VBF contributed by slope-normal motion, DPE is the dissipation of EPE, and RPE is the redistribution of EPE. Baroclinic instabilities get energy from the available potential energy of the background front. The energy transfer from mean available potential energy to EPE is through the HBF terms, which requires $\upsilon \u2032b\u2032\xaf>0$ given that HBF_{c} ≫ HBF_{n}.

By multiplying $u\xaf$ to the momentum equations in Eq. (D3) and taking the Reynolds averaging, the eddy kinetic energy (EKE) equation can be obtained as (the indicial notation is used here)

where SP is the shear production, DKE is the dissipation of EKE, RKE is the redistribution of EKE, and $sij=(1/2)[\u2061(\u2202ui\u2032/\u2202xj)/\u2061(\u2202uj\u2032/\u2202xi)]$ is the strain tensor. The VBF terms are the energy transfer between EKE and EPE. Growing instabilities require $w\u2032b\u2032\xaf>0$; VBF_{n} is an EKE source. Since $\upsilon \u2032b\u2032\xaf>0$, VBF_{c} is an EKE sink; it represents the energy transfer back to EPE due to overcoming the cross-slope gravity.

## REFERENCES

*J. Atmos. Sci.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Nat. Commun.*

*Tellus*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*J. Fluid Mech.*

*Annu. Rev. Fluid Mech.*

*J. Fluid Mech.*

*J. Geophys. Res. Oceans*

*J. Atmos. Sci.*

*Geophysical Fluid Dynamics*

*J. Mar. Res.*

*Tellus*

*J. Fluid Mech.*

*J. Atmos. Sci.*

*J. Atmos. Sci.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*Mar. Pollut. Bull.*

*Cont. Shelf Res.*

*J. Phys. Oceanogr.*

## Footnotes

Denotes content that is immediately available upon publication as open access.