Convective plumes in rotating systems

Abstract

Convective plumes emanating from fixed buoyant sources such as volcanoes, hot springs and oil spills are common in the atmosphere and the ocean. Most of what we know about their dynamics comes from scaling laws, laboratory experiments and numerical simulations. A plume grows laterally during its ascent mainly due to the process of turbulent entrainment of fluid from the environment into the plume. In an unstratified system, nothing hampers the vertical motion of the plume. By contrast, in a stratified system, as the plume rises, it reaches and overshoots the neutral buoyancy height – due to the non-zero momentum at that height. This rising fluid is then dense relative to the environment and slows down, ceases to rise and falls back to the height of the intrusion. For buoyant plumes occurring in the ocean or atmosphere, the rotation of the Earth adds an additional constraint via the conservation of angular momentum. In fact, the effect of rotation is still not well understood, and we addressed this issue in the study reported here. We looked for the steady states of an axisymmetric model in both the rotating and non-rotating cases. At the non-rotating limit, we isolated two regimes of convection depending on the buoyancy flux/momentum flux ratio at the base of the plume, in agreement with scaling laws. However, the inclusion of rotation in the model strongly affects these classical convection patterns: the lateral extension of the plume is confined at the intrusion level by the establishment of a geostrophic balance, and non-trivial swirl speed develops in and around the plume.

1 Introduction

Hydrothermal vents are often found in fracture areas at the sea floor. Each vent is a steady source of hot water and, sometimes, gas and other miscible and immiscible materials. Such a buoyant source generates a convective plume that rises and mixes with the surrounding water. The total equivalent heat flux generated by a single source is of the order of $10^{6}$  W (Converse, Holland & Edmond 1984; Carazzo, Kaminski & Tait 2008). Stein, Stein & Pelayo (1995) estimated a total heat flux, if all these sources were combined, of 32 TW. They argued that these sources can have a non-negligible effect on oceanic abyssal circulation. Scott, Marotzke & Adcroft (2001) also showed how hydrothermal vents affect the meridional overturning circulation. They demonstrated that the circulation that develops in response to these heat sources not only depends on the magnitude of the heat flux but is also strongly related to the stratification and the depth of the thermocline. When a plume rises, it entrains surrounding fluid, and the stratification is a key factor needed to arrest the vertical motion, as described by Morton, Taylor & Turner (1956) (henceforth MTT). Using simplified equations of the dynamics, MTT were able to obtain scaling laws to describe the shape of the plume and the height of the neutral level depending on the buoyancy flux and the Brunt–Väisälä frequency  $N$ .

While the original scaling laws were for plumes in the non-rotating case, Speer & Marshall (1995) added the effect of rotation to build a more complete set of scaling laws. The Coriolis parameter $f=2{\it\Omega}$ , where ${\it\Omega}$ is the rotation rate, introduces a new time scale that conflicts with the buoyancy frequency $N$ . Speer & Marshall (1995) circumvented this issue of multiple time scales by dissecting the evolution of the plume into a series of events. For an initial value problem where the buoyancy flux starts at $t=0$ , they identified three stages in the development of the plume.

1. (i) First, $t<N^{-1}$ : the rising plume. At this stage, the classical scaling laws apply. The system does not feel the effect of rotation.

2. (ii) Second, $N^{-1}<t<f^{-1}$ : the lines of constant angular momentum are deflected from their original positions. As a consequence, a cyclone forms at the base of the plume and an anticyclone forms at the intrusion level.

3. (iii) Finally, $t>f^{-1}$ : the plume shape adjusts to remain in geostrophic balance. Ultimately, the anticyclone at the intrusion level becomes baroclinically unstable and moves away from the plume according to its own dynamics (Helfrich & Battisti 1991; Fernando, Chen & Ayotte 1998).

At each stage of the plume development, a collection of length scales can be obtained to describe the base and the top of the plume, but the general shape is thought to remain unchanged during the evolution (Whitehead, Marshall & Hufford 1996).

Although several studies have addressed the question of the baroclinic vortex, fewer studies have focused on the shape of the plume between the source and the intrusion level, and the dynamical effect of the cyclone at the base of the plume (Julien et al. 1999; Yamamoto, Cenedese & Caulfield 2011). In fact, in the case of pure jets, the presence of swirl strongly affects the dynamics and creates new patterns of turbulence (Liang & Maxworthy 2005). These patterns are precisely what we investigate herein by addressing the following questions. Are the classical convection patterns modified by background rotation? Can we quantify the feedback of the swirl speed on the plume dynamics?

To elucidate the impact of rotation on plume dynamics, we built a model of intermediate complexity. This model is derived from the original Navier–Stokes equations and takes advantage of the radially symmetric property of the plume (see Fabregat et al. 2016a ). This model can be used in two different ways: either to reproduce the turbulent flow using a time integration or to look for the steady states of the system using an appropriate parameterization for turbulence. Our study explored the latter. Because the axisymmetric equations are an azimuthal average of the dynamics, we use this set of equations to understand the mean flow rather than the detailed turbulence behaviour. The strategy is to compute branches of solutions by means of a continuation technique, and thus to obtain the steady states of the equations for various values of the parameters.

The paper is organized as follows. Section 2 below is a review of the original MTT model used to describe a convective plume. We first recall the main hypothesis and then derive the equation for the top of the plume and neutral level in the non-rotating case. They serve as a basis of our study. In § 3, we introduce the axisymmetric model used in the remainder of the analysis and briefly describe the method of continuation. In § 4, we use this model at the non-rotating limit and compare it with scaling laws. We introduce rotation in § 5 to show how it affects the shape of the solution. We present results from a preliminary three-dimensional numerical experiment in § 6. We discuss our conclusions in § 7.

2 Scaling laws for convective plumes

A convective plume is defined here as a buoyant jet in which the buoyancy is supplied steadily from a point source (Emanuel 1994). It is usually difficult to describe in detail the dynamics of such plumes as the plume triggers high levels of turbulence. An elegant way to describe the shape of the plume is to use simplified equations of the dynamics. Here, we briefly review the classical scaling laws following MTT’s analysis. This approach relies on the following five hypotheses.

1. (1) The flow is steady and axisymmetric.

2. (2) The mean radial velocity $u_{e}$ at the edge of the plume and at any height $z$ above the source is proportional to the mean vertical velocity at the centre of the plume $w(r=0)$ (entrainment hypothesis),

   (2.1) $$\begin{eqnarray}u_{e}={\it\alpha}w(r=0),\end{eqnarray}$$
   where ${\it\alpha}$ is a proportionality constant, also known as the entrainment constant or dilution rate. The value of ${\it\alpha}$ is usually determined experimentally, and the precise value remains a matter for debate. However, commonly accepted values vary from ${\it\alpha}=0.054$ for jets to ${\it\alpha}=0.083$ for plumes (Turner 1986), and correspond, in fact, to a parameterization of the turbulence.

3. (3) The radial profiles of mean vertical velocity and mean buoyancy are similar at all heights.

4. (4) The vertical perturbation pressure gradient is small compared with the vertical buoyancy acceleration.

5. (5) The flow is Boussinesq.

By means of these five hypotheses, the system of equations describing the plume can be simplified but still yield an accurate representation of the dynamics in the plume.

2.1 Original equations

Because of hypotheses (1), (4) and (5), the momentum, continuity and buoyancy equations in cylindrical coordinates simplify to

(2.2a ) $$\begin{eqnarray}\displaystyle u\frac{\partial w}{\partial r}+w\frac{\partial w}{\partial z} & = & \displaystyle B,\end{eqnarray}$$

(2.2b ) $$\begin{eqnarray}\displaystyle \frac{\partial ru}{\partial r}+\frac{\partial rw}{\partial z} & = & \displaystyle 0,\end{eqnarray}$$

(2.2c ) $$\begin{eqnarray}\displaystyle u\frac{\partial B}{\partial r}+w\frac{\partial B}{\partial z} & = & \displaystyle 0,\end{eqnarray}$$

where $r$ and $z$ are the radial and vertical coordinates, $u$ and $w$ are the radial and vertical velocities, and $B$ is the buoyancy,

(2.3) $$\begin{eqnarray}B=g\frac{{\it\rho}}{{\it\rho}_{0}},\end{eqnarray}$$

where ${\it\rho}$ is the density anomaly with respect to a stable density profile ${\it\rho}_{1}(z)$ , ${\it\rho}_{0}$ is the constant background density (independent of $z$ ) and $g$ is the acceleration due to gravity. We integrate (2.2a – c ) at each height $z$ , over a horizontal disc of radius $R(z)$ . We use hypothesis (3) with the additional constraint that the profiles of any quantity inside the plume are normally distributed,

(2.4a ) $$\begin{eqnarray}\displaystyle w(r,z)=w^{\prime }(z)\exp (-r^{2}/R^{2}), & & \displaystyle\end{eqnarray}$$

(2.4b ) $$\begin{eqnarray}\displaystyle B(r,z)=B^{\prime }(z)\exp (-r^{2}/R^{2}). & & \displaystyle\end{eqnarray}$$

