2.1. Experiments

2.1.1. Apparatus

Velocity field measurements were conducted for steady, temperature-driven convection in water in an RBC cell of cross-sectional area $300\times 300$ mm$^2$, as shown in figure 1(a). The bottom hot copper plate, painted black to prevent reflection of light from its surface, was kept above a thin glass plate, which was above two aluminium plates that had a nichrome wire heater sandwiched between them. A constant heat flux was provided by a Variac connected to the wire heater. The heat flux was estimated by measuring the temperature difference $(T_1-T_2)$ across the glass plate in the heater assembly using T-type thermocouples. The sides and the bottom of the convection cell were insulated. The top cold plate, made of glass, was maintained at a constant temperature by circulating refrigerated water in a glass compartment above the plate. A water filled glass prism was fixed on top of this glass compartment to rectify the radial distortions in stereoscopic imaging. The vertical location of the cold plate assembly – which included the cold plate, the glass compartment and the glass prism – could be adjusted. The temperature difference between the bottom plate ($T_B$) and the top plate ($T_T$), $\Delta T=T_B-T_T$, was measured using thermocouples touching the hot copper plate at the bottom and the cold glass plate at the top. Measurements were made only after a steady state was attained after around three hours, when $\Delta T$ remained a constant, as monitored by a data logger (Agilent, model 34970A). The water layer height $H$ and the constant heat flux through the bottom copper plate were changed to vary the Rayleigh number $Ra=g\beta \Delta T H^3/\nu \alpha$, where $g$ is the gravitational acceleration, $\beta$ the coefficient of thermal expansion, $\nu$ the kinematic viscosity and $\alpha$ the thermal diffusivity. The values of $H$, $\Delta T$ and the corresponding estimated values of $Ra$ and Prandtl number ($Pr = \nu /\alpha$) are shown in table 1. Figure 1(b) shows an example of the pattern of line plumes observed on the hot surface, obtained by electrochemical dye visualisation (Baker Reference Baker1966) at $Ra= 8.26\times 10^6$ and $Pr= 6$ (Gunasegarane Reference Gunasegarane2015).

Figure 1.

(a) Schematic of the experimental set-up. (b) Top view of the line plumes at $Ra=8.26\times 10^6$ and $Pr=6$ detected by electrochemical dye visualisation (Gunasegarane Reference Gunasegarane2015). See Movie 1 available at in the supplementary movies for corresponding velocity field in RBC experiments.

Table 1.

Values of experimental parameters, length scales and PIV parameters. Physical properties of the fluid were estimated at $T_b$, where $T_b$ is the temperature of the bulk fluid. For PIV, one side of the square interrogation window $D_I$ was equal to 32 pixels and the overlap was equal to $50\,\%$. Vertical velocity at the centreline of a plume, estimated from the similarity solutions of Gebhart et al. (Reference Gebhart, Pera and Schorr1970), is shown as $w_{pc}$. Here $V_{sh}$ is the large-scale flow velocity given by (3.25).

2.2. Numerical simulations

The non-dimensional governing equations under the Boussinesq, incompressible approximation are

(2.1)\begin{gather} \boldsymbol{\nabla}^*\boldsymbol{\cdot} \boldsymbol{V}^*=0, \end{gather}

(2.2)\begin{gather} \frac{{\rm D}\boldsymbol{V}^*}{{\rm D}t^*} ={-}\boldsymbol{\nabla}^*P^*+ \left( \sqrt{\frac{Pr}{Ra}}\right) {\nabla}^*{^2} \boldsymbol{V}^* + T^* \boldsymbol{e_z}, \end{gather}

(2.3)\begin{gather} \frac{{\rm D} T^*}{{\rm D}t^*}=\frac{1}{ \sqrt{Ra\,Pr}}{\nabla}^*{^2} T^*, \end{gather}

where the lengths are normalised by $H$, time by $H/u_f$, $\boldsymbol {V}^*$ is the dimensionless velocity field obtained by normalising by the free-fall velocity $u_f=\sqrt {g\beta \Delta T H}$, $T^*=(T-T_C)/\Delta T$ and $\boldsymbol {e_z}$ is the unit vector in the $z$ direction. Here $T^*$ varies from a constant value of one at the bottom surface to zero at the top plate. The computations were done for a cylindrical domain with adiabatic sidewalls and an aspect ratio, $\varGamma = D/H$ of one, for $Ra = 2\times 10^7, 2\times 10^8$ and $Pr=1,\, 4.96, 10$.

The numerical schemes described in Verzicco & Camussi (Reference Verzicco and Camussi2003) are used for solving (2.1)–(2.3). In brief, we use a second-order finite difference scheme discretised on a staggered mesh in a cylindrical domain. The time integration is done using an explicit third-order Runge–Kutta scheme (Verzicco & Camussi Reference Verzicco and Camussi1997; Verzicco & Orlandi Reference Verzicco and Orlandi1996) for nonlinear terms, and implicit Crank–Nicholson scheme for the linear terms, as described in the fractional step method discussed in Kim & Moin (Reference Kim and Moin1985). The discretised equations are solved by a fractional-step procedure. The pressure equation is inverted using trigonometric expansion in the azimuthal direction and FISHPACK package (Swartzrauber Reference Swartzrauber1974) in the other two directions. The grid size $N_{\theta }\times N_r \times N_z$ used in the simulations is $513\times 129\times 385$. This ensures $33$ grid points in the thermal boundary layer for $(Ra,Pr)=(2\times 10^8,10)$, satisfying the criterion suggested by Shishkina et al. (Reference Shishkina, Stevens, Grossmann and Lohse2010). The maximum grid size in the vertical direction is $0.006$ and the Courant–Friedrichs–Lewy number was fixed at ${\rm CFL}=0.5$ with a maximum time step ${\rm d}t=0.0331$.

3.1. Theoretical velocity distributions

3.1.1. Plume

For a steady, 2D, laminar line plume, shown as the unhatched region in the schematic in figure 2(a), denoting the plume variables by the subscript $p$, the dimensional horizontal and vertical velocity profiles are (Gebhart et al. Reference Gebhart, Pera and Schorr1970)

