A Mathematical Study of Flow Dynamics and Stability of Bounded Cartesian Plumes

1. Introduction

When a binary fluid consisting of two components with different densities is cooled from below, the component with the higher freezing temperature solidifies first. The resulting solid crystals accumulate at the base of the container, forming a porous layer composed of solid crystals and the remaining lighter fluid [6]. This porous region is commonly referred to as a mushy layer [15]. If the crystallized material corresponds to the denser component, the crystals settle further and create a solid layer beneath the mushy region [4]. As solidification proceeds, the mushy layer increases in thickness and may eventually become unstable [25]. This instability can generate compositional plumes, where the lighter component rises from the mushy layer through the overlying melt [9]. These plumes are examples of directed fluid motion occurring within a surrounding fluid of different composition and physical properties. Such flows are known as compositional plumes (Al Mashrafi & Eltayeb, 2014a).

The mathematical study of compositional plume dynamics is essential for understanding the behavior of fluid mixtures and alloys subjected to thermal and pressure gradients. Such understanding has significant implications in a variety of practical applications, including industrial processes such as metal casting, geophysical phenomena such as solidification at the Earth’s inner core boundary and mantle convection, and environmental processes such as salt-finger formation and sea-ice dynamics.

In industrial applications, one of the challenges encountered in iron casting is the formation of defects known as freckles [21]. During the casting process, trapped gases escape through narrow channels within the molten metal. Upon solidification, these channels appear as thin dark streaks running along the iron bar. The presence of such imperfections can significantly reduce the mechanical strength and structural integrity of the final product.

In geophysical systems, the Earth’s outer core is believed to consist primarily of molten iron and nickel, together with smaller amounts of lighter elements such as sulfur, oxygen, hydrogen, and helium, while the inner core is largely solid iron and nickel [16]. Although temperature generally increases with depth within the Earth, it has been proposed that the inner core continues to grow through solidification occurring at the inner core boundary (ICB) because of extremely high pressures [2,24]. During this process, iron crystallizes while the lighter elements remain in liquid form, leading to the formation of a mushy layer. The dense iron crystals settle onto the ICB, contributing to the gradual growth of the inner core [18]. This mushy layer may become unstable [19], giving rise to narrow ascending columns of light fluid in the form of compositional plumes (Figure 1). These plumes are believed to interact with the Earth’s magnetic field and may play a role in sustaining the geodynamo responsible for generating the planet’s magnetic field [17].

Compositional plumes also occur in environmental settings, particularly in the phenomenon known as salt fingers (Figure 2). In polar and subpolar regions, situations can arise in which relatively warm and fresh water overlies colder and saltier water. Because heat diffuses more rapidly than salt, the upper water layer loses heat faster than it loses salinity [22]. Consequently, the surface water becomes colder while retaining a relatively high salt concentration, increasing its density. The denser fluid then sinks in the form of narrow compositional plumes, creating structures known as salt fingers [14]. This process contributes significantly to heat and salt transport in the oceans and therefore influences climate dynamics in cold regions.

The contrasting effects of compositional plumes in different applications, being detrimental in industrial casting processes while playing beneficial roles in geophysical and environmental systems, have motivated extensive research into their formation, evolution, and stability. Understanding the dynamics of these plumes remains an important objective in fluid mechanics and related scientific disciplines.

The dynamics of compositional plumes have been investigated both theoretically and experimentally. Experimental studies (see Figure 3) suggest that plume flow is generally stable [20], although some experiments have reported unstable behavior [5]. In addition, these studies showed that plume behavior depends on its position within the container, particularly whether it is located near a wall or away from it [13].

On the theoretical side, several studies have been undertaken to investigate the dynamics and stability of compositional plumes. The work of Eltayeb and Loper [9] was the first to examine the stability of compositional plumes. They investigated the stability of vertical interfaces across which an imposed compositional discontinuity existed in the presence of a vertically oriented stabilizing temperature gradient. They subsequently extended their analysis to more realistic plume configurations [10,11] . These studies consistently concluded that compositional plumes are inherently unstable. Furthermore, plume instability has been shown to persist in the presence of rotation [8], magnetic fields [12], and the combined effects of rotation and magnetic fields [7]. However, all these theoretical investigations were conducted in unbounded domains, and material diffusion was assumed to be negligible.

Comparison between theoretical predictions and experimental observations is complicated by differences in the fluid domain surrounding the plume. The present study therefore aims to investigate theoretically the influence of boundaries on plume dynamics. A simplified model is developed that incorporates the principal factors expected to affect the rise of a compositional plume in a bounded domain. The model is illustrated in Figure 4. A vertical column of fluid of finite thickness $2x_0$ rises through a surrounding fluid of different compositions and is bounded on either side by rigid vertical walls separated by a distance $d=a_1+a_2$ The coordinate system is chosen such that the plume interfaces are located at $x=x_0 , x=-x_0$ while the walls are situated at $x=a_1 , x=-a_2$

2. Formulation of the Model

We consider a two-component incompressible fluid consisting of a solvent (light material) with concentration $C$and temperature $T$. The two constituents are assumed to possess identical kinematic viscosity, $\nu$, and thermal diffusivity, $\kappa$, while molecular diffusion of the material component is neglected. Following Eltayeb and Loper [9], the fluid motion is described by the governing equations of momentum, mass conservation, heat transfer, concentration transport, and the equation of state.

These equations are

${\rho }_0\left[{\displaystyle \frac{\partial u}{\partial t}}+(u.\nabla )u\right]=-\nabla p+{\rho }_0\nu {\nabla }^{2}u-\rho g\widehat{\mathit{z}},$ (1)

$\nabla .u=0,$ (2)

${\displaystyle \frac{\partial T}{\partial t}}+u.\nabla T=\kappa {\nabla }^{2}T,$ (3)

${\displaystyle \frac{\partial C}{\partial t}}+u.\nabla C=0,$ (4)

${\displaystyle \frac{\rho }{{\rho }_0}}=1-\alpha \left(T-T_0\right)-\beta \left(C-C_0\right),$ (5)

Here, $\mathit{u}$denotes the velocity vector, $p$the pressure, $g$ the constant gravitational acceleration, $\widehat{\mathit{z}}$the unit vector directed vertically upward, and $t$the time. The parameters $\alpha$and $\beta$represent the coefficients of thermal and compositional expansion, respectively, while $\rho$ denotes the fluid density. Quantities with the subscript $0 : \left({\rho }_0,T_0,C_0\right)$refer to constant reference values. Throughout this study, the fluid is assumed to satisfy the Boussinesq approximation, whereby density variations are neglected in the governing momentum equations except in the buoyancy term, where they contribute to the gravitational force.

