Two-dimensional turbulence above topography: condensation transition and selection of minimum enstrophy solutions

Abstract

We consider two-dimensional flows above topography, revisiting the selective decay (or minimum enstrophy) hypothesis of Bretherton and Haidvogel. We derive a ‘condensed branch’ of solutions to the variational problem where a domain-scale condensate coexists with a flow at the (smaller) scale of the topography. The condensate arises through a supercritical bifurcation as the conserved energy of the initial condition exceeds a threshold value, a prediction that we quantitatively validate using direct numerical simulations. We then consider the forced–dissipative case, showing how weak forcing and dissipation select a single dissipative state out of the continuum of solutions to the energy-conserving system predicted by selective decay. As the forcing strength increases, the condensate arises through a supercritical bifurcation for topographic-scale forcing and through a subcritical bifurcation for domain-scale forcing, both predictions being quantitatively validated by direct numerical simulations. This method provides a way of determining the equilibrated state of forced–dissipative flows based on variational approaches to the associated energy-conserving system, such as the statistical mechanics of two-dimensional flows or selective decay.

1. Introduction

Rapidly rotating flows above topography are both a standard model for various large-scale oceanic flows and a fascinating playground to study how disorder disrupts the large-scale self-organization of two-dimensional flows. A first oceanic situation of interest is the Antarctic Circumpolar Current (ACC), where turbulent eddies are channelled by and somewhat locked to the topography, with important consequences for both the turbulent transport and the speed of the ACC itself (Straub 1993; Nadeau & Straub 2009; Abernathey & Cessi 2014; Nadeau & Ferrari 2015; Constantinou & Young 2017; Constantinou 2018; Constantinou & Hogg 2019; Barthel et al. 2022; Yung, Morrison & Hogg 2022). At a more local scale, the Lofoten vortex has attracted significant attention (Köhl 2007; Søiland & Rossby 2013; Trodahl et al. 2020; LaCasce, Palóczy & Trodahl 2024) because its direction of rotation is inconsistent with the natural tendency for flows to homogenize potential vorticity (PV) (Rhines & Young 1982).

Early on, the study of rotating flows above topography has been connected to the search for general organizing principles governing the evolution of quasi-two-dimensional flows. Bretherton & Haidvogel (1976) (BH in the following) applied ‘selective decay’ to such flows, an organizing principle borrowed from plasma physics (Kruskal & Kulsrud 1958; Taylor 1974; Montgomery, Turner & Vahala 1978). According to the selective decay principle, the forward-cascading invariant should be minimized for a given (conserved) value of the inverse-cascading invariant. Based on the inverse energy cascade and forward enstrophy cascade of two-dimensional flows, BH thus conjectured that the system evolves towards a state of minimum enstrophy, while conserving its initial energy. They report reasonable, although imperfect, agreement with their numerical simulations. A second, more recent variational approach is the statistical mechanics of two-dimensional flows (Miller 1990; Robert 1991; Robert & Sommeria 1991; Bouchet & Venaille 2012). The approach has been applied to several two-dimensional fluid systems, predicting various branches of solutions with different flow structures (Naso, Chavanis & Dubrulle 2011), together with ‘phase transitions’ between these various states (Chavanis & Sommeria 1996; Venaille & Bouchet 2009). In its most general form, the theory is based on the conservation of all the invariants of the system: energy, all the moments of the vorticity field, together with linear momentum, angular momentum or circulation, depending on the geometry of interest. However, for most practical computations the high-order (potential-)vorticity moments are discarded (e.g. Salmon, Holloway & Hendershott 1976; Venaille 2012) and the theory ends up being equivalent to selective decay (Carnevale & Frederiksen 1987), namely enstrophy minimization subject to the conservation of some first- and second-order invariants (see Brands et al. (1999) and Nycander & Lacasce (2004) for examples where all the invariants are simultaneously considered). While they provide invaluable insight, an important limitation of such variational approaches is that they apply to energy-conserving systems only, as opposed to forced–dissipative ones. Bridging this gap is a central motivation for the present study.

Recently, Siegelman & Young (2023) (SY in the following) revisited two-dimensional turbulence above topography, with a focus on the difference between low-energy flows and high-energy flows. In qualitative agreement with selective decay, enstrophy robustly decreases with time in their numerical simulations, while energy is conserved to a very good approximation. However, SY also notice that the flows attained in the long-time limit depart from the predictions of selective decay. For low initial energy, the numerical solutions exhibit intense cyclones above topographic mounts and anticyclones above topographic depressions that are not predicted by selective decay. Instead, such pinned vortices result from the tendency for anticyclones (cyclones) to migrate downslope (upslope), resulting in a vortex segregation mechanism that opposes PV homogenization (Carnevale, Kloosterziel & Van Heijst 1991; Köhl 2007; Trodahl et al. 2020; Solodoch, Stewart & McWilliams 2021). Above a threshold value of the initial energy, these vortices are no longer pinned to the topography, and instead they roam within the numerical domain.

