2. Mathematical formulation

This study considers a Newtonian fluid under the Boussinesq approximation, with density $\rho$ , kinematic viscosity $\nu$ and thermal diffusivity $\kappa$ . The fluid is bounded between two horizontal no-slip plates with a gap of width $d$ , where the top plate is usually considered as an ice–water interface and the bottom is perfectly insulating. By considering a monochromatic light spectrum, the system has an internal heat source $Q(z)$ ( $^{\circ} {\textrm C}$ s⁻¹), as sketched in figure 1. It is a simplified model of a short-wave radiation that is characterised by a single decay scale (Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019; Creyssels Reference Creyssels2020). According to the Beer–Lambert law, this type of heating is associated with a volumetric heat source that decays exponentially (Barnes & Hobbie Reference Barnes and Hobbie1960; Farmer Reference Farmer1975; Grace et al. Reference Grace, Stastna, Lamb and Scott2022, Reference Grace, Fogal and Stastna2023):

(2.1) \begin{equation} Q(z) = \frac {Q_0}{l}\exp \left (\frac {z-d}{l}\right ), \end{equation}

where $l$ is the attenuation length, which normally depends on radiation frequency and water clarity, and $Q_0$ is the magnitude of the thermal forcing at the top of the domain. Polychromatic light has a different effect: the heat flux is the sum of different wavelength bands that have different energy and attenuation depths (Mironov et al. Reference Mironov, Terzhevik, Kirillin, Jonas, Malm and Farmer2002; Bouffard et al. Reference Bouffard2019; Smirnov et al. Reference Smirnov, Smirnovsky, Zdorovennova, Zdorovennov, Efremova, Palshin and Bogdanov2023).

Figure 1.

Internally heated convection model with an internal heating source $Q$ , as expressed in equation (2.1). The decay scale of $Q$ is $l$ , and the $e$ -folding scale is located at a height $d - l$ . The model has an isothermal top layer and an insulating bottom layer.

Two additional models were also employed to investigate convection in lakes. When light penetrates the water column and reflects off the bottom, a reverse exponential decay happens from the bottom to the top. The intensity of the reflected light $I(z)$ at a depth $z$ (measured from the surface) can be modelled as similar to equation (2.1) but with a bottom heat flux $Q_{ bottom}$ . This model can well represent shallow aquatic systems or environments where the bottom strongly reflects light (sandy or icy bottoms) (Lepot et al. Reference Lepot, Aumaître and Gallet2018; Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019). If we consider that the bottom of the lake can absorb heat and re-inject it into the fluid, then we can regard it as a simple penetrative convection problem of bottom heating. In this case, there is a constant temperature difference between the upper and lower boundaries of the convection system (Veronis Reference Veronis1963).

The dimensionless governing equations can be obtained by introducing the layer height $d$ as the length scale and $d^2/\kappa$ as the thermal time scale. The temperature scale is $\Delta T = T_{M}-T_0$ , where $T_0={0}\,^\circ \textrm {C}$ is the reference temperature and $T_M={4}\,^\circ \textrm {C}$ . Thus $\theta = (T-T_0)/\Delta T$ is the dimensionless temperature and $\theta _{M} =(T_{M}-T_0)/\Delta T$ describes the effect under the maximal density with $\theta _{M} =1$ . Then the dimensionless governing equations for two-dimensional penetrative convection are as follows:

(2.2) \begin{equation} \nabla \cdot \boldsymbol{u}=0, \end{equation}

(2.3) \begin{equation} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\cdot \nabla \boldsymbol{u}=-\nabla p+Pr Ra(\theta -1)^2\boldsymbol{e}_z+Pr\nabla ^2\boldsymbol{u}, \end{equation}

(2.4) \begin{equation} \frac {\partial \theta }{\partial t}+\boldsymbol{u}\cdot \nabla \theta =\nabla ^2\theta + S(z), \end{equation}

where $ Ra := g\alpha (\Delta T)^2d^3/\nu \kappa$ is the Rayleigh number and $Pr:= \nu /\kappa$ is the Prandtl number, and

(2.5) \begin{equation} S(z) = \frac {M}{\delta }\exp \left (\frac {z-1}{\delta }\right ). \end{equation}

