Transitions of Spherical Thermohaline Circulation to Multiple Equilibria

TL;DR

This paper analyzes how thermohaline circulation in a spherical shell transitions to multiple equilibrium states, identifying conditions for different transition types and deriving exact formulas for transition numbers.

Contribution

It provides the first exact formulas for transition numbers in spherical thermohaline circulation and characterizes the nature of bifurcations depending on system parameters.

Findings

Transition type depends on the sign of the transition number.

Exact formulas for transition numbers are derived for specific spherical harmonic degrees.

Numerical results show different transition behaviors based on the Lewis number.

Abstract

The main aim of the paper is to investigate the transitions of the thermohaline circulation in a spherical shell in a parameter regime which only allows transitions to multiple equilibria. We find that the first transition is either continuous (Type-I) or drastic (Type-II) depending on the sign of the transition number. The transition number depends on the system parameters and $l_c$, which is the common degree of spherical harmonics of the first critical eigenmodes, and it can be written as a sum of terms describing the nonlinear interactions of various modes with the critical modes. We obtain the exact formulas of this transition number for $l_c=1$ and $l_c=2$ cases. Numerically, we find that the main contribution to the transition number is due to nonlinear interactions with modes having zero wave number and the contribution from the nonlinear interactions with higher frequency modes is negligible. In our numerical experiments we encountered both types of transition for $Le<1$ but only continuous transition for $Le>1$. In the continuous transition scenario, we rigorously prove that an attractor in the phase space bifurcates which is homeomorphic to the 2$l_c$ dimensional sphere and consists entirely of degenerate steady state solutions.

Tables

Table 1. Table 1. The numerical values of the nonlinear interaction terms for l c = 1 subscript 𝑙 𝑐 1 l_{c}=1 , r = 2 / π 𝑟 2 𝜋 r=2/\pi , Pr = 7.5 Pr 7.5 \text{Pr}=7.5 for different Le and R .

	$\text{R}=620$		$\text{R}=660$
Le	$d_{(1,1),(0,2)}$	$d_{(1,1),(2,2)}$	$d_{(1,1),(0,2)}$	$d_{(1,1),(2,2)}$
${10}^{-2}$	40.825	1.493	-23.53	-0.8537
${10}^{-1}$	1.955	0.074	0.275	0.0122
0.5	0.477	0.0196	0.424	0.0176
5	0.421	0.175	0.427	0.1177

Table 2. Table 2. The numerical values of the nonlinear interaction terms for l c = 2 subscript 𝑙 𝑐 2 l_{c}=2 , r = 2 ​ 3 / π 𝑟 2 3 𝜋 r=2\sqrt{3}/\pi , Pr = 7.5 Pr 7.5 \text{Pr}=7.5 for different Le and R .

	$\text{R}=620$			$\text{R}=660$
Le	$d_{(2,1),(0,2)}$	$d_{(2,1),(2,2)}$	$d_{(2,1),(4,2)}$	$d_{(2,1),(0,2)}$	$d_{(2,1),(2,2)}$	$d_{(2,1),(4,2)}$
${10}^{-2}$	40.825	8.359	0.592	-23.53	-4.797	-0.338
${10}^{-1}$	1.956	0.407	0.029	0.275	0.063	0.005
0.5	0.477	0.104	0.008	0.424	0.093	0.007
5	0.421	0.092	0.0069	0.427	0.093	0.007

Table 3. Table 3. Types of transitions with respect to R for l c = 1 subscript 𝑙 𝑐 1 l_{c}=1 , r = 2 / π 𝑟 2 𝜋 r=2/\pi , Pr = 7.5 Pr 7.5 \text{Pr}=7.5 for different Le values.

Le	Type-I transition	Type-II transition
${10}^{-2}$	if	if
${10}^{-1}$	if	if
$0.5$	if	if

Table 4. Table 4. Types of transitions with respect to R for l c = 2 subscript 𝑙 𝑐 2 l_{c}=2 , r = 2 ​ 3 / π 𝑟 2 3 𝜋 r=2\sqrt{3}/\pi , Pr = 7.5 Pr 7.5 \text{Pr}=7.5 for different Le values.

Le	Type-I transition	Type-II transition
${10}^{-2}$	if	if
${10}^{-1}$	if	if
$0.5$	if	if

Equations

$$K=\mathrm{sign}(1-\text{Le})\left[\frac{{\text{Le}}^2}{1-\text{Le}}(1+\frac{1}{\text{Pr}}){\sigma }_c-\tilde{\text{R}}\right],$$

$${\sigma }_c=\underset{l\in {\mathbb{Z}}^{+}}{\min }\frac{({\pi }^2+{\alpha }_l^2{)}^3}{{\alpha }_l^2}=\frac{({\pi }^2+{\alpha }_{l_c}^2{)}^3}{{\alpha }_{l_c}^2},$$

$${\alpha }_l^2=\frac{l(l+1)}{r^2},l\in {\mathbb{Z}}^{+},$$

$$\sigma =\text{R}-{\text{Le}}^{-1}\tilde{\text{R}},$$

$$q_1=\sum _{l=0}^{l_c}{d}_{(l_c,1),(2l,2)},$$