In the present study, we consider a spatial structure for the topography that differs from both BH and SY. While these authors include topography at all scales, we restrict attention to topography at a scale smaller than the extent of the fluid domain (strictly, although no scale separation is needed). Our goal is threefold. Firstly, we wish to continue the programme initiated by SY and further characterize the transition from the low-energy states to the high-energy states: How does the flow transition from a topographically locked state to the typical large-scale condensates that characterize weakly damped two-dimensional flows (Kraichnan 1975; Gallet & Young 2013; Laurie et al. 2014)? Secondly, we want to assess the predictive skill of the aforementioned variational approaches: Do they predict a sharp bifurcation or a smooth crossover between the two types of flows? Is the transition continuous or discontinuous? Can the theory quantitatively predict the strength of the large-scale condensate? Thirdly, we want to address the limitation of the variational approaches mentioned above, namely that they were designed for energy-conserving systems only, whereas natural and laboratory flows are typically forced and dissipative. While the conserved invariants of the system are the natural control parameters of the variational theories, such quantities are no longer conserved in the presence of forcing and dissipation. How does one make quantitative predictions in this context? How does one leverage the variational approaches to predict the equilibrated state of the forced–dissipative system?

We introduce the theoretical set-up in § 2, before focusing on the energy-conserving system in § 3. We predict a transition to condensation based on selective decay. In § 4, we report numerical simulations that quantitatively validate the predicted amplitude of the large-scale condensate. We turn to the forced–dissipative system in § 5. We show how weak forcing and damping select a single state out of the continuum of solutions to the energy-conserving system predicted by selective decay. We compare the predicted amplitude of the condensate with direct numerical simulations (DNS) of the forced–dissipative system, obtaining very good agreement for both topographic-scale forcing and domain-scale forcing. We conclude in § 6.

2. Rapidly rotating flows above topography

We consider a shallow-water flow above weak topography in a rapidly rotating frame, described within the $f$-plane quasi-geostrophic approximation. The horizontal velocity field ${\boldsymbol u}(x,y,t)$ stems from a streamfunction, ${\boldsymbol u} = - \boldsymbol {\nabla } \times [\psi (x,y,t) {\boldsymbol e}_z]$. The evolution of the system is governed by the material conservation of PV $q(x,y,t)$:

(2.1)$$\begin{gather} \partial_t q + J(\psi,q) = {\mathcal{F}}(x,y)+{\mathcal{D}} , \end{gather}$$

(2.2)$$\begin{gather}q= \Delta \psi + \eta(x,y) , \end{gather}$$

where $J(\,f,g)=(\partial _x\, f) (\partial _y g) - (\partial _y\, f)(\partial _x g)$ denotes the Jacobian operator, $\Delta$ denotes the Laplacian operator and $\eta (x,y)$ denotes the ‘topographic’ PV. For a layer of fluid of local depth $H+h(x,y)$, where $h(x,y) \ll H$ denotes the mean-zero fluctuations in bathymetry, the topographic PV is $\eta (x,y)=-f_0 h(x,y)/H$, with $f_0$ the Coriolis parameter. We consider (2.1) and (2.2) inside a doubly periodic domain, $(x,y)\in [0,2{\rm \pi} L]^2$, with $2{\rm \pi} L$ greater than the greatest wavelength of the topography (strictly, although no scale separation is required). We include forcing and dissipation on the right-hand side of (2.1). The forcing term ${\mathcal {F}}(x,y)$ corresponds to the wind-stress curl in a typical ocean setting (Salmon 1998). The dissipative term consists of linear ‘Ekman’ bottom drag (Pedlosky 1979) and hyperviscosity, ${\mathcal {D}}=-\kappa \Delta \psi - \nu \varDelta ^3 \psi$. In the following we consider both the ‘energy-conserving case’, characterized by ${\mathcal {F}}=0$, $\kappa =0$ and small hyperviscosity $\nu$, and the ‘forced–dissipative case’, characterized by ${\mathcal {F}} \neq 0$ and $\kappa \neq 0$.

3. Energy-conserving case: transition to condensation

When the right-hand side of (2.1) is set to zero, ${\mathcal {F}}= {\mathcal {D}}=0$, the system possesses two quadratic invariants. Denoting space average as $\langle \cdots \rangle$, the first invariant is the kinetic energy $E=\langle |\boldsymbol {\nabla } \psi |^2 \rangle$, where we omit the prefactor $1/2$ to alleviate the algebra. The second invariant is the enstrophy $Q=\langle q^2 \rangle$, omitting again the prefactor $1/2$. In a similar fashion to two-dimensional turbulence in the absence of topography, one expects the transient dynamics of the system to be characterized by an inverse cascade of energy together with a forward cascade of enstrophy. If a small hyperviscosity $\nu \neq 0$ is included, we expect $Q$ to be robustly dissipated once enstrophy has reached the small viscous dissipative scale. The Kolmogorov–Kraichnan dual cascade phenomenology further suggests that such ‘anomalous’ dissipation of enstrophy is maintained as $\nu \to 0$, whereas the energy $E$ is conserved in that limit. That enstrophy is strongly dissipated by weak small-scale damping while energy is approximately conserved is clearly illustrated by the recent DNS of SY. This intuition led BH to conjecture that the system eventually reaches the minimum enstrophy state associated with the conserved energy of the initial condition. This state is obtained by extremalizing the functional ${\mathcal {L}}\{ \psi \}=Q+\mu E$, where $\mu$ is a Lagrange multiplier ensuring energy conservation. This variational problem leads to

(3.1)\begin{equation} q = \Delta \psi + \eta = \mu \psi . \end{equation}

Note that the same relation can be deduced from the statistical mechanics of two-dimensional flows, if one retains only the invariants $E$ and $Q$. An argument in favour of selective decay as opposed to the statistical mechanics of conservative systems is that $Q$ consistently decreases during the transient phase observed in DNS.

