A barotropic model of eddy saturation

TL;DR

This paper demonstrates that eddy saturation, previously thought to require baroclinicity, can also occur in purely barotropic flows over topography, driven by the structure of geostrophic contours.

Contribution

It shows that eddy saturation is a fundamental property of barotropic dynamics and depends on the openness of geostrophic contours, expanding understanding beyond baroclinic models.

Findings

Eddy saturation occurs in barotropic flows with open geostrophic contours.

Closed geostrophic contours do not support eddy-saturated states.

Transient eddies result from barotropic-topographic instability.

Abstract

"Eddy saturation" refers to a regime in which the total volume transport of an oceanic current is insensitive to the wind stress strength. Baroclinicity is currently believed to be key to the development of an eddy-saturated state. In this paper, it is shown that eddy saturation can also occur in a purely barotropic flow over topography, without baroclinicity. Thus, eddy saturation is a fundamental property of barotropic dynamics above topography. It is demonstrated that the main factor controlling the appearance or not of eddy-saturated states in the barotropic setting is the structure of geostrophic contours, that is the contours of $f/H$ of the ratio of the Coriolis parameter to the ocean's depth. Eddy-saturated states occur when the geostrophic contours are open, that is when the geostrophic contours span the whole zonal extent of the domain. This minimal requirement for eddy-saturated states is demonstrated using numerical integrations of a single-layer quasi-geostrophic flow over two different topographies characterized by either open or closed geostrophic contours with parameter values loosely inspired by the Southern Ocean. In this setting, transient eddies are produced through a barotropic-topographic instability that occurs due tot the interaction of the large-scale zonal flow with the topography. Through the study of this barotropic-topographic instability insight is gained on how eddy-saturated states are established.

Equations

$$f_0+\beta y+\eta +{\nabla }^2\psi .$$

$${\partial }_t{\nabla }^2\psi +\mathsf{J}(\psi -Uy,{\nabla }^2\psi +\eta +\beta y)=-\mathsf{D}{\nabla }^2\psi ,$$

$${\partial }_{t}U=F-\mu U-⟨\psi {\partial }_x\eta ⟩,$$

$${\eta }^{\parallel }$$

$${\eta }^{\times }$$

$${\ell }_{\eta }^{\parallel }\approx {\ell }_{\eta }^{\times }=0.071L=55.03\text{km}.$$

$$\psi (x,y,t)=\bar{\psi }(x,y)+{\psi }^{\prime }(x,y,t),$$

$$\bar{U}⟨\bar{\psi }{\partial }_x\eta ⟩+\underset{\text{negligible}}{\underbrace{\overline{U^{\prime }⟨{\psi }^{\prime }{\partial }_x\eta ⟩}}}=2\mu \underset{\stackrel{\mathrm{def}}{=}\text{EKE}}{\underbrace{⟨\frac{1}{2}\overline{|\nabla \psi {|}^2}⟩}}+\begin{array}{c}\text{small hyperviscous}\\ \text{dissipation}\end{array}.$$

$$\mathbf{E}\stackrel{\mathrm{def}}{=}\left(\overline{(U^{\prime }+u^{\prime }){\zeta }^{\prime }},\overline{v^{\prime }{\zeta }^{\prime }}\right),$$

$$\mathrm{corr}({\eta }^{\parallel },\nabla \cdot \mathbf{E})\approx \mathrm{corr}({\partial }_x{\eta }^{\parallel },\nabla \cdot \mathbf{E})\approx 0.69,$$

$$\mathrm{corr}({\eta }^{\times },\nabla \cdot \mathbf{E})\approx 0.01\text{and}\mathrm{corr}({\partial }_x{\eta }^{\times },\nabla \cdot \mathbf{E})=-0.07.$$

$$\psi =\left[S(t)\sin (mx)+C(t)\cos (mx)\right]/m.$$

$$\frac{\mathrm{d}C}{\mathrm{d}t}$$

$$\frac{\mathrm{d}S}{\mathrm{d}t}$$

$$\frac{\mathrm{d}U}{\mathrm{d}t}$$

$$U=U^e,\psi =\underset{\stackrel{\mathrm{def}}{=}{\psi }^e}{\underbrace{\left[S^e\sin (mx)+C^e\cos (mx)\right]/m}}.$$

$$U(t)$$

$$\psi (x,y,t)$$

$$+[s(t)\sin (mx)+c(t)\cos (mx)]/m.$$

$$\frac{\mathrm{d}c}{\mathrm{d}t}$$

$$\frac{\mathrm{d}s}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\upsilon }{\mathrm{d}t}$$

$${\partial }_t{\nabla }^2\varphi =-\left[(U^e{\nabla }^2+\beta ){\partial }_x+\mu {\nabla }^2\right]\varphi -\{[S^e\cos (mx)$$

$$-C^e\sin (mx)]({\nabla }^2+m^2)+m{\eta }_0\sin (mx)\}{\partial }_y\varphi .$$

$$\varphi (x,y,t)={\mathrm{e}}^{\lambda t}{\mathrm{e}}^{\mathrm{i}{n}_{y}y}{\tilde{\varphi }}_{n_x}(x),$$

$${\tilde{\varphi }}_{n_x}\stackrel{\mathrm{def}}{=}{\mathrm{e}}^{\mathrm{i}{n}_{x}x}\sum _{M=-\infty }^{+\infty }{c}_M{\mathrm{e}}^{\mathrm{i}Mmx}.$$

$$\pm n_x,\pm (n_x\pm m),\pm (n_x\pm 2m),\dots .$$

$$\lambda {\nabla }_{n_y}^2{\tilde{\varphi }}_{n_x}={\text{L}}_{n_y}{\tilde{\varphi }}_{n_x},$$

$${\text{L}}_{n_y}$$

$$-C^e\sin (mx)]({\nabla }_{n_y}^2+m^2)+m{\eta }_0\sin (mx)\}.$$

$$\text{A}\left(\lambda \right)\left[\begin{array}{c}{c}_{-M_{\max }}\\ ⋮\\ c_{M_{\max }}\end{array}\right]=0,$$