For notational convenience, we omit the $^{\prime }$ from this point forward. If we use hypothesis (2) for the boundary condition on $u$ , then (2.2a – c ) become

(2.5a ) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}z}R^{2}w^{2} & = & \displaystyle 2R^{2}B,\end{eqnarray}$$

(2.5b ) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}z}R^{2}w & = & \displaystyle 2{\it\alpha}Rw,\end{eqnarray}$$

(2.5c ) $$\begin{eqnarray}\displaystyle \frac{\text{d}}{\text{d}z}R^{2}wB & = & \displaystyle -2R^{2}wN^{2},\end{eqnarray}$$

where

(2.6) $$\begin{eqnarray}N^{2}=-\frac{g}{{\it\rho}_{0}}\frac{\partial {\it\rho}_{1}}{\partial z}\end{eqnarray}$$

is the background buoyancy frequency.

To characterize the plume, Morton & Middleton (1973) introduced two quantities

(2.7) $$\begin{eqnarray}{\it\Gamma}=\frac{RB}{{\it\alpha}w^{2}}\end{eqnarray}$$

and

(2.8) $$\begin{eqnarray}{\it\Gamma}^{\prime }=\frac{B{\it\alpha}}{RN^{2}}.\end{eqnarray}$$

Hunt & Kaye (2005) associated ${\it\Gamma}(z=0)>1$ with lazy plumes and ${\it\Gamma}(z=0)<1$ with forced plumes. According to their definition, a lazy plume can be generated by slow upward emission of light fluid relative to the equivalent point source plume (relative excess of buoyancy). On the other hand, a forced plume can be generated by upward emission of light fluid with excess velocity relative to the point source plume (relative excess of momentum).

In an unstratified system ( $N^{2}=0$ ), equation (2.5c ) ensures that the buoyancy flux

(2.9) $$\begin{eqnarray}F=\frac{{\rm\pi}}{2}R^{2}wB\end{eqnarray}$$

through a horizontal section is constant at all heights $z$ . The solution of (2.5) is a simple scaling law for the vertical velocity and the radius of the plume,

(2.10a ) $$\begin{eqnarray}\displaystyle w & = & \displaystyle c_{1}F^{1/3}z^{-1/3},\end{eqnarray}$$

(2.10b ) $$\begin{eqnarray}\displaystyle R & = & \displaystyle c_{2}z,\end{eqnarray}$$

where $c_{1}$ and $c_{2}$ are constants of proportionality, which can be simply evaluated on assuming a particular value for the entrainment coefficient (see MTT). This result validates what would be obtained from straightforward dimensional analysis. For a point source release, the only important parameters are the buoyancy flux and the distance from the source. Given a buoyancy flux $F$ , one can compute the vertical velocity inside the plume at any height $z$ above this buoyancy source. Several experiments and numerical models have confirmed the validity of these scaling laws (Turner 1986; Woods 2010).

2.2 Dynamics in the stratified case

When the plume expands in a stratified medium, the dynamics is modified. As the plume rises from the source, it entrains relatively dense fluid and reaches a height where the plume is no longer buoyant, known as the level of neutral buoyancy, $H_{b}$ . The plume overshoots this level and falls back, around the rising flow, to the intrusion level, at which the plume expands laterally. The key parameters for the description of such a situation are the buoyancy flux at $z=0$ and the background stratification $N^{2}$ . If $N^{2}\neq 0$ , the scaling laws (2.10a,b ) are no longer solutions of (2.5a – c ), and only a numerical solution can be found. The buoyancy flux in (2.9) becomes a function of $z$ , and we denote the buoyancy flux at $z=0$ by $F_{0}\equiv F(z=0)$ . Henceforth, we use a subscript $0$ to denote the value of a variable at $z=0$ . MTT presented a generic solution of this system (see their figure 1). Near the buoyancy source, they recovered the solution of the unstratified system, and they demonstrated that the radius of the plume becomes infinite at the top of the plume. MTT defined two quantities to describe the plume:

1. (i) the top of the plume $H_{w}$ as the first height at which the vertical velocity vanishes (at $r=0$ ),

   (2.11) $$\begin{eqnarray}H_{w}=\left.z\right|_{w=0},\end{eqnarray}$$

2. (ii) the neutral level $H_{b}$ as the height at which the buoyancy vanishes (at $r=0$ ),

   (2.12) $$\begin{eqnarray}H_{b}=\left.z\right|_{B=0}.\end{eqnarray}$$

For an integrated model such as the one developed by MTT, $H_{w}$ and $H_{b}$ only describe the dynamics occurring at $r=0$ . In the axisymmetric model that we will introduce in § 3, these quantities may no longer be associated with the top of the plume and the neutral level.

We also define the overshoot as the difference ( $H_{w}-H_{b}$ ). In practice, as a plume overshoots the height $H_{b}$ and rises to $H_{w}$ , it entrains and mixes with (less dense) ambient fluid. The fluid at the top then falls back, again mixing with both ambient fluid and the rising core of the flow. The plume will only ever intrude at $H_{b}$ in the case that there is no mixing above the height $H_{b}$ ; in all other cases it will intrude at a height somewhere between $H_{b}$ and $H_{w}$ . We introduce a formal definition of the intrusion level in § 4. MTT also argued that $H_{w}$ and $H_{b}$ should obey the scaling law

(2.13) $$\begin{eqnarray}H_{s}=c\left(\frac{F_{0}}{N^{3}}\right)^{1/4},\end{eqnarray}$$

with $c$ a constant. A value of $c=3.8$ has often been proposed in order to match results obtained in laboratory experiments (Turner 1986).

From this point forward, we work with a non-dimensional form of (2.5a – c ) which is obtained on taking $L=500$  m as characteristic of the dominant length scale for oceanic convection and the corresponding characteristic velocity scale as $U=1$  m s $^{-1}$ . In order to assess the validity of this scaling law, we used a numerical solver to find the solutions of (2.5a )–(2.5c ) for several values of $F_{0}$ . It should be noted that we need $w_{0}\neq 0$ and $R_{0}\neq 0$ to obtain a non-zero solution (throughout we consider the virtual origin offset, Hunt & Kaye (2001), to be small compared with $L$ , and hence we neglect the effects of any virtual origin offset). In this example, we arbitrarily chose $w_{0}=1.0$ , $R_{0}=0.02$ , $N^{2}=5$ , and varied $B_{0}$ across several orders of magnitude. We take ${\it\alpha}=0.083$ (Turner 1986). In figure 1 we plot $H_{w}$ , $H_{b}$ and the scaling law (2.13) as a function $F_{0}$ in log–log space. There is good agreement between $H_{w}$ , $H_{b}$ and the scaling law for lazy plumes ( ${\it\Gamma}_{0}>1$ ): the slope of both $H_{w}$ and $H_{b}$ is 0.25. However, we observe a deviation for forced plumes ( ${\it\Gamma}_{0}<1$ ): $H_{w}$ is constant and $H_{b}$ varies linearly with $F_{0}$ (as indicated by the green dashed line). We discuss in detail the solutions of MTT’s model for forced plumes in appendix A. It should be noted that in our results (§ 4), the ratio ${\it\Gamma}_{0}/{\it\Gamma}_{0}^{\prime }=0.3$ is constant since we only vary $B_{0}$ . We repeated this experiment for various values of ${\it\Gamma}_{0}/{\it\Gamma}_{0}^{\prime }$ by varying the stratification, and we obtained a qualitatively similar curve. In the absence of complete field data and an apparent lack of experimental data, we are forced to use these scaling laws ((2.13) and appendix A) to validate our axisymmetric model in the following section.

Figure 1. The levels $H_{b}$ (green with circles) and $H_{w}$ (blue with crosses) as a function of the buoyancy flux $F_{0}$ in non-dimensional units. The black line is the scaling law (2.13). The blue dashed line (slope  $=$  0) is (A 11) and the green dashed line (slope  $=$  1) is (A 4).

3 The axisymmetric model

While MTT’s model and the scaling laws are extremely attractive because of their simplicity, we would like to study richer dynamics. We introduce a system of equations to study an axisymmetric plume, and abandon the entrainment hypothesis and the self-similarity hypothesis ((2) and (3) respectively). In fact, instead of parameterizing the entrainment, we use an advanced turbulence closure scheme. The ultimate goal is then to find the steady states of this set of equations.

3.1 Main equations

The non-hydrostatic Boussinesq equations for a radially symmetric flow on an $f$ -plane are

(3.1a ) $$\begin{eqnarray}\displaystyle \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}+w\frac{\partial u}{\partial z}-\frac{v^{2}}{r}-fv & = & \displaystyle -\frac{\partial p}{\partial r}+D_{u},\end{eqnarray}$$

(3.1b ) $$\begin{eqnarray}\displaystyle \frac{\partial v}{\partial t}+u\frac{\partial v}{\partial r}+w\frac{\partial v}{\partial z}+\frac{uv}{r}+fu & = & \displaystyle D_{v},\end{eqnarray}$$

(3.1c ) $$\begin{eqnarray}\displaystyle \frac{\partial w}{\partial t}+u\frac{\partial w}{\partial r}+w\frac{\partial w}{\partial z} & = & \displaystyle -\frac{\partial p}{\partial z}-\frac{{\it\rho}}{{\it\rho}_{0}}g+D_{w},\end{eqnarray}$$