$$\mathrm{sign}(q_{l_c})=\mathrm{sign}(d_{(l_c,1),(0,2)})=\frac{(1+\text{Le}){\alpha }_{l_c}^4}{16\pi \text{Le}({\text{R}}_0-\text{R})}({\text{R}}_1-\text{R}),$$

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=\nu \Delta \mathbf{u}-\frac{1}{{\rho }_0}(\nabla p+\rho g{\hat{e}}_z),$$

$$\frac{\partial T}{\partial t}+(\mathbf{u}\cdot \nabla )T={\kappa }_T\Delta T,$$

$$\frac{\partial S}{\partial t}+(\mathbf{u}\cdot \nabla )S={\kappa }_S\Delta S,$$

$$\text{div}\mathbf{u}=0,$$

$$\rho ={\rho }_0\left(1-a_T(T-T_0)+a_S(S-S_0)\right),$$

$${\mathbf{u}}_{\text{ss}}=0,$$

$$T_{\text{ss}}=T_0-(T_0-T_1)\frac{z}{h},$$

$$S_{\text{ss}}=S_0-(S_0-S_1)\frac{z}{h},$$

$$p_{\text{ss}}=p_0-g{\rho }_0\left(z+\frac{a_T}{2}(T_0-T_1)\frac{z^2}{h}-\frac{a_S}{2}(S_0-S_1)\frac{z^2}{h}\right).$$

$$\frac{\partial \mathbf{u}}{\partial t}=\text{Pr}(\Delta \mathbf{u}-\nabla p)+\text{Pr}(\text{R}T-\mathrm{sign}(S_0-S_1)\tilde{\text{R}}S){\hat{e}}_z-(\mathbf{u}\cdot \nabla )\mathbf{u},$$

$$\frac{\partial T}{\partial t}=\Delta T+w-(\mathbf{u}\cdot \nabla )T,$$

$$\frac{\partial S}{\partial t}=\text{Le}\Delta S+\mathrm{sign}(S_0-S_1)w-(\mathbf{u}\cdot \nabla )S,$$

$$\text{div}\mathbf{u}=0,$$

$$\text{R}=\frac{a_{T}g(T_0-T_1)h^3}{{\kappa }_T\nu },$$

$$\text{Pr}=\frac{\nu }{{\kappa }_T},$$

$$r=\frac{a}{h}.$$

$${\nabla }_{u}f=u\cdot \nabla f=\frac{1}{r}\left(u_{\theta }\frac{\partial f}{\partial \theta }+\frac{u_{\phi }}{\sin \theta }\frac{\partial f}{\partial \phi }\right),$$

$$\Delta f=\frac{1}{r^2\sin \theta }\left[\frac{\partial }{\partial \theta }\left(\sin \theta \frac{\partial f}{\partial \theta }\right)+\frac{1}{\sin \theta }\frac{{\partial }^{2}f}{\partial {\phi }^2}\right],$$

$${\nabla }_{u}v=\frac{1}{r}\left(u_{\theta }\frac{\partial v_{\theta }}{\partial \theta }+\frac{u_{\phi }}{\sin \theta }\frac{\partial v_{\theta }}{\partial \phi }-u_{\phi }{v}_{\phi }\cot \theta \right){\hat{e}}_{\theta }+\frac{1}{r}\left(u_{\theta }\frac{\partial v_{\phi }}{\partial \theta }+\frac{u_{\phi }}{\sin \theta }\frac{\partial v_{\phi }}{\partial \phi }+u_{\phi }{v}_{\theta }\cot \theta \right){\hat{e}}_{\phi },$$

$$\Delta u=\left(\Delta u_{\theta }-\frac{2\cos \theta }{r^2{\sin }^2\theta }\frac{\partial u_{\phi }}{\partial \phi }-\frac{u_{\theta }}{r^2{\sin }^2\theta }\right){\hat{e}}_{\theta }+\left(\Delta u_{\phi }+\frac{2\cos \theta }{r^2{\sin }^2\theta }\frac{\partial u_{\theta }}{\partial \phi }-\frac{u_{\phi }}{r^2{\sin }^2\theta }\right){\hat{e}}_{\phi },$$

$$u_t+{\nabla }_{u}u+w\frac{\partial u}{\partial z}+\text{Pr}\nabla p-\text{Pr}(\Delta +{\partial }_{zz})u=0,$$

$$w_t+{\nabla }_{u}w+w\frac{\partial w}{\partial z}+\text{Pr}\frac{\partial p}{\partial z}-\text{Pr}(\text{R}T-\mathrm{sign}(S_0-S_1)\tilde{\text{R}}S)-\text{Pr}(\Delta +{\partial }_{zz})w=0,$$

$$T_t+{\nabla }_{u}T+w\frac{\partial T}{\partial z}-w-(\Delta +{\partial }_{zz})T=0,$$

$$S_t+{\nabla }_{u}S+w\frac{\partial S}{\partial z}-\mathrm{sign}(S_0-S_1)w-\text{Le}(\Delta +{\partial }_{zz})S=0,$$

$$\text{div}u+\frac{\partial w}{\partial z}=0.$$

$$w=T=S=0,\frac{\partial u}{\partial z}=0,\text{at }z=0,1,$$