(3.1)\begin{gather} u_p(z,\eta_p)={-}4^{3/4}\frac{\nu}{z}Gr_z^{1/4}\left(\frac{3}{5}f_p(\eta_p)-\frac{2}{5}\eta_p f_p'(\eta_p) \right), \end{gather}

(3.2)\begin{gather} w_p(z,\eta_p)=4^{1/5}\left(\frac{g\beta Q_p}{C_p I}\right)^{2/5}\left(\frac{z}{\mu \rho}\right)^{1/5} f_p'(\eta_p). \end{gather}

Here, $f_p(\eta _p)$ is the dimensionless stream function obtained as the numerical solution of the self-similar plume equations; an approximation of $f_p$ is given in Appendix A (A1). In (3.1) and (3.2), a prime denotes differentiation with respect to the similarity variable $\eta _p$, defined as

(3.3)\begin{equation} \eta_p=\frac{x_p}{z}\left(\frac{Gr_z}{4}\right)^{1/4}, \end{equation}

where, $x_p$ is the horizontal distance from the plume centreline and $Gr_z=g\beta (T_0-T_\infty )z^3/\nu ^2$ is the Grashoff number based on the vertical distance $z$ from the hot plate. Here, $T_0$ and $T_\infty$ denote the temperatures at the plume centreline and the ambient, respectively. In (3.2), the heat flux per unit length of the line source is

(3.4)\begin{equation} Q_p=4^{3/4} (g \beta)^{1/4}\sqrt{\rho\mu} IC_pN^{5/4}, \end{equation}

where $N$ is a constant, $C_p$ is the specific heat at constant pressure and $I =\int _{-\infty }^\infty f_p'\phi \,{\rm d}\eta _p$, with $\phi =(T-T_\infty )/(T_0-T_\infty )$ and

(3.5)\begin{equation} T_0-T_\infty=Nz^{{-}3/5}. \end{equation}

Using (3.4) and (3.5) in (3.2) and (3.1), we obtain the dimensionless horizontal and vertical velocities within the plume as

(3.6)\begin{gather} u_p^*= \frac{u_p}{\nu\sqrt{Gr_z}/z} ={-}4^{3/4}Gr_z^{{-}1/4}\left(\frac{3}{5}f_p-\frac{2}{5}\eta_pf_p' \right), \end{gather}

(3.7)\begin{gather}w_p^*=\frac{w_p}{\nu\sqrt{Gr_z}/z}=2f_p', \end{gather}

where $f_p$ and $f_p'$ are given by (A1) and (A2).

Figure 2(b) shows the profiles of dimensionless horizontal and vertical velocities, (3.6) and (3.7), across a plume. The plume has an upward vertical velocity, which has a maximum at the centreline and reaches zero at its edge. Differentiating (3.7) with respect to $z$, we obtain the dimensionless vertical gradient of vertical velocity as

(3.8)\begin{equation} \frac{\partial w_p^*}{\partial z^*}= \frac{\partial w_p/\partial z}{\nu\sqrt{Gr_z}/z^2}=\frac{2}{5}\left( f_p'-2\eta_p f_p''\right). \end{equation}

From (3.8), we see that the plume region, shown as the unhatched region in figure 2(a), has a positive vertical spatial acceleration, $\partial w_p/\partial z\sim 1/z^{4/5}$; the positive vertical acceleration increases as we approach the hot plate. Simultaneously, owing to continuity, $\partial u_p/\partial x_p=-\partial w_p/\partial z\sim -1/z^{4/5}$ is a negative quantity. With $u_p$ being zero at the plume centreline, the horizontal velocities within the plumes are negative and decrease in magnitude as we approach the centre of the plumes from its edges. In summary, $\partial w_p/\partial z>0$ and ${\partial u_p/\partial x_p<0}$ in the unhatched plume region of figure 2.

3.1.2. Boundary layer

For a steady, 2D natural convection boundary layer, shown as the $45^\circ$ hatched regions in figure 2(a), the horizontal and vertical velocity profiles are (Pera & Gebhart Reference Pera and Gebhart1973)

(3.9)\begin{align} u_b&=5^{3/5}\frac{\nu}{x_b}Gr_x^{2/5}f_b'\quad{\rm and} \end{align}

(3.10)\begin{align} w_b&={-}\frac{\nu}{x_b}\left(\frac{Gr_x}{5}\right)^{1/5}\left(3f_b-2\eta_bf_b'\right), \end{align}

where the subscript $b$ denotes boundary layer variables. Here, $Gr_x=g\beta (T-T_\infty ) x_b^3/\nu ^2$ is the Grashoff number based on the horizontal distance $x_b$ within the boundary layer. Here $x_b = 0$ is the leading edge from where the boundary layer forms, and the boundary layer thickens as $x_b$ increases and we approach the plume. The dimensionless stream function $f_b(\eta _b)=\psi _b/(5\nu (Gr_x/5)^{1/5})$, an approximation of which is given by (A3), with $\psi _b$ being the stream function, and a prime denotes differentiation with respect to the boundary layer similarity variable,

(3.11)\begin{equation} \eta_b=\frac{z}{x_b}\left(\frac{Gr_x}{5}\right)^{1/5}. \end{equation}

Using (3.3), (3.11) can be rewritten in terms of $\eta _p$ as

(3.12)\begin{equation} \eta_b=\frac{1}{4^{1/10}5^{1/5}}\frac{Gr_z^{3/10}}{\eta_{p}^{2/5}}. \end{equation}

Rewriting (3.9) and (3.10) in dimensionless form using (3.12), we obtain the dimensionless horizontal and vertical velocities as

(3.13)$$\begin{gather} u_b^*=\frac{u_b}{\nu \sqrt{Gr_z}/z}=5^{3/5}4^{1/20}\frac{\eta_p^{1/5}}{Gr_z^{3/20}}f_b', \end{gather}$$

(3.14)$$\begin{gather}w_b^*=\frac{w_b}{\nu \sqrt{Gr_z}/z}=\frac{1}{5^{1/5}4^{1/10}}\frac{1}{\eta_p^{2/5}Gr_z^{1/5}}, \left(3f_b-2\eta_bf_b' \right), \end{gather}$$

where $f_b$ and $f_b'$ are given by (A3) and (A4).

The variations given by (3.13) and (3.14) along the boundary layer that feeds the plume at its base, as shown schematically in figure 2(a), are shown in figure 2(b). To plot (3.13) and (3.14) in figure 2(b), we write (A3) and (A4) in terms of $\eta _p$ using (3.12). Figure 2(b) shows that the horizontal velocity (towards the plume) inside the boundary layer region increases as we approach the plume, as also shown schematically in figure 2(a). From (3.9),

(3.15)\begin{equation} \frac{\partial u_b}{\partial x_b} =\frac{1}{5^{2/5}} \frac{\nu }{x_b^2} Gr_x^{2/5} \left(f_b'-2\eta_bf_b''\right) \end{equation}

is identified as a positive quantity, whose value decreases as $x_b^{-4/5}$ as we approach the plume from the leading edge of the boundary layer. This increase in horizontal velocity along the boundary layer occurs since the driving horizontal pressure gradient due to the density difference, $\Delta p \approx \rho g \beta \Delta T \delta _{nc}$, increases with increase in boundary layer thickness as we approach the plume. As discussed in § 3.1.1, inside the plume, the horizontal velocity decreases towards its centreline. Figure 2(b) shows that the boundary layer has a small vertical downward velocity far away from its leading edge, which increases to appreciable values closer to its leading edge; these boundary layers entrain ambient fluid into them. From (3.15), by continuity, $\partial w_b/\partial z\sim -x_b^{-4/5}$, a negative quantity. Hence, as shown schematically in the boundary layer regions in figure 2(a), for the regions within the boundary layers feeding the plumes on either side at the bottom of the plumes, $\partial u_b/\partial x_b >0$ and $\partial w_b/\partial z<0$.

In the ambient/bulk fluid region above the boundary layers, rising plumes cause flows due to the ambient fluid being entrained into them. Such entrainment flows are expected to have $\partial w/\partial z<0$ since the large-scale flow in the bulk is of an order higher in magnitude than the entrainment flows into the plumes and the boundary layers (Gunasegarane & Puthenveettil Reference Gunasegarane and Puthenveettil2014). Further, the kinematic blocking due to the presence of the bottom plate also slows down the downward vertical velocities as we approach the bottom plate. Hence, $\partial w/\partial z$ is expected to be negative in the regions shown as bulk in figure 2(a). In addition, the presence of plumes also results in increasing horizontal velocities in the bulk ($\partial u/\partial x>0$) due to the entrainment, as shown by the horizontal arrows in the bulk region in figure 2(a). We then see that the flow in and around a 2D laminar plume and a 2D laminar boundary layer feeding it at the bottom separates into two regions, viz. (i) the plumes where $\partial w/\partial z>0$ and $\partial u/\partial x<0$ and (ii) the bulk and the boundary layers where $\partial w/\partial z<0$ and $\partial u/\partial x>0$; these two regions are shown as the unhatched and the hatched regions, respectively, in figure 2(a). To separate the plumes from the boundary layers and the bulk in such flows, we then need a criterion that distinguishes the $\partial w/\partial z>0$, $\partial u/\partial x<0$ regions as plumes from the $\partial w/\partial z<0$, $\partial u/\partial x>0$ regions as bulk or boundary layers.

Extending the above considerations to a three-dimensional (3D) flow, the horizontal, 2D divergence of the horizontal velocity components, calculated in a plane parallel and close to the plate,

(3.16)\begin{equation} \boldsymbol{\nabla_H} \boldsymbol{\cdot} \boldsymbol{u}= \frac{\partial u}{\partial x} +\frac{\partial v}{\partial y}, \end{equation}

is the equivalent of $\partial u/\partial x$ in the similarity solutions, with the subscript $H$ indicating that the divergence operator is applied only on the horizontal ($x$ and $y$) components of the velocity field, namely $u$ and $v$. Since the flow is incompressible, $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u} = -\partial w/\partial z$; negative values of $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u}$ imply that such regions will have positive $\partial w/\partial z$ and vice versa. Further,

(3.17)\begin{equation} \boldsymbol{\nabla_H} \boldsymbol{\cdot} \boldsymbol{u} <0 \end{equation}

implies that the horizontal velocities reduce in magnitude along the direction of these velocities. As deduced in § 3.1.1, such regions are likely to be plumes. Similarly, regions with $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u} >0$ will have increasing horizontal velocities in the direction of these velocities as well as $\partial w/\partial z <0$; such regions are likely to be the ambient fluid or the boundary layer regions in-between plumes. We now analyse the nature of deformation of fluid elements within plumes and the associated boundary layers, in terms of the eigenvalues of the strain rate tensor, and arrive at the same kinematic condition (3.17) that would separate the plume regions from the boundary layer regions.

3.2. Deformation mechanism and strain rate tensor

Figure 3 shows the displacement and deformation of fluid elements as they pass from a local boundary layer to a plume, as observed in a vertical plane by using passive markers in the visualisation experiments of Weijermars (Reference Weijermars1988). Figure 2(a) shows a schematic representation of this observed evolution in shape of a spherical fluid element at three different locations along its path near the hot plate in RBC. The spherical fluid element is initially in the bulk (location 1), then moves into the boundary layer (location 2), and finally becomes a part of the plume (location 3). Balachandar (Reference Balachandar1992) has shown that the compressive strain rate orients along the vertical direction in the boundary layers. Then, in a horizontal plane inside the boundary layer, due to the horizontal shear forces present, the extensional strain rate is predominantly in the horizontal plane. Fluid elements get stretched along the extensional strain axes and show reluctance to be tilted by the vorticity field (Guala et al. Reference Guala, Lüthi, Liberzon, Tsinober and Kinzelbach2005). Hence, inside the boundary layer, the fluid element undergoes predominantly horizontal stretching, as shown in figures 2(a) and 3. At location 3 inside the plume (figure 2a), the extensional strain is vertical (Balachandar Reference Balachandar1992) due to the vertically accelerating plume, and the fluid element then gets stretched farther in the vertical direction, as could be noticed in figure 3. The extensional strain rate caused by this vertical stretching of the fluid element inside the plume, due to mass conservation, results in a compressive strain rate in the horizontal plane. Hence, in summary, the net horizontal strain rate transitions from being extensional to compressive as the fluid element travels from the boundary layer to the plume region.

Figure 3.

Side view of the observed deformations of circular blobs of markers near a hot surface at the bottom, as given by Weijermars (Reference Weijermars1988).

The extension or compression strain rates occurring at various points in a horizontal plane, which represent the deformation rates of fluid elements in that plane, can be described using the eigenvalues $\lambda _{1,2}$ of the 2D instantaneous strain rate tensor,

(3.18)\begin{equation} S= \begin{bmatrix} u_x & \frac{1}{2}(u_y+v_x) \\ \frac{1}{2}(u_y+v_x) & v_y \end{bmatrix}, \end{equation}

where $u$ and $v$ are the velocities in the $x$ and $y$ directions in the horizontal plane, and the subscripts $x$ and $y$ denote the gradients in these directions. The ordered eigenvalues of $S$ are

(3.19)\begin{gather} \lambda_1=\tfrac{1}{2}\left(\boldsymbol{\nabla_H} \boldsymbol{\cdot} \boldsymbol{u} +\gamma\right), \end{gather}

(3.20)\begin{gather}\lambda_2=\tfrac{1}{2}\left(\boldsymbol{\nabla_H} \boldsymbol{\cdot} \boldsymbol{u} -\gamma\right), \end{gather}

so that $\lambda _1>\lambda _2$; these eigenvalues provide the local principal strain rates at each location. In (3.19) and (3.20),

(3.21)\begin{equation} \gamma = \sqrt{\gamma_n^2+\gamma_s^2} \end{equation}

is the total strain rate in the horizontal plane, with $\gamma _n=u_x-v_y$ being the normal component of the strain rate, and $\gamma _s=u_y+v_x$ the shear component of the strain rate (Nolan, Mattia & Shane Reference Nolan, Mattia and Shane2020). The eigenvalues $\lambda _{1,2}$ take positive or negative values, corresponding to the extension or compression strain rates along the principal axes. The associated eigenvectors, denoted by $\boldsymbol {\xi _1}$ and $\boldsymbol {\xi _2}$, which are orthogonal to each other, give the lines of action of these pure extensional or compressive deformations of the fluid element at the corresponding location.

We define the dominant eigenvalue $\lambda _D$ at a spatial location as the eigenvalue that is larger in magnitude; $\lambda _D=\lambda _1$ if $|\lambda _1/\lambda _2|>1$ and $\lambda _D=\lambda _2$ if $|\lambda _1/\lambda _2|<1$. Since the boundary layers are dominated by extensional strains in the horizontal direction, $\lambda _D$ is expected to take positive values within the boundary layers. Similarly, since the vertical extensional strains inside the plumes result in compressive horizontal strains in a horizontal plane, $\lambda _D$ is expected to take negative values within the plumes. Hence, a spatial distribution of the dominant eigenvalue in a horizontal plane can distinguish plumes from boundary layers, with regions of

(3.22)\begin{equation} \lambda_D<0 \end{equation}

being the plumes, and those with

(3.23)\begin{equation} \lambda_D>0 \end{equation}

being the boundary layers. The relation between the signs of $\lambda _D$ and ${\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}}$ is depicted in table 2, based on the definitions in equations (3.19)–(3.20). With $\gamma >0$, it is straightforward to observe that $\lambda _D<0$ occurs when $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}<0$, also coinciding with $\lambda _2$ being the dominant eigenvalue. In other words, regions of negative horizontal divergence and negative $\lambda _D$ are identical, thus establishing a direct connection between the similarity-solution-based arguments in § 3.1 and the strain-rate-based arguments in the current subsection.

Table 2.

The signs of $\lambda _1$, $\lambda _2$, $\lambda _D$, and their relation to the sign of $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}$, as concluded based on (3.19) and (3.20).