(3.1d ) $$\begin{eqnarray}\displaystyle \frac{\partial T}{\partial t}+u\frac{\partial T}{\partial r}+w\frac{\partial T}{\partial z} & = & \displaystyle D_{T},\end{eqnarray}$$

(3.1e ) $$\begin{eqnarray}\displaystyle \frac{\partial ru}{\partial r}+\frac{\partial rw}{\partial z} & = & \displaystyle 0,\end{eqnarray}$$

where $T$ is the temperature, related to the density by the equation of state

(3.2) $$\begin{eqnarray}{\it\rho}={\it\rho}_{0}{\it\alpha}_{T}T,\end{eqnarray}$$

with ${\it\alpha}_{T}$ the thermal expansion coefficient. In (3.1a – d ) $D_{u,v,w,T}$ are dissipative terms described in § 3.2. The equation for $v$ , the swirl velocity, can be rewritten as

(3.3) $$\begin{eqnarray}\frac{\partial {\it\lambda}}{\partial t}+u\frac{\partial {\it\lambda}}{\partial r}+w\frac{\partial {\it\lambda}}{\partial z}=D_{{\it\lambda}},\end{eqnarray}$$

where

(3.4) $$\begin{eqnarray}{\it\lambda}=rv+\frac{f}{2}r^{2}\end{eqnarray}$$

is the angular momentum of the fluid. The continuity equation is solved by means of a stream function

(3.5a ) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\psi}}{\partial r} & = & \displaystyle rw,\end{eqnarray}$$

(3.5b ) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\psi}}{\partial z} & = & \displaystyle -ru.\end{eqnarray}$$

We also define the vorticity

(3.6) $$\begin{eqnarray}{\it\zeta}=\frac{\partial w}{\partial r}-\frac{\partial u}{\partial z},\end{eqnarray}$$

or in terms of the stream function

(3.7) $$\begin{eqnarray}{\it\zeta}=\frac{1}{r}\left(\frac{\partial ^{2}{\it\psi}}{\partial r^{2}}+\frac{\partial ^{2}{\it\psi}}{\partial z^{2}}-\frac{1}{r}\frac{\partial {\it\psi}}{\partial r}\right).\end{eqnarray}$$

Using these definitions, we can rewrite the original system as

(3.8a ) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\zeta}}{\partial t} & = & \displaystyle -J\left({\it\psi},\frac{{\it\zeta}}{r}\right)-\frac{1}{r^{3}}\frac{\partial {\it\lambda}^{2}}{\partial z}+g{\it\alpha}_{T}\frac{\partial T}{\partial r}+D_{{\it\zeta}},\end{eqnarray}$$

(3.8b ) $$\begin{eqnarray}\displaystyle \frac{\partial {\it\lambda}}{\partial t} & = & \displaystyle -\frac{1}{r}J({\it\psi},{\it\lambda})+D_{{\it\lambda}},\end{eqnarray}$$

(3.8c ) $$\begin{eqnarray}\displaystyle \frac{\partial T}{\partial t} & = & \displaystyle -\frac{1}{r}J({\it\psi},T)+D_{T},\end{eqnarray}$$

with $J(a,b)=(\partial a/\partial r)(\partial b/\partial z)-(\partial a/\partial z)(\partial b/\partial r)$ . We will consider these equations in a cylinder of radial extent $L_{r}$ and vertical extent $L_{z}$ . The boundary conditions for this system are given in table 1.

Table 1. Boundary conditions for the system of (3.8).

We set ${\it\psi}={\it\psi}_{bc}$ and ${\it\zeta}=0$ at $z=0$ , $z=L_{z}$ and $r=0$ , which is equivalent to free-slip boundary conditions. The constant ${\it\psi}_{bc}$ is adjusted to have a non-zero injection velocity through the bottom boundary,

(3.9) $$\begin{eqnarray}w=w_{0}\exp \left(-\frac{r^{2}}{R_{0}^{2}}\right),\end{eqnarray}$$

and no flux at $z=L_{z}$ and $r=0$ . At $r=L_{r}$ , we consider an open boundary with zero vertical velocity: $\partial _{n}{\it\psi}=0$ (where $\partial _{n}$ is the derivative in the direction normal to the boundary).

For the swirl speed, we have $v=0$ at $r=0$ (by construction), and we opt for the free-slip boundary condition: $\partial _{n}v=0$ on all other boundaries. This translates into ${\it\lambda}=0$ at $r=0$ , $\partial _{n}{\it\lambda}=0$ at the bottom, and top, and the mixed boundary condition $\partial _{n}{\it\lambda}=({\it\lambda}+fr^{2}/2)/r$ at $r=L_{r}$ . Interestingly, the effect of rotation only appears via this boundary condition. However, the joint effect of this boundary condition and the diffusion operator $D_{{\it\lambda}}$ in fact implies that the angular momentum field at rest in the domain is ${\it\lambda}=fr^{2}/2$ .

We use two different boundary conditions for the temperature at the top and bottom in order to maintain a stratification $N^{2}$ . At the top we set the temperature to $T=T^{top}$ (Dirichlet), and at the bottom we use a no-flux boundary condition.

The buoyancy flux is applied via a relaxation at the first grid point $(T-T^{forc})/{\it\tau}$ , with ${\it\tau}$ an arbitrarily small time constant and $T^{forc}$ a fixed reference temperature profile. This choice is necessary to avoid a modulation of the buoyancy flux due to the turbulence scheme (which caries a variable thermal diffusivity). The prescribed temperature profile is

(3.10) $$\begin{eqnarray}T^{forc}=T_{0}\exp \left(-\frac{r^{2}}{R_{0}^{2}}\right),\end{eqnarray}$$

with $T_{0}$ a constant and $R_{0}$ a characteristic length for the radial extension of the prescribed profile. We also apply a no-flux boundary condition at $r=0$ and $r=L_{r}$ ( $\partial _{n}T=0$ ). The buoyancy flux driving the plume is given by

(3.11) $$\begin{eqnarray}F_{0}=\frac{{\rm\pi}}{2}w_{0}R_{0}^{2}g{\it\alpha}_{T}T_{0}+2{\rm\pi}g{\it\alpha}_{T}{\rm\Delta}z\int _{0}^{R_{0}}\frac{1}{{\it\tau}}(T-T^{forc})r\,\text{d}r,\end{eqnarray}$$

where the first term is the classical buoyancy flux (see (2.9)) and the second term is a pseudo-diffusive flux, with ${\rm\Delta}z$ the vertical extension of the first grid point.

3.2 Turbulence closure

The idea behind a steady-state axisymmetric model of a convective plume is to let the turbulence closure represent the mean effect of all the eddies present in the real plume. In this study, we use a simplified version of the Smagorinsky (1963) closure scheme, used in Rotunno & Emanuel (1987) to study tropical cyclones. The dissipative terms in (3.1a – d ) are computed via a Reynolds decomposition of the axisymmetric system: each variable $X$ is written as the sum of an azimuthally averaged component $\overline{X}$ plus a departure $X^{\prime }$ . The resulting dissipative terms are

(3.12a ) $$\begin{eqnarray}\displaystyle D_{u} & = & \displaystyle \frac{1}{r}\frac{\partial r\overline{u^{\prime }u^{\prime }}}{\partial r}+\frac{\partial \overline{u^{\prime }w^{\prime }}}{\partial z}-\frac{\overline{v^{\prime }v^{\prime }}}{r},\end{eqnarray}$$

(3.12b ) $$\begin{eqnarray}\displaystyle D_{v} & = & \displaystyle \frac{1}{r^{2}}\frac{\partial r^{2}\overline{u^{\prime }v^{\prime }}}{\partial r}+\frac{\partial \overline{v^{\prime }w^{\prime }}}{\partial z},\end{eqnarray}$$

(3.12c ) $$\begin{eqnarray}\displaystyle D_{w} & = & \displaystyle \frac{1}{r}\frac{\partial r\overline{u^{\prime }w^{\prime }}}{\partial r}+\frac{\partial \overline{w^{\prime }w^{\prime }}}{\partial z},\end{eqnarray}$$

(3.12d ) $$\begin{eqnarray}\displaystyle D_{T} & = & \displaystyle \frac{1}{r}\frac{\partial r\overline{u^{\prime }T^{\prime }}}{\partial r}+\frac{\partial \overline{w^{\prime }T^{\prime }}}{\partial z}.\end{eqnarray}$$

The Reynolds stresses are computed with a linear eddy viscosity model,

(3.13a ) $$\begin{eqnarray}\displaystyle \overline{u^{\prime }u^{\prime }} & = & \displaystyle 2{\it\nu}\frac{\partial u}{\partial r},\end{eqnarray}$$

(3.13b ) $$\begin{eqnarray}\displaystyle \overline{u^{\prime }w^{\prime }} & = & \displaystyle {\it\nu}\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right),\end{eqnarray}$$

(3.13c ) $$\begin{eqnarray}\displaystyle \overline{v^{\prime }v^{\prime }} & = & \displaystyle 2{\it\nu}\frac{u}{r},\end{eqnarray}$$