Full text

††journal: jpo

1 Introduction

The Southern Ocean, and in particular the Antarctic Circumpolar Current (ACC), are key elements of the climate system. The ACC is driven by a combination of strong westerly winds and buoyancy forcing. Straub (1993) advanced the remarkable hypothesis that the equilibrated ACC zonal transport should be insensitive to the strength of the wind stress forcing. This insensitivity was later verified in eddy-resolving ocean models of the Southern Ocean and is now referred to as eddy saturation (Hallberg and Gnanadesikan 2001; Tansley and Marshall 2001; Hallberg and Gnanadesikan 2006; Hogg et al. 2008; Nadeau and Straub 2009; Farneti et al. 2010; Hogg 2010; Meredith et al. 2012; Nadeau and Straub 2012; Morisson and Hogg 2013; Munday et al. 2013; Nadeau et al. 2013; Abernathey and Cessi 2014; Farneti et al. 2015; Nadeau and Ferrari 2015; Howard et al. 2015; Marshall et al. 2017). Some indications of eddy saturation are seen in observations (Böning et al. 2008; Firing et al. 2011).

There is evidence that the strength of the westerly winds over the Southern Ocean, which force the ACC, are increasing (Thompson and Solomon 2002; Marshall 2003; Yang et al. 2007; Swart and Fyfe 2012). Recently, Hogg et al. (2015) using satellite altimetry data identified that along with the strengthening of the westerlies comes also a linear trend of the Southern Ocean surface eddy kinetic energy (EKE) while the ACC zonal transport remains insensitive or even has decreased. This way, Hogg et al. (2015) concluded that the ACC is in an eddy-saturated state. Given the strengthening of the Southern Ocean westerly winds over the last decades, and the potential enhanced strengthening under global warming forcing (Bracegirdle et al. 2013), the question that naturally arises is how will the ACC transport respond? Thus, understanding the mechanisms behind eddy saturation is particularly relevant. This paper attempts to shed some more insight in the mechanism underlying eddy saturation.

Initially the explanation for eddy saturation given by Straub (1993) relied on baroclinic processes and on the existence channel walls. In the following detailed models of Nadeau and Straub (2009, 2012), Nadeau et al. (2013), and Nadeau and Ferrari (2015) the arguments for explaining eddy saturation were barotropic in heart. Specifically, Nadeau and Ferrari (2015) argued that the circulation can be decomposed to a circumpolar mode and a gyre mode (with the latter not contributing to the total transport). Nadeau and Ferrari (2015) showed that the wind stress curl spins barotropic gyres and, furthermore, that an increase of wind stress spins up the gyres while leaves the circumpolar mode intact; thus the transport remains insensitive. Still, though, all those explanations relied in baroclinic instability for producing transient eddies to transfer the momentum down from the surface of the ocean to the bottom. On the other hand, Marshall et al. (2017) and Mak et al. (2017) recently showed that eddy saturation can emerge as a result of an eddy flux parametrization that was introduced by Marshall et al. (2012). In agreement with Straub (1993), Marshall et al. (2017) also relate the production of EKE to the vertical shear of the zonal mean flow. Again, baroclinic instability is identified as the main source of EKE. In this paper, we show that eddy saturation can be observed without baroclinic eddies, without channel walls and without any wind stress curl, i.e. without any gyres.

According to Johnson and Bryden (1989) different density layers are coupled via interfacial form stress that transfers momentum downwards from the sea-surface to the bottom. But at the bottom it is topographic form stress that transfers momentum from the ocean to the solid earth (Munk and Palmén 1951). Note that only the standing eddies result in time-mean topographic form stress. Thus, it seems reasonable that topography should play a dominant role in understanding what sets up the total vertically-integrated transport. There is a consensus that topography acts as the main sink in the zonal momentum balance (in agreement with Munk and Palmén (1951)). But, can the topography also have an active role in setting up the momentum balance, e.g., by shaping the standing eddy field and its associated form stress? Abernathey and Cessi (2014) argued in favor of the latter: they showed that isolated topographic features result in localized absolute baroclinic instability (i.e. baroclinic eddy growth in situ above the topographic features) and an associated almost-barotropic standing wave pattern. Transient eddies are suppressed away from the topography. Furthermore, relative to the flat-bottom case, the thermocline is shallower and isopycnal slopes are smaller. Thus, the usual arguments assuming that the isopycnal slopes are so steep as to be marginal with respect to the flat-bottomed baroclinic instability cannot explain the ACC equilibration (and thus neither eddy saturation) in this model configuration with localized topography. In addition, Thompson and Naveira Garabato (2014), and more recently Youngs et al. (2017), further emphasized the role of the standing eddies (or standing meanders) in setting up the momentum balance and ACC transport. In this paper, we emphasize the role of the standing eddies and their form stress in determining transport in an eddy-saturated regime within a barotropic setting.

Lately, there has been increasing evidence arguing for the importance of the barotropic processes in determining the ACC transport. Ward and Hogg (2011) studied the ACC equilibration using a multi-layer primitive equation wind-driven model with an ACC-type configuration starting from rest. Following turn-on of the wind a strong barotropic current forms within several days that is able to transfer most of the imparted momentum to the bottom; only after several years does the momentum start being transferred vertically via interfacial form stress. The fast response is that a bottom pressure signal arises a few days after turn-on, and the associated topographic form stress couples the ocean to the solid earth. Subsequently, for about ten years, and contrary to the statistical equilibrium scenario described by Johnson and Bryden (1989), interfacial from stress transfers momentum vertically from the bottom upwards. At equilibrium both eastward momentum is transfered from surface downwards and also westward momentum from bottom upwards. On the other hand, studies using in situ velocity measurements, satellite altimetry, and output from the Southern Ocean State Estimate (SOSE) also argue in favor of the importance of the bottom-velocity component of the ACC transport, which comes about from the barotropic component of the flow111Traditionally, the baroclinic component of the flow is obtained through the thermal-wind relationship after assuming zero flow at the bottom. Thus, any transport due to the bottom flow is not thought as part of the “baroclinic component” of the ACC transport. (Rintoul et al. 2014; Peña Molino et al. 2014; Masich et al. 2015; Donohue et al. 2016). For example, using measurement from the cDrake experiment Donohue et al. (2016) estimated that the bottom-velocity component of the ACC transport accounts for about 25% of the total transport.