Figure 4 shows the spatial distribution of $\lambda _D$ and the corresponding eigenvector (figure 4a,c), and the regions with negative values of the dimensionless horizontal divergence,

(3.24)\begin{equation} \zeta=\frac{\boldsymbol{\nabla_H} \boldsymbol{\cdot} \boldsymbol{u}}{V_{sh}/\bar{\lambda}}, \end{equation}

overlaid with the horizontal velocity field (figure 4b,d). Here, $V_{sh}$ is the large-scale flow velocity, estimated from the Reynolds number based on $V_{sh}$ as (Gunasegarane & Puthenveettil Reference Gunasegarane and Puthenveettil2014)

(3.25)\begin{equation} Re_{sh}=\frac{V_{sh}H}{\nu}=0.404Ra^{4/9}Pr^{{-}2/3}; \end{equation}

the values of $V_{sh}$ for the present range of $Ra$ are listed in table 1. In (3.24),

(3.26)\begin{equation} \bar{\lambda}=C_1Z_w Pr^{n_1} \end{equation}

is the mean spacing between the line plumes, where

(3.27)\begin{equation} Z_w=\left(\frac{\nu\alpha}{g\beta \Delta T_w}\right)^{1/3}=H\left(\frac{2}{Ra}\right)^{1/3} \end{equation}

is a length scale near the plate, $C_1=47.5$ and $n_1=0.1$ (Puthenveettil et al. Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011). The distributions in figure 4 are in a central, horizontal 73.32 $\times$ 61.32 mm$^2$ region, at a height of 1.5 mm from the hot surface in the experimental set-up shown in figure 1(a). The eigenvalues and eigenvectors were calculated from the PIV-based velocity field measurements in the plane within the Prandtl–Blasius and natural convection boundary layer thicknesses. Figure 4(a,b) and figure 4(c,d) correspond to ($Ra,Pr$)=($3.3\times 10^6$, 5.21) and ($1.21\times 10^9, 5.09$), respectively. As expected, the shape and extent of regions with negative $\zeta$ in figures 4(b) and 4(d) coincide with those having negative $\lambda _D$ in figures 4(a) and 4(c).

Figure 4.

(a) Negative dominant eigenvalues $\lambda _D$ of the 2D strain rate tensor $S$ (3.18), overlaid with the corresponding eigenvectors for $(Ra,Pr)=(3.33\times 10^6, 5.21)$; (b) corresponding distribution of the dimensionless horizontal divergence field, overlaid with the horizontal velocity vector field in a horizontal plane at a height of 1.5 mm from the hot plate. Panels (c,d) are the same as (a,b), but for $(Ra, Pr)= (1.21\times 10^9, 5.09)$. An arrow length of 1 mm is equal to the magnitude of horizontal velocity of 0.67 mm s$^{-1}$ in (b) while that in (d) is equal to 3.33 mm s$^{-1}$. All the plots in this figure are based on 2D3C PIV measurements. See Movies 2, 3 and 4 (see supplementary movies) showing different motions of plumes at three different $Ra$.

In figure 4(a), the domain splits into slender, line-like regions that have negative $\lambda _D$, surrounded by regions of positive $\lambda _D$. Note that the eigenvectors corresponding to $\lambda _D$ within the negative $\lambda _D$ regions in figure 4(a) are mostly perpendicular to the longitudinal boundaries of these regions. Owing to the orthogonality between the eigenvectors of $S$, the eigenvectors corresponding to the non-dominant eigenvalues within the negative $\lambda _D$ regions mostly align along the length of these regions (not shown here). The negative $\zeta$ regions in figure 4(b) have velocity vectors pointed towards them, as would be expected for plumes, since feeding of plumes by the boundary layers on their sides is predominant in the absence of strong large-scale flows at low $Ra$ (Gunasegarane & Puthenveettil Reference Gunasegarane and Puthenveettil2014).

A distribution of negative $\lambda _D$ in slender, line-like regions is also observed in figure 4(c) at a higher Rayleigh number of $Ra=1.2\times 10^9$. Figures 4(a) and 4(c) show that an increase in $Ra$ decreases the mean spacing between the negative $\lambda _D$ regions, as expected for line plumes from (3.26) and (3.27). On comparing with the velocity vector field shown in figure 4(d), the line-like regions in figure 4(c) could be seen to be more aligned with the shear near the hot plate, compared with that in figure 4(a). Further, the orientation of the eigenvectors corresponding to $\lambda _D$, shown in figure 4(c), is seen to be more tilted along the shear direction; in contrast, for lower $Ra$, the eigenvectors are more perpendicular to the length of the $\lambda _D<0$ regions in figure 4(a). A corresponding tilt of the eigenvectors of the non-dominant eigenvalue away from the longitudinal orientation of the $\lambda _D<0$ regions thus occurs at the higher $Ra$ (not shown here).