An appealing aspect of enstrophy minimization is that fields satisfying (3.1) are steady solutions to (2.1) with ${\mathcal {F}}={\mathcal {D}}=0$. Additionally, the solution to (3.1) corresponding to the absolute enstrophy minimum is stable according to Arnol'd's stability criterion (Arnol'd 1978). The following subsections aim at computing this absolute minimizer.

3.1. The Bretherton–Haidvogel branch

The standard approach to solving (3.1) is the one proposed by BH. They assume that $\mu >-k_0^2$, where $k_0=1/L$ denotes the minimum (gravest) wavenumber of the square domain. Introducing the Fourier decompositions of $\eta$ and $\psi$ as

(3.2a,b)\begin{equation} \eta(x,y) = \sum_{{\boldsymbol k}L \in \mathbb{Z}^2} \hat{\eta}_{\boldsymbol k} \, {\rm e}^{{\rm i} {\boldsymbol k}\ \boldsymbol{\cdot}\ \boldsymbol{x}} , \quad \psi(x,y,t) = \sum_{{\boldsymbol k}L \in \mathbb{Z}^2} \hat{\psi}_{\boldsymbol k}(t) \, {\rm e}^{{\rm i} {\boldsymbol k}\ \boldsymbol{\cdot}\ \boldsymbol{x}} ,\end{equation}

equation (3.1) has a unique solution obtained in spectral space:

(3.3)\begin{equation} \hat{\psi}_{\boldsymbol k} =\frac{\hat{\eta}_{\boldsymbol k}}{\mu+k^2}.\end{equation}

From (3.3) the energy $E$ and the enstrophy $Q$ can be expressed in terms of $\mu$. This leads to a parametric curve for $Q$ versus $E$, and one thus predicts (among other things) the final enstrophy in terms of the initial energy of the system.

3.2. Monoscale topography

Although not a necessary assumption, restricting attention to monoscale topography provides an illustrative example. We thus consider that the topography is peaked around a single wavenumber $k_\eta$, such that $\Delta \eta = - k_\eta ^2 \eta$. For such topography (3.3) leads to

(3.4a–c)\begin{equation} \psi = \frac{\eta}{\mu+k_\eta^2} , \quad E = \frac{k_\eta^2}{(\mu+k_\eta^2)^2} \langle \eta^2 \rangle, \quad Q = \frac{\mu^2}{(\mu+k_\eta^2)^2} \langle \eta^2 \rangle . \end{equation}

It proves convenient to introduce the energy scale $E_\eta =\langle \eta ^2 \rangle /k_\eta ^2$ associated with the topography (following SY, $E_\eta$ also corresponds to the energy of a flow with homogenized PV; see discussion in § 6). Denoting the dimensionless energy as $\tilde {E}=E/E_\eta$, one can eliminate $\mu$ to express $Q$ in terms of $E$:

(3.5)\begin{equation} \frac{Q}{\langle \eta^2 \rangle} = (\sqrt{\tilde{E}}-1 )^2 .\end{equation}

3.3. The condensed branch

When the system has a size greater than the largest scale of the topography (when $k_0=1/L< k_\eta$ for monoscale topography) we argue that the BH branch above is not always the absolute minimizer of the enstrophy. Indeed, consider (3.1) with $\mu =-k_0^2$:

(3.6)\begin{equation} \Delta \psi + k_0^2 \psi =- \eta . \end{equation}

The operator on the left-hand side is no longer invertible and the solution is (restricting attention to monoscale topography, the generalization to an arbitrary spectrum for the topography being straightforward)

(3.7)$$\begin{gather} \psi(x,y) = \psi_\eta(x,y) + A \, \psi_0(x,y) , \end{gather}$$

(3.8)$$\begin{gather}\psi_\eta(x,y) = \frac{\eta(x,y)}{k_\eta^2-k_0^2} . \end{gather}$$

In (3.7), $\psi _0(x,y)$ is a harmonic function at the gravest scale of the domain (that is, $\Delta \psi _0 =- k_0^2 \psi _0$) normalized such that $\langle \psi _0^2 \rangle =1$, and $A \geq 0$ denotes the amplitude of this large-scale ‘condensate’. Terms $A$ and $\psi _0(x,y)$ are related to the two gravest Fourier coefficients of $\psi$ through $A \, \psi _0(x,y) = \hat {\psi }_{(1/L,0)} \, {\rm e}^{{\rm i} x/L} + \hat {\psi }_{(0,1/L)} \, {\rm e}^{{\rm i} y/L}+ \text {c.c.}$, that is,

(3.9)\begin{equation} A =\sqrt{2 |\hat{\psi}_{(1/L,0)}|^2 + 2 |\hat{\psi}_{(0,1/L)}|^2} .\end{equation}

The amplitude $A$ of the large-scale condensate can be expressed in terms of the conserved energy $E$ of the initial condition:

(3.10)\begin{equation} E=A^2 \langle |\boldsymbol{\nabla} \psi_0|^2 \rangle + \langle |\boldsymbol{\nabla} \psi_\eta|^2 \rangle , \end{equation}

that is,

(3.11)\begin{equation} A= \sqrt{\frac{E- \langle |\boldsymbol{\nabla} \psi_\eta|^2 \rangle}{\langle |\boldsymbol{\nabla} \psi_0|^2 \rangle}} = \frac{1}{k_0} \sqrt{E- \langle |\boldsymbol{\nabla} \psi_\eta|^2 \rangle} .\end{equation}

For monoscale topography the expression of $\langle |\boldsymbol {\nabla } \psi _\eta |^2 \rangle$ reads

(3.12)\begin{equation} \frac{\langle |\boldsymbol{\nabla} \psi_\eta|^2 \rangle}{E_\eta} = \frac{\tilde{k}_\eta^4}{(\tilde{k}_\eta^2-1)^2} , \end{equation}

where $\tilde {k}_\eta =k_\eta /k_0$, and after substitution in (3.11) one obtains the following expression for the dimensionless amplitude of the condensate $\tilde {A}=k_0 A/\sqrt {E_\eta }$:

(3.13)\begin{equation} \tilde{A}= \sqrt{\tilde{E} - \tilde{E}_c} , \quad \text{with } \tilde{E}_c = \frac{\tilde{k}_\eta^4}{(\tilde{k}_\eta^2-1)^2}. \end{equation}

This prediction takes a particularly simple form in the limit of scale separation between the domain size and the wavelength of the topography, $\tilde {k}_\eta \gg 1$, which leads to $\tilde {E}_c \simeq 1$:

(3.14)\begin{equation} \tilde{A}= \sqrt{\tilde{E} - 1} . \end{equation}

As illustrated by the plot in figure 1, equation (3.14) displays a transition to large-scale condensation above a critical value $E_\eta$ of the energy of the initial condition. An interpretation of this result is that the topography plays the role of a reservoir that can host a small-scale flow of energy up to $E_\eta$. Any energy beyond $E_\eta$ spills out of this reservoir and is subject to inverse energy transfers. This excess energy eventually ends up accumulating in a large-scale condensate, hence the non-zero value of $A$.

Figure 1. The amplitude $\tilde {A}$ of the condensate as a function of the (conserved) energy of the initial condition $\tilde {E}$. The theory predicts the existence of a condensed branch for $\tilde {E}>1$ in the scale separation limit (slightly larger threshold for finite $\tilde {k}_\eta$). The DNS are performed using $\tilde {k}_\eta =12$. The data points agree well with the predicted amplitude, for the two initial conditions introduced in § 4. Some departure from the theory remains in the immediate vicinity of the threshold for condensation, as a consequence of the extra vortices pinned to the topography.

The enstrophy associated with the condensed branch is

(3.15)\begin{equation} \frac{Q}{\langle \eta^2\rangle}=\frac{1}{\tilde{k}_\eta^2} \left(\tilde{E} - 1 - \frac{1}{\tilde{k}_\eta^2-1} \right) , \end{equation}

valid only when $A$ exists, that is, for $\tilde {E}>\tilde {E}_c$. In the scale separation limit this expression reduces to

(3.16)\begin{equation} \frac{Q}{\langle \eta^2\rangle}=\frac{1}{\tilde{k}_\eta^2} (\tilde{E} - 1 ) , \end{equation}

valid for $\tilde {E}>1$.

In figure 2 we show the enstrophy $Q$ as a function of the energy $E$ for both the BH branch and the condensed branch. One can check that the enstrophies $Q(E)$ evaluated on the two branches (3.5) and (3.16) are equal and tangent to one another at the critical value $\tilde {E}=\tilde {E}_c$. For $\tilde {E}<\tilde {E}_c$ the absolute enstrophy minimizer is the BH branch, while for $\tilde {E}>\tilde {E}_c$ the absolute minimizer corresponds to the condensed branch. Absolute minimization of enstrophy thus predicts a transition to condensation as $E$ exceeds $E_c$, the amplitude of the condensate being given by (3.13). These predictions are particularly compact in the limit of scale separation, where the dimensionless threshold for condensation reduces to $\tilde {E}_c=1$ and the amplitude of the condensate is given by (3.14).

Figure 2. Enstrophy $Q$ in the (quasi-)stationary state as a function of the (conserved) energy of the initial condition $\tilde {E}$, using the same symbols as in figure 1. Also shown are the results of the minimization: BH branch in blue and condensed branch in red, for $\tilde {k}_\eta =12$. The absolute minimum (over the branches) is shown as a black dashed line. The DNS data are close to the BH branch for $\tilde {E} \leq 1$ and close to the condensed branch for $\tilde {E} \geq 1$. The inset is a zoom on the data for $\tilde {E} \geq 1$. The departure from the absolute minimum near the condensation threshold is due to the extra vortices pinned to the topography.

4. Numerical simulations

To test the predictions above, we perform DNS of (2.1) and (2.2) by time-stepping $\psi$ using a pseudo-spectral solver. The resolution is $1024^2$ before de-aliasing following the $2/3$ rule. We first consider the enstrophy-decaying energy-conserving case, ${\mathcal {F}}=0$ and $\kappa =0$ with a small value of the hyperviscosity $\nu k_0^3/ \sqrt {E_\eta }=2\times 10^{-10}$. The topography is approximately monoscale and centred around the dimensionless wavenumber $\tilde {k}_\eta =12$. Specifically, we pick random values for the Fourier coefficients of all the modes such that $12-1/4 \leq kL \leq 12+1/4$. All the numerical runs are performed using the same realization of the topography, displayed in figure 3(a).

Figure 3. For the topography shown in (a), a numerical run initialized with energy $\tilde {E}=0.743$ displays no condensation (b), whereas a run initialized with $\tilde {E}=1.44$ settles in a condensed state in the long-time limit (c).

We consider two types of initial conditions.

1. (i) The ‘$\psi$ initial condition’: we set some initial complex Fourier amplitude at wavenumber $(k_x L,k_y L)=(0,2)$ for $\psi$ to match the target initial value of $\tilde {E}$.

2. (ii) The ‘$\phi$ initial condition’: we set some initial complex Fourier amplitude at wavenumber $(k_x L,k_y L)=(0,2)$ for a field $\phi (x,y)$, before initializing $\psi$ as $\psi (x,y,t=0)=\psi _\eta (x,y)+\phi (x,y)$.

Both types of initial conditions correspond to initial perturbations at rather large scale. Starting from the $\psi$ initial condition, the flow spontaneously develops smaller structures through interaction with the topography. By contrast, the $\phi$ initial condition readily includes the small-scale flow structure predicted by enstrophy minimization. This initial condition only allows one to achieve an initial energy $\tilde {E}\geq 1$. By contrast, one can achieve any positive initial value for the energy $\tilde {E}$ using the $\psi$ initial condition.

As discussed by SY, in such simulations the enstrophy rapidly decays from its initial value. By contrast, the energy decreases extremely slowly and is conserved to a very good approximation. We wish to characterize the (quasi-)stationary state obtained once the transient has subsided. In figure 3 we show snapshots obtained at the end time of two simulations differing in the value of the initial energy. The streamfunction $\psi$ agrees qualitatively with the theory: the simulation with $\tilde {E}<1$ displays no large-scale condensation, whereas the simulation with $\tilde {E}>1$ displays a strong large-scale condensate. By contrast, in line with previous studies the PV $q$ exhibits strong isolated vortices that are not predicted by enstrophy minimization (Zhao, Chieusse-Gérard & Flierl 2019; Siegelman & Young 2023; LaCasce et al. 2024). When condensation arises, however, these vortices appear to be rather sparsely distributed above the faint background PV structure of the condensed branch, given by (3.1) with $\mu =-k_0^2$. There is thus hope for large-scale condensation to be quantitatively predicted by enstrophy minimization despite the emergence of the isolated vortices.

To characterize the condensation transition quantitatively, we extract from the DNS the amplitude $A$ as defined in (3.9). In figure 1 we show the resulting dimensionless condensate amplitude $\tilde {A}$ as a function of the dimensionless initial energy $\tilde {E}$. The numerical data are in very good agreement with the theoretical prediction from enstrophy minimization. The $\phi$ initial condition leads to particularly good agreement because this initial condition is close to the enstrophy minimizer, thus preventing the emergence of strong isolated PV vortices. The more generic $\psi$ initial condition also leads to very good quantitative agreement. There is some departure from the theoretical prediction in the immediate vicinity of the bifurcation threshold only, as a consequence of the isolated vortices pinned to the topography.

In figure 2, we plot the equilibrated enstrophy $Q$ extracted from the DNS. Above the transition $Q$ is significantly lower than the enstrophy of the BH branch, because the latter is not the absolute minimizer. The equilibrated $Q$ lies extremely close to the prediction of the condensed branch for DNS initialized using the $\phi$ initial condition. It is slightly greater than (but close to) this absolute minimum for the $\psi$ initial condition, as a result of the extra isolated vortices. Once again, this slight departure from the absolute minimum of $Q$ has only a modest impact on the value of the condensate amplitude $A$, which closely follows the theoretical prediction.

5. Forced–dissipative case

The natural objection to the variational approaches to computing the statistically steady states of fluid flows is that these approaches are usually based on conservative dynamics, while real flows are forced and dissipative. For instance, the statistical mechanics of two-dimensional flows assumes that energy and enstrophy are conserved and keep their initial values, while DNS point to a transient phase with a robust decrease in enstrophy. In that respect, selective decay (enstrophy minimization) is one step closer to realistic flows, allowing for such a net decrease in enstrophy as the system evolves towards the minimizer, in line with the forward enstrophy cascade of two-dimensional flows. However, selective decay also assumes that energy is conserved and set by the initial condition, whereas natural and laboratory flows typically receive energy from some external forcing and dissipate it through a combination of damping mechanisms. To make contact between such forced–dissipative systems and the variational framework discussed above, we now consider (2.1) with non-zero forcing ${\mathcal {F}}$ injecting energy into the system, and non-zero friction coefficient $\kappa$ extracting energy from the system (we also retain the hyperviscous coefficient). The goal is at a fundamental level: more than describing realistic flows of geophysical interest, we wish to establish a direct quantitative connection between the variational approach described above and forced–dissipative systems. We thus illustrate how weak forcing and dissipation select a unique state out of the continuum of solutions to the energy-conserving problem predicted by selective decay.

5.1. Selection of minimum enstrophy states

Consider (2.1) with weak forcing and damping:

(5.1)$$\begin{gather} \partial_t q + J(\psi,q) =\epsilon[ {\mathcal{F}}^{(1)}(x,y)+{\mathcal{D}}^{(1)} ] , \end{gather}$$

(5.2)$$\begin{gather}q= \Delta \psi + \eta, \end{gather}$$

where $\epsilon \ll 1$ is a small parameter, ${\mathcal {D}}^{(1)}=-\kappa ^{(1)} \Delta \psi - \nu ^{(1)} \varDelta ^3 \psi$ and the topography $\eta$ is ${O}(1)$. Expanding the fields in powers of $\epsilon$:

(5.3)$$\begin{gather} q(x,y,t) = q^{(0)}(x,y,T)+\epsilon q^{(1)}(x,y,T) + \cdots , \end{gather}$$

(5.4)$$\begin{gather}\psi(x,y,t) = \psi^{(0)}(x,y,T)+\epsilon \psi^{(1)}(x,y,T) + \cdots , \end{gather}$$

where we have introduced the slow time variable $T=\epsilon t$, that is, $\partial _t \to \partial _t +\epsilon \partial _T$.

Collecting terms of ${O}(\epsilon ^0)$ leads to the set of equations discussed in the preceding section:

(5.5a,b)\begin{equation} \partial_t q^{(0)} + J(\psi^{(0)},q^{(0)}) = 0 , \quad q^{(0)}= \Delta \psi^{(0)} + \eta . \end{equation}

We consider a solution to these equations that lies on the condensed branch, allowing the amplitude of the condensate to vary slowly in time:

(5.6)\begin{equation} \psi^{(0)}(x,y,T) = \psi_\eta(x,y) + A(T) \psi_0(x,y) . \end{equation}

Following the traditional approach, the slow evolution equation for the amplitude $A(T)$ can be obtained from a solvability condition at next order. An alternative and more enlightening way of obtaining this amplitude equation consists of writing the energy evolution equation to lowest order. Upon multiplying equation (5.1) with $\psi$ before averaging over $x$ and $y$ and neglecting terms of order $\epsilon ^2$, one obtains

(5.7)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}T} \left( \frac{\tilde{A}^2}{2} \right) \simeq-\frac{ \langle {\mathcal{F}}^{(1)} \eta \rangle }{\langle \eta^2 \rangle} - \kappa^{(1)} - \nu^{(1)} k_\eta^4 -\frac{\langle {\mathcal{F}}^{(1)} \psi_0 \rangle }{k_0 \sqrt{E_\eta}} \tilde{A} - \kappa^{(1)} \tilde{A}^2 , \end{equation}