Here, $\delta :=l/d$ and $M := dQ_0/(\Delta T\kappa )$ are dimensionless parameters that quantify attenuation length and incoming thermal forcing, respectively. Once radiation penetrates through the ice, the value of $Q_0$ is about $7.25\times 10^{-6}$ $^\circ \textrm {C}$ $\cdot$ m s⁻¹ (Kirillin et al. 2012; Ulloa et al. Reference Ulloa, Winters, Wüest and Bouffard2019; Grace et al. Reference Grace, Stastna, Lamb and Scott2022). Assuming that $d=1$ m and taking freshwater as an example ( $\kappa =10^{-7}\,{\textrm m}^2\, {\textrm s}{^-}{^1}$ ), we estimate $M=18.125$ . We chose the range of $\delta$ from $0.1$ to $0.8$ to examine the impacts of both smaller and larger attenuation lengths. By integrating the heat flux $S(z)$ , the total heat input to the system is

(2.6) \begin{equation} E = \int _0^1\frac {M}{\delta }\exp \left (\frac {z-1}{\delta }\right ){\textrm d}z = M\left (1-\exp \left (-\frac {1}{\delta }\right )\right ). \end{equation}

Clearly, the total energy input increases when $M$ is larger. For given $M$ , the total energy input decreases as $\delta$ increases.

The boundary conditions for velocity and temperature are

(2.7) \begin{equation} \boldsymbol{u}|_{z=0,1} = 0, \end{equation}

(2.8) \begin{equation} \partial _z\theta |_{z=0}=0,\,\,\theta |_{z=1}=0. \end{equation}

To investigate the heat transport efficiency of the system, we start from the definition of heat flux:

(2.9) \begin{equation} q := \frac {\lambda (T_{bottom}-T_{top})}{d}, \end{equation}

where $\lambda$ is the thermal conductivity. Applying the boundary conditions, we have a fixed temperature at the top layer, $T_{top}= {0}\, ^\circ \textrm {C}$ . In our configuration, the Nusselt number is defined as $q_{cond}/q_{conv}$ ( $q_{cond}$ and $q_{conv}$ are the conductive and convective heat fluxes, respectively) to reflect the influence of the intensity of convection on the heat transport, and calculated as (Goluskin Reference Goluskin2016; Kazemi et al. Reference Kazemi, Ostilla-Mónico and Goluskin2022)

(2.10) \begin{equation} Nu = \frac {\theta _{st}(0)}{\bar {\theta }(0)}, \end{equation}

where $\theta _{st}$ is the dimensionless base state solution to the temperature equation with $\boldsymbol{u}=\boldsymbol{0}$ . The detailed expression can be found in equation (3.2). The overbar denotes an average over the horizontal surfaces and long enough time with

(2.11) \begin{equation} \bar {(\cdot )}= \frac {1}{\tau L}\int _0^{\tau }{\textrm d}t\int _0^L\,(\cdot )\,{\textrm d}x, \end{equation}

where $\tau$ can be arbitrary for steady solutions, and for turbulent flows it is chosen as four thermal times in simulations when the flow is fully developed.

3. Base state and instability

Firstly, we consider a static base state:

(3.1) \begin{equation} \boldsymbol{u} = \textbf {0},\quad \frac {\partial ^2\theta _{st}}{\partial z^2}+S(z)=0. \end{equation}

Thus, static state temperature $\theta _{st}$ is

(3.2) \begin{equation} \begin{aligned} \theta _{st} & = -M\delta \exp \left (\frac {z-1}{\delta }\right )+C_1z+C_2,\\[6pt] C_1 := & M\exp \left (-\frac {1}{\delta }\right ),\quad C_2 := M\delta -M\exp \left (-\frac {1}{\delta }\right ). \end{aligned} \end{equation}