$$H=\{(u,w,T,S)\in L^2(\Omega {)}^5\mid \text{div}u=w=T=S=0\text{ at }z=0,1\},$$

$$H_1=\{(u,w,T,S)\in H^2(\Omega {)}^5\mid \text{div}u=w=T=S=0,\frac{\partial u}{\partial z}=0,\text{ at }z=0,1\}.$$

$$⟨{\Psi }_1,{\Psi }_2⟩={\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^1(u_1\cdot \overline{u_2}+w_1\overline{w_2}+T_1\overline{T_2}+S_1\overline{S_2})r^2\sin \theta dzd\phi d\theta .$$

$$L(\Psi )=\mathcal{P}\left[\begin{array}{c}\text{Pr}(\Delta +\frac{{\partial }^2}{\partial z^2})u\\ \text{Pr}(\Delta +\frac{{\partial }^2}{\partial z^2})w+\text{Pr}(\text{R}T-\mathrm{sign}(S_0-S_1)\tilde{\text{R}}S)\\ (\Delta +\frac{{\partial }^2}{\partial z^2})T+w\\ \text{Le}(\Delta +\frac{{\partial }^2}{\partial z^2})S+\mathrm{sign}(S_0-S_1)w\end{array}\right],\Psi =(u,w,T,S)\in H_1,$$

$$G({\Psi }_1,{\Psi }_2)=-\mathcal{P}\left[\begin{array}{c}{\nabla }_{u_1}{u}_2+w_1\frac{\partial u_2}{\partial z}\\ {\nabla }_{u_1}{w}_2+w_1\frac{\partial w_2}{\partial z}\\ {\nabla }_{u_1}{T}_2+w_1\frac{\partial T_2}{\partial z}\\ {\nabla }_{u_1}{S}_2+w_1\frac{\partial S_2}{\partial z}\end{array}\right],{\Psi }_i=(u_i,w_i,T_i,S_i)\in H_1.$$

$$\frac{\partial \Psi }{\partial t}=L\Psi +G(\Psi ),\Psi (0)={\Psi }_0,$$

$$\text{Pr}\left((\Delta +{\partial }_{zz})u-\nabla p\right)=\beta u,$$

Full text

1. Introduction

An important source of climate low frequency variability is the so called thermohaline circulation (THC). The underlying mechanism of THC is well established. THC is essentially driven by the temperature and freshwater fluxes in the ocean-atmosphere interface which in turn produces density gradients. These gradients are much sharper in the vertical direction compared to the horizontal directions and therefore are associated with an overturning.

There are indications that the Atlantic circulation has varied in the past [2]. Taking into account the enormous effects the ocean circulation has on the climate, the sensitivity, stability and transitions of the large scale ocean circulation became an important issue in climate research [12]. The studies of the THC using a hierarchy of ocean models, starting with a very simple box model [16] shows that the presence of heat and salt, with their different influences on the density field, may lead to different stable steady flow patterns.

This paper arises out of a research program to generate rigorous mathematical results on climate variability developed from the viewpoint of dynamical transitions [10, 5, 6]. The basic philosophy of dynamic transition theory is to search for the full set of transition states, giving a complete characterization of stability and transition. The set of transition states which is often represented by a local attractor may lie near or away from the basic state. This theory has recently been successfully applied to many branches of nonlinear sciences, see [19, 13, 4, 18] among others.

One of the main focus of the dynamic transition theory is the identification of transition states and the classification of dissipative systems into three distinct transition types, namely continuous (Type-I), catastrophic (Type-II) and random (Type-III) transitions which describe the nature of the state of the system as the control parameter crosses a critical threshold. The transition states stay in a close neighborhood of the basic state in the case of continuous transition, and lie outside of a neighborhood of the basic state in the case of a catastrophic transition. For a random type transition, a neighborhood of the basic state consists of two disjoint open regions with a continuous transition taking place in one region while a catastrophic transition occurring in the other one.

We briefly recall several recent papers which investigate the problem in simpler settings than the one studied in this paper. The paper [1] considers the THC in a 2D rectangular box and considers the oscillatory transitions and finds that both continuous and drastic transitions are possible. In [9], the authors investigate the dynamic transitions of THC in a 3D rectangular enclosure where they show that THC exhibits transitions to either multiple equilibria or to time periodic solutions and that there are parameter regimes leading to continuous or catastrophic transitions. Finally, in [17], the authors consider the dynamic transitions in a spherical shell but only for the pure thermal convection case without salinity gradients. The transition type in that case is well-known to be continuous and the transition is described by an attractor bifurcation [11]. Hence in [17], the authors aim to describe the structure of the local bifurcated attractor.

The main objective of this paper is to carry out the dynamic transition analysis of the THC problem in the important case of a spherical shell domain. For simplicity, the spatial domain is considered as the product of a two dimensional sphere $S_{a}^{2}$ with radius $a$ and an interval $(0,h)$ where $h$ is is the height of the fluid layer. This approach is mainly motivated by the fact that the aspect ratio for the large scale atmosphere and ocean is small; see among others [7, 14]. The main challenge in carrying the analysis of [9] to a spherical shell domain is due to the detailed calculations of nonlinear interactions of spherical harmonics (both scalar and vectorial) that need to performed to compute the reduction of the infinite dimensional system to a finite one.

It is known [9] that the THC system exhibits first transitions to both multiple equilibria and to spatio-temporal oscillations depending on the parameter

[TABLE]

where Pr represents the Prandtl number, Le is the Lewis number, $\tilde{\text{R}}$ is the saline Rayleigh number. Here $\sigma_{c}$ and $l_{c}\in\mathbb{Z}^{+}$ are defined by

[TABLE]

which exist thanks to the convexity of the function $x\to\frac{(\pi^{2}+x)^{3}}{x},\,x>0$. In (2), the wave number $\alpha_{l}$ is defined as

[TABLE]

with $r=a/h$ denoting the aspect ratio of the sphere. The well known result of Rayleigh Bénard convection is that the minimum of $\sigma_{c}$ is achieved when $\alpha_{l_{c}}=\frac{\pi}{\sqrt{2}}$ in which case $\sigma_{c}=\frac{27}{4}\pi^{4}$, see [3].

In this paper, we restrict ourselves to the parameter regime $K>0$ which only allows transitions to multiple equilibria and will address the case of transitions to spatio-temporal oscillations ($K<0$ regime) in another paper. In the case $K>0$, the first dynamic transition of the system occurs as the control parameter $\sigma$ defined by

[TABLE]

crosses the critical threshold $\sigma_{c}$, leading to multiple equilibria. Here R is the thermal Rayleigh number.

Our main results are as follows. We find that the transition of the THC problem at $\sigma=\sigma_{c}$ for $l_{c}=1$ or $l_{c}=2$ is either Type-I or Type-II depending on the sign of the transition number $q_{l_{c}}$. We present the exact formulas for the transition number $q_{l_{c}}$ which depend on the system parameters for $l_{c}=1$ and $l_{c}=2$ cases in (30)–(31). In these cases, the transition number $q_{l_{c}}$ is given by the sum

[TABLE]

where $d_{(l_{c},1),(2l,2)}$ is a number determined by the nonlinear interactions of the critical modes with all modes corresponding to the spherical harmonics of degree $2l$, $0\leq l\leq l_{c}$ and vertical wave number $2$ (there are $2l+1$ of such modes). In the continuous transition scenario, an attractor in the phase space bifurcates on $\sigma>\sigma_{c}$ which is a $2l_{c}$-dimensional homological sphere. We show that this sphere is indeed homeomorphic to $S^{2l_{c}}$, the 2$l_{c}$ dimensional sphere and consists entirely of degenerate steady state solutions.

Our numerical explorations in the parameter space suggest that

[TABLE]

when R is away from $\text{R}_{1}=\frac{\sigma_{c}}{1-\text{Le}^{2}}$, indicating that the nonlinear interactions of the critical modes with the zero wave number modes determine the type of transition and that the contribution from higher frequency modes to the transition number is negligible. Also our numerical simulations suggest the following. For $\text{Le}<1$, both continuous and drastic transitions are possible and a continuous transition is preferred for $\text{R}<\text{R}_{\ast}$ while a drastic transition is preferred for $\text{R}_{\ast}<\text{R}<\text{R}_{0}=\frac{\text{Le}+\text{Pr}}{(1-\text{Le})\text{Pr}}\sigma_{c}$ with $R_{\ast}\approx R_{1}$. For $\text{Le}>1$ flows, the transition type does not change and is always Type-I for all R.

The paper is organized as follows. In Section 2, the mathematical setting of the problem is introduced. In Section 3, the linear stability analysis of the main equations is summarized. Section 4 details our main theorem and its proof. Finally in Section 5, we explore the transition number numerically.

2. Governing Equations and the Functional Setting

As mentioned in the Introduction, we consider a spherical shell $S_{a}^{2}\times(0,h)$ as the spatial domain for the motion of the large scale ocean. Here $S_{a}^{2}$ represents the 2D sphere with radius $a$ and $h$ denotes the height of the fluid layer. The governing equations are the familiar Boussinesq equations (see [14, 9] among others):

[TABLE]

where ${\bf u}$ is the velocity, $T$ is the temperature, $S$ is the salinity, $\hat{e}_{z}$ is the unit vector in the z-direction, $\nu$, $\kappa_{T}$, $\kappa_{S}$, $g$ are all positive constants denoting the kinematic diffusivity, the thermal diffusivity, the saline diffusivity and the gravitational constant, respectively. The fluid density $\rho$ is given by the linear equation of state

[TABLE]

where $a_{T}$ and $a_{S}$ are assumed to be positive constants and $\rho_{0}$ is the density at the lower surface. The constants $T_{0}$, $S_{0}$ represent the fixed temperature and salinity at the lower boundary $z=0$, whereas we denote the fixed temperature and salinity at the upper boundary $z=h$ by $T_{1}$ and $S_{1}$.

In (5), the opposing effects of temperature and salinity are evident. Both $T_{0}>T_{1}$ (heated from below) and $S_{0}<S_{1}$ (salted from above) are destabilizing mechanisms. We will treat the general case which allows one or both of these conditions to be satisfied. In other words, a competition between a stabilizing and a destabilizing mechanism is allowed.

The trivial steady state solution to the problem (4)–(5) is given by

[TABLE]

We nondimensionalize the equations (4) exactly in the same way as in [9], and obtain

[TABLE]

where $w$ is the vertical velocity in the direction of $\hat{e}_{z}$, R is the thermal Rayleigh number, $\tilde{\text{R}}$ is the saline Rayleigh number, Pr is the Prandtl number and Le is the Lewis number defined as

[TABLE]

In (7), the unknowns now represent deviations from the steady state solutions given by (6). Also the nondimensional spatial domain is $\Omega=S_{r}^{2}\times(0,1)$ where the aspect ratio is

[TABLE]

Let ${\bf u}=(u,w)$ be the 3D velocity vector where $u=u_{\theta}\hat{e}_{\theta}+u_{\varphi}\hat{e}_{\varphi}$ is the 2D horizontal velocity field. Hereafter, $\nabla$, $\nabla_{u}$, div and $\Delta$ will denote both the scalar and vectorial differential operators in the horizontal direction, that is on the sphere $S^{2}_{r}$, given by

[TABLE]

where $f$ is a scalar, $u=u_{\theta}\hat{e}_{\theta}+u_{\varphi}\hat{e}_{\varphi}$ and $v=v_{\theta}\hat{e}_{\theta}+v_{\varphi}\hat{e}_{\varphi}$ are 2D vectors.

With the above notations, (7) can be written as

[TABLE]

In this study, we consider the equations (9) supplemented with the free-slip boundary conditions

[TABLE]

and note that our analysis can be expanded to other types of boundary conditions as well.

The below functional spaces are needed to recast the equations (9) and (10) in an abstract form.

[TABLE]

We recall that the inner product in $H$ is defined for vectors $\Psi_{i}=(u_{i},w_{i},T_{i},S_{i})$ as

[TABLE]

Let the linear operator $L:H_{1}\rightarrow H$ be defined by

[TABLE]

and the bilinear operator $G:H_{1}\times H_{1}\rightarrow H$ be defined by

[TABLE]

Here $\mathcal{P}:L^{2}(\Omega)^{5}\rightarrow H$ is the Leray projection.

Now the problem can be cast as an abstract ODE as

[TABLE]

where $G(\Psi)=G(\Psi,\Psi).$

3. Linear Stability

We first consider the eigenvalue problem for the linearized equations of (9)

[TABLE]

supplemented with the boundary conditions with (10).

The z-independent solutions of (12)–(14) clearly satisfy $w=T=S=0$. In this case, taking the curl of (12) and using the incompressibility condition is equivalent to

[TABLE]

It is easy to obtain the solutions of this equation and they are given by (20).

For the general case, we use the separation of variables in the form

[TABLE]

With (15), the incompressibility equation (14) is equivalent to the Helmholtz equation on the sphere

[TABLE]

which has solutions $f=Y_{lm}(\theta,\varphi)$ only when $\alpha^{2}=\alpha_{l}^{2}=\frac{l(l+1)}{r^{2}}$, where $Y_{lm}$ are the spherical harmonics, $l,m\in\mathbb{Z}$, $l\geq 1$ and $-l\leq m\leq l$.

Eliminating the pressure term $p$ by taking $\nabla$(13)-$\frac{\partial}{\partial z}$(12), using (15)–(16) and denoting $D=\frac{d}{dz}$, the equations (12)–(14) become

[TABLE]

[TABLE]

[TABLE]

From (10), the boundary conditions of the above system of ODEs are

[TABLE]

If $l=0$, then $u=0$, $w=0$ and the equations (17)–(19) reduce to

[TABLE]

whose solutions are given in (21).

On the other hand, if $l\neq 0$, plugging

[TABLE]

into (17), (18) and (19), we find the coefficients $\Theta_{lmn}$ and $\Phi_{lmn}$ given by (23) and the compatibility condition (22) for the existence of eigenvalues.

3.1. Eigenpairs

We will denote the eigenvalues of the linearized operator by $\beta_{ln}$ and their corresponding eigenvectors by $\Psi_{lmn}=(u,w,T,S)$ where $l,m,n\in\mathbb{Z}$, $l\geq 0$, $-l\leq m\leq l$, $n\geq 0$. By the discussion in the previous section, we now summarize all the eigenpairs.

If $l\neq 0$ and $n=0$, then

[TABLE]

where

[TABLE]

If $l=0$ and $n\neq 0$, then

[TABLE]

If $l\neq 0$ and $n\neq 0$, there are three distinct eigenvalues which we order as $\operatorname{Re}(\beta^{1}_{ln})\geq\operatorname{Re}(\beta^{2}_{ln})\geq\operatorname{Re}(\beta^{3}_{ln})$ corresponding to the three distinct solutions of

[TABLE]

where

[TABLE]

As the relation (22) is independent of $m$, to each $\beta_{ln}^{k}$ ($l\neq 0$, $n\neq 0$, $k\in\{1,2,3\}$) corresponds $2l+1$ eigenfunctions $\Psi_{lmn}^{k}$, $m=-l,\dots,l$ given by

[TABLE]

By the symmetry $\overline{Y}_{lm}=(-1)^{m}Y_{l-m}$ of the spherical harmonics, it follows that if $\overline{\beta}_{ln}^{k_{1}}=\beta_{ln}^{k_{2}}$ then

[TABLE]

In particular when $\beta_{ln}^{k}\in\mathbb{R}$, we have $\overline{\Psi}_{lmn}^{k}=(-1)^{m}\Psi_{l-mn}^{k}$.

3.2. Adjoint Problem

For our analysis, we also require the explicit form of the adjoint eigenvectors. The adjoint problem of (12)-(14) can be written as

[TABLE]

In the case $n=0$, $l\neq 0$ and in the case $l=0$, $n\neq 0$, the adjoint eigenfunctions are as given in (20) and (21) respectively. When $l\neq 0$ and $n\neq 0$, the eigenfunctions are now given by

[TABLE]

3.3. Principle of Exchange of Stabilities

A crucial parameter in the first transition is the critical wave integer $l_{c}$ defined in (2) which depends only on the aspect ratio $r$. This dependence can be easily obtained explicitly from (2) by solving $\sigma_{c}(l)=\sigma_{c}(l+1)$ and the following lemma holds.

Lemma 1**.**

Define

[TABLE]

Then the critical integer $l_{c}$ defined implicitly by (2) is

[TABLE]

Numerical values of the first three $r_{l}$, computed from (25), are

[TABLE]

We remark here that when the aspect ratio is critical, i.e. $r=r_{l}$ for some $l\in\mathbb{Z}^{+}$, the eigenvalues $\beta_{l_{c}1}^{1}$ and $\beta_{l_{c}+1,1}^{1}$ with wave numbers $\alpha_{l_{c}}$ and $\alpha_{l_{c}+1}$ will become critical simultaneously. This type of transition was studied for the pure 2D Bénard convection in a rectangular container in [15] leading to competition among the pure $l_{c}$-modes, the pure $l_{c}+1$-modes and the mixed modes (superpositions of $l_{c}$ and $l_{c}+1$ modes).

The dispersion relation (22) is exactly the same as in [9] in which the thermohaline circulation problem is considered in a rectangular domain. Thus the following principle of exchange of stabilities follows directly from the corresponding one in [9].

Theorem 1**.**

Assume that $\text{Le}\neq 1$ and that the aspect ratio $r$ defined by (8) is not critical, i.e. $r\neq r_{l}$ for any $l\in\mathbb{Z}^{+}$, with $r_{l}$ as defined by (25). Consider $K$, $\sigma$, $\sigma_{c}$ and the integer $l_{c}$ defined by (1), (3), (2) and (26) respectively.

1. (1)

If $K>0$ then $\sigma_{c}$ is the first critical Rayleigh number and the first critical eigenvalue $\beta^{1}_{l_{c}1}$ has multiplicity $2l_{c}+1$. Moreover the condition

[TABLE]

is satisfied. 2. (2)

If $K<0$, let

[TABLE]

Then

[TABLE]

is the first critical Rayleigh number. Moreover the first two critical eigenvalues satisfy $\beta^{1}_{l_{c}1}=\overline{\beta^{2}_{l_{c}1}}$ and each has multiplicity $2l_{c}+1$ satisfying

[TABLE]

Since the sign of the parameter $K$ defined by (1) plays a crucial role in the transition, let us investigate it in some detail. For this let us define

[TABLE]

Then

[TABLE]

If $\text{Le}>1$, it is clear that $\text{R}_{0}<0$ and $K>0$. From these observations, we find that

[TABLE]

at the criticality $\sigma=\sigma_{c}$ or $\eta=\eta_{c}$, see also Figure 1.

4. Main Theorem and its Proof

In this section we describe the transitions of the system (11) to multiple equilibria. For simplicity, we will consider the cases for which the critical integer $l_{c}$ defined by (2) is $1$ or $2$. The reason for this assumption is that the reduction to the center manifold, although still computable, becomes increasingly difficult to compute as $l_{c}$ is increased.

We recall that according to Lemma 2, the aspect ratio ratio $r=a/h$ must be in the range $0<r<r_{1}$ for the $l_{c}=1$ case and $r_{1}<r<r_{2}$ for the $l_{c}=2$ case.

Now let us define the following coefficients

[TABLE]

[TABLE]

[TABLE]

[TABLE]

In the above notation, $d_{(l_{c},1),(l,n)}$ is the coefficient arising from the nonlinear interactions of the critical modes, i.e. eigenmodes corresponding to the $\beta^{1}_{l_{c}1}$ eigenvalue with all the modes corresponding to $\beta_{ln}$ eigenvalue.

Theorem 2**.**

Assume the conditions of Theorem 1 and that $l_{c}=1$ or $l_{c}=2$. Consider $q_{l_{c}}$ defined by (30) and (31). If $\sigma<\sigma_{c}$ then the basic solution $\psi=0\in H$ is locally asymptotically stable.

• •

If $q_{l_{c}}>0$ then there is a continuous (Type-I) transition at $\sigma=\sigma_{c}$ and an attractor $\Sigma_{\sigma}$ bifurcates on $\sigma>\sigma_{c}$ which is homeomorphic to the 2$l_{c}$* dimensional sphere. $\Sigma_{\sigma}$ consists of degenerate steady states and has the following approximation*

[TABLE]

Moreover there is an open neighborhood $U\in H$ of [math] such that $\Sigma_{\sigma}$ attracts $U\setminus\Gamma_{\sigma}$ where $\Gamma_{\sigma}$ is the stable set of [math] of codimension $2l_{c}+1$.

• •

If $q_{l_{c}}<0$ then there is a drastic (Type-II) transition at $\sigma=\sigma_{c}$ and a repeller bifurcates on $\sigma<\sigma_{c}$ and no steady state solution bifurcates on $\sigma>\sigma_{c}$.

4.1. Proof of the Main Theorem

The proof lies mainly on the reduction of the equation (11) onto the center manifold which is tangent to the critical eigenspace

[TABLE]

near the criticality $\sigma=\sigma_{c}$ and a careful analysis of the reduced equations. The condition $x_{-m}=(-1)^{m}\overline{x_{m}}$ is required in light of (24) so that $E_{1}$ consists of real valued functions.

Let

[TABLE]

Now a crucial element of the reduction procedure is the determination of the approximation of the center manifold function which captures the local dynamics near the criticality $\sigma=\sigma_{c}$. We will make use of the following approximation of the center manifold, see [10],

[TABLE]

in our analysis. Here $P_{2}:H\rightarrow E_{2}:=\{\Psi\in H_{1}\mid\langle\Psi,\Psi^{1\ast}_{l_{c},m,1}\rangle=0\text{ for all }\lvert m\rvert\leq l_{c}\}$ is the canonical projection onto the orthogonal complement, $\mathcal{L}$ is the restriction of the linear operator $L$ in (11) onto $E_{2}$ and

[TABLE]

where $\lvert x\rvert^{2}=\sum_{m=-l_{c}}^{l_{c}}\lvert x_{m}\rvert^{2}$.

According to (32), the center manifold is determined by the following nonlinear interactions of the critical modes

[TABLE]

where

[TABLE]

for $\Psi_{j}=(u_{j},w_{j},T_{j},S_{j})$, $j=1,2$.

Due to (33), the integral of triple product of spherical harmonics

[TABLE]

plays a crucial role in the determination of the above nonlinear interactions. The coefficients $c_{l_{1},l_{2},l_{3},m_{1},m_{2},m_{3}}$ are related to Clebsch-Gordan coefficients and are zero if

[TABLE]

or

[TABLE]

Since the critical modes always have wave index 1 in the z-direction, the product in (33) will vanish unless $n=0$ or $n=2$. As a summary of the above remarks, we obtain the following expansion of the center manifold function

[TABLE]

where by (32),

[TABLE]

Next we carry out the computation of the center manifold coefficients in (35) for $l_{c}=1$ and $l_{c}=2$ separately. In both cases, we find

[TABLE]

for all $1\leq l\leq 2l_{c}$ and $\lvert m\rvert\leq l_{c}$.

From (34) and (36), we may write the center manifold for $l_{c}=1$ as

[TABLE]

where,

[TABLE]

and $f_{l}^{k}$, $c_{l}^{k}$ are given by (28).

From (34) and (36), we may write the center manifold for $l_{c}=2$ as

[TABLE]

Here

[TABLE]

and $f_{l}^{k}$, $c_{l}^{k}$ are given by (28).

Now, we are ready to write down the reduced equations. We plug in $\Psi=\Phi_{c}+\Phi_{cm}$ into (11) and project the resulting equation onto $E_{1}$ to obtain

[TABLE]

Since $\int_{0}^{1}s_{1}(z)s_{2}(z)s_{3}(z)dz=0$ where $s_{i}(z)=\cos\pi z$ or $s_{i}(z)=\sin\pi z$, one gets right away that

[TABLE]

Also as $\Phi_{cm}=O(\lvert x\rvert^{2})$, one gets $G(\Phi_{cm},\Phi_{cm})=O(\lvert x\rvert^{4})$ and hence the equations (37) become

[TABLE]

for each $m\in\{-l_{c},\dots,l_{c}\}$.

We compactly present the results of tedious computations, both manual and by symbolic computation software, of the nonlinear terms in (38) for the cases $l_{c}=1$ and $l_{c}=2$ separately below. For the $l_{c}=1$ case, the reduced equations become

[TABLE]

where $q_{1}$ is defined as in (30).

For the $l_{c}=2$ case, the reduced equations are

[TABLE]

where $q_{2}$ is defined as in (31).

Since $x_{-m}x_{m}=(-1)^{m}\frac{1}{2}(\lvert x_{-m}\rvert^{2}+\lvert x_{m}\rvert^{2})$, letting

[TABLE]

the reduced equations (39) and (40) can be written compactly as

[TABLE]

By taking the product of (41) with $\overline{x_{m}}$ and summing from $m=-l_{c}$ to $m=l_{c}$, we derive

[TABLE]

Since the critical eigenspace is $2l_{c}+1$ dimensional which is odd, using Krasnosel’skii Theorem (see Theorem 1.10 in [8]) the existence of a bifurcated nontrivial steady state solution of the main equations at $\sigma=\sigma_{c}$ can be shown exactly as in [17].

In the case $q_{l_{c}}>0$ and $\sigma>\sigma_{c}$ case, from (42), it follows that $\lvert x\rvert^{2}\to\frac{\beta^{1}_{l_{c}1}(\sigma)}{q_{l_{c}}}$ as $t\to\infty$ for all sufficiently small initial conditions. By the attractor bifurcation theorem in [8], the bifurcated attractor $\Sigma_{\sigma}$ is homeomorphic to the $2l_{c}$ dimensional sphere $S^{2l_{c}}$. Moreover since the main equations possesses $S^{2l_{c}}$-symmetry, this steady state solution will generate a $S^{2l_{c}}$ set of steady states. Thus $\Sigma_{\sigma}$ consists solely of steady state solutions which are all degenerate since the Jacobian determinant of the right hand side of (41) vanishes at the criticality $\beta^{1}_{l_{c}1}(\sigma)=0$ for those steady state solutions. That proves the assertions of our main theorem.

5. Numerical Computations of the Transition Number

As shown in our main theorem, the type of first transition to multiple equilibria depends on a nondimensional parameter $q_{l_{c}}$, $l_{c}\in\mathbb{Z}^{+}$ which in turn depends on five system parameters: the aspect ratio $r$, the Lewis number Le, the Prandtl number Pr, the thermal Rayleigh number R, the saline Rayleigh number $\tilde{\text{R}}$. These parameters are related by the equation $\text{R}-\text{Le}^{-1}\tilde{\text{R}}=\sigma_{c}$, where $\sigma_{c}$ is given in (2) at the onset of transition, leaving four degrees of freedom in the determination of $q_{l_{c}}$. Hence, in the following discussion we fix $\tilde{\text{R}}$ by the choice of other four parameters. Also the parameter regime we are interested in our main theorem is the region $K>0$ where $K$ is given by (1). According to (27) this corresponds to either $\text{Le}>1$ regime or $\text{Le}<1$ and $\text{R}<\text{R}_{0}=\dfrac{\text{Le}+\text{Pr}}{(1-\text{Le})\text{Pr}}\sigma_{c}$ regime, see Figure 1.

Our main theorems for $l_{c}=1$ and $l_{c}=2$ cases prove that the type of transition is governed by the parameters $q_{1}=d_{(1,1),(0,2)}+d_{(1,1),(2,2)}$ and $q_{2}=d_{(2,1),(0,2)}+d_{(2,1),(2,2)}+d_{(2,1),(4,2)}$, respectively. We recall that $d_{(l_{c},1),(l,n)}$ is the coefficient arising from the nonlinear interactions of the critical modes with all the modes corresponding to the $\beta_{ln}$ eigenvalue.

We first present the numerical values of these nonlinear interaction terms when $\text{Pr}=7.5$, $\text{Le}=10^{-2}$, $\text{Le}=10^{-1}$, $\text{Le}=0.5$, $\text{Le}=5$ and $\text{R}=620$, $\text{R}=640$ for the case $l_{c}=1$ (with $r=2/\pi$) in Table 1 and $l_{c}=2$ (with $r=2\sqrt{3}/\pi$) in Table 2.

Our numerical results in Table 1 and Table 2 strongly indicate that

[TABLE]

Thus the sign of $d_{(l_{c},1),(0,2)}$ given by (29) plays a crucial role in the determination of the transition to multiple equilibria. So we investigate the sign of $d_{(l_{c},1),(0,2)}$ in detail now. Recalling that $\text{R}_{1}=\frac{\sigma_{c}}{1-\text{Le}^{2}}$, we observe that $\text{R}_{0}>\text{R}_{1}>0$ if $\text{Le}<1$ and $\text{R}_{1}<0$ if $\text{Le}>1$. Thus we obtain, see also Figure 2,

[TABLE]

The above remarks suggest that $\operatorname{sign}(q_{1})=\operatorname{sign}(d_{(1,1),(0,2)})$ and $\operatorname{sign}(q_{2})=\operatorname{sign}(d_{(2,1),(0,2)})$ at least when R is away from $\text{R}_{1}=\frac{\sigma_{c}}{1-\text{Le}^{2}}$. Hence the contribution to the transition number $q_{l_{c}}$ of the nonlinear interactions of the critical modes with the higher frequency modes is negligible compared to the interactions with zero wavenumber modes.

For $\text{Le}=10^{-2},$ $\text{Le}=10^{-1}$, $\text{Le}=0.5$ and in both $l_{c}=1$ and $l_{c}=2$ cases, we find that the transition is Type-I ($q_{l_{c}}>0$) if $\text{R}<\text{R}_{\ast}$ and Type-II ($q_{l_{c}}<0$) if $\text{R}_{\ast}<\text{R}<\text{R}_{0}$ as shown in Table 3, Table 4 and Figure 3 with $R_{\ast}\approx R_{1}=\frac{\sigma_{c}}{1-\text{Le}^{2}}$. The numerical values of $\text{R}_{1}$ when $\sigma_{c}=\frac{27}{4}\pi^{4}$ is $\text{R}_{1}=657.577$ when $\text{Le}=10^{-2}$, $\text{R}_{1}=664.153$ when $\text{Le}=10^{-1}$ and $\text{R}_{1}=876.682$ when $\text{Le}=0.5$. For $\text{Le}=5$, $q_{l_{c}}$ does not change sign and the transition is Type-I ($q_{l_{c}}>0$) for all R, see Figure 3.