Under the Boussinesq approximation, density variations are neglected throughout the governing equations except in the buoyancy term of the momentum equation, where they contribute to the gravitational force. The system admits a hydrostatic equilibrium state. The equations (1) - (5) allow a hydrostatic balance governed by

${\displaystyle \frac{d{p}_h}{dz}}+\rho g=0, {\displaystyle \frac{d^2{T}_h}{d{z}^2}}=0, {\mathit{u}}_h=0, C_h=C_0.$ (6)

Motivated by experimental studies of plumes originating from mushy layers, the temperature distribution is prescribed as

$T_h=T_0+\gamma z,$ (7)

where $\gamma$is a positive constant. This temperature profile increases with height, producing a thermally stable stratification. Consequently, any instability that develops is expected to arise primarily from compositional effects rather than thermal convection.

To express the governing equations in dimensionless form, appropriate characteristic scales are introduced. Since the focus is on instabilities associated with the ascent of buoyant fluid, the maximum concentration of the light material in the basic state, $\tilde{C}$, is selected as the concentration scale. To preserve the influences of both thermal and compositional stratification, the characteristic temperature scale is chosen such that the thermal and compositional buoyancy contributions remain comparable. Then the unit of temperature is

$T_u={\displaystyle \frac{\beta }{\alpha }}\widetilde{ C}$ (8)

Because the observed buoyant plume is narrow, viscous forces must remain of the same order as the buoyancy forces. This requirement leads to the selection of characteristic length and velocity scales based on the classical salt-finger length scale [9,23]. These scales are given by

$L={\left({\displaystyle \frac{\nu \kappa }{\alpha \gamma g}}\right)}^{{\displaystyle \frac{1}{4}}}, U=\beta \tilde{C}{\left({\displaystyle \frac{g\kappa }{\alpha \gamma \nu }}\right)}^{{\displaystyle \frac{1}{2}}}$ (9)

The associated convective time scale is then adopted as the relevant timescale for analyzing the growth of instabilities [1].

$t_c={\displaystyle \frac{L}{U}}$ (10)

A corresponding pressure scale is also introduced as

$\tilde{p}={\rho }_0\beta \tilde{C}{\left({\displaystyle \frac{\nu g^3\kappa }{\alpha \gamma }}\right)}^{{\displaystyle \frac{1}{4}}}$ (11)

Then the governing equations (1) - (5) are subsequently recast in nondimensional form:

$R\left[{\displaystyle \frac{\partial u}{\partial t}}+(u.\nabla )u\right]=-\nabla \left(p+{\displaystyle \frac{z}{\beta \tilde{C}}}\right)+{\nabla }^{2}u+\left(T-T_0+C-C_0\right)\widehat{\mathit{z}},$ (12)

$\nabla .u=0,$ (13)

$R\sigma \left[{\displaystyle \frac{\partial T}{\partial t}}+u.\nabla T\right]={\nabla }^{2}T,$ (14)

${\displaystyle \frac{\partial C}{\partial t}}+u.\nabla C=0,$ (15)

The resulting dimensionless system contains two principal parameters: the Grashoff number and the Prandtl number. The Grashoff number measures the ratio of buoyancy forces generated by concentration differences between the plume and the ambient fluid to viscous forces, whereas the Prandtl number represents the ratio of momentum diffusivity to thermal diffusivity. These two dimensionless numbers are written in the form

$R={\displaystyle \frac{UL}{\nu }}, \sigma ={\displaystyle \frac{\nu }{\kappa }}.$ (16)

A Cartesian coordinate system is adopted throughout the study to describe the flow field, with the 𝑧 -axis oriented vertically upward and the 𝑥 - and 𝑦 -axes lying in the horizontal plane. The dependent variables are decomposed into hydrostatic, basic-state, and perturbation components. The hydrostatic contribution accounts for the equilibrium pressure and temperature distributions, while the basic-state variables depend only on the horizontal coordinate 𝑥, representing a steady vertical plume. Small-amplitude perturbations are then introduced into this basic state to investigate its stability characteristics. Then the variables of the system (12) − (15) take the form

$u(x,y,z,t)=0+\bar{w}\left(x\right)\widehat{\mathit{z}}+\epsilon u^{†}(x,y,z,t),$ (17)

$C(x,y,z,t)=C_0+\bar{C}\left(x\right)+\epsilon C^{†}(x,y,z,t),$ (18)

$p(x,y,z,t)=p_h+\bar{p}\left(x\right)+\epsilon p^{†}(x,y,z,t),$ (19)

$T(x,y,z,t)=T_h+\bar{T}\left(x\right)+\epsilon T^{†}(x,y,z,t),$ (20)

such that the variables with subscript h have hydrostatic contribution and given by

$T_h=T_0+{\displaystyle \frac{(z-z_0)}{\sigma R}},{ p}_h=p_0-{\displaystyle \frac{\left(z-z_0\right)}{\beta \tilde{C}}}+{\displaystyle \frac{{\left(z-z_0\right)}^2}{2\sigma R}},$ (21)

in which $z_0$ is a reference value.

The boundary conditions are specified as follows. First, the heat flux, momentum flux, and all flow variables except concentration are continuous across the plume interfaces. Second, the bounding walls are rigid and maintained at the hydrostatic temperature distribution. Third, the plume interfaces are treated as material surfaces. Since material diffusion is assumed negligible, no transfer of light material occurs across these interfaces, consistent with the analysis of Al Mashrafi and Eltayeb [1].

Substituting the decomposed variables (17) - (20) into the governing equations (12) - (15) and collecting terms independent of the perturbation amplitude yields the equations governing the basic state. These equations are given by

$-{\displaystyle \frac{d\bar{p}}{dx}}\widehat{\mathit{x}}+\left({\displaystyle \frac{d^2\bar{w}}{d{x}^2}}+\bar{C}+\bar{T}\right)\widehat{\mathit{z}}=0,$ (22)

${\displaystyle \frac{d^2\bar{T}}{d{x}^2}}=\bar{w}.$ (23)

Retaining terms proportional to the perturbation amplitude of order $\epsilon$produces the linearized perturbation equations. These equations form the basis of the stability analysis that are used to determine the behavior and growth of disturbances within the plume. These perturbation equations are

$R\left[{\displaystyle \frac{\partial u^{†}}{\partial t}}+\bar{w}\widehat{\mathit{z}}.\nabla u^{†}+\left(u^{†}.\nabla \bar{w}\right)\widehat{\mathit{z}}\right]=-\nabla p^{†}+{\nabla }^2{u}^{†}+\left(T^{†}+C^{†}\right)\widehat{\mathit{z}},$ (24)