The water body beneath the ice is heated via radiation which is the only form of energy input of the system. Figure 2(a) represents the base state of the temperature for varied attenuation depths and the corresponding density profiles are displayed in figure 2(b). As the density is quadratically dependent on temperature, the water body is stratified into a stably stratified layer (SSL), whose temperature ranges from 0 to $\theta _M$ , and an unstably stratified layer (USL). When most radiation is absorbed near the ice–water interface, characterised by a small attenuation depth $\delta$ , fluid deep in the water body cannot be efficiently heated and maintains a uniform temperature, forming a neutrally stratified layer (NSL). For a system with larger $\delta$ , the depth of the NSL becomes smaller and the USL is thicker. However, when radiation penetrates even deeper, the magnitude of radiation absorbed by the water body, characterised by $M/\delta$ , is lower. Therefore, the system reaches a maximum bottom temperature $\theta _{st}(0) \approx 5.40$ and a minimum density $\rho _{st}(0)$ at $\delta \approx 0.56$ when ${\textrm d}\theta _{st}(0)/{\textrm d}\delta = 0$ .

Figure 2.

(a) The static state temperature profiles $\theta _{st}$ with respect to the height $z$ for varied attenuation ratios $\delta$ . (b) The corresponding static state density profiles $\rho _{st}$ for different $\delta$ . The horizontal dashed line is the dividing line between the SSL and the USL.

We add small perturbations to the base state, and then the velocity field and temperature field become

(3.3) \begin{equation} \boldsymbol{u} = \boldsymbol{0}+\tilde {\boldsymbol{u}},\quad \theta = \theta _{st}+\tilde {\theta }. \end{equation}

For two-dimensional incompressible flow, we can introduce a streamfunction $\psi$ defined as $\tilde {\boldsymbol{u}} = \nabla \times (\psi \boldsymbol{e}_y)$ such that

(3.4) \begin{equation} \tilde {u} = -\partial _z\psi, \quad \tilde {w}=\partial _x\psi . \end{equation}

Therefore, the linearised equations for small perturbations are given by

(3.5) \begin{equation} \frac {\partial }{\partial t}(\nabla ^2\psi ) = 2Pr Ra(\theta _{st}-1)\frac {\partial \tilde {\theta }}{\partial x}+Pr\nabla ^2(\nabla ^2\psi ), \end{equation}

(3.6) \begin{equation} \frac {\partial \tilde {\theta }}{\partial t}+\frac {\partial \psi }{\partial x}\frac {\partial \theta _{st}}{\partial z}=\nabla ^2\tilde {\theta }, \end{equation}

and the corresponding boundary conditions become

(3.7) \begin{equation} \begin{aligned} & \psi |_{z=0,1} = \partial _z\psi |_{z=0,1} = 0,\\[6pt] & \partial _z\tilde {\theta }|_{z=0}= 0, \,\,\tilde {\theta }|_{z=1}=0. \end{aligned} \end{equation}

In a standard way, we conduct the normal mode analysis:

(3.8) \begin{equation} \tilde {\psi } = \hat {\psi }(z){\textrm e}^{{ i}kx+\sigma t}, \quad \tilde {\theta } = \hat {\theta }(z){\textrm e}^{{ i}kx+\sigma t}, \end{equation}

where $k$ is the wavenumber and $\sigma$ is the growth rate. Then the governing equation is transformed into an eigenvalue problem of $\sigma$ .

Figure 3.

(a) Marginal stability curve in the wavenumber ( $k$ )–Rayleigh number ( $Ra$ ) plane under different attenuation ratios $\delta$ from $0.1$ to $0.8$ , with a fixed thermal forcing $M =18.1250$ . (b) The surface plot shows how the critical Rayleigh number $Ra_c$ is affected by $M$ and $\delta$ , for $15\leqslant M\leqslant 30$ and $0.1\leqslant \delta \leqslant 1$ .

The relationship between the critical Rayleigh numbers $Ra_c$ and wavenumbers $k$ is shown in figure 3(a). The dashed line directly indicates that the relationship between $Ra_c$ and critical wavenumber $k_c$ is not monotonic. Within our tested range, $Ra_c$ has a minimum point around $\delta = 0.6$ (as shown in the inset) with $Ra_c \sim O(10)$ , which is two decades smaller than that in classical RBC: $Ra_c \sim O(10^3)$ (Chandrasekhar Reference Chandrasekhar1961; Cross & Hohenberg Reference Cross and Hohenberg1993). Intuitively, we expect to observe a more stable flow if the input energy is reduced. However, when $\delta$ increases up to 0.6, the thickness of the USL increases even though the energy input decreases. Thus, the system can be more unstable and the critical Rayleigh number becomes smaller. As $\delta$ becomes larger ( $\delta \gt 0.6$ ), further reduction in the input energy causes a thinner unstable layer, leading to a more stable system. Thus, we observe that the critical Rayleigh number increases when $\delta \gt 0.6$ . Furthermore, we took a transect through the plane presented in figure 3(b). Results are shown in figure 16 in Appendix A. It is clear that the system would be more unstable as the total input energy becomes higher.