where, for brevity, we restrict attention to the scale separation limit $\tilde {k}_\eta \gg 1$. After dividing by $\kappa ^{(1)}$, (5.7) can be recast in terms of the original (non-expanded) variables as

(5.8)\begin{equation} \frac{1}{2 \kappa}\frac{\mathrm{d}}{\mathrm{d}t} \tilde{A}^2 = Gr_\eta-1-\frac{\nu k_\eta^4}{\kappa} - \tilde{A}^2 +Gr_0 \tilde{A} . \end{equation}

For body-forced flows the dimensionless strength of the forcing is usually quantified by a Grashof number (Doering & Gibbon 1995). In (5.8) we have thus introduced a topographic-scale Grashof number $Gr_\eta$ and a domain-scale Grashof number $Gr_0$, defined, respectively, as

(5.9a,b)\begin{equation} Gr_\eta =- \frac{\langle {\mathcal{F}} \eta\rangle }{\kappa \langle \eta^2 \rangle } \quad \text{and} \quad Gr_0 =- \frac{\langle{\mathcal{F}} \psi_0 \rangle }{k_0 \kappa \sqrt{E_\eta}} . \end{equation}

The forcing thus enters the evolution equation (5.8) for the amplitude of the condensate in two different ways: through its correlation with the smaller-scale topography and through its correlation with the domain-scale condensate. In the following we consider these two situations independently.