(3.13d ) $$\begin{eqnarray}\displaystyle \overline{u^{\prime }v^{\prime }} & = & \displaystyle {\it\nu}r\frac{\partial }{\partial r}\left(\frac{v}{r}\right),\end{eqnarray}$$

(3.13e ) $$\begin{eqnarray}\displaystyle \overline{v^{\prime }w^{\prime }} & = & \displaystyle {\it\nu}\frac{\partial v}{\partial z},\end{eqnarray}$$

(3.13f ) $$\begin{eqnarray}\displaystyle \overline{w^{\prime }w^{\prime }} & = & \displaystyle 2{\it\nu}\frac{\partial w}{\partial z},\end{eqnarray}$$

(3.13g ) $$\begin{eqnarray}\displaystyle \overline{u^{\prime }T^{\prime }} & = & \displaystyle {\it\nu}\frac{\partial T}{\partial r},\end{eqnarray}$$

(3.13h ) $$\begin{eqnarray}\displaystyle \overline{w^{\prime }T^{\prime }} & = & \displaystyle {\it\nu}\frac{\partial T}{\partial z},\end{eqnarray}$$

where we have neglected the isotropic part of the stress tensor (Kundu & Cohen 2008) and used the same parameter ${\it\nu}$ for thermal diffusivity and momentum kinematic viscosity. Henceforth, we also neglect all spatial derivatives of ${\it\nu}$ : these terms mimic the effect of an additional advection operator. We verified that the effect of this operator is negligible compared with the actual advective terms. We use the Smagorinsky (1963) formulation for  ${\it\nu}$ ,

(3.14) $$\begin{eqnarray}{\it\nu}=l_{s}^{2}S+{\it\nu}_{0},\end{eqnarray}$$

where ${\it\nu}_{0}$ is a background kinematic viscosity, $S$ is the deformation given by

(3.15) $$\begin{eqnarray}S^{2}=2\left[\left(\frac{\partial u}{\partial r}\right)^{2}+\left(\frac{u}{r}\right)^{2}+\left(\frac{\partial w}{\partial z}\right)^{2}\right]+\left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}\right)^{2}+\left(\frac{\partial v}{\partial r}-\frac{v}{r}\right)^{2}+\left(\frac{\partial v}{\partial z}\right)^{2}\end{eqnarray}$$

and $l_{s}$ is the Smagorinsky length scale, usually chosen to be of the order of the grid size to represent subgrid-scale processes. In our case, this turbulence closure is more than a subgrid-scale model. Since it parameterizes all eddy activity, $l_{s}$ should then be equal to the typical size of eddies. However, since we have no a priori value for $l_{s}$ , we will use MTT’s model to calibrate $l_{s}$ . The final expressions of the dissipative terms are

(3.16a ) $$\begin{eqnarray}\displaystyle D_{{\it\zeta}} & = & \displaystyle {\it\nu}\left(\frac{\partial ^{2}{\it\zeta}}{\partial r^{2}}+\frac{\partial ^{2}{\it\zeta}}{\partial z^{2}}+\frac{1}{r}\frac{\partial {\it\zeta}}{\partial r}-\frac{{\it\zeta}}{r^{2}}\right),\end{eqnarray}$$

(3.16b ) $$\begin{eqnarray}\displaystyle D_{{\it\lambda}} & = & \displaystyle {\it\nu}\left(\frac{\partial ^{2}{\it\lambda}}{\partial r^{2}}+\frac{\partial ^{2}{\it\lambda}}{\partial z^{2}}-\frac{1}{r}\frac{\partial {\it\lambda}}{\partial r}\right),\end{eqnarray}$$

(3.16c ) $$\begin{eqnarray}\displaystyle D_{T} & = & \displaystyle {\it\nu}\left(\frac{\partial ^{2}T}{\partial r^{2}}+\frac{\partial ^{2}T}{\partial z^{2}}+\frac{1}{r}\frac{\partial T}{\partial r}\right).\end{eqnarray}$$

3.3 Non-dimensional parameters

We make the system of equations non-dimensional using a reference length $L$ , a reference velocity $U$ and a reference temperature $T$ , and use typical values relevant to oceanic convection: $L=500$  m, $U=1$  m s $^{-1}$ and $T=1\,^{\circ }\text{C}$ . From this point forward, all variables are in their non-dimensional form. The non-dimensional parameters that are relevant to this system of equations are the Rossby number

(3.17) $$\begin{eqnarray}Ro=\frac{U}{fL}\end{eqnarray}$$

and the Froude number

(3.18) $$\begin{eqnarray}F_{r}=\frac{U}{NL},\end{eqnarray}$$

which is proportional to $\sqrt{{\it\Gamma}^{\prime }/{\it\Gamma}}$ . For oceanic convection, we have $Ro=20$ and $F_{r}=15$ (Speer & Marshall 1995). In the absence of stratification (e.g. well mixed shallow areas), $F_{r}\rightarrow \infty$ . The plume is lighter than its surroundings and rises freely (cf. the scaling law in § 2) until it eventually encounters the upper boundary. Since the interaction of the plume with the boundary is beyond the scope of this study, we use the Froude number (stratification) as a way to stop the upward motion. In the numerical formulation, the background Reynolds number $Re=UL/{\it\nu}_{0}$ is only meaningful in the absence of any motion, so we use $l_{s}$ to characterize the level of turbulence. We now review the methodology used to analyse this set of equations.

3.4 Methodology: bifurcation analysis

We use a bifurcation analysis to find the steady states of the system of (3.8) as a function of the external parameters. We briefly recall the main steps of this methodology. We first write (3.8) in a compact form,

(3.19) $$\begin{eqnarray}\unicode[STIX]{x1D648}\frac{\partial \boldsymbol{x}}{\partial t}=\boldsymbol{F}(\boldsymbol{x},p),\end{eqnarray}$$

with $\boldsymbol{x}\equiv ({\it\psi},{\it\lambda},T)$ the state vector, $\boldsymbol{F}$ a map of $\mathbb{R}^{n}\times \mathbb{R}\rightarrow \mathbb{R}^{n}$ (right-hand side of (3.8)), $p$ an adjustable parameter and $\unicode[STIX]{x1D648}$ a weight matrix (stream function–vorticity transformation matrix; (3.7)). We define the Jacobian

(3.20a ) $$\begin{eqnarray}\displaystyle \text{J}(\boldsymbol{x}_{0},p_{0}) & = & \displaystyle \left.\frac{\partial \boldsymbol{F}(\boldsymbol{x},p)}{\partial \boldsymbol{x}}\right|_{\boldsymbol{x}_{0},p_{0}},\end{eqnarray}$$

(3.20b ) $$\begin{eqnarray}\displaystyle \boldsymbol{k}(\boldsymbol{x}_{0},p_{0}) & = & \displaystyle \left.\frac{\partial \boldsymbol{F}(\boldsymbol{x},p)}{\partial p}\right|_{\boldsymbol{x}_{0},p_{0}}.\end{eqnarray}$$

We look for the steady-state solutions such that

(3.21) $$\begin{eqnarray}\boldsymbol{F}(\boldsymbol{x},p)=0.\end{eqnarray}$$

This steady-state search is performed using the pseudo-arclength continuation method (Keller 1977), briefly described as follows. The implicit function theorem states that the existence of a solution $\boldsymbol{x}_{p_{0}}$ at $p_{0}$ , together with differentiability of $\boldsymbol{F}$ and non-singularity of the Jacobian $\text{J}$ , will imply the existence of $\boldsymbol{x}_{p}$ for $p$ near $p_{0}$ , and also imply that $\boldsymbol{x}_{p}$ is a smooth function of  $p$ .

In practice, $\text{J}$ might be singular for certain values of $p$ . The singularity disappears if we introduce an additional variable $s$ and solve the augmented system

(3.22a ) $$\begin{eqnarray}\displaystyle \boldsymbol{F}(\boldsymbol{x}(s),p(s)) & = & \displaystyle 0,\end{eqnarray}$$

(3.22b ) $$\begin{eqnarray}\displaystyle n(\boldsymbol{x}(s),p(s),s) & = & \displaystyle 0.\end{eqnarray}$$

We find a solution of this system with Newton iterations. The second equation is the normalization of the tangent vector with respect to the parameter

(3.23) $$\begin{eqnarray}{|\boldsymbol{x}^{\prime }(s)|}^{2}+{|p^{\prime }(s)|}^{2}=1,\end{eqnarray}$$

where the prime $(\text{}^{\prime })$ denotes the derivative with respect to the auxiliary variable $s$ . We discretize this equation as

(3.24) $$\begin{eqnarray}\frac{\boldsymbol{x}(s_{1})-\boldsymbol{x}(s_{0})}{{\rm\Delta}s}\boldsymbol{\cdot }\boldsymbol{x}^{\prime }(s_{0})+\frac{p(s_{1})-p(s_{0})}{{\rm\Delta}s}p^{\prime }(s_{0})=1,\end{eqnarray}$$

with ${\rm\Delta}s=s_{1}-s_{0}$ . Hence, one can write (3.22b ) as