$\nabla .u^{†}=0,$ (25)

$\sigma R\left[{\displaystyle \frac{\partial T^{†}}{\partial t}}+\bar{w}{\displaystyle \frac{\partial T^{†}}{\partial z}}+u^{†}.\nabla \bar{T}\right]+u^{†}.\widehat{\mathit{z}}={\nabla }^2{T}^{†},$ (26)

${\displaystyle \frac{\partial C^{†}}{\partial t}}+\bar{w}{\displaystyle \frac{\partial C^{†}}{\partial z}}+u^{†}.\nabla \bar{C}=0.$ (27)

3. Basic State Solution

The governing equation (15) is identically satisfied by the basic state, allowing the concentration distribution to be specified freely. In this study, a top-hat concentration profile is adopted.

$\bar{C}\left(x\right)=\{\begin{array}{c}1, \left|x\right|\le x_0\\ 0, -a_2\le x<-x_0,x_0<x\le a_1\end{array}.$ (28)

By introducing the expression

$F\left(x\right)=\bar{T}\left(x\right)-i\bar{w}\left(x\right),$ (29)

the governing equations (22) and (23) reduce to

$\bar{p}=0, {\displaystyle \frac{d^{2}F}{d{x}^2}}-iF=i\bar{C},$ (30)

subject to the boundary conditions

$\begin{array}{c}\left(i\right) F, {\displaystyle \frac{dF}{dx}} \mathrm{a}\mathrm{r}\mathrm{e} \mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{i}\mathrm{n}\mathrm{u}\mathrm{o}\mathrm{u}\mathrm{s} \mathrm{a}\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s} x=\pm x_0 \\ \left(\mathrm{i}\mathrm{i}\right) F=0 \mathrm{a}\mathrm{t} x=a_1,-a_2\end{array}\}.$ (31)

The solution of the system (28) – (31) is given by