5.2. Topographic-scale forcing: continuous transition to condensation

Consider first the situation where there is no forcing at large scale, $\langle {\mathcal {F}} \psi _0 \rangle =0$ and therefore $Gr_0=0$. Equation (5.8) reduces to

(5.10)\begin{equation} \frac{1}{2 \kappa}\frac{\mathrm{d}}{\mathrm{d}t} \tilde{A}^2 = Gr_\eta-Gr_\eta^{(c)}- \tilde{A}^2 , \end{equation}

where we have introduced a critical Grashof number $Gr_\eta ^{(c)}=1+{\nu k_\eta ^4}/{\kappa }$.

For $Gr_\eta < Gr_\eta ^{(c)}$, $\tilde {A}^2$ decreases and becomes negative according to this equation. The ansatz (5.6) is then invalid and instead the system ends up on the non-condensed BH branch. By contrast, a steady solution $\tilde {A} \geq 0$ exists for $Gr_\eta \geq Gr_\eta ^{(c)}$. We thus predict the following value for the equilibrated amplitude of the condensate:

(5.11)$$\begin{gather} \tilde{A} = 0 \quad \text{for } Gr_\eta< Gr_\eta^{(c)} , \end{gather}$$

(5.12)$$\begin{gather}\tilde{A} = \sqrt{Gr_\eta-Gr_\eta^{(c)}} \quad \text{for } Gr_\eta \geq Gr_\eta^{(c)} . \end{gather}$$