Constantinou and Young (2017) discussed the role of standing eddies using a simple barotropic model forced by an imposed steady wind stress and retarded by a combination of bottom drag and topographic form stress (Hart 1979; Davey 1980; Holloway 1987; Carnevale and Frederiksen 1987). Constantinou and Young (2017) used a random monoscale topography and argued that a critical requirement for eddy saturation is that the ratio of the planetary potential vorticity (PV) to the topographic PV is large enough so that the geostrophic contours (that is the contours of the planetary plus the topographic contribution to PV) are “open” in the zonal direction. Here, we demonstrate that what matters for eddy saturation is not the actual value of the ratio of the planetary potential vorticity (PV) over the topographic PV, but rather the structure of the geostrophic contours.

The main goal of the present paper is to provide insight on how eddy saturation is established in this barotropic setting. We do that by comparing numerical results using two simple sinusoidal topographies. We show that this barotropic model without baroclinic instability can exhibit impressive eddy saturation provided that the geostrophic contours are open (see section 3). We do not claim here that baroclinic processes are not important for setting up the momentum balance in the Southern Ocean. Instead, we emphasize the role of barotropic dynamics and the fact that we can still observe eddy saturation without baroclinicity. In this barotropic model, transient eddies arise as an instability due to the interaction of the large-scale zonal flow with the topography. The simple topography used here allows us to study in detail the barotropic–topographic instability that gives rise to transient eddies and thus provide insight on how eddy-saturated states appear (see section 5). We show, in this way, that topography does not only have a passive role in setting the zonal momentum balance acting as the main sink of zonal momentum, but it can also have an active role by producing transient eddies that shape the standing eddy field and its associated form stress.

2 Setup

Consider the quasi-geostrophic dynamics of a barotropic fluid of depth $H-h(x,y)$ on a beta-plane. The fluid velocity consisting of a large-scale domain-averaged zonal flow, $U(t)$, along the zonal $x$ direction plus smaller-scale eddies with velocity $(u,v)$. The eddy component of the flow is derived via an eddy streamfunction $\psi(x,y,t)$ through $(u,v)=(-\partial_{y}\psi,\partial_{x}\psi)$; the streamfunction of the total flow is $-U(t)y+\psi(x,y,t)$. The relative vorticity of the flow is $\nabla^{2}\psi$, where $\nabla^{2}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\partial_{x}^{2}+\partial_{y}^{2}$, and the quasi-geostrophic potential vorticity (QGPV) of the flow is

[TABLE]

In (1), $f_{0}$ is the Coriolis parameter at the center of the domain, $\beta$ is the meridional PV gradient and $\eta(x,y)\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}f_{0}h(x,y)/H$ is the topographic contribution to QGPV or simply the topographic PV. The QGPV and the large-scale flow evolve through:

[TABLE]

(Hart 1979; Davey 1980; Holloway 1987; Carnevale and Frederiksen 1987). In (2a) $\mathsf{J}$ is the Jacobian, $\mathsf{J}(a,b)\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}(\partial_{x}a)(\partial_{y}b)-(\partial_{y}a)(\partial_{x}b)$, and $\mathsf{D}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\mu+\nu_{4}\nabla^{4}$ is the dissipation operator, which includes linear Ekman drag $\mu$ and hyperviscosity $\nu_{4}$ used for numerical stability. On the right of (2b), $\langle\psi\partial_{x}\eta\rangle$ is the topographic form stress, with angle brackets denoting an average over the $(x,y)$-domain. The large-scale flow $U(t)$ in (2b) is forced by the constant $F=\tau/(\varrho_{0}H)$, where $\tau$ is the uniform surface zonal wind stress and $\varrho_{0}$ is the reference density of the fluid, while is being retarded by a combination of bottom drag and topographic form stress. Being constant through the domain, the wind stress has no curl. The domain is periodic in both the zonal and the meridional directions, with size $2\pi L\times 2\pi L$. Following Bretherton and Karweit (1975), we have in mind a mid-ocean region that is smaller than ocean basins but much larger than the length scale of ocean macroturbulence. The role of hyperviscosity is limited only to the removal of small-scale enstrophy. Thus the hyperviscosity has a very small effect on larger scales and energy dissipation is mainly due to drag $\mu$.

The model formulated in (2) is the simplest model that can be used to investigate beta-plane turbulence above topography driven by a large-scale zonal wind stress applied at the surface of the fluid. It has been used in the past for studying the interaction of zonal flows with topography (Hart 1979; Davey 1980; Holloway 1987; Carnevale and Frederiksen 1987) and recently by Constantinou and Young (2017) for studying the geostrophic flow regimes above random monoscale topography.

Inspired by the Southern Ocean we take $L=775\;\textrm{km}$, $H=4\,\textrm{km}$, $\varrho_{0}=1035\;\textrm{kg}\;\textrm{m}^{-3}$, $f_{0}=-1.26\times 10^{-4}\;\textrm{s}^{-1}$ and $\beta=1.14\times 10^{-11}\;\textrm{m}^{-1}\textrm{s}^{-1}$. Also, we take $h_{\rm rms}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\sqrt{\langle h^{2}\rangle}=200\;\textrm{m}$, which implies that $\eta_{\mathrm{rms}}=6.3\times 10^{-6}\;\textrm{s}^{-1}$. For Ekman drag we use $\mu=6.3\times 10^{-8}\;\textrm{s}^{-1}\approx(180\;\textrm{day})^{-1}$ (Arbic and Flierl 2004).

Results are reported using two sinusoidal topographies:

[TABLE]

Both topographies in (3) are monoscale, that is they are characterized by a single length-scale \ell_{\eta}=\eta_{\mathrm{rms}}\big{/}|\bm{\nabla}\eta|_{\rm rms}, which is:

[TABLE]

Thus, for both topographies the ratio of planetary PV to topographic PV is $\beta/|\bm{\nabla}\eta|_{\rm rms}=0.1$.

An important factor controlling the behavior of the flow is the geostrophic contours, that is the level sets of $\beta y+\eta(x,y)$. Figure 1 shows the structure of geostrophic contours for the two topographies $\eta^{\parallel}$ and $\eta^{\times}$ using the parameters given above. Here, we distinguish between “open” and “closed” geostrophic contours: open geostrophic contours span the full zonal extent of the domain; see, for example, the two contours marked with thick black curves in each panel of figure 1. For any nonzero $\beta$, all geostrophic contours for topography $\eta^{\parallel}$ are open while for topography $\eta^{\times}$ there exist both closed as well as open geostrophic contours. However, for topography $\eta^{\times}$ open geostrophic contours, such as the one shown in figure 1(b), are only found in the vicinity of narrow channels that span the horizontal extent of the domain snaking around local maxima and minima of $\beta y+\eta^{\times}$. We say, therefore, that $\eta^{\times}$ is characterized by closed geostrophic contours.

It is useful to decompose the eddy streamfunction $\psi$ into time-mean “standing eddies” with streamfunction $\bar{\psi}$, and residual “transient eddies” $\psi^{\prime}$, so that:

[TABLE]

where the time-mean is $\bar{\psi}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\lim_{T\to\infty}\frac{1}{T}\int_{t_{0}}^{t_{0}+T}\!\!\!\psi(x,y,t^{\prime})\,{\rm d}t^{\prime}$, with $t_{0}$ the time the flow needs to reach a statistically steady state. Similarly, all fields can be decomposed into time-mean and transient components e.g., $U(t)=\bar{U}+U^{\prime}(t)$. We are interested on how the time-mean large-scale flow $\bar{U}$, which is directly related to the time-mean total zonal transport, depends on wind stress forcing $F$.

All solutions presented in this paper employ $512^{2}$ grid points; this resolution allows about 6 grid points within $\ell_{\eta}$. We have verified that results remain unchanged with double the resolution. The hyperviscocity coefficient is set to $\nu_{4}=2.27\times 10^{9}\;\textrm{m}^{4}\,\textrm{s}^{-1}$. System (2) is evolved using the Exponential Time Differencing 4th-order Runge-Kutta time-stepping scheme (Cox and Matthews 2002; Kassam and Trefethen 2005).222A Matlab code used for solving (2) is available at the github repository: https://github.com/navidcy/QG_ACC_1layer. Time-averaged quantities are calculated by averaging the fields over the interval $30\leq\mu t\leq 60$.

3 Results

3.1 Variation of the time-mean large-scale flow with wind stress forcing

Figure 2 shows how the time-mean large-scale flow $\bar{U}$ varies with wind stress forcing $F$ for the two topographies. A generic feature of model (2), for both topographies used here (as well as other more complex and/or multiscale topographies) is the following: as wind stress strength varies from weak to stronger values the flow transitions from a regime with relatively small time-mean large-scale flow $\bar{U}$ to a regime with large $\bar{U}$. Specifically, for very high wind stress forcing values the time-mean large-scale flow becomes $\bar{U}\approx F/\mu$, which is the value for $\bar{U}$ in the flat-bottom case (cf. (2b) with $\eta=0$). The flat-bottom transport value is marked with the dashed-dotted line in both panels of figure 2. These two flow regimes are referred to as the lower-branch and the upper-branch (Davey 1980), and they are indicated in figure 2. In both branches, the flow reaches a steady state without any transients after an initial adjustment on the time scale of $O(\mu^{-1})$. For intermediate wind stress forcing values the flow develops transients and becomes turbulent.

In the remainder of this subsection we describe in detail the qualitative features of the flow on the lower and upper branches, as well as the transition from one branch to the other for the two topographies.

For weak wind stress forcing values (i.e. on the lower branch) the equilibrated solutions of (2) are time-independent without any transient eddies. These steady states have a large-scale flow $U$ and an associated stationary eddy flow field $\psi$ that both vary linearly with wind stress forcing $F$. As a result, EKE varies quadratically with wind stress, i.e., as $F^{2}$. Figure 3(a) shows how EKE varies with wind stress forcing for topography $\eta^{\parallel}$ and confirms the $F^{2}$ dependence on the lower branch (similar behavior is found for topography $\eta^{\times}$; not shown). As the wind stress increases beyond a certain value the lower-branch steady states undergo an instability and develop transients. Even though the topography used here imposes a length-scale $\ell_{\eta}$, the transient eddies can have scales much smaller than $\ell_{\eta}$. The onset of this instability is roughly at $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\approx 7\times 10^{-3}$ for both topographies. (The stability of the lower-branch states for topography $\eta^{\parallel}$ is studied in detail in section 5.)

For wind stress forcing values beyond the threshold of the lower-branch instability the flow above the two topographies are qualitatively very different: topography $\eta^{\parallel}$ shows eddy saturation while topography $\eta^{\times}$ does not (see figure 2).

For topography $\eta^{\parallel}$, as wind stress continues to increase beyond the onset of the transient eddy instability, the time-mean large-scale flow $\bar{U}$ ceases to grow linearly with $F$; instead $\bar{U}$ grows at a much slower rate. In fact, $\bar{U}$ shows only a very weak dependence on wind stress strength in the shaded region of figure 2(a): $\bar{U}$ increases only about 4-fold in the course of a 60-fold wind stress increase from $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})=7.0\times 10^{-3}$ up to $4.1\times 10^{-1}$. This regime is identified as the eddy saturation regime.333Note, that Nadeau and Ferrari (2015) found a similar weak dependence of the transport to wind stress strength in a series of experiments they performed without wind stress curl, i.e. similar to the forcing we used here. In their experiments with wind stress curl the eddy saturation was more pronounced.