The surface plot presented in figure 3(b) provides a detailed view of how $Ra_c$ is influenced by the two dimensionless parameters $\delta$ and $M$ . We chose the ranges of $15\leqslant M\leqslant 30$ and $0.1\leqslant \delta \leqslant 1$ to examine the influence of incoming heat thermal heat flux and forcing length scales on the linear and nonlinear dynamics. If the attenuation depth is fixed, an increase in the incoming heat flux ( $M$ ) would lead to a higher temperature difference within the USL. This suggests that the system retains more heat, resulting in a lower density in the unstable stratified layer, eventually leading to a lower $Ra_c$ and making the system unstable easily.

4. Weakly nonlinear stability

To analyse the convective flow, it is crucial to incorporate the nonlinear terms into the system. Therefore, we conduct a weakly nonlinear analysis of the system, focusing on the dynamics near $Ra_c$ . We aim to identify the specific type of local bifurcation it undergoes. To investigate weakly nonlinear solutions, we use multi-scale expansion and determine solvability conditions.

The governing equations (2.3)–(2.4) are restated by the streamfunction and temperature:

(4.1) \begin{equation} \frac {\partial \nabla ^2\psi }{\partial t}+J(\psi, \nabla ^2\psi ) = 2RaPr(\theta _{st}-1+\tilde {\theta })\frac {\partial \tilde {\theta }}{\partial x}+Pr\nabla ^4\psi, \end{equation}

(4.2) \begin{equation} \frac {\partial \tilde {\theta }}{\partial t}+J(\psi, \tilde {\theta }) = -\frac {\partial \psi }{\partial x}\frac {\partial \theta _{st}}{\partial z}+\nabla ^2\tilde {\theta }, \end{equation}

where $J(\alpha, \beta ) = ({\partial \alpha }/{\partial x}) ({\partial \beta }/{\partial z})- ({\partial \alpha }/{\partial z}) ({\partial \beta }/{\partial x})$ and the boundary conditions are with reference to equation (3.7).

When $Ra$ is slightly larger than $Ra_c$ , there is an unstable mode that becomes dominant through linear theory, whereas other modes decay exponentially in time. However, the unstable mode in linear theory cannot reveal the full dynamics of the problem as the nonlinear term would be significant when the perturbation grows. We invoke the multi-scale analysis and define a small parameter $\epsilon$ :

(4.3) \begin{equation} \epsilon ^2r = Ra-Ra_c,\quad |Ra-Ra_c|\ll 1, \end{equation}

where $r$ is a rescaled parameter of $r\sim O(1)$ and the slow time scale $\tau$ is defined as

(4.4) \begin{equation} \tau = \epsilon ^2t. \end{equation}

Thus, this leads us to the assumption that the weakly nonlinear solution is of the form

(4.5) \begin{equation} \psi (x, z, \tau ) := \epsilon ^\alpha \hat {\psi }_1(z){\textrm e}^{{ i}kx}A(\tau ) + \epsilon ^{2\alpha }\left [\hat {\psi }_{20}(z)+\hat {\psi }_{22}(z){\textrm e}^{{ i}2kx}\right ]A^2(\tau )+\cdots +\text{c.c.} \end{equation}

(4.6) \begin{equation} \tilde {\theta }(x, z, \tau ) := \epsilon ^\alpha \hat {\theta }_1(z){\textrm e}^{{ i}kx}A(\tau ) + \epsilon ^{2\alpha }\left [\hat {\theta }_{20}(z)+\hat {\theta }_{22}(z){\textrm e}^{{ i}2kx}\right ]A^2(\tau )+\cdots + \text{c.c.} \end{equation}

The system of equations possesses a unique set of solutions if and only if it satisfies a solvability condition, which yields the following Landau-type equation:

(4.7) \begin{equation} \frac {{\textrm d}A}{{\textrm d}\tau } = \lambda _1rA+\lambda _2|A|^2A, \end{equation}