$F\left(x\right)=\{\begin{array}{c}-A\mathit{sin}h(k{a}_1\left)\mathit{sin}h[k\right(x+a_2\left)\right]; -a_2\le x<-x_0\\ -A\mathit{sin}h(k{a}_1\left)\mathit{sin}h[k\right(x+a_2\left)\right]+\mathit{cos}h[k(x+x_0)]-1; -x_0\le x\le x_0\\ A\mathit{sin}h(k{a}_2\left)\mathit{sin}h[k\right(x-a_1\left)\right];{ x}_0<x\le a_1\end{array},$ (32)

where $A$ and $k$ are given by

$A={\displaystyle \frac{2\mathit{sin}h(k{x}_0)}{\mathit{sin}h(kd)}}, k={\displaystyle \frac{1}{\sqrt{2}}}\left(1+i\right), d=a_1+a_2.$ (33)

The basic-state solution was computed numerically, and representative results are presented in Figure 5 and Figure 6 [1]. Figure 5 demonstrates the influence of the sidewalls on the plume structure. For a plume of fixed thickness, decreasing the separation between the sidewalls intensifies the plume-induced flow. When the plume is positioned closer to one wall than the other, the downward flow outside the plume becomes weaker in the region adjacent to the nearest wall and stronger on the opposite side. As a result, the flow field loses its symmetry.

Figure 6 highlights the effect of plume thickness and wall proximity on the flow characteristics. The results indicate that the sidewalls exert a significant influence on the solution. When the plume is located at the midpoint between the walls, the profiles remain symmetric. However, this symmetry is progressively disrupted as the plume approaches one of the sidewalls. In addition, the oscillatory behaviour of the velocity field can produce regions of negative velocity, corresponding to downward motion within the plume when the plume thickness is sufficiently large. This phenomenon alters the overall transport of material and affects the net flux carried by the plume.

4. The Stability Analysis

To examine the linear stability of the basic-state solution given by (32) , we use the perturbation equations (24) - (27). The interface located at the plane $x=x_0$is assumed to be subjected to a small-amplitude harmonic disturbance, as illustrated in Figure 7, of the form

$x=x_0+\epsilon \mathit{exp}\left(\mathsf{\Omega }t+i(my-nz)\right)+c.c.,$ (34)

where $m$and $n$denote the horizontal and vertical wavenumbers, respectively. The symbol c.c. represents the complex conjugate, while $\mathsf{\Omega }$is a constant amplitude parameter that may conveniently be written as

$\mathsf{\Omega }={\mathsf{\Omega }}_{\mathrm{r}}+\mathrm{i}{\mathsf{\Omega }}_{\mathrm{i}}.$. (35)

This representation enables the disturbance field to be expressed as a normal-mode solution, thereby facilitating the analysis of the growth or decay of perturbations and the determination of the stability characteristics of the basic flow.

The real component ${\mathsf{\Omega }}_{\mathrm{r}}$controls the temporal evolution of the disturbance amplitude and, consequently, determines its stability characteristics. When the real part remains negative for all possible combinations of wavenumbers $m$ and $n$, the plume is considered stable. In contrast, the plume becomes unstable if at least one combination of wavenumbers $m$ and $n$produces a positive value. A neutrally stable state exists when the real part ${\mathsf{\Omega }}_r$is zero for all wavenumber values.

If the most amplified mode is associated with both wavenumbers $m$ and $n$ being non-zero, the disturbance is classified as a three-dimensional (3D) mode. When $m=0$, the disturbance is referred to as a two-dimensional (2D) mode. The situation in which $n=0$ and $m\ne 0$ is not observed. The imaginary component ${\mathsf{\Omega }}_i$, on the other hand, governs the propagation of the disturbance by determining its phase velocities. The vertical phase speed $u_z$ and the horizontal phase speed $u_h$ are defined as follows.

$u_z={\displaystyle \frac{{\mathsf{\Omega }}_i}{n}}, u_h={\displaystyle \frac{{\mathsf{\Omega }}_i}{m}}={\displaystyle \frac{n{u}_z}{m}}.$ (36)

It should be noted that $u_h$is defined only when $m\ne 0$ .

The disturbance (34) propagates through the fluid, influencing the second interface and inducing perturbations in the system variables. As a result, the interface located at $x=-x_0$can be expressed in the form

$x=-x_0+\epsilon {\eta }_1\mathit{exp}\left(\mathsf{\Omega }t+i(my-nz)\right)+c.c.,$ (37)

where ${\eta }_1$represents the amplitude of the interface displacement, as determined by the boundary conditions of the problem, while the perturbation variables can be expressed in the form

$\left\{u^{†},C^{†},T^{†},p^{†}\right\}=\left\{-inu,nmv,w, C,T,-inp\right\}\mathit{exp}\left(\mathsf{\Omega }t+i(my-nz)\right)+c.c.,$ (38)

where the factors $-in$$nm$$-in$ are introduced into the variables $u$$v$ and $p$respectively, for notational convenience.

Substituting the expressions given by (39) into Equations (24) – (27) yields the following system of ordinary differential equations in the variable$x$

$Du-m^{2}v+w=0,$ (40)

$\mathsf{\Delta }u-Dp=R\bar{\mathsf{\Omega }}u,$ (41)

$\Delta v-p=R\bar{\mathsf{\Omega }}v,$ (42)

$\mathsf{\Delta }w+T+C+n^{2}p=R\left(\bar{\mathsf{\Omega }}w-inuD\bar{w}\right),$ (43)

$\mathsf{\Delta }T-w=\sigma R\left(\bar{\mathsf{\Omega }}T-inuD\bar{T}\right),$ (44)

$\bar{\mathsf{\Omega }}C=0,$ (45)

where

$b^2=m^2+n^2, D\equiv {\displaystyle \frac{d}{dx}}, \mathsf{\Delta }\equiv D^2-b^2, \bar{\mathsf{\Omega }}=\mathsf{\Omega }-in\bar{w}.$ (46)

It is convenient to derive the following three equations

$\mathsf{\Delta }\varsigma =R\left(-invD\bar{w}+\bar{\mathsf{\Omega }}\varsigma \right),$ (47)

$n^{2}u=m^2\varsigma -D\left(w+p\right)-R\bar{\mathsf{\Omega }}u,$ (48)

$\mathsf{\Delta }p-T-C=2inRuD\bar{w},$ (49)

The first equation (47) is obtained by differentiating (50) with respect to $x$and subsequently subtracting (51). The second equation (52) is derived by differentiating (53) once and then subtracting (54). The third equation (55) is obtained by applying the operator $\mathsf{\Delta }$ to (56) and making use of Equations (40) - (57). The variable $\varsigma$, which is associated with the vertical component of the vorticity, is defined by

$\varsigma =\left(Dv-u\right)={\displaystyle \frac{1}{nm}}\widehat{\mathit{z}}.curlu.$ (58)

It follows from this equation (59) that the perturbation concentration vanishes everywhere within the fluid domain. Consequently

$C=0$ (60)

The system is subject to the following boundary conditions

$D\left(w+p\right),v,w,T,p,Dv,DT are continuous across x=\pm x_0,$ (61)

$⟨Dw\left(x_0\right)⟩=⟨\bar{C}\left(x_0\right)⟩, ⟨Dw(-x_0)⟩={\eta }_1⟨\bar{C}(-x_0)⟩,$ (62)

$-inu\left(x_0\right)=\mathsf{\Omega }-in\bar{w}\left(x_0\right), -inu(-x_0)=\left[\mathsf{\Omega }-in\bar{w}(-x_0)\right]{\eta }_1.$ (63)

$D\left(w+p\right)=v=w=T=\varsigma =0 at x=a_1,x=-a_2,$ (64)

where $⟨f⟩$ denotes the operator defined by

$⟨f\left(\alpha \right)⟩=f\left({\alpha }^{-}\right)-f\left({\alpha }^{+}\right),$ (65)

representing the discontinuity, or jump, across $x=\alpha$. It should be noted that the governing equation (66)has been used to replace the continuity condition on $u$ together with its vanishing at the boundaries by the condition $D\left(w+p\right)$, since $\varsigma$ is likewise continuous across the interfaces and vanishes at the boundaries. The jump conditions (67) express the requirement that the interfaces behave as material surfaces.

Earlier studies on compositional plumes have shown that the flow is unstable when the Grashoff number, $R,$is small [10]. Since the Grashoff number measures the strength of the plume through the maximum amplitude of the basic concentration distribution, instability at small $R$ implies instability throughout the entire range of possible $R$ values. Accordingly, the variables of the system and $\mathsf{\Omega }$ may be expanded in powers of the small parameter $R$in the form

$f\left(x,y,z,t\right)={\sum }_{s=0}^{\infty }{f}_s\left(x,y,z,t\right)R^s, \mathsf{\Omega }={\sum }_{s=1}^{\infty }{\mathsf{\Omega }}_s{R}^{s-1},R<<1,$ (68)

where $f\left(x,y,z,t\right)$stands for any of the perturbation variables $u$, $v$, $w$, $p$ or $T$.

Upon substituting the expansion (69) into Equations (40) - (44), (47) - (49) and the associated boundary conditions (61) - (64) , and requiring the coefficients of $R^s$ ($s$ = 0 , 1 , 2 , ….) to vanish independently, a sequence of ordinary differential equation systems is obtained. Solving these systems successively provides the expression for the growth rate. To leading order, the stability characteristics of the interfaces are determined by the two systems corresponding to $n=0$(hereafter referred to as Problem 0) and $n=1$(hereafter referred to as Problem 1).

4.1. Problem 0

The coefficients of ${ R}^0$in the system of R⁰ in the system of equations (42) − (44), (47) − (49) are given as follows

$D{u}_0-m^2{v}_0+w_0=0,$ (70)

$\mathsf{\Delta }{w}_0+T_0+n^2{p}_0= 0,$ (71)

$\mathsf{\Delta }{T}_0-w_0= 0,$ (72)

$\mathsf{\Delta }{\varsigma }_0=0,$ (73)

$n^2{u}_0=m^2{\varsigma }_0-D\left(w_0+p_0\right),$ (74)

$\mathsf{\Delta }{p}_0-T_0= 0,$ (75)

subject to the following boundary conditions

$D(p_0+w_0)=w_0={\varsigma }_0=T_0=0 at x=a_1,x=-a_2,$ (76)

${\varsigma }_0,w_0,T_0,p_0,D{T}_0,D\left(p_0+w_0\right) are continuous across x=\pm x_0,$ (77)

$⟨D{w}_0\left(x_0\right)⟩=⟨\bar{C}\left(x_0\right)⟩, ⟨D{w}_0(-x_0)⟩={\eta }_1⟨\bar{C}(-x_0)⟩,$ (78)

$-in{u}_0\left(x_0\right)={\mathsf{\Omega }}_1-in\bar{w}\left(x_0\right), -in{u}_0(-x_0)=\left[{\mathsf{\Omega }}_1-in\bar{w}(-x_0)\right]{\eta }_1.$ (79)

The system, (70) - (75) subject to the prescribed boundary conditions, can be solved directly using the elimination method to determine the unknown variables. Subsequently, applying the boundary conditions (79) yields expressions for the growth rate ${\mathsf{\Omega }}_1$ and the interface displacement ${\eta }_1$ [1]. This results in the following expression

${\mathsf{\Omega }}_1={\displaystyle \frac{in}{2}}\left(-S_1\pm \sqrt{{\mathsf{\Delta }}_p}\right), {\mathsf{\Delta }}_p=S_1^2-4S_2,$ (80)

${\eta }_1={\displaystyle \frac{-N_{j-}}{\left(\frac{{\Omega }_1}{in}\right)-\bar{w}(-x_0)+M_{j+}}},$ (81)

where

$S_1=N_{j+}+M_{j+}-\bar{w}\left(x_0\right)-\bar{w}(-x_0),$ (82)

$S_2=\left\{N_{j+}-\bar{w}\left(x_0\right)\right\}\left\{M_{j+}-\bar{w}(-x_0)\right\}-N_{j-}{M}_{j-},$ (83)

$N_{j\pm }={\sum }_{j=1}^3{\displaystyle \frac{-{\lambda }_j{F}_j}{\mathit{sin}h({\lambda }_{j}d)}}\mathit{sin}h\{{\lambda }_j(x_0-a_1)\left\}\mathit{cos}h\{{\lambda }_j\right(x_0\pm a_2)\},$ (84)

$M_{j\pm }={\sum }_{j=1}^3{\displaystyle \frac{-{\lambda }_j{F}_j}{\mathit{sin}h({\lambda }_{j}d)}}\mathit{sin}h\{{\lambda }_j(x_0-a_2)\left\}\mathit{cos}h\{{\lambda }_j\right(x_0\pm a_1)\},$ (85)

with

$F_j={\displaystyle \frac{{\mu }_j^2}{{\lambda }_j(2{\mu }_j+3n^2)}},$ (86)

with ${\mu }_j$ $\left(j=\mathrm{1,2},3\right)$ representing the three roots of the cubic equation

${\mu }_j^3+{\mu }_j+n^2=0, where {\lambda }_j=\sqrt{{\mu }_j+b^2}.$ (87)

The nature of the roots of the cubic equation (87) ensures that the relevant expressions $N_{j\pm }$ and $M_{j\pm }$ are real. Consequently, the quantities $S_1$ and $S_2$ are real , implying that the discriminant ${\mathsf{\Delta }}_p$ appearing in the equation (80) is also real. When the discriminant is positive, the equation admits two distinct real roots; when it vanishes, the roots are real and identical. In both cases, ${\mathsf{\Omega }}_1$is purely imaginary, indicating that the system is neutrally stable. Therefore, it is necessary to proceed to the next-order approximation (Problem 1). Here, ${\mathsf{\Omega }}_1$is given by

${\Omega }_1^{\left(a\right)}=i{\displaystyle \frac{n}{2}}\left(-S_1+\sqrt{{\Delta }_p}\right),{\mathsf{\Omega }}_1^{\left(b\right)}=i{\displaystyle \frac{n}{2}}\left(-S_1-\sqrt{{\Delta }_p}\right),$ (88)

where the superscripts $a$ and $b$ are used to distinguish between the two modes.

When the discriminant is negative, the equation possesses a pair of complex-conjugate roots. Consequently, the corresponding values of ${\mathsf{\Omega }}_1$ have real parts with opposite signs. The root associated with the positive real part characterizes the unstable mode of the system. The two roots can be expressed as follows:

${\mathsf{\Omega }}_1^{\left(a\right)}={\displaystyle \frac{n}{2}}\left\{-\sqrt{-{\mathsf{\Delta }}_p}-i{\mathrm{S}}_1\right\},{ \mathsf{\Omega }}_1^{\left(b\right)}={\displaystyle \frac{n}{2}}\left\{\sqrt{-{\mathsf{\Delta }}_p}-i{S}_1\right\}.$ (89)

Using the expression (81), the amplitudes ${\eta }_1^{\left(a\right)}$ and ${\eta }_1^{\left(b\right)}$ can be expressed as

${\eta }_1^{(k)}={\displaystyle \frac{-N_{j-}}{\left({\Omega }_{\hspace{0.17em}1}^{\hspace{0.17em}(k)}/in\right)-\hspace{0.17em}\hspace{0.17em}\overline{w}(-x_0)+M_{j+}}}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}k=\{\hspace{0.17em}a\hspace{0.17em},\hspace{0.17em}b\hspace{0.17em}\}.$ (90)