In summary, the regions of negative $\lambda _D$ at the smaller $Ra$ of $3.33\times 10^6$ have incoming flows with fluid elements that predominantly get compressed in directions perpendicular to the length of these regions. Such a flow and deformation is expected in a rising plume when the fluid elements get stretched vertically, resulting in dominant compression perpendicular to these plumes in a horizontal plane. These negative $\lambda _D$ regions tend to align along the dominant shear direction at higher $Ra$; plumes are known to align along the direction of large-scale flow (Puthenveettil & Arakeri Reference Puthenveettil and Arakeri2005; Shevkar et al. Reference Shevkar, Gunasegarane, Mohanan and Puthenveettil2019). At higher $Ra$, the mean spacing between these regions is smaller, resulting in the length of these regions being larger, similar to the trend observed for the length of plume regions. In the presence of shear at higher $Ra$, the compressive deformation inside these regions occurs more in the direction of shear. Such an effect would occur at larger $Ra$ since the external large-scale flow would result in a velocity field within the plume which is inclined to the longitudinal axes of the plumes (Shevkar et al. Reference Shevkar, Gunasegarane, Mohanan and Puthenveettil2019). We hence expect these regions of negative $\lambda _D$ (or negative $\zeta$) to be the plumes, while the regions with positive $\lambda _D$ (or positive $\zeta$) to be the boundary layers. Such an understanding of plume regions also imply that the boundaries of plume regions would occur at $|\lambda _1/\lambda _2|=1$, where both the extensional and the compressional strains are equal in magnitude. We now look at the connection between these regions of negative dominant eigenvalues (plumes) and LCS.

3.3. LCS and plumes

Lagrangian coherent structures is a dynamical systems-based framework to identify underlying robust features that organise tracer patterns resulting from finite time stirring in any fluid flow (Haller Reference Haller2015). While the relevance of stable and unstable manifolds for long-term advective stirring has long been known in steady, periodic and quasiperiodic flows (Beigie, Leonard & Wiggins Reference Beigie, Leonard and Wiggins1994), rigorous definitions have emerged for equivalent features that are relevant for finite-time stirring in flows of arbitrary time dependence also (Haller & Yuan Reference Haller and Yuan2000). Numerical detection of these LCS has often invoked the finite-time Lyapunov exponent (FTLE), defined as (Haller Reference Haller2002; Shadden, Lekien & Marsden Reference Shadden, Lekien and Marsden2005)

(3.28)\begin{equation} \sigma_{t_0}^{t_i}(\boldsymbol{x})=\frac{1}{|t_i|}\ln \sqrt{\lambda_{m}}, \end{equation}

where $\lambda _{m}$ is the maximum eigenvalue of the Cauchy–Green strain tensor,

(3.29)\begin{equation} \Delta=({\boldsymbol{\nabla}} \phi_{t_0}^{t_0+t_i}(\boldsymbol{x}))^{\rm T} {\boldsymbol{\nabla}} \phi_{t_0}^{t_0+t_i}(\boldsymbol{x}), \end{equation}

with the flow map $\phi _{t_0}^{t_0+t_i}$: $\boldsymbol {x}(t_0)\longmapsto \boldsymbol {x}(t_0+t_i)$ calculated from the advected positions of tracers from $t_0$ to $t_0+t_i$, where $t_i$ is the finite time interval of integration and superscript ${\rm T}$ indicates the transpose. Physically, FTLE is a measure of the maximum separation, after a time $t_i$, between the advected positions of infinitesimally close initial conditions around $\boldsymbol {x}$ at $t_0$. The direction associated with the maximally separating (over time $t_i$) initial conditions at $t_0$ is given by $\boldsymbol {\xi _m}$, the eigenvector corresponding to the eigenvalue $\lambda _m$.

Ridges of the backward ($t_i<0$) and forward ($t_i>0$) time FTLE fields have previously been computed as candidates to represent attracting and repelling LCS (together referred to as hyperbolic LCS), respectively, including in 2D (Mathur et al. Reference Mathur, Haller, Peacock, Ruppert, Jori and Swinney2007) and 3D (Green, Rowley & Haller Reference Green, Rowley and Haller2007) turbulent flows. However, since the maximum separation of initial conditions is not necessarily normal to the FTLE ridges, all FTLE ridges may not represent hyperbolic LCS (Haller Reference Haller2002). More recent approaches have identified hyperbolic LCS as maximally attracting or repelling strain lines, which, by definition, are normal to the direction $\boldsymbol {\xi _m}$ along which maximal separation of initial conditions occurs (Farazmand & Haller Reference Farazmand and Haller2012b). Attracting LCS, by definition, represent distinguished material lines along which fluid blobs get stretched. Patterns in dye release experiments therefore tend to resemble attracting LCS (Shadden, Dabiri & Marsden Reference Shadden, Dabiri and Marsden2006; Jimreeves et al. Reference Jimreeves, Mathur, Govardhan and Arakeri2018). In the turbulent ocean, even diffusive tracers such as the temperature have been shown to spatially organise, on sufficiently large scales, around the attracting LCS computed using optimally chosen integration times (Balasuriya, Ouellette & Rypina Reference Balasuriya, Ouellette and Rypina2018; Mathur et al. Reference Mathur, Jimreeves, Sharma and Agarwal2019). In RBC, Bandaru et al. (Reference Bandaru, Kolchinskaya, Padberg-Gehle and Schumacher2015) reported that, on the hot plate, thermal plumes seem to form along with attracting LCS in the skin friction field.

The choice of integration time $t_i$ for a given flow depends on the time scales over which one is interested in understanding advective stirring and the resulting tracer patterns. In addition, constraints from the flow field dataset also impose restrictions on the value of $t_i$. For example, (i) the velocity field (from experiments or simulations) of an unsteady flow may be available only for a finite period of time, or (ii) the velocity field may be available only over a finite spatial extent, in which case numerically advected tracer particles would leave the domain of measurements in finite time. In the current study, the velocity field measurements are in a horizontal plane, and tracer particles tend to leave this plane in a few $\Delta t$, where $\Delta t\sim 0.1$ s is the temporal resolution of our measurements (see table 1). We therefore seek LCS in the limit of $t_i\to 0$, i.e. the instantaneous limit of LCS. Physically, this limit would reveal structures that drive tracer patterns in short time scales (Mattia & Haller Reference Mattia and Haller2016).

Thermal plumes entrain fluid from their longitudinal edges, thus motivating us to define plumes as hyperbolic attracting LCS. To elucidate the attracting nature of plumes even at short time scales, the horizontal deformation of a fluid parcel in the vicinity of a plume, as evidenced in a direct numerical simulation at $(Ra, Pr)=(2\times 10^8, 1)$, captured using the Lagrangian particle tracking code of Sharma & Sameen (Reference Sharma and Sameen2019), is shown in figure 5. The fluid parcel, initiated close to the edge of a plume, is instantly pulled into the plume region and is also observed to get stretched in a way so as to align with the plume orientation in the horizontal plane. In summary, in a horizontal plane, plumes are sought as instantaneous limits of attracting LCS, and we proceed to show that regions of negative horizontal divergence are equivalent to regions in the vicinity of sufficiently strong attracting LCS in the instantaneous limit.

Figure 5.

Evolution of a fluid parcel in the vicinity of a plume, as evidenced in a numerical simulation, captured using the Lagrangian particle tracking code of Sharma & Sameen (Reference Sharma and Sameen2019), at ($Ra, Pr) = (2\times 10^8,1$). (a) Fluid particles located outside a plume at the initial time ${t^*}=0$, and (b) the advected fluid particle locations at a later time of ${t^*}=0.6$. In both (a,b), the background colour indicates the temperature distribution, with the red region being the plume central region.

For sufficiently small integration times $t_i$, considering the 2D motion within a horizontal plane, the Cauchy–Green strain tensor (3.29) can be written as $\Delta = I + 2t_iS$, where $I$ is the $2\times 2$ identity matrix and $S$ (3.18) is the Eulerian rate of strain tensor (Nolan et al. Reference Nolan, Mattia and Shane2020). As a consequence, the forward and backward time FTLE fields (3.28) are given by

(3.30a,b)\begin{equation} \lim_{t_i\to 0_+}\sigma_{t_0}^{t_i}=\lambda_{1},\quad \lim_{t_i\to 0_-}\sigma_{t_0}^{t_i}={-}\lambda_{2}, \end{equation}

where $\lambda _1$ and $\lambda _2$ are the eigenvalues of $S$, as given in (3.19)–(3.20). Therefore, ridges in the $-\lambda _2$ field are candidates to represent attracting LCS in the instantaneous limit. Furthermore, the eigenvectors corresponding to the forward and backward time FTLE are the same as the eigenvectors of $S$ corresponding to $\lambda _1$ and $\lambda _2$, respectively.