(3.25) $$\begin{eqnarray}n(\boldsymbol{x}(s),p(s),s)=(\boldsymbol{x}(s_{1})-\boldsymbol{x}(s_{0}))\boldsymbol{\cdot }\boldsymbol{x}^{\prime }(s_{0})+(p(s_{1})-p(s_{0}))p^{\prime }(s_{0})-{\rm\Delta}s=0.\end{eqnarray}$$

Given $(\boldsymbol{x}_{0},p_{0},s_{0})$ , a solution of (3.22), we use a predictor–corrector scheme to find a nearby solution $(\boldsymbol{x}_{1},p_{1},s_{1})$ . The predictor is computed using a Euler scheme,

(3.26) $$\begin{eqnarray}(\boldsymbol{x}_{1},p_{1},s_{1})=(x_{0},p_{0},s_{0})+{\rm\Delta}s(\boldsymbol{x}_{0}^{\prime },p_{0}^{\prime },1),\end{eqnarray}$$

and with $x_{0}^{\prime }$ and $p_{0}^{\prime }$ , which are solutions of

(3.27a ) $$\begin{eqnarray}\displaystyle \text{J}(\boldsymbol{x}_{0},p_{0})\boldsymbol{x}_{0}^{\prime }+\boldsymbol{k}(\boldsymbol{x}_{0},p_{0})p_{0}^{\prime } & = & \displaystyle 0,\end{eqnarray}$$

(3.27b ) $$\begin{eqnarray}\displaystyle |\boldsymbol{x}_{0}^{\prime }|^{2}+|p_{0}^{\prime }|^{2} & = & \displaystyle 1.\end{eqnarray}$$

3.5 Numerical formulation

The system of (3.8) plus boundary conditions (table 1) is discretized with finite differences on a regular A-grid (stream function, vorticity and temperature are collocated on the same point). The lateral and vertical extensions of the domain are $L_{r}=1$ , $L_{z}=2$ in the experiments without rotation, and $L_{r}=2$ , $L_{z}=2$ when rotation is present. The grid step is ${\rm\Delta}r={\rm\Delta}z=0.005$ in all cases. The Jacobian operator (advection term) is discretized with an energy–enstrophy conserving form (Arakawa & Lamb 1977). We carried out the bifurcation analysis with the large-scale library LOCA (Salinger et al. 2002), which is part of the trilinos software (Heroux et al. 2005). We also performed several time integrations of the model to examine the evolution of a passive tracer. In this case, the vorticity/stream function matrix (3.7) was inverted using an LU decomposition (Amestoy et al. 2006). We used a second-order Adams–Bashforth time stepping with a time step of ${\rm\Delta}t=10^{-4}$ .

4 Non-rotating limit ( $\text{}f=0$ )

We begin our analysis with the solution of (3.8a – c ) at the non-rotating limit ( $f=0$ ). In this case, ${\it\lambda}=v=0$ is a solution of (3.8b ). No forcing term can increase the swirl speed. The remaining vorticity equation is strictly two-dimensional (2D), and the buoyancy is advected by a 2D velocity field. We use this configuration as a starting point of our study. We describe the shape of the steady state and compare the solution with the solution of MTT’s model (as reasoned in § 2.2).

4.1 General behaviour

In figure 2, we plot a typical steady-state solution of the model; more precisely, we plot the non-dimensional stream function ${\it\psi}$ and temperature in an azimuthal section $(r,z)$ . The destabilizing source buoyancy flux (3.11) is obtained with $R_{0}=0.02$ , $T_{0}=200$ and $w_{0}=1$ , such that $F_{0}=0.1$ . The buoyancy flux destabilizes the flow from below and induces a local vertical velocity. Between $z=0$ and $z=0.6$ , the plume entrains the surrounding fluid (horizontal stream lines coming from the far edge). The core of the plume overshoots above the neutral level $H_{b}$ and reaches the height $H_{w}$ , where the plume is negatively buoyant. It falls back on the side (negative vertical velocity), and the flow exits through the far edge (open) boundary at the intrusion level. We define the intrusion level as

(4.1) $$\begin{eqnarray}H_{i}=\arg \!\max (u_{r\infty }(z)),\end{eqnarray}$$

the height at which the outward radial velocity is maximum away from the core of the plume (where $w=\partial {\it\psi}/\partial r\simeq 0$ ). In our model, we measure $H_{i}$ at the open boundary $r=L_{r}$ .

Figure 2. (a) Stream function ${\it\psi}$ and (b) temperature $T$ of a typical steady state for $F_{0}=0.1$ . The levels $H_{w}$ and $H_{b}$ (green dotted lines) are computed in the axisymmetric model at the first grid point and $r=0$ respectively (see next section for details). The level $H_{i}$ (green dotted line) is computed at $r=1$ . The contour interval is $5\times 10^{-4}$ in (a) and $1.0$ in (b). Negative contours are plotted with dashed lines. The thick red line in (b) is a contour level of a passive tracer released at the source of the plume and advected with the flow shown in (a). The passive tracer reaches $r=1$ at $t=28$ . The thick blue line in (a) is $R(z)$ computed with the original MTT model (2.5).

At $z=H_{i}$ , the stream lines are exactly horizontal, which means that the plume expands to infinity. We changed the size of the domain to confirm that the dynamics of the plume was not significantly impacted by the presence of the boundary (not shown).

The turbulence scheme plays a crucial role in the establishment of the plume pattern. The Smagorinsky scheme naturally increases the thermal diffusivity and kinematic viscosity (diffusion of momentum) in places where the vorticity and shear are high. In the present case, we used $l_{s}=0.01$ (see next section). The deformation $S$ (see (3.15)), which controls the diffusivity locally, is maximum at the base and at the top of the plume. The value of the background Reynolds number (5000 in the present and subsequent experiments) has a negligible effect on the shape of the plume.

We obtain similar results (stream function and temperature profiles) when we increase the buoyancy flux, while we remain in the lazy plume regime, as long as the top of the plume does not interact with the upper boundary. In the forced plume scenario (not shown), the plume rises from the bottom boundary with negligible buoyancy flux, reaches its maximum height ( $H_{w}$ ) and falls back to (or very close to) the bottom boundary.

The red line in figure 2 corresponds to a contour level of a passive tracer steadily released at the base of the plume and advected by the steady-state velocity field. We stopped the time integration when the tracer reached the far edge boundary so that the tracer envelope revealed the structure of an intermediate stage (if we carry the time integration further, we reach the steady state where the red lines flatten outside the core of the plume and closely follow the temperature lines). We recognize the classical shape of a convective plume with a rising cone, an overshoot and an intrusion level. For reference, plotted within figure 2 is the plume width $R(z)$ (blue line) computed using MTT’s equations (2.5) and with the same parameters ( $R_{0}$ , $B_{0}$ and $w_{0}$ ). It should be noted that in our derivation of MTT’s model (§ 2), we used $R(z)$ as a characteristic width of the plume (cf. (2.4)) and not as a measure of the boundary of the plume. We are confident that the top of the plume is better represented in the present axisymmetric model than in MTT’s model because we are able to compute the downward vertical velocity in the overshoot region (which is impossible in MTT’s model).

Figure 3. The levels $H_{w}$ (blue with crosses) and $H_{b}$ (green with circles) in the axisymmetric model (solid lines) and in MTT’s model (dashed lines). The solid black line is $F_{0}^{1/4}$ . The level $H_{i}$ in the axisymmetric model (red line) is slightly above, but almost indistinguishable from, the green line ( $H_{b}$ ). The buoyancy flux $F_{0}$ varies via the bottom boundary condition $T_{0}$ . The parameters at the inlet are $w_{0}=0.1$ and $R_{0}=0.02$ .

4.2 Comparison with MTT’s model

Since the axisymmetric model involves more physics than the scaling laws, it is interesting to compare the two models. According to MTT’s model, $H_{w}$ is expected to vary with $F_{0}^{1/4}$ for lazy plumes (2.13), and is expected to be constant for forced plumes (A 10)–(A 11). To verify these power laws, we vary the buoyancy flux and measure $H_{w}$ and $H_{b}$ in the model. We compute $H_{w}$ using the vertical velocity calculated at the first grid point in the radial direction to avoid the singularity at $r=0$ . We compute $H_{b}$ with the temperature at $r=0$ , calculated with the first-order no-flux boundary condition

(4.2) $$\begin{eqnarray}T(r=0)={\textstyle \frac{4}{3}}T(r={\rm\Delta}r)-{\textstyle \frac{1}{3}}T(r=2{\rm\Delta}r).\end{eqnarray}$$