where $\lambda _1$ and $\lambda _2$ are constants (see Appendix B). The Landau equation (4.7) defines the amplitude function and is used to determine whether the stability is supercritical or subcritical. The equilibrium amplitude can be solved as

(4.8) \begin{equation} A_{eq} = \sqrt {-\frac {\operatorname {Re}(\lambda _1)r}{\operatorname {Re}(\lambda _2)}}. \end{equation}

Figure 4.

Comparison of the weakly nonlinear analysis with a full branch continuation prediction at the neighbourhood of the critical Rayleigh number $Ra_c$ for $\delta = 0.6$ and different $M$ . Values of $\sqrt {\langle \boldsymbol{u},\boldsymbol{u}\rangle }$ are taken over the whole domain $\varOmega = [-\pi /k,\pi /k]\times [0,1]$ with equation (B24) defined in Appendix B. The dashed lines are the results from full branch continuation and the solid lines are for weakly nonlinear analysis.

Notice that the supercriticality or the subcriticality of the bifurcation depends on the sign of $\operatorname {Re}(\lambda _1)$ and $\operatorname {Re}(\lambda _2)$ . If $\operatorname {Re}(\lambda _1)\cdot \operatorname {Re}(\lambda _2)\lt 0$ , the system will reach the stable steady solutions when $r\gt 0$ . On the contrary, if they have the same sign, non-trivial solutions exist when $r$ is negative in the system for $Ra\lt Ra_c$ and is normally known as a subcritical bifurcation (Schmid & Henningson Reference Schmid and Henningson2002). By examining different factors, we determined that the instability of this configuration is always supercritical for all $\delta$ and $M$ considered. As shown for results obtained with weakly nonlinear analysis (figure 4), there is a bifurcation branch at $\sqrt {\langle \boldsymbol{u},\boldsymbol{u}\rangle }=0$ that is stable for $Ra\lt Ra_c$ and unstable for $Ra\gt Ra_c$ . At $Ra_c$ , a stable branch of solutions emerges. This means that the finite-amplitude equilibrium state would be achieved if the static state loses stability. An obvious finding is that with an increase of $M$ , the greater the temperature difference in the USL, the stronger the buoyancy, and therefore the greater the kinetic energy.

Figure 5.

Schematic diagram of mean temperature profiles $\bar {\theta }$ of base state and convection. The red line represents the base state and the black line convection. The vertical dashed line gives the position of $\theta =\theta _M$ .

Figure 6.

The flow field and temperature field predicted by full branch continuation (a) and weakly nonlinear analysis (b) techniques at point $\diamond$ . The colour bar marks the contours of the temperature field.

Compared to the Veronis model (Veronis Reference Veronis1963), the involvement of the internal heating source and a perfect insulating bottom layer in our model leads to different interesting dynamics. The Veronis model argued that with a finite-amplitude perturbation, the mixing would increase the thickness of the unstably stratified layer (Moore & Weiss Reference Moore and Weiss1973), resulting in loss of instability at smaller $Ra\lt Ra_c$ . Such convection is sustainable via the nonlinear mechanism. Yet in our model, the nonlinear mechanism forms a homogeneous layer near the bottom. In this case, the USL becomes thinner (as shown in figure 5) and thus enhances stability. Therefore, convection always appears via a supercritical pitchfork bifurcation.

By comparing the results obtained from the weakly nonlinear analysis with those derived from a fully nonlinear branch, the weakly nonlinear solutions are in agreement with the fully nonlinear solution in the neighbourhood of critical states and start to diverge for higher Rayleigh numbers away from the critical point (see figure 4). Such divergence is originated from the truncation error in weakly nonlinear methods, which would be larger for stronger convection intensity, i.e. larger $||u||_E$ , induced from more intense radiation, i.e. larger $M$ . Naturally, this divergence results from the weakly nonlinear analysis using a three-mode expansion as opposed to the full branch continuation curve which uses 60 Fourier modes. This low-order approximation serves as a valuable method for comparing the behaviour of the system under various conditions in the neighbourhood of critical points. Figure 6 presents a comparison between steady solutions derived from the two methods, showing both flow fields and temperature fields at the specific point $\diamond$ (marked in figure 4). It has been observed that there is one pair of convective cells exhibiting highly similar convection patterns.