The eddy saturation regime in figure 2(a) terminates at $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\approx 4.1\times 10^{-1}$ by a discontinuous transition to the upper-branch, marked with a dashed vertical line. The upper branch steady solutions are characterized by much larger values of time-mean large-scale flow $\bar{U}$ compared to ones in the eddy saturation regime or on the lower-branch. This discontinuous transition has been coined drag crisis (Constantinou and Young 2017). The location of the drag crisis exhibits hysteresis (Charney and DeVore 1979). As we increase wind stress forcing the drag crisis and the transition to the upper branch occurs at $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\approx 4.1\times 10^{-1}$. However, if we initiate our model with the upper-branch solutions and start decreasing the wind stress forcing, we can find solutions on the upper branch down to $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\approx 2.44\times 10^{-2}$: see figure 2(a). We note that to obtain the upper-branch solutions we have to be very delicate in initiating the system appropriately. Thus, although the upper-branch solutions shown in figure 2(a) are linearly stable, their basin of attraction is smaller than those of the co-existing eddy-saturated solutions. The extent of the region of multiple stable equilibria is explained by studying the stability of the upper-branch solution (see section 5).

On the other hand, the case with topography with closed geostrophic contours $\eta^{\times}$ in figure 2(b), does not exhibit any drag crisis nor any multiple equilibria. Also, the simulations using topography $\eta^{\times}$ do not show any sign of eddy saturation; there is only a slight, barely noticeable, decrease from the linear growth of the time-mean large-scale flow $\bar{U}$ in the range $8.0\times 10^{-4}\leq F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\leq 3.0\times 10^{-3}$, and for those wind stress values the flow does not have any transient eddies. Transient eddies appear for wind stress forcing $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})>4.2\times 10^{-3}$. For wind stress forcing values beyond the onset of transient eddies the flow transitions from the lower to the upper branch in a continuous manner between $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})\approx 8.0\times 10^{-3}$ and $4.2\times 10^{-2}$ and at a rate much faster than linear.

3.2 The flow regimes

As described in the previous subsection, we distinguish three qualitatively different flow regimes for each topography. There exist, for both topographies, a lower-branch flow regime for weak wind stress forcing and an upper-branch flow regime for strong wind stress forcing. These flow regimes consist of steady flows without any transient eddies. In between the lower- and upper-branch flow regimes there exists a regime in which the flow has a transient component: the “eddy saturation regime” for topography $\eta^{\parallel}$ case and the “transient regime” for topography $\eta^{\times}$. These intermediate regimes are characterized by flows that feature strong transients and are turbulent, especially the eddy saturated case. Moreover, in these regimes the flow shows energy in a wide range of spatial scales. Figures 5 and 5 show the energy spectra and a snapshot of the flow fields for a typical representative case of each of the three flow regimes (the cases presented in figures 5 and 5 are the ones marked with a square, a triangle and a diamond in panels (a) and (b) of figure 2).

The eddy field $\psi$ in the eddy saturated regime is characterized as two-dimensional turbulence. The wind stress $F$ directly drives the large-scale flow $U$ in (2b) and, in turn, $U$ drives the eddy field through the term $U\partial_{x}\eta$ in (2a). This leads to a forward transfer of energy from the largest possible scale to the eddies on the length scale of the topography $\ell_{\eta}$ and then nonlinearity transfers energy from length scale $\ell_{\eta}$ to all other scales. This turbulent regime is anisotropic; for example, for the case marked with a triangle in figure 2(a) and also shown in figure 5(c), the transient eddy velocities are related by $\overline{\langle u^{\prime 2}\rangle}\approx 0.57\overline{\langle v^{\prime 2}\rangle}$. The turbulent eddy flow exhibits regimes that resemble the energy and enstrophy inertial ranges. For length scales smaller than $\ell_{\eta}$, the energy spectrum even shows a slope close to $-3$ resembling the enstrophy inertial range as predicted by homogeneous isotropic turbulence arguments. Characteristically, for the open geostrophic contours case (figure 5) the eddy saturated solution shows energy content in spatial scales other than $\ell_{\eta}$ that is at least eight orders of magnitude larger compared to the steady lower- and upper-branch solutions. Similar behavior is seen for the closed geostrophic contours case of figure 5. The difference between these two cases is that for topography $\eta^{\times}$ the lower- and upper-branch solutions are not monochromatic, i.e., they do not show energy content only at scale $\ell_{\eta}$: the nonlinearity due to non-vanishing Jacobian term $\mathsf{J}(\psi,\eta^{\times})$ induces energy to cascade to scales smaller than scale $\ell_{\eta}$.

3.3 Further characteristics of the eddy saturation regime for topography $\eta^{\parallel}$

It has been noted that in the eddy-saturation regime, even though the total ACC transport varies very little with wind stress, the domain-averaged EKE is approximately linearly related to the wind stress (e.g., see Hallberg and Gnanadesikan (2001, 2006)). Figure 3(a) shows how this is reflected in this barotropic model: indeed, in the eddy saturation regime the EKE varies with wind stress forcing at a rate much slower than quadratic. It is also apparent from figure 3(a) that strong transient eddies characterize the eddy saturation regime; this can be seen by the diminishing of the standing wave EKE (sEKE), that is the EKE that results from the standing waves alone.

The barotropic model (2) predicts that for a fully eddy-saturated state, i.e. when $\bar{U}$ does not vary at all with $F$, the EKE should vary linearly with wind stress forcing. If the flow is completely saturated, then from (2b) we get that the form stress $\langle\bar{\psi}\partial_{x}\eta\rangle$ varies linearly with $F$. With that in hand, we can understand why EKE varies linearly with wind stress from the energy power integral, i.e. \overline{\langle-\psi\times\text{\eqref{eq:z_NL}}\rangle} gives:

[TABLE]