In the absence of sidewall boundaries, the two modes correspond to configurations in which the plume interfaces oscillate either in phase, resulting in a sinuous mode, or out of phase, resulting in a varicose mode (see Figure 8). In both cases, ${\mathsf{\Omega }}_1$is purely imaginary, implying that the disturbances are neutrally stable at the present order of approximation. The introduction of sidewalls breaks this symmetry, except in the special case where the plume is positioned equidistantly between the boundaries. However, when the discriminant is real, both modes remain neutrally stable and propagate with distinct phase speeds while retaining either an in-phase or an out-of-phase structure. By contrast, when the discriminant is negative, the plume becomes unstable and ${\eta }_1^{\left(k\right)}$assumes complex values. Consequently, a phase difference arises between the two interfaces, such that they exhibit neither an in-phase nor an out-of-phase relationship.

It follows from the governing equations (89) that, if one of the two modes exhibits a positive growth rate, the other must necessarily possess a negative growth rate. Consequently, at most one of the modes can be unstable. Since $n>0$, the MS mode is identified as the unstable mode. Figure 9 illustrates the behaviour of the MS mode in the wavenumber plane for representative values of $x_0,a_1,a_2$. It is evident that, as the plume thickness decreases, the region of the wavenumber plane corresponding to instability expands, although the maximum growth rate is reduced. Conversely, as the plume becomes thicker, the unstable region contracts and eventually vanishes when $x_0$ exceeds $0.6$. The instability is characterized by a two-dimensional mode, and the vertical wavenumber decreases with increasing $x_0$.

in the $\left(m,n\right)$plane for $d=10$. Results are shown for $(x_0,a_2)$ = (a) $\left(\mathrm{0.1,0.25}\right)$ , (b) $\left(\mathrm{0.2,0.4}\right)$ , (c) $\left(\mathrm{0.3,0.5}\right)$ , and (d) $\left(\mathrm{0.4,0.55}\right)$. Instability occurs within the contour, while the flow remains neutrally stable outside the curve ${\mathsf{\Omega }}_1=0$.

The preferred instability mode is identified as the one corresponding to the maximum growth rate (76) as a function of the wavenumbers $m$ and $n$for fixed values of the parameters $x_0$ , $a_1$ and $a_2$. In other words, we solve the governing differential equations and determine the mode that exhibits the largest growth rate:

${\displaystyle \frac{\partial }{\partial m}}{\mathsf{\Omega }}_1^{\left(k\right)}\left(x_0,a_1,a_2\right)=0, {\displaystyle \frac{\partial }{\partial n}}{\mathsf{\Omega }}_1^{\left(k\right)}(x_0,a_1,a_2)=0.$ (91)

Maximization of the growth rate (89) reveals that instability occurs only for the MS mode when $a_2-x_0\le 0.25$. Furthermore, the unstable mode is two-dimensional $\left(m_c=0\right)$and propagates vertically upwards $(U_c>0)$.

A representative set of profiles for the preferred mode is presented in Figure 10. The growth rate increases from zero at the boundary, reaches a maximum at a certain distance from the wall, $a_0 ,$and subsequently decreases. This distance $a_0$ depends on both the plume thickness and the parameters $a_2, d$. The vertical wavenumber, $n_c$, decreases as the plume moves away from the wall when the plume is sufficiently thin; however, for relatively thick plumes exhibit the opposite behavior and increases with distance from the wall.

The instability region in the $\left(x_0,a_2-x_0\right)$plane is illustrated in Figure 11. It can be observed that the unstable region expands with increasing plume thickness, reaches a maximum extent, and then gradually shrinks as the plume thickness approaches zero. In general, the plume remains neutrally stable throughout the domain, except in the case of very thin plumes located close to the wall, where instability may arise. We next investigate the stability characteristics of Problem 1.

plane, corresponding to growth rates of order $O\left(1\right)$. The region marked $N_1$ represents neutral stability, while the region marked $U_1$indicates the onset of instability. In all unstable cases, the preferred mode is the modified sinuous (MS) mode; the modified varicose (MV) mode is found to be stable.

4.2. Problem 1

The coefficients associated with the perturbation terms of $R^1$ in the governing equations yield the following set of equations

$\mathsf{\Delta }{u}_1= M_u,$ (92)

$\mathsf{\Delta }{v}_1-p_1= M_V,$ (93)

$\mathsf{\Delta }{w}_1+T_1+n^2{p}_1= M_w,$ (94)

$\mathsf{\Delta }{T}_1-w_1= M_T,$ (95)

$\mathsf{\Delta }{p}_1-T_1= M_p,$ (96)

$n^2{u}_1=m^2{\varsigma }_1-D\left(w_1+p_1\right)-{\bar{\mathsf{\Omega }}}_1{u}_0,$ (97)

$\Delta {\varsigma }_1=-im\left(w_0-n^2{v}_0\right)D\bar{w},$ (98)

where we have defined the expressions

$\begin{array}{c}{M}_u=D{p}_1 +\left({\mathsf{\Omega }}_1-in\bar{w}\right)u_0, M_v= \left({\mathsf{\Omega }}_1-in\bar{w}\right)v_0\\ M_w= \left({\mathsf{\Omega }}_1-in\bar{w}\right)w_0-in{u}_{0}D\bar{w}, M_p=2in{u}_{0}D\overline{w }\\ M_T= \sigma \left({\mathsf{\Omega }}_1-in\bar{w}\right)T_0-in\sigma u_{0}D\bar{T} \end{array}\}.$ (99)

The governing equations are subject to the following boundary conditions

$v_1=w_1=T_1={\varsigma }_1=D(p_1+w_1)=0 at x=a_1,x=-a_2,$ (100)

$v_1,w_1,T_1,p_1,{\varsigma }_1,D{\varsigma }_{1}D{v}_1,D{T}_1,D{w}_1,D\left(p_1+w_1\right) are continuous across x=\pm x_0,$ (101)

${\mathsf{\Omega }}_2=-in{u}_1\left(x_0\right),{ \mathsf{\Omega }}_2{\eta }_1=-in{u}_1(-x_0)$ (102)

The governing equations, together with the boundary conditions (92) - (102), can be solved by deriving the solvability condition for the resulting nonhomogeneous system[1]. This analysis yields an expression for the growth rate${\mathsf{\Omega }}_2$which may be written in the form

${\mathsf{\Omega }}_2^{\left(k\right)}={\displaystyle \frac{-1}{1+{\left({\eta }_1^{\left(k\right)}\right)}^2}}\left\{{\mathsf{\Omega }}_{21}^{\left(k\right)}+{\mathsf{\Omega }}_{22}^{\left(k\right)}+{\mathsf{\Omega }}_{23}^{\left(k\right)}+{\mathsf{\Omega }}_{24}^{\left(k\right)}+{\displaystyle \frac{i}{n}}{\widehat{g}}^{\left(k\right)}\right\}.$  (103)

In the above expression, $k$takes the modes $MS$ or $MV$, corresponding to the sinuous and varicose modes, respectively. The quantities ${\mathsf{\Omega }}_{21}^{\left(k\right)}$ , ${\mathsf{\Omega }}_{22}^{\left(k\right)}$ , ${\mathsf{\Omega }}_{23}^{\left(k\right)}$ , ${\mathsf{\Omega }}_{24}^{\left(k\right)}$ and ${\widehat{g}}^{\left(k\right)}$ are given by the following integrals:

${\mathsf{\Omega }}_{21}^{\left(k\right)}=-i{\mathsf{\Omega }}_1^{\left(k\right)}{\int }_{-a_2}^{a_1}\left\{n{u}_0^{\left(k\right)}{G}^{\left(k\right)}-{\displaystyle \frac{1}{n}}{w}_0^{\left(k\right)}{H}^{\left(k\right)}\right\}dx,$ (104)

${\mathsf{\Omega }}_{22}^{\left(k\right)}={\displaystyle \frac{i{\Omega }_1^{\left(k\right)}}{n}}{\sum }_{j=1}^3{C}_j{\int }_{-a_2}^{a_1}\left\{w_0^{\left(k\right)}+\sigma {\mu }_j{T}_0^{\left(k\right)}\right\}{H}_j^{\left(k\right)}dx,$ (105)

${\mathsf{\Omega }}_{23}^{\left(k\right)}={\int }_{-a_2}^{a_1}\left(-u_0^{\left(k\right)}\left\{n^2\bar{w}{G}^{\left(k\right)}+H^{\left(k\right)}D\bar{w}\right\}+\bar{w}{w}_0^{\left(k\right)}{H}^{\left(k\right)}\right)dx,$ (106)