Figure 7.

Variation of equilibrium amplitude $|A_{eq}|$ as a function of Rayleigh number for different values of $Pr$ at $\delta =0.3$ .

Furthermore, our model can be extended to other materials. For example, the density of tellurium–cadmium–mercury alloys ( $Pr \approx 0.1{-} 0.5 \lt 1$ ) also quadratically depends on temperature, similar to water. Thus, we considered a wide range of $Pr$ to study its influence on the nonlinear solutions. However, most of our results are for $Pr=10$ , i.e. for freshwater. We also determined that the bifurcation behaviours are supercritical for different Prandtl numbers. Figure 7 illustrates the equilibrium amplitude ( $A_{eq}$ ) away from the critical point. With higher $Pr$ , the magnitude of $|A_{eq}|$ increases. We assume that when the thermal diffusivity remains constant, a higher $Pr$ corresponds to a more viscous fluid. Thus, it means that the system needs more energy to drive fluid motion and results in a larger equilibrium amplitude $|A_{eq}|$ . The observation that $|A_{eq}|$ becomes invariant for $Pr \gt 10$ indicates that beyond a certain threshold, further increases in $Pr$ do not significantly affect the equilibrium amplitude. This indicates that the influence of the Prandtl number diminishes, as inertial effects become negligible when $Pr$ is large.

5. Fully nonlinear dynamics

When $Ra\gt Ra_c$ , the flow becomes unstable and convection cells will form. We seek steady convection cell solutions using the Newton method and a pseudo-arclength continuation technique to track the branches. The details of this method can be found in Appendix C. To obtain non-trivial steady solutions, we added the critical mode to the base state as the initial guess of the Newton iteration. Moving slightly away from the critical point indicated by the linear stability, convergence happens quickly when the amplitude of the eigenvector is well chosen. The amplitude can be chosen from the weakly nonlinear analysis. We also used DNS and simulated 10 thermal times to understand the behaviour of the system, and the snapshots are taken at $t_{thermal} = 10$ s. A Fourier-spectral method was employed to discretise the spatial domain and a finite-difference method was used with non-uniform grids that cluster near the walls (Chebyshev points). The time discretisation was accomplished with a second-order Adams–Bashforth–Crank–Nicolson method. The numerical method is adapted based on Ding & Wu (Reference Ding and Wu2021) and Ding & Ouyang (Reference Ding and Ouyang2023).

The bifurcation diagram is plotted in the $Ra{-}Nu$ plane for varied values of $\delta$ . We chose the fixed domain width $[-\pi /k_c,\pi /k_c]$ , where $k_c$ is the critical wavenumber of instability. Near the critical bifurcation points, figure 8 aligned with the results gained from the weakly nonlinear analysis that the instability is supercritical for all $\delta$ . The bifurcation of the system is calculated up to $Ra \sim O(10^6)$ , alongside two-dimensional DNS results, where the average Nusselt number is obtained after convection is fully developed. When $\delta =0.8$ , figure 9(a) shows the mean temperature profiles $\bar {\theta }$ for two-dimensional exact coherent states for $Pr =10$ and different $Ra$ . As $Ra$ increases and convection strengthens, the dimensionless temperature decreases, and the interior of the system strongly mixes and becomes closer to isothermal. When convection is sufficiently strong, the maximum value of $\bar {\theta }$ shifts from the bottom boundary to a location inside the layer. Figure 9(b,c) illustrates the flow fields and the temperature fields of the steady solution and unsteady motions for $Ra=10^6$ . The DNS results in figure 9(c) show that the unsteady motions cause the fluid to mix strongly and become closer to isothermal. The temperature becomes uniform within the mixing region. Plumes arise from the USL upper boundary layer and vigorously blend the fluid in the remaining area. Together with the numerical results shown in figure 8, the Nusselt number of the unsteady motion (figure 9 c) is $Nu\approx 4$ , whereas the steady solution (figure 9 b) has $Nu\approx 3.5$ . Thus, the unsteady motion transports heat more efficiently.

Figure 8.