Seeking regions of negative $\lambda _2$ as regions containing ridges in the $-\lambda _2$ field, we investigate the spatial distribution of $-\lambda _2$ in figure 6. Figure 6(a) shows the spatial distribution of $-\lambda _2$ estimated from the experimental measurements at $(Ra,Pr) = (3.33\times 10^6,5.21)$, the same case as in figure 4(a). Superimposed on the $-\lambda _2$ field are the corresponding eigenvectors $\boldsymbol {\xi _2}$, which seem to be orthogonal to the ridges in most of the domain. This orthogonality indicates that ridges of the $-\lambda _2$ field are strain lines, and being local maximisers of the normal attraction rate (ridges are locally maximising features of the scalar field), they represent robust attracting LCS. When attracting LCS occur in the vicinity of repelling LCS, fluid parcels experience a combination of attracting and repelling behaviour, much like in the vicinity of saddle-type fixed points. Despite the presence of attracting LCS, repelling behaviour may dominate depending on the strength of the repelling LCS in the neighbourhood. Hence, to identify only the attraction-dominated regions as plume-like, we eliminate regions where $\lambda _1/(-\lambda _2)$ is larger than unity from the regions identified using $\lambda _2<0$, as shown later in figure 7(a). It is worth recalling from (3.30a,b) that $\lambda _1$ represents the forward time FTLE, and is hence a measure of the repelling behaviour in the region. As shown in figure 6(b), some regions of large and positive $\lambda _1$ indeed overlap with regions of positive $-\lambda _2$. For example, the region at around ($-15$, 25) mm where $\lambda _1$ is large and positive in figure 6(b) is not plume-like even though $-\lambda _2$ is also positive there. Moreover, the positive $\lambda _1$ region around ($-15$, 25) mm in figure 6(b) has the corresponding eigenvectors being orthogonal to the ridges in the $\lambda _1$ field, indicating the occurrence of strong repelling LCS in the region. Figure 7, as discussed later in this subsection, indeed confirms the occurrence of repelling LCS in this region.

Figure 6.

(a) Spatial distribution of $-\lambda _2$, overlaid with the corresponding eigenvectors $\boldsymbol {\xi _2}$ for ($Ra$, $Pr$) = (3.33 $\times 10^6$, 5.21); (b) corresponding distribution of $\lambda _1$, overlaid with the corresponding eigenvectors $\boldsymbol {\xi _1}$ in a horizontal plane at a height of 1.5 mm from the hot plate. Panels (c,d) are the same as (a,b), but for ($Ra$, $Pr$) = ($1.21\times 10^9$, 5.09). All the plots in this figure are based on experimental 2D3C PIV measurements.

Figure 7.

Attracting and repelling LCS in the instantaneous limit (integration time $t_i=0$) based on 2D3C velocity measurements in the experiment at ($Ra, Pr)=(3.33\times 10^6,5.21)$. (a) Attracting LCS (blue lines) overlaid over $\lambda _D$ field (regions of positive $\lambda _D$ are shown in white), (b) attracting (blue) and repelling (red) LCS plotted together.

In summary, in figure 6(a,b), regions of negative $\lambda _2$, with $-\lambda _2/\lambda _1<0$ or $-\lambda _2/\lambda _1>1$, indicate regions of attracting LCS with no repelling LCS that are stronger, and are hence plume-like. Based on the various scenarios shown in table 2, plumes can now be identified as regions where $\lambda _2$ is dominant and negative, which coincide with regions of $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}<0$. It is insightful to recall from the study of Mensa et al. (Reference Mensa, Özgökmen, Poje and Imberger2015) that horizontal divergence and tracer concentration fields were found to be highly correlated in a model ocean mixed layer with weak wind forcing. Our analysis provides a complementary perspective that the horizontal divergence field is highly correlated with regions of sufficiently strong attracting LCS. Attracting LCS, by their very definition, significantly influence tracer concentration fields, thus establishing that our discussion is consistent with the conclusions of Mensa et al. (Reference Mensa, Özgökmen, Poje and Imberger2015).

In figure 6(c,d), we reproduce figures 6(a) and 6(b) for a larger $Ra$ of $1.21\times 10^9$, which corresponds to a larger shear imposed on the regions close to the hot plate due to a stronger large-scale flow at the higher $Ra$. Regions of positive $-\lambda _2$ extract plume-like regions (figure 6c), but with the corresponding eigenvectors being somewhat non-orthogonal to the ridges of the $-\lambda _2$ field. We attribute this non-orthogonality to the presence of shear, which in some extreme scenarios can even introduce non-hyperbolic FTLE ridges (Haller Reference Haller2002). Regions inside the plumes, however, are known to have relatively low shear (Shevkar et al. Reference Shevkar, Gunasegarane, Mohanan and Puthenveettil2019; Blass et al. Reference Blass, Zhu, Verzicco, Lohse and Stevens2020), owing to which ridges in the $-\lambda _2$ field are likely to represent attracting LCS even at high $Ra$. In summary, attracting LCS in regions where repelling LCS are either absent or sufficiently weak are interpreted as the centrelines of plume regions. Finally, the $\lambda _1$ field in the regions of positive $-\lambda _2$ shown in figure 6(d) show that stronger repelling LCS can occur in the regions of attracting LCS, with such regions not considered plume-like.

We now extract variational theory-based attracting and repelling LCS (Farazmand & Haller Reference Farazmand and Haller2012a), based on the instantaneous limits of forward-time and backward-time FTLEs ($\lambda _1$ and $-\lambda _2$, respectively) and the corresponding eigenvectors, using the MATLAB-based software discussed in Onu, Huhn & Haller (Reference Onu, Huhn and Haller2015). Specifically, attracting LCS are obtained as trajectories of the differential equation

(3.31)\begin{equation} \frac{{\rm d}\boldsymbol{r}}{{\rm d}t}={\pm} \boldsymbol{\xi_1}(\boldsymbol{r}), \end{equation}

where $\boldsymbol {r}=(x,y)$ is the position vector and $\xi _1$ is the eigenvector corresponding to $\lambda _1$. These trajectories are initialised from local maxima of $-\lambda _2$, that lie in $\lambda _D<0$ regions, and are stopped based on the specification of a maximum length. Correspondingly, repelling LCS are obtained as trajectories of ${{\rm d}\boldsymbol {r}}/{{\rm d}t}=\pm \boldsymbol {\xi _2}(\boldsymbol {r})$, with initialisation points being from local maxima of $\lambda _1$ that lie in $\lambda _D>0$ regions. More details can be found in Onu et al. (Reference Onu, Huhn and Haller2015).

For the experimental measurements corresponding to figure 6(a), figure 7(a) shows the attracting LCS plotted on top of the $\lambda _D$ field, for which only the negative valued regions are shown in colour. While all the $\lambda _D<0$ regions have attracting LCS passing through them, attracting LCS also pass through regions where $\lambda _D>0$. Noticeable examples of such regions are marked as A, B, C, D, E and F in figure 7(a). Figure 7(b) shows that these regions have corresponding repelling LCS passing in their vicinity, with the repelling behaviour being dominant since $\lambda _D = \lambda _1>0$. It is also noteworthy that regions of intersection between attracting and repelling LCS are saddle-like and hence less likely to be plume regions.

Finally, to highlight the ability of $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}$ to pick up regions of negative $\lambda _D$, we plot the correlation between $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}$ and $\lambda _D$ from the entire spatial domain, averaged over 1399 snapshots, as a function of $Ra$ in figure 8. Values close to unity indicate that $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}$ can be a reliable indicator of regions with negative $\lambda _D$. Furthermore, the similarly large correlation between $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}$ and $\lambda _2$ suggests that regions with negative $\lambda _2$ and positive dominant $\lambda _1$ represent a relatively small fraction of the spatial domain at all $Ra$.

Figure 8.

The temporal mean values of the correlation coefficient between ◆ (blue), $\boldsymbol {\nabla }_H \boldsymbol {\cdot } \boldsymbol {u}$ values and the dominant eigenvalues and ■ (red), $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ and $\lambda _2$ values for the present range of $Ra$ in experiments. The correlation coefficient between arrays $P$ and $Q$ is defined as $r=\sum _{m} \sum _n (P_{mn}-\langle P\rangle _A)(Q_{mn}-\langle Q\rangle _A)/\sqrt {(\sum _{m} \sum _n (P_{mn}-\langle P\rangle _A)^2)(\sum _{m} \sum _n (Q_{mn}-\langle Q\rangle _A)^2)}$, where $\langle P\rangle _A$ and $\langle Q\rangle _A$ indicate spatial means.

3.4. Spatial distribution of eigenvalues

Based on (3.19)–(3.21), we define the dimensionless second eigenvalue as

(3.32)\begin{equation} \kappa=2\lambda_2/\gamma, \end{equation}