${\mathsf{\Omega }}_{24}^{\left(k\right)}={\sum }_{j=1}^3-C_j{\int }_{-a_2}^{a_1}{H}_j^{\left(k\right)}\left[-\bar{w}\left\{w_0^{\left(k\right)}+\sigma {\mu }_j{T}_0^{\left(k\right)}\right\}+\left\{\left(2{\mu }_j^2+1\right)D\bar{w}-\sigma {\mu }_{j}D\bar{T}\right\}{u}_0^{\left(k\right)}\right]dx,$ (107)

${\widehat{g}}^{\left(k\right)}=g^{\left(k\right)}(-a_2)+g^{\left(k\right)}\left(a_1\right),$ (108)

where $G^{\left(k\right)}$, $H^{\left(k\right)}$, $H_j^{\left(k\right)}$, $C_j$, $g^{\left(k\right)}(-a_2)$, $g^{\left(k\right)}\left(a_1\right)$ are specified as follows

$G^{\left(k\right)}\left(x\right)={\displaystyle \frac{1}{b}}\{\begin{array}{c}-F_{1b}^{\left(k\right)}\mathit{sin}h[b(x+a_2)]; -a_2\le x<-x_0\\ \left(-F_{1b}^{\left(k\right)}\mathit{sin}h[b(x+a_2)]-{\eta }_1^{\left(k\right)}\mathit{sin}h[b(x+x_0\left)\right]\right); -x_0\le x\le x_0\\ F_{2b}^{\left(k\right)}\mathit{sin}h[b(x-a_1)];{ x}_0<x\le a_1\end{array},$ (109)

$H^{\left(k\right)}\left(x\right)=\{\begin{array}{c}-F_{1b}^{\left(k\right)}\mathit{cos}h[b(x+a_2)]; -a_2\le x<-x_0\\ -F_{1b}^{\left(k\right)}\mathit{cos}h[b(x+a_2)]-{\eta }_1^{\left(k\right)}\mathit{cos}h[b(x+x_0\left)\right]; -x_0\le x\le x_0\\ F_{2b}^{\left(k\right)}\mathit{cos}h[b(x-a_1)];{ x}_0<x\le a_1\end{array},$ (110)

$H_j^{\left(k\right)}\left(x\right)=\{\begin{array}{c}-F_{1j}^{\left(k\right)}\mathit{cos}h[{\lambda }_j(x+a_2)]; -a_2\le x<-x_0\\ -F_{1j}^{\left(k\right)}\mathit{cos}h[{\lambda }_j(x+a_2)]-{\eta }_1^{\left(k\right)}\mathit{cos}h[{\lambda }_j(x+x_0\left)\right]; -x_0\le x\le x_0\\ F_{2j}^{\left(k\right)}\mathit{cos}h[{\lambda }_j(x-a_1)];{ x}_0<x\le a_1 \end{array},$ (111)

$C_j={\displaystyle \frac{-n^2}{3n^2+2{\mu }_j}},$ (112)

$g^{\left(k\right)}\left(a_1\right)={\sum }_{j=1}^3-C_j{\mu }_j{F}_{2j}^{\left(k\right)}\left(D{T}_1^{\left(k\right)}\left(a_1\right)+{\mu }_j{p}_1^{\left(k\right)}\left(a_1\right)\right),$ (113)

$g^{\left(k\right)}(-a_2)={\sum }_{j=1}^3-C_j{\mu }_j{F}_{1j}^{\left(k\right)}\left(D{T}_1^{\left(k\right)}(-a_2)+{\mu }_j{p}_1^{\left(k\right)}(-a_2)\right),$ (114)

and $F_{1b}^{\left(k\right)}$, $F_{2b}^{\left(k\right)}$, $F_{1j}^{\left(k\right)}$ and $F_{2j}^{\left(k\right)}$ are

$F_{1b}^{\left(k\right)}={\displaystyle \frac{-1}{\mathit{sin}h(bd)}}\left\{{\eta }_1^{\left(k\right)}\mathit{sin}h(b(a_1+x_0))+\mathit{sin}h(b(a_1-x_0\left)\right)\right\},$ (115)

$F_{2b}^{\left(k\right)}={\displaystyle \frac{-1}{\mathit{sin}h(bd)}}\left\{{\eta }_1^{\left(k\right)}\mathit{sin}h(b(a_2-x_0))+\mathit{sin}h(b(a_2+x_0\left)\right)\right\},$ (116)

$F_{1j}^{\left(k\right)}={\displaystyle \frac{-1}{\mathit{sin}h({\lambda }_{j}d)}}\left\{{\eta }_1^{\left(k\right)}\mathit{sin}h({\lambda }_j(a_1+x_0))+\mathit{sin}h({\lambda }_j(a_1-x_0\left)\right)\right\},$ (117)

$F_{2j}^{\left(k\right)}={\displaystyle \frac{-1}{\mathit{sin}h({\lambda }_{j}d)}}\left\{{\eta }_1^{\left(k\right)}\mathit{sin}h({\lambda }_j(a_2-x_0))+\mathit{sin}h({\lambda }_j(a_2+x_0\left)\right)\right\}.$ (118)

The growth rate ${\mathsf{\Omega }}_2^{\left(k\right)}$ defined by (103), can be decomposed into five distinct contributions, each reflecting the influence of specific system parameters and variables. The terms ${\mathsf{\Omega }}_{21}^{\left(k\right)}$ and ${\mathsf{\Omega }}_{23}^{\left(k\right)}$ arise from interactions between the perturbation variables and the growth rate of the leading-order system, whereas ${\mathsf{\Omega }}_{23}^{\left(k\right)}$ and ${\mathsf{\Omega }}_{24}^{\left(k\right)}$ result from interactions between the leading-order variables and the basicstate quantities. While ${\mathsf{\Omega }}_{21}^{\left(k\right)}$ and ${\mathsf{\Omega }}_{23}^{\left(k\right)}$ are independent of the Prandtl number,${\mathsf{\Omega }}_{22}^{\left(k\right)}$ and ${\mathsf{\Omega }}_{24}^{\left(k\right)}$ depend linearly on $\sigma$ The term ${\widehat{g}}^{\left(k\right)}$ represents the contribution induced by the presence of boundaries. Since the leading-order variables are real and the first-order variables are imaginary, the resulting growth rate ${\mathsf{\Omega }}_2^{\left(k\right)}$ is purely real and therefore governs the amplification or decay of the disturbance.

Numerical computations of the growth rate (103) indicate that the plume is unstable, with disturbances growing at a rate of$O\left(R\right)$. For a given location in parameter space $\left(x_0,a_2,\sigma \right)$, the preferred mode, corresponding to the maximum growth rate, may be either the modified sinuous (MS) mode or the modified varicose (MV) mode, depending on the complex interplay among the governing parameters. When one parameter is varied while the remaining two are held fixed, a transition in the preferred mode may occur, with the dominant instability changing from MS to MV, or vice versa, once a critical parameter value is reached.

Furthermore, variations in a parameter can cause a mode of a given type, whether MS or MV, to undergo a transition between two-dimensional and three-dimensional structures, or conversely, as the parameter passes through a critical threshold. This behavior arises because the growth-rate expression (103) may possess multiple local maxima. As the parameter increases, the initially dominant maximum may decrease while a secondary maximum increases, until a crossover point is reached at which the latter exceeds the former and becomes the preferred mode. Figure 12 illustrates a representative example of this behavior for a selected set of parameter values.

, $\sigma$ = 10 , $d=10$, and $a_2=3$ for (a), (b) and $a_2=5$ for (c), (d). Note that the MS mode is preferred for $a_2=3$ and the MV mode is preferred when $a_2=5$.

Figure 13 illustrates the dependence of the preferred instability mode on the Prandtl number, $\sigma$, and facilitates comparison with the limiting case in which sidewalls are absent. For low values of $\sigma$, the modified sinuous (MS) mode is the preferred mode of instability, whereas the modified varicose (MV) mode becomes dominant as $\sigma$increases. This behavior is consistent with the unbounded configuration reported by Eltayeb and Loper [10]. The critical Prandtl number, ${\sigma }_0$ , at which the preferred mode switches from MS to MV depends on the distance between the plume and the nearest sidewall. As the wall approaches the plume, ${\sigma }_0$increases, indicating that the presence of sidewall boundaries suppresses the MV mode relative to the MS mode. The boundaries also exert a stabilizing influence on the plume, as evidenced by the reduction in the magnitude of the growth rate with decreasing wall separation $d$. It is noteworthy that, regardless of the values of the governing parameters $d$ and $\sigma$, the preferred MS mode remains three-dimensional, whereas the preferred MV mode remains two-dimensional when the plume is located midway between the two sidewalls.

, for $x_0=2$, when the plume is positioned midway between the two sidewalls $\left(a_1=a_2\right)$. Curves (i) and (ii) correspond to two different sidewall spacings: (i) $d=10$and (ii) $d=20$. The solid curve represents the modified sinuous (MS) mode, whereas the dashed curve represents the modified varicose (MV) mode.

In contrast to the unbounded plume, which is unstable of $O\left(R\right)$, the bounded Cartesian plume exhibits two distinct instability regimes, characterized by growth rates of orders $O\left(1\right)$and$O\left(R\right)$. The region associated with the larger growth rate, $O\left(1\right)$, occupies only a small portion of the parameter space and is strongly influenced by the spacing of the sidewalls. The corresponding regime diagram is presented in Figure 14 for $d=10$. Instability with growth rate $O\left(1\right)$occurs only when the plume is relatively thin, with a thickness not exceeding approximately one-half of the salt-finger length scale, and when its distance from the nearest sidewall is less than about 0.25. Furthermore, as the plume approaches the wall, the growth rate decreases, indicating that the stabilizing influence of the boundary becomes increasingly significant in the immediate vicinity of the sidewall.

plane for $d=10$. In panel (a), the regions labeled $U_1$ and $U_2$ correspond to instabilities with growth rates of order $O\left(1\right)$and $O\left(R\right)$, respectively. The shaded area lies outside the physically admissible domain, since $x_0$ cannot exceeds $a_2$. Panel (b) presents an enlarged view of the region $a_2-x_0\le 0.25$ and $x_0\le 0.6$ from panel (a).

The preferred mode of instability is associated with plume interfaces represented by (34) and (37). The amplitude at the interface $x=x_0$is fixed at unity, whereas the amplitude at $x=-x_0$governed by the ratio ${\eta }_1$, whose value is determined by the parameters defining the preferred mode for given values of $x_0,a_2,\sigma ,d$. Sample interface profiles for selected preferred modes are displayed in Figure 15. It is observed that the interfaces approach each other closely at periodically spaced positions along the plume, suggesting the possibility of plume breakup and the subsequent formation of discrete blob-like structures.

at $d=10$. For visual clarity, both profiles are displayed using the same magnification factor $\epsilon$($=0.1)$. (a) $x_0=0.1$, $a_2=0.2$, $n_c=2.08, m_c=0$ and (b) $x_0=0.5$ , $a_2=0.6$, $n_c=0.63,m_c=0$. Notice that the interfaces come into close proximity at periodically recurring positions along the domain.

5. Conclusions

The mathematical model describing the dynamics of bounded Cartesian plumes has been investigated. This study investigates the dynamics of compositional plumes in fluids of infinite extent. To gain insight into the influence of lateral boundaries on plume stability and evolution, a simplified Cartesian geometry is employed. The corresponding unbounded configuration was previously examined by Eltayeb and Loper [9,10,11]. This formulation facilitates a direct comparison between plumes evolving in bounded and unbounded domains, thereby highlighting the effects introduced by the presence of boundaries.

The effect of boundaries on the Eltayeb-Loper model is investigated by introducing two vertical sidewalls on either side of the ascending buoyant fluid column. The inclusion of these boundaries introduces two additional dimensionless parameters into the problem: the distance between the sidewalls, $d$, and the distance from the plume to the nearest sidewall, $a_2$. In the absence of sidewalls, the basic state is symmetric with respect to the $x-$axis, measured normal to the plume interfaces. Consequently, the perturbations can be classified into two independent symmetry classes: the even mode, commonly referred to as the varicose (V) mode, and the odd mode, known as the sinuous (S) mode. The introduction of sidewalls generally breaks this symmetry whenever the plume is not positioned equidistantly between the boundaries. Nevertheless, the two classes of disturbances persist and are subsequently identified as the modified varicose (MV) mode and the modified sinuous (MS) mode.

The presence of sidewalls alters the basic state of a bounded Cartesian plume, with its behavior strongly influenced by the distance between the plume and the nearest boundary. As a result, the stability characteristics of the bounded plume differ significantly from those of its unbounded counterpart. The stability analysis reveals the emergence of a new instability mode when the plume is located sufficiently close to a sidewall. Unlike the unbounded case, where the growth rate scales as $O\left(R\right)$the order of magnitude of the growth rate in the bounded configuration depends on both the plume thickness and its proximity to the nearest wall. When the plume is thin and situated close to a boundary, the instability develops with a growth rate of $O\left(1\right)$. For other combinations of plume thickness and distance from the nearest sidewall, the plume remains unstable with a growth rate of the same order as that observed in the absence of boundaries. Nevertheless, even in this regime, the magnitude of the growth rate decreases as the separation between the sidewalls becomes smaller.