Among the three possible ways to increase the buoyancy flux, we chose to increase the prescribed temperature at the boundary, $T_{0}$ , the two other ways being either to increase the vertical velocity at the inlet, $w_{0}$ , or to increase the radius, $R_{0}$ . In figure 3, we plot $H_{w}$ , $H_{b}$ and $H_{i}$ as a function of the buoyancy flux $F_{0}$ with solid lines (axisymmetric model) and dashed lines (MTT). We only plot $H_{w}$ and $H_{b}$ that are higher than $2{\rm\Delta}z$ and lower than $L_{z}$ , where ${\rm\Delta}z$ is the vertical resolution in the axisymmetric model. The black line is given as a reference ((2.13)). We recover the two regimes (forced and lazy plume), and the transition occurs at the same value of ${\it\Gamma}_{0}=1$ . For the chosen value of the Smagorinsky length scale ( $l_{s}=0.01$ ), the power law matches the prediction of the scaling law (slope of 0.25) for ${\it\Gamma}_{0}>1$ . This slope increases for lower values of $l_{s}$ and decreases for higher values of $l_{s}$ . This good correspondence between MTT’s model and the present model could not be achieved without the Smagorinsky turbulent closure. We were not able to prove the linear relationship between $H_{b}$ and $F_{0}$ for ${\it\Gamma}_{0}\ll 1$ (A 4) due to a lack of vertical resolution in our model.

Further, as pointed out by a reviewer, many of the differences between the curves from MTT and from the axisymmetric model come down to the fact that the MTT model does not allow for any interaction between dense fluid falling from $H_{w}$ . For example, the MTT model overpredicts the top of the plume for ${\it\Gamma}_{0}<1$ , presumably because in a real plume the excess momentum carries the flow far above the intrusion height and the plume mixes with falling fluid which would lower the top of the plume. In addition, the mixing above $H_{b}$ of lighter ambient fluid means that the intrusion height is higher than $H_{b}$ , which is what we obtain in figure 3.

For ${\it\Gamma}_{0}>1$ , the agreement between our axisymmetric model and the scaling laws (MTT and appendix A) is far better. This is to be expected since for these plumes the flow has little momentum at $H_{b}$ and therefore only rises slightly higher; thereby the mixing above $H_{b}$ (not parameterized in the scaling laws) is minimal.

Due to its computational cost, we could not establish the stability of every single steady state in the usual way (using an eigenvalue solver). However, we performed several time integrations of the system with the steady state as the initial condition. All of the steady states that we tested (either forced or lazy plumes) were stable.

Moreover, the open boundary at the far edge is a crucial element to get the right behaviour of the model. It is still possible to perform a bifurcation analysis in a closed domain, but the recirculation near the far edge boundary produces non-trivial deviation from the scaling law (not shown).

To summarize, we have described the dynamics of a convective plume in the absence of rotation. Typical flow patterns correspond to a rapid ascent near the centre with an eventual overshoot. The radial flow is directed towards the centre of the plume near its base and outward at the intrusion level. This investigation serves as a reference case with which we will compare solutions in the presence of rotation.

5 The impact of rotation

When $f\neq 0$ , the system has richer dynamics. At rest (quiescent ambient in solid-body rotation), the angular momentum is ${\it\lambda}=0.5fr^{2}$ . If we apply a flow pattern similar to the one presented in the previous section, we can expect a convergence of angular momentum near the centre at the bottom and a divergence at the intrusion level. This corresponds to cyclonic and anticyclonic flow respectively (for $f>0$ ). This non-trivial swirl speed pattern will in turn affect the azimuthal vorticity (see (3.8)).

5.1 Increasing the buoyancy flux

We compute the steady-state solutions in a similar configuration to the reference case (cf. figure 2), but with $Ro=20$ . As before, we increase the buoyancy flux $F_{0}$ and look for the steady-state solutions. The curve showing $H_{w}$ as a function of $F_{0}$ is almost indistinguishable in the rotating axisymmetric and non-rotating axisymmetric cases (figure 3). The level $H_{b}$ is also identical in the forced plume regime; in the lazy plume regime, the curve $H_{b}(F_{0})$ has a slope of 0.25 but is closer to $H_{w}$ (not shown).

Figure 4. Steady state for $Ro=20$ and $F_{0}=0.1$ . (a) Stream function, ${\it\psi}$ , (b) temperature (black), $T$ , envelope of the plume from the time integration of a passive tracer (thick red), (c) angular momentum, ${\it\lambda}$ , (d) swirl velocity (shaded area is negative), $v$ . The passive tracer reaches $r=1$ at $t=55\;(\simeq 3f^{-1})$ .

We plot in figure 4 the stream function, temperature, angular momentum and swirl velocity (simply derived from (3.4)) of the steady state computed with $F_{0}=0.1$ and $Ro=20$ . This computation was performed in a larger domain ( $L_{z}=2$ and $L_{r}=2$ ) to minimize the effect of the far-field open boundary. Some tests with even wider domains obtained similar solutions. The plots are truncated in the radial direction for visualization purposes.

Figure 4(c) illustrates how the isolines of angular momentum are deflected under the action of the plume; an undisturbed profile would exhibit vertical lines ( ${\it\lambda}=0.5fr^{2}$ ). The presence of the buoyant source squeezes these lines towards the centre ( $r=0$ ) at the bottom. In other words, it corresponds to the entrainment of ambient fluid coming from larger radial locations, i.e. with a relative excess of angular momentum when it is in the plume. This convergence of angular momentum corresponds to a maximum swirl speed of $v=0.2$ (as shown in figure 4 d), a value that is not negligible compared with the vertical velocity at the source ( $w_{0}=1$ ). At the intrusion level, the deflection of the isolines of angular momentum is less pronounced than at the base of the plume. At this level, the order of magnitude for the swirl speed is $v\simeq -0.002$ (shaded area in figure 4 d). In this experiment, the Rossby number controls the magnitude of the swirl speed. A less obvious consequence of rotation is that the stream function is squeezed between the neutral level and the top of the plume so that the vertical thickness of the plume at the intrusion level is smaller in the rotating case. This feature is also visible when comparing the two envelopes (red lines in figures 4 b and 2 b).

A striking feature in the rotating case is that the circulation in the $(r,z)$ -plane generated by the plume is laterally confined: the far-field stream lines are no longer horizontal as they were in the non-rotating case. The tilt of these stream lines indicates a downward flow just outside of the core of the plume and below the intrusion level. This downward flow is accentuated in the immediate vicinity of the central core. The signature of this effect is the bending of the temperature isolines (figure 4 b). This structure causes the radial distribution of velocity and buoyancy to be drastically altered from the non-rotating case for which the time-averaged radial profiles are approximately Gaussian (hypothesis (3)). Since the radial profiles are not self-similar in the rotating case, one can never expect to find a power-law scaling solution from a similarity solution. The horizontal spread at the intrusion level is confined laterally as a consequence of the establishment of a geostrophic balance. That is, the radial pressure gradient, which in the non-rotating case was inducing an outward radial flow at the intrusion level, is now balanced by the Coriolis force, with the effect of shutting off the radial flow. Such balance is responsible for confining the fluid laterally. The relevant quantity to characterize the horizontal spread of the plume is the first Rossby radius of deformation, which for a linear stratification is simply

(5.1) $$\begin{eqnarray}R_{d}=\frac{NH}{{\rm\pi}f}\simeq \frac{(NF_{0})^{1/4}}{f},\end{eqnarray}$$

see for example Gill (1982). In (5.1), we used (2.13) as a scaling for the characteristic height $H$ and $c\simeq {\rm\pi}$ . For the plume plotted in figure 4, $R_{d}=13.7$ , such that the aspect ratio $H_{w}/R_{d}\ll 1$ , i.e. the vertical extent of the plume is far less than the Rossby deformation radius. Since this aspect ratio is small, we can still evaluate $H_{i}$ (the height of the intrusion level). We will not be able to do so when $H_{w}/R_{d}\lesssim 1$ (as is the case in figure 5), because there is no clear radial distance at which $\partial {\it\psi}/\partial r=0$ for all $z$ in the computational domain (see figure 5).

The lateral extension decreases when the rotation rate increases. This is shown in figure 5, where a similar steady-state calculation was carried with $Ro=0.4$ . In this configuration, $R_{d}=0.27$ , and hence the time evolution of the passive tracer obeys different dynamics in this case. At $t=55$ ( $\simeq 137f^{-1}$ ), the passive tracer uniformly fills the first envelope (red solid line in figure 5 b). The envelope then expands in the radial and vertical directions (above the top of the plume, at the intrusion level and outside of the core of the plume) via diffusive processes. This situation differs from the previous non-rotating cases where the passive tracer was confined at the intrusion level by advective processes.

It should be noted that even if the plume is laterally confined, there is still a mass flux at the far edge boundary, balancing the mass flux at the inlet to ensure mass conservation (not visible in figure 4 because this flow is weak compared with the circulation in the plume). As in the non-rotating case, we verified that these steady states are stable through time integrations of the model for several values of the buoyancy flux.

Figure 5. The same as figure 4(a,b) but with $Ro=0.4$ . For this configuration, $R_{d}=0.27$ and is marked with a dotted green line. The envelope of the passive tracer is shown at two different time intervals, $t=55$ ( $\simeq \,137f^{-1}$ ) (solid red line) and $t=1450$ ( $\simeq \,3625f^{-1}$ ) (dashed red line).

5.2 Multiple steady states