This prediction corresponds to a supercritical transition to condensation as the dimensionless forcing strength $Gr_\eta$ increases. The bifurcation threshold reduces to the compact expression $Gr_\eta ^{(c)} = 1$ in the limit where hyperviscosity damps energy at a negligible rate as compared with friction, $\nu k_\eta ^4\ll \kappa$.

As discussed above, the forcing enters the amplitude equation (5.10) through its correlation with the topography. We first illustrate this situation numerically using a forcing that is perfectly correlated with the topography:

(5.13)\begin{equation} {\mathcal{F}}(x,y)=-Gr_\eta \kappa \eta(x,y). \end{equation}

We consider the same realization of the topography as in the preceding sections, together with damping coefficients $\kappa /(k_0 \sqrt {E_\eta })=0.0021$ and $\nu k_0^3/ \sqrt {E_\eta }=2.8 \times 10^{-8}$. The corresponding critical Grashof number is $Gr_\eta ^{(c)} \simeq 1.28$. We perform a suite of numerical simulations for various values of $Gr_\eta$, extracting the amplitude of the condensate once the transient has subsided. We show $\tilde {A}$ as a function of $Gr_{\eta}$ in figure 4, together with the theoretical prediction (5.12). The agreement is very good, both for the threshold value $Gr_\eta ^{(c)}$ and for the amplitude of the condensate above threshold. Together with the simple forcing described above, we have considered a second forcing structure where ${\mathcal {F}}(x,y)$ is proportional to $\eta (x,y)+\xi (x,y)$, where $\xi (x,y)$ is a random field with the same statistics as $\eta (x,y)$ but corresponding to another realization (the realizations of $\eta$ and $\xi$ are the same for the entire suite of simulations). Even though the resulting forcing ${\mathcal {F}}(x,y)$ is less correlated to the topography than in the previous suite of simulations, the condensate amplitude is again well predicted by (5.12), as shown in figure 4.