representing twice the ratio of the corresponding extensional or compressional strain rate – as the case may be based on whether $\lambda _2$ is positive or negative – with the total strain rate. Instantaneous snapshots of the spatial distribution of $\kappa$ are shown in figures 9(a) and 9(b) for $Ra=3.33\times 10^6$ and $Ra=1.2\times 10^9$, respectively. Regions with positive $\lambda _1$ and $\lambda _2$ (case 1 in table 2) are the boundary layer regions where strong local impingement occurs, and the fluid elements get stretched in both the horizontal principal directions; such regions are shown in maroon in figure 9. Regions with positive $\lambda _1$ and negative $\lambda _2$ with $\lambda _1$ being dominant (case 2 in table 2) are the boundary layer regions with little or no impingement, where the extensional strain dominates over the compressional strain; such regions are shown in red in figure 9. In the plume regions where the compressional strain dominates over the extensional strain, $\lambda _1$ is positive and $\lambda _2$ is negative with $\lambda _2$ being dominant (case 3 in table 2); such regions are shown in green/yellow in figure 9. Regions with negative $\lambda _1$ and $\lambda _2$ (case 4 in table 2) seem to occur at the centre of the plumes, such regions are shown in blue/cyan in figure 9. At low $Ra$ (figure 9a), a converging (in all directions) type of flow is observed in these regions, possibly due to the larger vertical velocity at the centre of the plumes; fluid elements would then get compressed in all directions in the horizontal plane due to the large vertical stretching by the rising plume. At higher $Ra$ (figure 9b) these plume centre regions with a converging flow are smaller in extent, possibly due to the stronger effect of shear from the large-scale flow. In table 2, cases 1 and 4 correspond to impingement portions of the boundary layer and central plume regions, respectively. Regions in the vicinity of boundary layer – plume boundaries are more faithfully captured by cases 2 and 3; case 2 being on the boundary layer side and case 3 being on the plume side.

Figure 9.

Spatial distribution of the dimensionless $\lambda _2$ (3.32), with the corresponding horizontal velocity vector field overlaid at (a) $Ra=3.33\times 10^6$ and (b) $Ra=1.21\times 10^9$.

The shear fields created due to LSF can be seen from the directions and magnitudes of horizontal velocity vectors in figures 9(a) and 9(b), as well as in figures 4(b) and 4(d), at $Ra=3.33\times 10^6$ and $Ra=1.21\times 10^9$, respectively. The scales of vector lengths are given in the caption of figure 4. In RBC, inside a box, large-scale flow is known to get oriented along one of the diagonals of a box and may change its orientation along the other diagonal with time. At lower $Ra$, due to the elliptic shape of LSF path (Niemela & Sreenivasan Reference Niemela and Sreenivasan2003), only a small near-plate area comes in direct contact with the LSF. As seen in figure 9(a), at lower $Ra$, the vectors are of smaller magnitudes and are randomly oriented, and so are the plumes; implying that the effect of shear due to the LSF on the near-plate region is minimal. Increase in $Ra$ increases the large-scale flow strength (3.25). Also, with increasing $Ra$, the shape of the LSF path gets changed from elliptic to increasingly square-like (Niemela & Sreenivasan Reference Niemela and Sreenivasan2003), which in-effect causes larger and larger near-wall areas to get affected by the LSF. At high $Ra$, due to this increased LSF strength, and the square-like nature of the LSF path along the diagonal of the box, a unidirectional shear field, with plumes oriented in the direction of shear, is created in the centre of a horizontal plane close to the hot plate (Shevkar et al. Reference Shevkar, Gunasegarane, Mohanan and Puthenveettil2019), as shown in figure 9(b). At any $Ra$, the effect of LSF is felt the most at the central region of the bottom plate, and the effect decreases as we move away from the centre towards the sidewalls.

For various eigenvalue combinations, we now categorise the flow field into four different regions based on the value of $\kappa$. It is evident from (3.20) that $\kappa =-1$ is satisfied when $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}=0$. As a result, regions with $\kappa \geq -1$ correspond to boundary layer regions, while regions with $\kappa <-1$ correspond to plume regions. As discussed earlier in the section, cases 1 and 2, separated by the $\lambda _2 = 0$ boundary, further divide the boundary layer regions into impingement ($\lambda _1>\lambda _2>0$) and non-impingement ($\lambda _1>0>\lambda _2$ with $\lambda _1$ being dominant) portions. In terms of $\kappa$, impingement and non-impingement portions of the boundary layer regions are given by $\kappa >0$ and $-1<\kappa <0$, respectively. Similarly, the plume regions given by $\kappa <-1$ are divided further as convergent-in-all-directions ($\lambda _2<\lambda _1<0$) and remaining ($\lambda _2<0<\lambda _1$ with $\lambda _2$ being dominant) portions. In terms of $\kappa$, $\kappa <-2$ and $-2<\kappa <-1$ correspond to convergent-in-all-directions plume and remaining plume regions, respectively.

The aforementioned classification based on the values of $\kappa$ is shown using different colours in figure 9. While the impingement and non-impingement boundary layer regions are shown in maroon and red, respectively, the convergent-in-all-directions and remaining plume regions are shown in blue/cyan and green/yellow, respectively. Based on similar arguments, it is also possible to perform the classification based on the values of non-dimensionalised $\lambda _1$, given by $\chi = 2\lambda _1/\gamma$. Specifically, the boundary layer regions correspond to $\chi >2$ (impingement) and $1<\chi <2$ (non-impingement), while the plume regions correspond to $\chi <0$ (convergent-in-all-directions) and $0<\chi <1$ (remaining plume regions).

In summary, the plume regions have $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u}<0$, which can be alternatively expressed as $\kappa <-1$ or $\chi <1$. Since these regions have $\lambda _D=\lambda _2$, which is negative, these regions have the dominant compression more than half the total strain, i.e. $|\lambda _2|>\gamma /2$. Within such plume regions, case 3 regions tend to occur on the edges of the plumes where shear strain, and hence the total strain, could be significant so that the dominant compressional strain $\lambda _2$ in these regions is in the range $\gamma /2<\lambda _2<\gamma$. These are the regions of the plumes that are affected by the shear due to the large-scale flow; we expect the extent of these regions to increase with $Ra$. At the convergent-in-all-directions portions of the plumes, where shear effects are not significant, we find case 4 regions where the dominant compressive strain $\lambda _2>\gamma$, the total strain.

The boundary layer regions have $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u}>0$, alternately expressed as $\chi >1$ or $\kappa >-1$. These regions have $\lambda _D=\lambda _1$, which is positive. Then, these regions have dominant extensional strain more than half the total strain, $\lambda _1>\gamma /2$. These boundary layer regions have within them case 1 regions, which occur in the impingement regions, and have dominant extensional strain $\lambda _1$ more than the total strain $\gamma$. These regions are surrounded by case 2 regions, which occur at the edge of these boundary layer regions, and have non-negligible shear strain, and hence total strain $\gamma$, so that the dominant extensional strain satisfies $\gamma >\lambda _1>\gamma /2$. These are the regions of the boundary layers that are affected by the external shear. Hence, it appears that the effects of external shear due to the large-scale flow are more strongly felt at the edges of the boundary layers and plumes, with the beginning of the boundary layer regions and the centre of plumes being relatively less affected by the large-scale flow.

The present criterion for detecting plumes, $\boldsymbol {\nabla _H} \boldsymbol {\cdot } \boldsymbol {u}<0$, which we showed to be equivalent to $\lambda _D<0$ in § 3.2 and equivalent to $\kappa =2\lambda _2/\gamma <-1$ in this subsection, picks up the plume regions along with their edges. In addition, the central plume regions are detected by local minima of $-\lambda _2$ or $\kappa$, where shear effects are minimal. Similarly, attracting LCS in the instantaneous limit occur in regions where shear is not dominant, and can hence be expected to pass through the middle of the $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}<0$ regions. While $\boldsymbol {\nabla _H}\boldsymbol {\cdot }\boldsymbol {u}=0$ identifies the plume boundaries, and hence allows us to estimate the widths of the plume structures, the LCS perspective only allows for us to extract lines that pass through the centre of the plume regions. Additional criteria are required for the LCS-based approach to detect the plume boundaries.

4.1. Verification with the theoretical estimate of length and area of plumes

The length of plumes in an area $A$, as given by Puthenveettil et al. (Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011), is

(4.1)\begin{equation} L_p=A/\bar{\lambda}, \end{equation}

where $\bar {\lambda }$ is given by (3.26). Substituting (3.26) in (4.1), we obtain

(4.2)\begin{equation} L_p=\frac{1}{2^{1/3}C_1Pr^{n_1}}\frac{A}{H}Ra^{1/3}, \end{equation}

thus giving $L_p$ as a function of $Ra$ and $Pr$ (Puthenveettil et al. Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011). Equation (4.2) has been verified with measurements of plume lengths – obtained by representing the plume regions in areas $A$ with short linear segments – over six decades of $Ra$, three decades of $Pr$ and a wide range of $A$ by Puthenveettil et al. (Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011), using planforms from experimental and numerical studies. We now compare the total length of negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions in an area $A$ with equation (4.2) to quantitatively verify whether the negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions are plumes. The total length of the negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions in various areas $A$, from our experimental and computational velocity fields, were measured by representing these regions with short linear segments, as done in Puthenveettil et al. (Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011). The sum of the length of these segments gives the total length of the negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions. Figure 10(a) shows the measured total length of the negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions, along with (4.2) for the range of $Ra$ and $Pr$ considered in the present study. The error bars show the range of data from measurements at different instants at the same $Ra$. The excellent agreement between (4.2) and the measured lengths of negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions indicates that these regions are indeed the plumes.