If we consider the energy transfers amongst the various flow components, we deduce that $\overline{U^{\prime}\langle\psi^{\prime}\partial_{x}\eta\rangle}$ should be less than $\bar{U}\langle\bar{\psi}\partial_{x}\eta\rangle$ (Constantinou and Young 2017). Moreover, the numerical solutions of (2) suggest that in the eddy saturation regime the ratio of $\overline{U^{\prime}\langle\psi^{\prime}\partial_{x}\eta\rangle}$ to $\bar{U}\langle\bar{\psi}\partial_{x}\eta\rangle$ is less than about $10^{-3}$ and, therefore, negligible. Thus, (4) implies that EKE is proportional to form stress $\langle\bar{\psi}\partial_{x}\eta\rangle$, i.e. EKE varies linearly with $F$.

Results both from quasi-geostrophic models (Hogg and Blundell 2006; Nadeau and Straub 2012; Nadeau and Ferrari 2015) and also primitive-equation models (Marshall et al. 2017) demonstrated the somehow counterintuitive fact that in the eddy saturation regime the total transport increases with increasing bottom drag. Hogg and Blundell (2006) suggested that this effect comes about because increasing the bottom drag decreases the strength of the transient eddies that are responsible for transferring momentum to the bottom where topographic form stress act; thus the strength of the zonal current increases. Nadeau and Straub (2012) on the other hand, argued that increase of the bottom drag damps the gyre circulation and its associated form stress; thus the zonal current increases to provide the drag needed to balance the momentum imparted by the wind stress.

In the barotropic model studied here, we also show that transport increases with increasing bottom drag in the eddy saturation regime. Figure 3(b) shows how the time-mean large-scale flow $\bar{U}$ varies when wind stress forcing is kept fixed at $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})=1.4\times 10^{-1}$ while bottom drag $\mu$ varies: we find that $\bar{U}$ increases with $\mu$ (up to a point). Here the forcing is a steady mean zonal wind stress with no curl and thus large gyres do not form. Therefore, our results line-up more with the explanation given by (Hogg and Blundell 2006) without, however, precluding the explanation by Nadeau and Straub (2012). Within model (2) large bottom drag damps the eddy field $\psi$, thereby decreasing the form stress. Hence, increased bottom drag strengthens the time-mean large-scale flow $\bar{U}$, which is necessary to balance the wind stress forcing $F$ in the eddy saturation regime (cf. (2b)). If we increase bottom drag $\mu$ beyond a certain threshold the flow becomes laminar and transitions from the eddy-saturation regime to the lower-branch solution (for the case in figure 3(b) this transition occurs at $\mu/\eta_{\mathrm{rms}}=3\times 10^{-2}$). On the lower-branch the time-mean large-scale flow $\bar{U}$ decreases with increasing $\mu$.

4 The effect of transient eddies to the standing wave

The eddy PV fluxes

[TABLE]

appears as forcing in the standing wave field $\bar{\psi}$. To elucidate how the zonal momentum balance is shaped by the transient eddies for the two different topographies, we compute the divergence of the eddy PV fluxes, $\bm{\nabla}\bm{\cdot}\bm{E}$. Figure 6 compares the structure of the eddy PV flux divergence (shading) with the actual standing wave $\bar{\psi}$ (contours), for the two example cases shown in figures 5(c) and 5(c). There is an apparent qualitative difference between the two cases. For the open geostrophic contours case in panel (a) the eddy PV fluxes are dominant over the whole domain and, to a large extent, they are aligned with the standing wave streamlines. This is not at all the case for the topography with closed geostrophic contours in panel (b). In model (2) we can verify, by considering the energy budgets among components (Constantinou and Young 2017), that the only source of transient eddy energy is the energy conversion from the standing-wave component through the term $\left\langle\bar{\psi}\,\bm{\nabla}\bm{\cdot}\bm{E}\right\rangle$. Thus, the stronger the transient eddies are the more stronger the correlation between eddy PV flux divergence and time-mean streamlines should be. In section 5 we show that for the case shown in figure 6(a) there is a strong transient-eddy instability that explains this almost-perfect alignment between eddy PV fluxes and standing wave streamlines.

More importantly, however, notice how for the case with $\cos{(14x/L)}$-topography shown in figure 6(a) the eddy PV fluxes are offset from the topography contours by 1/8-th wavelength (see figure 6(c)). This offset in the eddy forcing induces a standing flow that projects onto the $\sin{(14x/L)}$ rather than the $\cos{(14x/L)}$, i.e., projects onto $\partial_{x}\eta^{\parallel}$ rather than $\eta^{\parallel}$. Since the time-mean topographic form stress is $\langle\bar{\psi}\partial_{x}\eta\rangle$, it is apparent that only the $\sin$-component of the standing wave can produce non-zero topographic form stress. Thus, the mis-alignment between the topography and the transient-eddy forcing, induces a standing-eddy field which contributes to the form stress.

In conclusion, the transient eddies have the ability to shape the standing-wave field so that it projects more onto $\partial_{x}\eta$ and, therefore, enhancing the time-mean topographic form stress. This shaping of the standing-wave field by the transient eddies is prominent in figure 6(a) with topography $\eta^{\parallel}$ but not in the case with $\eta^{\times}$. Specifically, for the case shown in figure 6(a) the off-set by 1/8-th wavelength implies that $\bm{\nabla}\bm{\cdot}\bm{E}$ projects equally444A mis-alignment of 1/8-th wavelength implies that both correlations in (6a) are $1/\sqrt{2}$. to both $\eta$ and $\partial_{x}\eta$ (i.e. onto the $\cos$ and $\sin$). Indeed,

[TABLE]

5 The stability of the lower- and upper-branch solutions for topography $\eta^{\parallel}$

For topography of the form $\eta=\eta_{0}\cos{(mx)}$ (as is (3a)) we can understand both the transitions among the three flow regimes, as well as the role the transient eddies play on eddy saturation by studying a 3-dimensional invariant manifold of model (2).

It is easy to verify in this case that there exists a 3-dimensional invariant manifold by assuming that the eddy streamfunction is:

[TABLE]

Inserting (7) into (2) the Jacobian term $\mathsf{J}(\psi,\eta)$ vanishes identically, and thus (2) reduce to:

[TABLE]

(after also ignoring $\nu_{4}$ in $\mathsf{D}$). Steady equilibrium solutions of (8) can be found in the form:

[TABLE]