Figure 4. Condensate amplitude as a function of the forcing expressed in terms of the Grashof number. The DNS data points are in good agreement with the theoretical prediction (5.12) associated with the selection of minimum enstrophy solutions (dashed line).

5.3. Large-scale forcing: discontinuous transition to condensation

Consider now the situation where the forcing entering the amplitude equation (5.8) is at large scale only, $Gr_\eta =0$ and $Gr_0\neq 0$. For simplicity, we also consider that the hyperviscous contribution to this amplitude equation is negligible, $\nu k_\eta ^4 \ll \kappa$. Equation (5.8) reduces to

(5.14)\begin{equation} \frac{1}{2 \kappa}\frac{\mathrm{d}}{\mathrm{d}t} \tilde{A}^2 = Gr_0 \tilde{A} -1 - \tilde{A}^2, \end{equation}

with steady solutions

(5.15)\begin{equation} \tilde{A}_\pm= \tfrac{1}{2} \left( Gr_0 \pm \sqrt{Gr_0^2-4} \right) .\end{equation}

These solutions exist only for $Gr_0 \geq 2$, while for $Gr_0<2$ the system lies on the uncondensed BH branch, with $\tilde {A}=0$. One can furthermore check that $A_-$ is unstable according to (5.14), while $A_+$ is a stable branch of solutions. We thus predict a subcritical bifurcation to large-scale condensation as $Gr_0$ exceeds the threshold value $Gr_0^{(c)}=2$. It is interesting to notice that the same supercritical bifurcation of the energy-conserving system induces either a supercritical or a subcritical bifurcation in the weakly forced–dissipative system, depending on the structure of the forcing.

We investigated this situation numerically with the same resolution and the same realization of the topography as before. The friction coefficient is still $\kappa /(k_0 \sqrt {E_\eta })=0.0021$, while we employ the much smaller value of the hyperviscosity $\nu k_0^3/ \sqrt {E_\eta }=2\times 10^{-10}$, which ensures that $\nu k_\eta ^4 \ll \kappa$. The forcing is

(5.16)\begin{equation} {\mathcal{F}}(\kern 0.09em y)=-\sqrt{2 E_\eta} k_0 \kappa Gr_0 \cos(\kern 0.09em y) ,\end{equation}

associated with a structure $\psi _0=\sqrt {2}\cos y$ for the large-scale flow. To investigate possible bistability, the initial condition is either very low-amplitude fluctuations or the end state of a previous simulation.

The amplitude of the condensate is shown as a function of $Gr_0$ in figure 5. Once again, the agreement with the theory is excellent at finite distance from the instability threshold, with slight discrepancies near threshold. In particular, the amplitude of the condensate is very well captured by the upper branch in (5.15) for $Gr_0 \geq 2$. We observe only a very narrow region of bistability, possibly because the uncondensed branch has other directions of instability, besides direct condensation following (5.14). Additionally, the region of stability of the condensed branch extends somewhat below $Gr_0=2$, down to $Gr_0\simeq 1.8$. This $10\,\%$ correction in threshold value may be a consequence of the extra vortices pinned to the topography. As for the previous cases considered above, we conclude that the amplitude of the condensate is accurately predicted by the theory, except in the immediate vicinity of the instability threshold where the pinned vortices impact the dynamics of the system.

Figure 5. Condensate amplitude as a function of the large-scale forcing amplitude expressed in terms of $Gr_0$. The DNS data points are in good agreement with the theoretical prediction $A_+$ in (5.15) associated with the selection of minimum enstrophy solutions (dashed line).

5.4. Structure of the condensate

We have mainly discussed the amplitude of the condensate, setting aside its degenerate structure $\psi _0(x,y)$. The latter depends on the relative magnitudes of the complex Fourier coefficients $\hat {\psi }_{(1/L,0)}$ and $\hat {\psi }_{(0,1/L)}$, as well as on their phases. In simulations of the energy-conserving system, the precise structure of the condensate (see e.g. figure 3) seems to be determined predominantly by the precise realization of the topography and, to a lesser extent, by the initial condition. When driving the system with the domain-scale forcing (5.16) together with moderately small friction, the condensate always adopts the large-scale structure of the forcing, $\psi _0(x,y)=\sqrt {2} \cos (y)$, and the theoretical predictions are accurately satisfied. However, if $\epsilon$ on the right-hand side of (5.1) is made extremely small (very weak forcing and dissipation), there seems to be a competition between the forcing, which favours the structure $\psi _0(x,y)=\sqrt {2} \cos (y)$, and the precise realization of the topography, which favours an alternative structure, such as the one visible in figure 3. The prediction (5.15) for the condensate amplitude does not hold if the system adopts this alternative condensate structure. An in-depth characterization of such phase dynamics in the very-weak-drag regime would deserve a dedicated study that goes beyond the scope of the present article. One way to lift the degeneracy of the condensate structure consists of including weak domain-scale topography in a perturbative fashion, as described in Appendix A. Then a single condensate structure $\psi _0(x,y)$ proportional to the domain-scale topography arises, the condensate amplitude being still determined by forcing and dissipation (or by the initial energy for the energy-conserving system).