Figure 10.

(a) Variation of the measured total length of the negative $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}$ regions with $Ra$. The error bars show the ranges of values obtained from 10 different instants. Here $\blacktriangle$ (red) are experiments at the $Ra$ and $Pr$ shown in table 1. Computations at $Ra=2\times 10^7$ and $2\times 10^8$ for $\circ$ (blue), $Pr=1$; $\diamond$ (magenta), $Pr=10$; $\square$ (green), $Pr=4.96$; and $--$ (4.2). (b) Variation of the measured area fraction of the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ regions with $Ra$ with $n=0.09$. The error bars show one standard deviation of measurements from various instants. Symbols denote the same quantities as in (a), $--$ (4.9).

We now obtain a theoretical estimate of the area fraction $A_p/A$ of the plumes, where $A_p$ is the area of the plumes in an area $A$, so as to compare the measured values of the area fractions of the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ regions with this theoretical estimate; this serves as another verification of the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ criterion. Plumes are a result of instability of the natural convection boundary layers on either side of the plumes. Due to this, the mean thickness of the plumes in a horizontal plane ($t_{pc}$) at the height the boundary layers turn into plumes is known to scale as the boundary layer thickness (Van der Poel et al. Reference Van der Poel, Verzicco, Grossmann and Lohse2015). Then,

(4.3)\begin{equation} t_{pc}\sim 2 \delta_{vc}, \end{equation}

twice the critical velocity boundary layer thickness, $\delta _{vc}$. For any $Pr$, the total area of the plumes in an area $A$ is then

(4.4)\begin{equation} A_p\sim t_{pc} L_p, \end{equation}

where $L_p$ is given by (4.2). From (4.3) and (4.4),

(4.5)\begin{equation} \frac{A_p}{A} =C(Pr)2\delta_{vc}\frac{L_p}{A}, \end{equation}

where $C(Pr)$ is an unknown function of $Pr$ to take into account the effects of $Pr$ on the instability of the boundary layers. Using the expressions for the critical thermal boundary layer thickness $\delta _c\sim 4\bar {\lambda }/(Gr_\lambda Pr)^{1/5}$ and the ratio $\delta _c/\delta _{vc}\sim 1/\sqrt {Pr}$, given as (A11) and (A17), respectively, in Puthenveettil et al. (Reference Puthenveettil, Gunasegarane, Agrawal, Schmeling, Bosbach and Arakeri2011),

(4.6)\begin{equation} \delta_{vc}\sim \frac{4\bar{\lambda} Pr^{2/5}}{Gr_\lambda^{1/5}}. \end{equation}

Here, $Gr_\lambda =g\beta \Delta T_w \overline {\lambda }^3/\nu ^2$ is the Grashoff number based on the mean plume spacing $\bar {\lambda }$. Substituting (3.26) for $\bar {\lambda }$ into (4.6), we obtain

(4.7)\begin{equation} \delta_{vc}\sim \left(\frac{C_1}{2}\right)^{2/5}Pr^{0.54}Z_w, \end{equation}

where $Z_w$ is the near wall length scale given by (3.27). Substituting (4.7) and (4.2) in (4.5), we obtain the theoretical estimate of the area fraction of plumes as

(4.8)\begin{equation} \frac{A_p}{A}=0.15\,C(Pr)\,Pr^{0.44}. \end{equation}

Interestingly, (4.8) – an important new result of the present work – shows that the area fraction of plumes near the hot surface in turbulent RBC is independent of $Ra$, and is only a function of $Pr$. Even though the total length of plumes ($L_p$) increases as $Ra^{1/3}$ by (4.2), their corresponding mean thickness ($2\delta _{vc}$) decreases as $Ra^{-1/3}$ (see (4.7) and (3.27)) so that the area fraction of plumes $A_p/A$ becomes independent of $Ra$.

Figure 10(b) shows the measured areas of the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ regions, normalised by the total area, from the experimental and the computational velocity fields for $1\le Pr\le 10$ and $5\times 10^5\le Ra \le 1 \times 10^9$. The measured area fractions shown in the figure are averaged from those obtained at various instants at a given ($Ra$, $Pr$), with the error bars showing one standard deviation range of these instantaneous values. The measured mean area fractions of the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ regions collapse on to the relation

(4.9)\begin{equation} \frac{A_p}{A}= 0.55\,Pr^{{-}0.09}, \end{equation}

shown by the dashed line in the figure. The measured mean area fractions are also independent of $Ra$, in agreement with the theoretical relation (4.8), thereby validating the $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ criterion. Comparison of (4.9) with (4.8) implies that $C(Pr)=3.67Pr^{-0.53}$ captures the $Pr$ effects on the instability of the boundary layer.

4.2. Comparison with other criteria for detecting plumes

In this subsection, we compare the plumes extracted by the ${\nabla ^{*}}_H\boldsymbol {\cdot }\boldsymbol {u^{*}}<0$ criterion, when applied to the velocity fields obtained from our numerical computations, with those extracted by the criterion in (1.1) (Gastine et al. Reference Gastine, Wicht and Aurnou2015) and the criterion