Both the lower- and upper-branch solutions for the case with topography $\eta^{\parallel}$ are exactly such steady solutions in the form of (9); see panels (b) and (d) of figure 5.

Starting with the work by Hart (1979), the stability of solutions (9) to perturbations that lie within this invariant manifold has been extensively studied (Hart 1979; Pedlosky 1981; Källén 1982; Rambaldi and Flierl 1983; Yoden 1985). Here, we study the stability of (9) to any general perturbation that may or may not lie within the invariant manifold. Charney and Flierl (1980) studied the stability of the inviscid and unforced version of (12), i.e., with $\mu=F=0$ and assuming the large-scale mean flow $U^{e}$ to be a given external parameter of the problem. Here, we perform the stability of the forced–dissipative problem. Thus, we use $U^{e}$ that results from a steady solutions of (2) (or, for $\eta^{\parallel}$, a steady solution of (8)). We investigate how this instability may elucidate why eddy saturation is observed.

5.1 Stability calculation

Assume perturbations about the steady state (9):

[TABLE]

In (10) we write the perturbation streamfunction as a sum of a perturbation that lies within the 3-state invariant manifold $[s\,\sin{(mx)}+c\,\cos{(mx)}]/m$ and a perturbation $\phi$ that does not project on the invariant manifold. That is, $\phi$ satisfies: $\langle\phi\,{\rm e}^{{\rm i}mx}\rangle=0$.

By inserting (10) into (2) we get the linearized equations for the perturbations. When projected on the invariant manifold the perturbation equations read:

[TABLE]

while the orthogonal complement of the perturbation equations is:

[TABLE]

Eigenanalysis of (11) determines the stability of steady state (9) to perturbations within the invariant manifold; eigenanalysis of (12) determines the stability of (9) to perturbations outside the manifold. Details for the eigenanalysis of (12) are provided in the Appendix.

5.2 Stability results

We can understand the various flow transitions for topography $\eta^{\parallel}$ by studying the stability of steady states (9). The region of stability of (9), both within and outside the low-dimensional manifold (i.e. with respect to both (11) and (12)), is marked with the thick semitransparent curve in figure 2(a). The onset of this barotropic–topographic instability explains: (i) the appearance of transient eddies as wind stress is increased and (ii) the termination of the upper-branch solution as wind stress is decreased.

Figure 7(a) shows the large-scale flow $U^{e}$ of the steady states (9) as a function of wind stress forcing $F$, together with their stability as predicted by (11) and (12). Also shown in figure 7(a) are the numerical results of (2) for comparison. It is clear that the transition from the lower branch to the turbulent regime, in which eddy saturation is observed, is triggered by the instability due to (12) rather than by (11). This instability of the lower branch first occurs at $F/(\ell_{\eta}\eta_{\mathrm{rms}}^{2})=7.15\times 10^{-3}$. Figure 7(c) shows the most unstable eigenfunction just beyond the marginal point of instability. This eigenfunction has small-scale meridional structure; the meridional wavenumber of the eigenfunction shown in figure 7(c) is $n_{y}=2.5m$. Thus, the instability of the lower branch due to (11) introduces $y$-dependence in the flow; see figure 7(c).

Next, we would like to investigate how this instability contributes to the occurrence of eddy saturation. What is interesting here is how the instability maximum growth rate varies with the large-scale flow $U^{e}$. Figure 7(b) shows that in the regime which we find eddy saturation the instability maximal growth rate due to (11) increases dramatically with $U$; growth rates increase by a factor of $10^{3}$ when $U$ increases only by a factor of 15. The transient eddy source, thus, increases significantly for only small changes in $U$. In this way, and by analogy with arguments for baroclinic eddy saturation, larger wind stress forcing implies the need for stronger eddies and this leads to eddy saturation. Furthermore, this suggests that the flow adjusts to a state close to marginal stability for this barotropic–topographic instability, similarly to the baroclinic marginal stability argument first invoked by Stone (1978).

The stability calculations presented in this section are based on that the topography does not depend on $y$. However, the same flow transitions occur for flows with open geostrophic contours above complex topographies. In that case, finding the lower- and upper-branch equilibrium steady states is much more painful since some of the Jacobian terms, for example $\mathsf{J}(\psi^{e},\eta)$, are not identically zero. Therefore, the lower- and upper-branch solutions do not lie within any low-dimensional manifold. Methods for obtaining such equilibria were developed by Tung and Rosenthal (1985).

6 Discussion

The results reported here demonstrate that eddy saturation can occur even without baroclinicity. A bare barotropic setting is capable of producing an eddy saturation regime in which the transport is insensitive to wind stress forcing. The model shows additionally some other features of eddy saturation that were previously observed in more elaborate ocean models: EKE varies close to linearly with wind stress forcing (instead of quadratically) and also transport increases with increasing bottom drag.

In this barotropic model with no baroclinic eddies, the flow relies on the barotropic–topographic instability to produce transient eddies. For the $y$-independent topography with open geostrophic contours $\eta^{\parallel}$, we extended the stability analysis of Hart (1979) to encompass any general flow perturbations. Thus, we managed to identify the various flow transitions shown in figure 2(a) with the onset of an instability. The stability calculations presented section 5 crucially depend on that the topography $\eta^{\parallel}$ does not depend on $y$. However, in the context of this barotropic model, we expect that the same flow transitions above more complex topography $\eta(x,y)$ with open geostrophic contours would result from a similar barotropic–topographic instability. More importantly, by studying the barotropic–topographic instability we showed that for the case with eddy saturation the transient eddy source (i.e. the barotropic–topographic instability) had a particular property. The growth rate of the instability increases dramatically with small changes of the large-scale flow. This explains why eddy saturation occurs: an increase of the wind stress strength requires the generation of more transient eddies to carry this excess of momentum to the solid Earth through form stress; but only a slight change in the transport can give rise to the transient eddies needed to balance the momentum. This is analogous with arguments for baroclinic eddy saturation.