The model was initially calibrated using $H_{w}$ and $H_{b}$ of MTT’s model in the non-rotating case. We now explore the robustness of the steady states when the system is slightly pushed away from this state of reference. This step is in fact a crucial aspect of dynamical system theory for two reasons. (i) The initial calibration of $l_{s}$ does not guarantee that the axisymmetric model correctly captures all the details of the plume. In fact, in Fabregat et al. (2016a ), we highlighted some differences between this axisymmetric model and a three-dimensional (3D) model, such as the mixing rate at the base of the plume. (ii) A real-world turbulent system will often make excursions in parameter space. The reasons for these excursions might be non-uniformity of the turbulence level, non-stationarity of the buoyancy source or inhomogeneity in the medium. Therefore, it is important to understand how the steady states are perturbed when some parameters of the model are changed. In this section, we first investigate the effect of a variation of $l_{s}$ , the Smagorinsky length scale, on the structure of the steady states, and then introduce another bifurcation parameter $C_{b}$ .

Figure 6. Bifurcation diagram ( $\overline{w}$ as a function of $l_{s}$ ) in the non-rotating (black curve with stars) and rotating cases (red curve with crosses). Here, $S_{1}$ and $S_{2}$ mark the positions of the two saddle-node bifurcations. All of these steady states are stable except between $S_{1}$ and $S_{2}$ .

Figure 6 shows the bifurcation diagram of the height-averaged vertical velocity $\overline{w}$ at the centre of the plume when $l_{s}$ varies. This figure reads from right ( $l_{s}=0.01$ : reference case) to left. The black curve is the bifurcation diagram in the absence of rotation. In this case, $\overline{w}$ increases when $l_{s}$ decreases. A similar curve is obtained on decreasing ${\it\alpha}$ in MTT’s model (not shown), suggesting that increasing $l_{s}$ is, in some sense, equivalent to reducing the entrainment and mixing of the plume. The global shape of the plume is the same for all the points of this curve, and both $H_{w}$ and $H_{b}$ increase when $l_{s}$ decreases.

The red curve is a bifurcation diagram in the presence of rotation ( $Ro=20$ ). For $0.004<l_{s}<0.01$ , this curve follows the reference case without rotation. In this regime, the core of the plume is identical in both rotating and non-rotating cases, but as we saw in § 5.1 the far field is different. For $l_{s}<0.004$ , $\overline{w}$ decreases when $l_{s}$ decreases. Two saddle-node bifurcations occur: $S_{1}$ near $l_{s}=0.002$ and $S_{2}$ near $l_{s}=0.004$ . The branch past $S_{2}$ has the unconventional result $\overline{w}<0$ , implying that the average vertical velocity above the buoyant source is negative. We call this branch the centrifugal branch for reasons that will become obvious at the end of this section. In the range $0.002<l_{s}<0.004$ , we have three coexisting steady states and we can expect a hysteresis-type behaviour. The positions of the two saddle-node bifurcations are functions of the other parameters of the problem, such as the buoyancy flux and the Coriolis parameter. The exact positions of the two saddle-node bifurcations depend on the choice of the continuation parameter. However, what really matters in this analysis is the presence of two radically different solutions.

To prove the robustness of these two states, we repeat our bifurcation analysis with another bifurcation parameter, $C_{b}$ . Now, we keep $l_{s}$ constant when it appears in the stream function and temperature equations and we use

(5.2) $$\begin{eqnarray}l_{s}^{{\it\lambda}}=\left(1-\frac{C_{b}}{2}(1-\tanh (40(z-10{\rm\Delta}z)))\right)l_{s}\end{eqnarray}$$

for the angular momentum equation. The effect of this profile is to reduce $l_{s}$ near the injection site for the angular momentum only. We arbitrarily chose the cutoff to be located at 10 grid points above the source; other choices give similar results. Variation of this parameter has a weak impact on the shape of the stream function and temperature in the non-rotating case (with changes of $H_{w}$ of less than 1 % when $C_{b}$ varies between 0 and 1). The bifurcation diagram obtained when $C_{b}$ varies from 0 ( $l_{s}^{{\it\lambda}}=l_{s}$ ) to 1 is shown in figure 7. We recover the same structure for the steady states: two radically different solutions ( $\overline{w}>0$ and $\overline{w}<0$ ) separated by two saddle-node bifurcations.

Figure 7. Bifurcation diagram ( $\overline{w}$ as a function of $C_{b}$ ) for the rotating case only. In contrast to figure 6, this figure reads from left to right.

We analyse the structure of this new solution in figure 8, where we plot the steady state obtained at $l_{s}=0.003$ . In this figure, the upward motion is no longer centred on the axis of symmetry, but is deflected as a consequence of an outward radial velocity. Near the centre, we now have downward vertical velocities, so that the dynamics in the core of the plume is drastically different from in the non-rotating case. This observation demonstrates the limit of the validity of the entrainment hypothesis for plumes in rotating systems. This downward flow is particularly visible in figure 8(b), where the maximum temperature at a given height is shifted from the centreline. The height of the intrusion level remains comparable to the non-rotating case (not shown), but the radial velocity at the intrusion level is now as intense as the vertical velocity near the source. The quantities $H_{w}$ and $H_{b}$ are no longer adequate to describe the plume since they reflect the dynamics at $r=0$ . In this case, we define $H_{w}^{\prime }$ as the maximum height reached by the stream line ${\it\psi}=0$ . This solution no longer exhibits an overshoot as there is no pronounced downward vertical velocity outside the core of the plume near $z=H_{w}^{\prime }$ . The time integration (red curve in figure 8) confirms that the plume does not overshoot the intrusion level.

The presence of multiple steady states can also be highlighted with time integrations of the model. We start from an initial state of rest and let the system reach an equilibrium. After the initial transient phase, the model quickly settles into a state that matches the non-rotating limit (black curve in figure 6). In fact, the angular momentum is close to its initial condition and the effect of rotation is unimportant. As time integration progresses, the system starts to feel the effect of rotation, and the plume and its surroundings start to swirl. For time integration performed with $C_{b}=1$ (or $l_{s}<0.002$ ), we observe that, after a sudden jump, the system reaches the centrifugal branch with negative mean vertical velocity at the centre.

Figure 8. Typical steady state in the centrifugal branch ( $l_{s}=0.003$ ). (a) Stream function, ${\it\psi}$ , (b) temperature (black), $T$ , and envelope of passive tracer obtained with a time integration (red). The passive tracer reaches $r=1$ at $t=78$ ( $\simeq \,4f^{-1}$ ). This plot is very similar to the solution found at $C_{b}=1$ .

From the previous analysis, it is obvious that this new branch of solution is a consequence of rotation, and we provide a physical explanation for it. As we observed, a convection pattern typically gathers angular momentum at the centre and at the bottom, which translates to a non-zero swirl speed. For high values of this swirl speed, the core of the plume is in cyclostrophic balance: the radial pressure gradient balances the centrifugal force, (Smith 1980)

(5.3) $$\begin{eqnarray}\frac{\partial p}{\partial r}=\frac{v^{2}}{r}.\end{eqnarray}$$

We integrate this equation between $r=0$ and the edge of the core of the plume $r^{\prime }$ , where the dynamic pressure field is negligible. Differentiating the result with respect to $z$ gives

(5.4) $$\begin{eqnarray}\left.\frac{\partial p}{\partial z}\right|_{r=0}=-\frac{\partial }{\partial z}\int _{r=0}^{r^{\prime }}\frac{v^{2}}{r}\,\text{d}r.\end{eqnarray}$$

According to this equation, an adverse pressure gradient occurs when the swirl speed decays or spreads with height. In our model, the parameter $C_{b}$ (or $l_{s}$ ) controls the sharpness of the gradient of the angular momentum ${\it\lambda}$ and thus directly affects the intensity of the adverse pressure gradient. We showed that, past a given threshold (the saddle-node bifurcation $S_{1}$ ), the adverse pressure gradient along the jet axis is strong enough to create a stagnation point, also known as vortex breakdown (see, e.g., Billant, Chomaz & Huerre 1998), near the base of the plume. In this axisymmetric model, the steady states are stable when the stagnation point is either near the intrusion level or near the base of the plume. The intermediate cases where the stagnation point is at midheight (steady states between $S_{1}$ and $S_{2}$ ) are unstable according to figure 6.

It should be noted that the flow reversal (as shown in figure 8) is constrained by the axisymmetric nature of the system. Such an axisymmetric pattern would naturally be unstable to perturbations in the azimuthal direction. As a result, a real plume would generally precess around the central axis; this centrifugal mode could then be seen as a time mean of such a precessing plume.

5.3 Co-dimension-two bifurcations

According to the previous bifurcation diagram (figure 6), in the absence of rotation, there is only one steady state, but when rotation is present, multiple steady states may coexist. Our goal is now to clarify which type of plume (lazy or forced) can reach the centrifugal mode, depending on the Rossby number. This will result in a ‘phase diagram’ in the $({\it\Gamma}_{0},Ro)$ space. To do so, we now use the classical results of bifurcation theory (Kuznetsov 2004) to identify the origin of the two saddle-node bifurcations.