The bifurcation diagram in the $Ra$ – $Nu$ plane representing exact coherent solutions (ECS) and two-dimensional DNS results for different attenuation depths $\delta$ . The Prandtl number $Pr=10$ and the dimensionless heat flux $M=18.125$ . The vertical dashed line indicates the location of $Ra = 2\times 10^4$ .

Figure 9.

(a) The mean temperature field $\bar {\theta }$ for $\delta =0.8$ . The red solid line is that of the base state and the red dashed line is that of DNS at $Ra=10^6$ . The others, from right to left, are from the steady-state solutions for $Ra = (0.01, 0.1, 1, 10, 100) \times 10^4$ . Flow fields and temperature fields from (b) two-dimensional steady solution and (c) DNS at $Ra=10^6$ .

Figure 10.

The flow field and temperature field at $Ra=20\,000$ and $Pr=10$ , $\delta = 0.3$ . (a) The upper branch solution, (b) the middle branch solution and (c) the lower branch solution. (d) The horizontal mean temperature profile $\bar {\theta }$ for the three cases. (e) The corresponding density profiles for the three cases.

Returning to figure 8, it is obvious that the system would experience a saddle-node bifurcation when $\delta$ is small and the dashed lines indicate that they are unstable. An interesting observation is that there are multiple steady solutions for the same Rayleigh number. Taking $\delta =0.3$ for example, there are three steady solutions at $Ra = 20\,000$ and figure 10(a–c) illustrates the flow fields and temperature fields for each solution. There is a pair of same-sized convection cells in the system at point (a); however, at points (b) and (c), the shape of convection cells gradually changed and several small-sized convection cells appeared, which impeded the effective heat transport with smaller $Nu$ . Our results suggest that the form of a pair of convection cells does not exist if $Ra$ exceeds a certain value, and a saddle-node bifurcation takes place and connects to the bifurcation branch of the solutions with multiple pairs of convection cells. By analysing the corresponding mean temperature and density profiles (as shown in figure 10 d,e), we determined that at case (c), the USL becomes thicker and the maximum temperature becomes higher. This suggests that multiple convection cells started to form at the upper part of the USL but with opposite movements; thus, they will impede the rise of warm flow and start to form another convection layer. In this case, the original convection area is compressed and the heat flow is accumulated, which reduces the heat transfer efficiency.

Furthermore, we compare the steady solutions and turbulent snapshots for $\delta =0.3$ and $Ra =10^6$ . The steady solution shown in figure 11(a) presented that plumes gradually formed a convection mixing layer at the upper layer of the USL. In this case, the fluid has divided into two convection layers with a decreasing staircase temperature profile with respect to heights (as shown in figure 11 c). There is a higher temperature in the bottom and the staircase structure can be usually observed in later winter ice-covered lakes (Jansen et al. Reference Jansen2021). However, the unsteady results gained from DNS presented in figure 11(b) showed that the convection is more vigorous and unsteady, and the temperature profile is almost isothermal within the interior. These results suggest that the upper convection layer produced by secondary stratification hinders the effective mixing of heat in the system.

Figure 11.

The flow field and temperature field for (a) exact coherent solution and (b) DNS result for $Ra=10^6$ at $\delta =0.3$ with other parameters $M=18.125$ and $Pr=10$ . (c) The mean temperature $\bar {\theta }$ and density $\bar {\rho }$ profiles. The blue line is the exact coherent solution (ECS) and the red line is the DNS result.

Figure 12.

(a) The bifurcation diagram in the $Ra$ – $Nu$ plane with different values of $M$ at $\delta =0.6$ and $Pr=10$ (ECS, exact coherent solution). (b) The horizontal mean temperature $\bar {\theta }$ for various values of $M$ at $Ra=10^4$ and $\delta =0.6$ .

Figure 12(a) illustrates the bifurcation diagram at the same attenuation depth $\delta =0.6$ with different dimensionless incoming heat flux $M$ . It is obvious that higher $M$ provides higher energy to the USL, and it is easier for the system to lose stability with a smaller critical value $Ra_c$ . For various values of $M$ , mean temperature profiles at $Ra = 10^4$ are plotted in figure 12(b). The profiles reveal that at a fixed attenuation depth $\delta$ , the position of maximum temperature is about the same and a higher incoming heat flux leads to a thinner SSL. Due to the larger thickness of the USL ( $d_{USL}$ ) and the higher temperature difference ( $\Delta T_{ USL}$ ), the characteristic Rayleigh number that is defined in terms of the USL