The factor which controls the appearance or not of eddy saturated states in a barotropic setting is whether or not the geostrophic contours are open. This was demonstrated here using two simple sinusoidal topographies that bare this distinction. However, the bare existence of the eddy saturation in this simple setting depends on the presence of a transient-eddy source (i.e. of an instability). The range of wind stress forcing values that eddy saturation could potentially occur is limited by the range of wind stress forcing values for which the barotropic–topographic instability occurs. And the latter crucially depends on the height of the topography. For the flat-bottom case with $h=0$ all geostrophic contours are open but there is no eddy saturation since there is no instability. Figure 8 shows the steady states (9) together with their stability (similarly to figure 7(a)) for various topography heights. The range of wind stress forcing values that we have instability scales roughly as $h_{\rm rms}^{2}$.

We note, that the flow characteristics that were described in section 3 are not quirks of the simple sinusoidal topographies (3). The eddy saturation regime, the drag crisis, and multiple equilibria that are present here when geostrophic contours are open, have been all also recently found to exist in this model for a flow above a random monoscale topography with open geostrophic contours (Constantinou and Young 2017) and also above a multi-scale topography that has a $k^{-2}$ power spectrum (not reported).

In conclusion, the results presented here emphasized that barotropic processes might play a role in shaping the zonal momentum balance in the ACC and in producing eddy saturation. In particular, barotropic processes are responsible at least for the component of the transport that is related to the bottom flow velocity of the ACC. (This bottom-velocity component of the ACC transport accounts for about 25% of the total transport according to Donohue et al. (2016).) Even though the ACC is strongly affected by baroclinic processes, our results here reinforce the increasing evidence arguing for the importance of the barotropic mechanisms in determining the ACC transport (Ward and Hogg 2011; Thompson and Naveira Garabato 2014; Youngs et al. 2017). Additionally, the results of this paper argue that topography does not only have a passive role in the momentum budget acting merely as a sink for zonal momentum, but it can also have an active role by producing transient eddies that shape the standing eddy field and its form stress.

Acknowledgements.

I am mostly grateful to William Young for his support and his insightful comments. Discussions with Anand Gnanadesikan, Petros Ioannou, Spencer Jones, Sean Haney, Cesar Rocha, Andrew Thompson, Gregory Wagner, and Till Wagner are greatly acknowledged. Also, I thank Louis-Philippe Nadeau and David Straub for their constructive review comments. This project was supported by the NOAA Climate and Global Change Postdoctoral Fellowship Program, administered by UCAR’s Cooperative Programs for the Advancement of Earth System Sciences. \appendixtitleStability of steady states (9) to perturbations outside the low-dimensional manifold For the eigenanalysis of (12), first note that if $\phi$ is independent of $y$ then (12) implies stability. Therefore, a necessary condition instability of (12) is $\partial_{y}\phi\neq 0$. In this case, also note that:

1. (i)

(12) is homogeneous in $y$, i.e., if $\phi(x,y,t)$ is a solution so is $\phi(x,y+a,t)$ for any $a$, and 2. (ii)

the coefficients of (12) remain unchanged under $x\mapsto x+2\pi\kappa/m$ for any integer $\kappa$, i.e., if $\phi(x,y,t)$ is a solution so is $\phi(x+2\pi\kappa/m,y,t)$.

Remark (i) implies that the linear operator in (12) commutes with the translation operator in $y$ and therefore they share the same eigenfunctions. Thus, the eigensolution has $y$-dependence $\propto{\rm e}^{{\rm i}n_{y}y}$. Remark (ii) implies that the linear operator in (12) commutes with the appropriate translation operator by $2\pi\kappa/m$ in $x$, and thus the eigensolution has an $x$-dependence in the form of a Bloch wavefunction.555For a more detailed discussion on Bloch eigenfunctions the reader is referred, for example, to the textbook by Ziman (1979) or to the Appendix D in the thesis of Constantinou (2015). Therefore, we search for the eigensolutions of (12) as:

[TABLE]

with

[TABLE]

In (13) it is understood that the eigenvalue $\lambda$ depends on both $n_{x}$ and $n_{y}$. In (14), $n_{x}$ is the zonal Bloch wavenumber and it is restricted to take values $|n_{x}|\leq m/2$ (in condensed matter physics this is referred to as the first Brillouin zone). To see why $n_{x}$ is restricted, consider what (14) (or equivalently remark (ii)) implies for the zonal spectral power of $\phi$. Eigenfunction $\phi$ cannot have spectral power in any arbitrary zonal wavenumber, but instead it can only have power at zonal wavenumbers:

[TABLE]

From (15) we can verify that Bloch wavenumbers $n_{x}$ and $\pm(n_{x}-m)$ are completely equivalent and thus we only need to study the stability for $|n_{x}|\leq m/2$ (cf. Ziman (1979); Constantinou (2015)).

With the ansantz (13) the stability of (12) reduces to a generalized eigenvalue problem for each meridional wavenumber $n_{y}$:

[TABLE]

where $\nabla^{2}_{n_{y}}\stackrel{{\scriptstyle{\scriptscriptstyle\mathrm{def}}}}{{=}}\partial^{2}_{x}-n_{y}^{2}$ and

[TABLE]

The problem is further reduced if we insert (14) into (16). Then, the eigenproblem (16) reduces into an infinite system of equations whose unknowns are the coefficients $c_{M}$ and the eigenvalue $\lambda$. Truncating the sum in (14) up to $|M|\leq M_{\max}$ leaves us with $2M_{\max}+1$ equations for $2M_{\max}+2$ unknowns: $c_{-M_{\max}},c_{-M_{\max}+1},\dots,c_{M_{\max}}$ and $\lambda$. These equations are compactly written as:

[TABLE]

where A is a $(2M_{\max}+1)\times(2M_{\max}+1)$ matrix. A non-trivial solution of (18) exists only when $\det[\textsf{A}\left(\lambda\right)]=0$, and this last condition determines the eigenvalues $\lambda$ that correspond to the eigenfunction with wavenumbers $n_{x}$ and $n_{y}$ (Lorenz 1972; Gill 1974; Charney and Flierl 1980).