(4.10)\begin{equation} w^*T^{*'} > 0, w^{*}>0. \end{equation}

Here, in addition to $w^{*}T^{*'}>0$ mentioned by Schumacher (Reference Schumacher2009), we include $w^{*}>0$ since we are only interested in rising plumes close to the bottom plate. The criterion (1.1), which we use with $c=0.25$, is a temperature based criterion while the criterion (4.10) is a flux-based criterion. The plumes extracted by the three criteria (3.17), (1.1) and (4.10) in a horizontal plane at a height within the natural convection boundary layer, are compared for $Pr=1,\, 4.98$ and $10$ at $Ra=2\times 10^7$. The plumes from these criteria are assessed vis-à-vis the dimensionless positive vertical velocity ($w^{*}>0$) field and the dimensionless temperature field $T^*$ in the same plane. A theoretical comparison of the present $\boldsymbol {\nabla }_H\boldsymbol {\cdot }\boldsymbol {u}<0$ criterion with the $w>0$ condition using similarity solutions described in § 3.1 is given in Appendix B; a similar comparison with the flux criterion, using 2D boundary layer equations, is given in Appendix C. We also show in Appendix D that application of methods for identification of coherent structures in shear turbulence (Hunt et al. Reference Hunt, Wray and Moin1988; Jeong & Hussain Reference Jeong and Hussain1995; Zhou et al. Reference Zhou, Adrian, Balachandar and Kendall1999) fail to detect the complete plume structure.

Figure 11(c–e) show the plumes detected by the above three criteria, along with the dimensionless temperature field in figure 11(a) and the $w^{*}>0$ regions in figure 11(b), all at $z/H=0.012$ for $Ra=2\times 10^7$ and $Pr=1$. Overall, all the three criteria detect similar regions, which have higher temperatures and $w^{*}>0$. However, the horizontal divergence criterion (3.17) detects high temperature upwellings such as the connecting plumes at D and E (see figure 11d), isolated initiation of plumes at F and upwellings in-between merging plumes at G (see figure 11e); these commonly occurring regions are all excluded by both the criteria (1.1) and (4.10). In addition, there are a few regions, such as A and B in figure 11(b) and C in figure 11(c), which have temperatures quite higher than the mean temperature, but have $w^{*}<0$. Such regions are seen to occur in-between two merging plumes. A purely temperature-based criterion, such as (1.1), erroneously identifies these regions as plumes, as seen in figure 11(d). The kinematic condition $\boldsymbol {\nabla }^{*}_H\boldsymbol {\cdot }\boldsymbol {u^{*}}<0$ excludes such regions since these are not regions where there is an upward acceleration. A flux-based criterion, such as (4.10), also excludes such regions (see figure 11e) since the flux is downward due to $w^{*}<0$ in such regions. Further, (1.1) and (4.10), in general, select narrower plume regions, missing out the upwellings which occur at the edges of upwelling regions that have higher than mean temperatures, as seen in the specific case of the plume at H in figure 11(e).

Figure 11.

Plumes detected by various criteria at $Ra=2\times 10^7$ at three $Pr$: (a–e), $Pr=1$, $z/H= 0.012$; (f–j), $Pr=4.96$, $z/H=0.029$; (k–o), $Pr=10$, $z/H=0.043$. (a,f,k) The temperature field; (b,g,l), the positive vertical velocity field; (c,h,m), $\boldsymbol {\nabla ^{*}_H} \boldsymbol {\cdot } \boldsymbol {u^{*}}<0$; (d,i,n), the criterion by Gastine et al. (Reference Gastine, Wicht and Aurnou2015) (1.1); (e,j,o), the criterion by Schumacher (Reference Schumacher2009) with $w^{*}>0$ (4.10).

In figure 12, we compare the evolution of the instantaneous area fractions of plumes ($A_{pi}/A$) with time in a horizontal plane at $z/H=0.012$ for $(Ra, Pr)=(2\times 10^7,1)$ using the above three criteria (3.17), (1.1) and (4.10), along with the $w^{*}>0$ regions and the $T^{*}>0.75$ regions. A theoretical estimate for the corresponding time averaged area fraction of plumes ($A_p/A$) was obtained in (4.8). In the present calculations of $A_{pi}/A$, a circular area within $r^*=0.352$ is taken as $A$ to avoid the influence of the sidewall. The overall trend of $A_{pi}/A$ with time is the same for all these criteria, even though there are differences in the magnitudes and the details at specific time instants, such as at $t^{*}=7, 28, 51$ and $61$. The area fractions estimated from $\boldsymbol {\nabla }^{*}_H\boldsymbol {\cdot }\boldsymbol {u^{*}}<0$ follow the trend of the $w^{*}>0$ regions quite closely, but is always slightly lower in magnitude. This slightly lower area of plumes estimated from (3.17), compared with that from $w^{*}>0$, can also be seen in figure 11(c,b), especially in the larger plume area at F and E in figure 11(b), compared with the corresponding regions in figure 11(c). The edges of the $w^{*}>0$ regions, which do not have buoyancy driven vertical acceleration, are, possibly correctly, not considered as plumes by the horizontal divergence criterion. Further, as shown in (B4a–d), (3.17) implies that, in addition to $w=0$, the edge of the plume must also have $\partial w/\partial z =0, \partial w/\partial x =0$ and $\partial u/\partial x=0$; these more stringent conditions seem to shift the plume boundary calculated from (3.17) slightly inward into the upwelling regions, so as to result in slightly lower plume areas compared with that with a $w>0$ condition. It must also be remembered that a $w>0$ criterion would wrongly identify all the upward flow regions in the up-flow region of the large-scale flow as plumes.

Figure 12.

Time series of instantaneous plume area fractions detected by $T^{*}\ge 0.75$ (in black), $w^{*}>0$ (in red), $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot } \boldsymbol {u^{*}} >0$ (in green), (1.1) (in magenta) and  (4.10) (in blue) at $Ra=2\times 10^7$ and $Pr=1$. The inset shows the $T^{*}>0.75$ regions in green, yellow and red colour at $t^{*}=51$, with the $w^{*}=0$ boundary shown by the dashed line, while the regions where $T^{*}<0.75$ and $w^{*}<0$ are shown in light blue.

The criteria (1.1) and (4.10) closely follow the trend of the area evolution of the $T^{*}\ge 0.75$ regions, albeit with noticeably smaller areas at all instants. Similarly, the areas captured by (1.1) and (4.10) are much lower than those from (3.17) and $w^{*}>0$. These three criteria, viz. $T^{*}\ge 0.75$, (1.1) and (4.10), however, show substantial deviations from the velocity based criteria, viz. (3.17) and $w^{*}>0$, at some specific time instants, such as at $t^{*}=7, 28, 51$ and $61$. The plume areas estimated by $T^{*}\ge 0.75$ at these times show a sudden rise to values approximately 40 % more than the mean area fractions; similar trends are shown by the (1.1) and the (4.10) criteria, even though their mean area fractions are lower. We examine one such instant at $t^{*}=51$ in the inset in figure 12. The inset figure shows that there are large regions where $T^{*}> 0.75$, which, however, have $w^{*}$ slightly less than zero, possibly due to reasons like the hotter rising fluid bouncing back down after hitting the top plate, or due to the stronger downward entrainment of hotter fluid in-between two nearby rising strong plumes. The spike in $A_{pi}/A$ shown by the $T^{*}\ge 0.75$ condition and by the criteria (1.1) and (4.10) are due to such regions being erroneously classified as plumes.

At a higher Prandtl number, $Pr=4.96$, the structure is less intricate in figure 11(g) compared with that in figure 11(b), possibly due to the increase in viscous effects. Comparing figures 11(g) and 11(b), we see that the upwellings at $Pr=4.96$ are mostly of the same thickness as at $Pr=1$. The plumes, however, have thinner higher temperature regions, as could be noticed by comparing figures 11(f) and 11(a). A temperature criterion such as $T^{*}>0.75$ at this $Pr$ would then detect much smaller area of plumes, as seen in the temporal evolution of plume areas at $Pr=4.96$ in figure 13(a). The $T^{*}\ge 0.75$ criterion, which gave $A_{pi}/A\approx 0.5$ for $Pr=1$ in figure 12, now detects a much smaller area of plumes at $Pr=4.96$, with $A_{pi}/A\approx 0.1$ in figure 13(a). The fraction of areas of velocity upwellings at $Pr=4.96$, however, remains at approximately $0.45$ in figure 13(a), around the same value as at $Pr=1$ in figure 12.

Figure 13.

(a) The time series of instantaneous plume area fractions at $Ra=2\times 10^7$ and $Pr=4.96$ detected by $T^{*}\ge 0.75$ (in black), $w^{*}>0$ (in red), $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot } \boldsymbol {u^{*}} >0$ (in green), (1.1) (in magenta) and  (4.10) (in blue). (b) Correlation coefficients between $T^{*}$ and $\boldsymbol {\nabla ^{*}_H} \boldsymbol {\cdot } \boldsymbol {u^{*}}$ (in black), $T^{*}$ and $w^{*}$ (in red) and $w^{*}$ and $\boldsymbol {\nabla ^{*}_H} \boldsymbol {\cdot } \boldsymbol {u^{*}}$ (in green). Circles denote $Ra=2\times 10^7$ and triangles denote $Ra=2\times 10^8$. The error bars are plotted using the ranges of values obtained in five instants.

Figure 11(h) shows that the criterion (3.17) captures all these upwellings at this higher $Pr$ with their actual thicknesses. The $w^{*}>0$ areas and the $\boldsymbol {\nabla }^{*}_H\boldsymbol {\cdot }\boldsymbol {u^{*}}<0$ areas coincide at all instants at this higher $Pr$, as seen in figure 13(a). The divergence criterion, however, seems to slightly underestimate the downward flow, non-plume regions in-between large upwellings, as seen at I in figure 11(g); such regions are, however, only a small fraction of the plume area. On the contrary, the criteria (1.1) and (4.10), capture approximately $40\,\%$ fewer upwelling areas, as seen in figure 13(a); the reason being the much thinner regions characterised as plumes by these criteria at $Pr=4.96$ (see figure 11i,j vis-à-vis figure 11g). These two criteria seem to be able to capture only the thermal boundary layer regions of the plumes, which are thinner compared with the velocity boundary layers at $Pr>1$. The two criteria (1.1) and (4.10) miss large upwell regions such as an entire plume at L in figure 11(i). The two criteria also miss out upwellings on the outer peripheries of higher temperature regions, such regions such as J and K in figure 11(h).

A further increase in $Pr$ reduces the intricacy of the plume patterns further; see the upwelling patterns shown in figure 11(l) at $Pr=10$ and those in figure 11(g,b) at $Pr=4.96$ and $Pr=1$, respectively. The thickness of the upwellings in these figures remains approximately the same, while the high temperature regions in figure 11(k) become even thinner than in figure 11(f). The $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot } \boldsymbol {u^{*}}<0$ criterion identifies all the upwelling regions as plumes, while the criteria (1.1) and (4.10) mostly identify only the thinner, higher temperature regions as plumes, often missing entire large upwelling regions, such as the plume at M in figure 11(l). However, the divergence criterion (3.17) seems to slightly overestimate the plume regions with downwelling regions near the edges of the upwellings, especially when they are closed polygons, such as at N and O in figure 11(m), also being identified as plumes; these are, however, quite small areas compared with the overall plume areas.

This increasing difference between the areas of upwelling regions and the areas of higher temperature regions with an increase in $Pr$ could be seen in figure 13(b), which shows the correlation coefficient between $T^{*}$ and $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot } \boldsymbol {u^{*}}<0$ for various $Ra$ and $Pr$. As discussed earlier, with increase in $Pr$, the higher temperature region areas, compared with the upwelling region areas, shrink in their thickness (see figures 11f, 11g vis-à-vis figures 11a, 11b, and figures 11k, 11l vis-à-vis figures 11f, 11g). This occurs since the thermal boundary layers become thinner compared with the viscous boundary layer with increase in $Pr$. Correspondingly, the correlation between $T^{*}$ and $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot }\boldsymbol {u^{*}}$ in figure 13(b) decreases from around $0.85$ to $0.7$ with increase in $Pr$ from $1$ to $10$, showing that temperature-based criteria will detect increasingly smaller areas as plumes, compared with the upwelling regions, with increasing $Pr$. Due to the same reason, a decrease in correlation between $T^{*}$ and $w^{*}$ with increase in $Pr$ is also seen in figure 13(b). In contrast, the correlation between $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot }\boldsymbol {u^{*}}$ and $w^{*}$ is high and varies little (between $1$ and $0.95$) with $Pr$, showing the capability of the $\boldsymbol {\nabla ^{*}_H}\boldsymbol {\cdot }\boldsymbol {u^{*}}<0$ criterion to pick up the upwellings as plumes at all $Pr\ge 1$.