(5.1) \begin{equation} Ra^* = \frac {\alpha g\Delta T_{ USL}^3d_{ USL}^3}{\nu \kappa } \end{equation}

Figure 13.

(a) The mean temperature $\bar {\theta }$ and (b) density $\bar {\rho }$ profiles. The horizontal dashed lines are the dividing lines of each layer. The other parameters are $\delta =0.6$ and $M=10$ .

becomes larger, and the system becomes more efficient in transporting heat. An interesting observation is that $Nu$ becomes invariant for $M=10$ even though $Ra$ is raised. We plotted the mean temperature and density profiles (figure 13). We observed that the thickness of the USL and the temperature at the bottom layer do not change with increasing $Ra$ (as shown in figure 13). Thus, the Nusselt number is basically the same in the range of $Ra\in [2\times 10^5,10^6 ]$ . In addition, we have explored the $Nu$ scaling, $Nu\sim Ra^\gamma$ , within the range of $Ra\in [10^4, 10^6 ]$ for $M=25$ , and obtained the exponent $\gamma \approx 0.12$ by data fitting. However, when $M$ is small, the scaling law is unclear, and we may be able to get a clear scaling law when $Ra$ is larger.

Figure 14.

(a) The influence of $Pr$ is shown on the bifurcation diagram in the $Ra$ – $Nu$ plane. The blue lines are for $\delta = 0.3$ and the red lines are for $\delta =0.8$ . (b) The Nusselt number $Nu$ versus the Prandtl number $Pr$ at $Ra = 5000$ (the vertical dashed line in a) with exact steady solutions (ECS) and DNS solutions. Blue patterns are for $\delta =0.3$ and red ones are for $\delta =0.8$ .

Figure 15.

The flow field and temperature field from steady-state analysis at (a) $Pr = 0.1$ , (b) $Pr=1$ , (c) $Pr=10$ and (d) $Pr=100$ . The other parameters are $\delta =0.3$ , $Ra=5000$ . (e) The horizontal averaged temperature profile $\bar {\theta }$ for various values of $Pr$ . (f) The root mean square of the velocity $\sqrt {\langle \boldsymbol{u},\boldsymbol{u}\rangle }$ for different $Pr$ with $\delta =0.3$ and $0.8$ and $Ra=5000$ .

The influence of the Prandtl number is examined over three decades, $10^{-1}\leqslant Pr\leqslant 10^2$ , and results are illustrated in figure 14. When $Ra$ is small, the data in figure 14(a) clearly show that $Nu$ increases as $Pr$ increases. However, when $Ra$ is large, it is more reasonable to compare the heat transport of optimal steady-state or turbulent flow (Ding & Wu Reference Ding and Wu2021), and we leave this for future study. It is obvious that the behaviour of the system near $Ra_c$ is always supercritical, which is aligned with the outcomes from weakly nonlinear analysis (as shown in figure 7). When $Pr$ exceeds a certain threshold, the flow pattern in the whole field becomes similar; for example, saddle-node bifurcations also appear when $Pr=100$ , which is similar to the case with $Pr = 10$ for $\delta =0.3$ . Figure 14(b) illustrates the influence of $Pr$ on $Nu$ for $Ra=5000$ . We observed that DNS results are aligned with exact coherent solutions and $Nu$ slightly increases when $Pr\lt 10$ , then becomes invariant when $Pr\gt 10$ . Figure 15 displays the variations in fluid flow patterns and temperature fields with different $Pr$ when $Ra=5000$ . It is obvious that the system has a cooler bottom temperature when the Prandtl number increases (as shown in figure 15 e). When $Pr$ becomes larger, $\sqrt {\langle \boldsymbol{u},\boldsymbol{u}\rangle }$ increases, which makes the system strongly mixed and results in a higher $Nu$ . However, beyond a threshold, the flow pattern seems not to change and $\sqrt {\langle \boldsymbol{u},\boldsymbol{u}\rangle }$ decreases slightly (as shown in figure 15 f). The value of $Nu$ slightly increases then becomes insensitive to $Pr$ , implying that a weaker flow can also transport more heat, which is somewhat counterintuitive.