6. Conclusion

We have characterized the transition to large-scale condensation for two-dimensional flows over smaller-scale topography, revisiting the selective decay or minimum enstrophy hypothesis of BH. We have computed a condensed branch of solutions to the minimization problem. This branch of solutions is the absolute enstrophy minimizer above some threshold value of the (conserved) energy. As this control parameter increases, selective decay predicts a continuous transition to large-scale condensation. This bifurcation adds to the various ‘phase transitions’ of quasi-two-dimensional fluid systems predicted based on the statistical mechanics of two-dimensional flows, such as the transition in flow structure as the aspect ratio of a two-dimensional rectangular fluid domain is varied (Chavanis & Sommeria 1996; Bouchet & Venaille 2012). The present bifurcation is also very similar to the emergence of a Rossby wave in the enstrophy minimizing state of a channel flow on the $\beta$ plane (Young 1987), the connection between the two problems being apparent if $\beta$ is interpreted as a uniform topographic slope.

For simplicity, we have illustrated the present theory using monoscale topography, but the approach is easily extended to more complex topography. As compared with the recent study of flows above topography by SY, the main qualitative difference is probably that SY include large topographic variations at the domain scale. Authors SY report a transition from isolated pinned vortices to freely roaming vortices as the initial energy exceeds a threshold value $E_\sharp =\langle | \boldsymbol {\nabla } \varDelta ^{-1} \{ \eta \} |^2 \rangle$, omitting the prefactor $1/2$ for consistency with our definition for $E$ (SY identify $E_\sharp$ as the energy of a flow with homogeneous PV $q=0$, which for monoscale topography reduces to the energy scale $E_\eta$ that we introduced on dimensional grounds). By contrast, in the absence of topography at the domain scale we observe a transition to condensation as the initial energy exceeds a threshold value $E_c=\langle |\boldsymbol {\nabla } \psi _\eta |^2 \rangle$. Interestingly, while $E_c \neq E_\sharp$ for arbitrary topography, $E_c$ and $E_\sharp$ reduce to the same expression in the limit of separation between the topographic scale and the domain scale. This indicates that the study of SY and the present study may correspond to two manifestations of the same transition. An appealing aspect of the present set-up is that selective decay provides quantitative predictions for the amplitude of the condensate above threshold.

Indeed, another goal of the present study is to assess the predictive skill of such a variational approach through comparison with DNS. The DNS of the energy-conserving system show that selective decay quantitatively predicts the amplitude of the condensate. For initial conditions tailored to be close to the absolute minimizer, the amplitude of the condensate is in excellent agreement with the theory. For more generic initial conditions isolated vortices arise, polluting the PV field $q$. The vortices affect the very near-threshold behaviour of the condensate amplitude only, the agreement remaining very good at larger distance from threshold. That the vortices predominantly affect the near-threshold behaviour probably originates from the fact that they become more scarce as condensation proceeds. Indeed, the isolated vortices are located at the local extrema of the total streamfunction, as discussed in SY. In the absence of large-scale condensation, the streamfunction resembles the topography, with many such extrema. By contrast, at finite distance above threshold the condensate dominates the total streamfunction, and the latter has fewer local extrema to host isolated vortices (see figure 3).

The above-mentioned dependence on initial conditions is perhaps a complication of the energy-conserving system. After all, natural and laboratory flows are typically forced and dissipative, with energy supplied by an external source and removed by some combination of damping mechanisms. The energy of the equilibrated state is then an emergent property of the system, as opposed to the initial energy of the system being a control parameter of the variational approaches to conservative flows. As such, the variational approaches to energy-conserving systems cannot predict the strength of forced–dissipative flows. With the goal of bridging the gap between these two situations, we have introduced a perturbative expansion that illustrates how weak forcing and dissipation select the equilibrated state of the flow out of the continuum of states predicted by selective decay (or by the statistical mechanics of two-dimensional flows). We have restricted attention to steady forcing for simplicity, but the approach probably holds for slowly varying forcing.

This selection mechanism is very similar to – and was inspired by – the selection of solitary wave solutions to the complex Ginzburg–Landau equation by weak dissipation (Fauve & Thual 1990). The theoretical predictions for the resulting amplitude of the condensate agree quantitatively with DNS of the forced–dissipative system. Interestingly, the same continuous transition to condensation for the energy-conserving system turns into a discontinuous transition when the system is forced at large scale, while it remains continuous when the system is forced at the smaller topographic scale.

This method provides an avenue for predicting the equilibrated state of forced–dissipative flows based on variational approaches initially designed for conservative systems (selective decay or statistical mechanics). Its robustness and predictive skill remain to be further assessed through application to various systems of interest, be they numerical simulations, laboratory experiments or observational data.