There are two co-dimension-two bifurcations that could involve two saddle nodes in a 2D-parameter space: the cusp and the Bogdanov–Takens (BT) bifurcations. The cusp bifurcation corresponds to a merging of two saddle nodes, in which case, no other branch exists. On the contrary, the BT bifurcation involves the collision of two branches. The normal form of the BT bifurcation has been extensively studied, and we know that it involves two Hopf bifurcations on either side of the saddle nodes (see Kuznetsov 2004). In our case, let us add the dimension ${\it\Gamma}_{0}$ to the bifurcation diagram, so that figure 6 is just a slice in this new 3D space $(l_{s},\overline{w},{\it\Gamma}_{0})$ . We tracked the two saddle nodes ( $S_{1}$ and $S_{2}$ ) when decreasing ${\it\Gamma}_{0}$ (from lazy plume to forced plume), and noticed that $S_{1}$ and $S_{2}$ merged in a cusp bifurcation at ${\it\Gamma}_{0}={\it\Gamma}_{0}^{c}$ . We see that, for this value of the Rossby number, there is only one steady state for ${\it\Gamma}_{0}<{\it\Gamma}_{0}^{c}$ . We repeat this operation for different values of the Rossby number in order to obtain a curve of ${\it\Gamma}_{0}^{c}$ as a function of $Ro$ (plotted in figure 9). This figure indicates the region in parameter space ( $Ro,{\it\Gamma}_{0}$ ) where multiple steady states may be found. This plot confirms that only lazy plumes may evolve in the centrifugal mode, and provides the order of magnitude of the threshold ${\it\Gamma}_{0}$ . It should be noted, however, that the quantitative results obtained in this section might be biased due to the imperfection of the turbulence closure scheme and would require more investigation (see the suggestions in the conclusion).

Figure 9. The position of the cusp bifurcation, ${\it\Gamma}_{0}^{c}$ , as a function of $Ro$ . The grey area is the region where multiple steady states coexist.

6 Illustration in 3D dynamics

We now investigate whether these solutions are artefacts of the axisymmetric model or whether they are also found in realistic configurations. The point is not to attempt to find precisely the same patterns in the axisymmetric and three-dimensional models since we know to expect differences (see § 5.2). Rather, it is to see whether we find elements in the 3D configuration that are reminiscent of the axisymmetric model, such as the vertical velocity profiles in the core of the plume. To address this question, we use a non-hydrostatic ocean model: the MITgcm (Marshall et al. 1997) in a 3D configuration. Since this study was initially motivated by a desire to understand the fate of oil after a deep-water blowout, we choose a configuration that mimics just such a situation. The lateral extension of the domain is 7 km, with periodic boundary conditions and a horizontal resolution of 20 m. The vertical extension is 1000 m, with a resolution of 5 m. At this resolution, we do not expect to fully resolve turbulence, but our goal is merely to have an overview of a 3D plume in a rotating environment. A buoyant tracer that represents the gas fraction is injected at the bottom of the model by means of simple relaxation on four grid cells in the middle of the domain. This method ensures that we are creating a lazy plume ( ${\it\Gamma}_{0}$ is infinite since the injection velocity is zero). The density of this tracer is 1 kg m $^{-3}$ , and it is relaxed to a target concentration of 0.01 with a time scale of 5 s such that the buoyancy flux $F_{0}=50$  m $^{4}$  s $^{-1}$ . This value is computed a posteriori by measuring the flux at the top of the grid cell where a buoyancy flux is applied. The initial vertical stratification is $10$  K km $^{-1}$ , which corresponds to $N^{2}=2\times 10^{-5}$  s $^{-2}$ . The Coriolis parameter is $f=10^{-4}$  s $^{-1}$ . The kinematic viscosity and thermal diffusivity coefficients are set to 0.05 m $^{2}$  s $^{-1}$ , which corresponds to $Re\simeq 10^{6}$ . We do not use any other turbulence closure since we expect to partially resolve the turbulent processes in this 3D configuration. Using these parameters, (2.13) predicts an intrusion level of 584 m.

Figure 10. Vertical cross-section in the middle of the domain showing two time averages of the vertical velocity for two distinct periods: (a) early stage ( $0<t<2f^{-1}$ ) and (b) mature stage ( $4f^{-1}<t<8f^{-1}$ ), shown in dimensional units (m for the axis and m s $^{-1}$ for the colour bar). To highlight negative velocities, the colour scales do not extend over all positive values.

In figure 10, we plot the mean vertical velocity in the early stages of the development of the plume (a) and once the plume is mature (b). We interrupt the numerical simulation after $t\sim 10f^{-1}$ because after that time baroclinic instability will affect the dynamics of the plume, which is beyond the scope of this study (see Helfrich & Battisti 1991). In figure 10(a), the plume penetrates almost to the surface and corresponds to an overshoot profile. Indeed, the intrusion level (identified as the height at which the tracer extends laterally) is approximately 400 m (36 % lower than the prediction). Weak negative vertical velocities are also visible on the edge of the core of the plume. Once the plume is in a mature stage (figure 10 b), the mean vertical velocity evolves in a ‘V-shape’ profile. The vertical penetration is deflected away from the central axis. On this axis, we observe weak downward velocities. This situation is clearly reminiscent of the steady states on the centrifugal branch of the axisymmetric model. The shapes of the mean swirl speed and angular momentum patterns are also similar to the axisymmetric case (not shown). We state again that this 3D numerical simulation should not be regarded as a validation of the 2D axisymmetric model but rather as a preliminary experiment suggesting the need for a systematic 2D/3D comparison.

7 Conclusions

While convective plumes are often described either by scaling laws or by large eddy simulations (LES), it was our intent to bridge the gap between these two extremes by proposing a model of intermediate complexity to see what new information we gain from this midway step. To describe a convective plume, we used an axisymmetric model where turbulent processes are parameterized by the Smagorinsky (1963) closure scheme. We also take advantage of the continuation method to compute the steady states of the model as a function of the parameters (source forcing ${\it\Gamma}_{0}$ , Smagorinsky length scale $l_{s}$ and Rossby number $Ro$ ). When the system is not rotating ( $Ro\rightarrow \infty$ ), in the absence of sufficient observational or experimental data, we can calibrate our axisymmetric model using the well-known MTT model. We separated the forced plume regime, for which the height of the intrusion level is independent of the buoyancy flux, and the lazy plume regime, for which the height of the intrusion layer varies with $F_{0}^{1/4}$ .

We then introduced the effect of rotation in the system. The intrusion layer is thinner in the rotating case while the core of the plume remains largely unchanged. The swirl speed resulting from the displacement of angular momentum lines is maximum near the base of the plume. As expected, we observe a cyclonic swirl speed near the base of the plume and an anticyclonic swirl speed at the intrusion level. The plume is laterally confined at the intrusion level due to the establishment of a geostrophic balance.

We then explored the robustness of these patterns and investigated the effect of a change of the Smagorinsky length scale $l_{s}$ . For lower values of $l_{s}$ (less turbulent plumes), a new branch of solution emerges. On this branch, the convection occurs on the outer part of the plume as a result of the centrifugal force. This convection pattern is characterized by a downward velocity above the buoyancy source and the absence of overshoot. At the intrusion level, the outward radial velocity is more intense than in the non-rotating case. We also found a region in the parameter space where the classical and centrifugal branches coexist. The coexistence of these two states is important because it has an impact on the shape of the plume and can potentially affect the subsequent development of a baroclinic instability. Last, we searched for a co-dimension-two bifurcation which could indicate whether another branch could be present. We found that the two saddle-node bifurcations merge in a cusp bifurcation at a critical value of the buoyancy flux parameter ${\it\Gamma}_{0}^{c}$ . We repeated this operation for several values of the Rossby number and plotted a phase diagram in the $(Ro,{\it\Gamma}_{0})$ space, indicating where multiple steady states were to be expected.

Of course, further investigation with 3D models and laboratory experiments is needed to get the correct quantitative details (e.g. exact positions of the bifurcation points). These future quantitative laboratory of numerical experiments will have to overcome the imperfections of the turbulence closure scheme that are inevitable in a 2D axisymmetric model.

Finally, we investigated the robustness of these solutions and showed that some elements survive in a preliminary 3D model (where turbulence was partly resolved and no longer parameterized). This was illustrated with an example of oceanic convection. At the initial stage, when the system does not feel the effect of rotation, we obtained the ‘classic’ solution. As time goes on, and angular momentum is advected by the plume, a swirl speed develops and the convection pattern evolves into a centrifugal solution. We observe downward vertical velocities above the source and a tilted convection.

Although scaling laws and MTT’s model provide a very attractive description of the plume in the non-rotating case, they have not been designed to describe the convective pattern in the rotating case. Using a model of intermediate complexity, we were able to highlight how a plume deforms in the presence of rotation. The next step is to systematically compare this axisymmetric model with 3D LES or direct numerical simulations (see, for example, the study of Craske & van Reeuwijk 2015, and Fabregat et al. 2016b ) in order to refine the description of the turbulence in the plume. We also wish to make an extensive comparison of these convective patterns with laboratory experiments. A precise understanding of the plume dynamics could help in adjusting the parameterization of convection in ocean models. This could, in turn, improve our knowledge of the abyssal circulation in the ocean.