2.1. Problem formulation

We consider a three-dimensional (3-D) porous domain with uniform porosity $\phi$ and permeability $K$ . The medium is fully saturated with fluid, whose density depends on the temperature field, $T^*$ . The local gradients of fluid density within the system drive the flow. Figure 1 shows a sketch of the domain; here, two different flow configurations, discussed in detail in the following, are described.

We indicate by $x^*$ and $y^*$ the horizontal directions and by $z^*$ the vertical direction (perpendicular to the walls) along which the gravitational acceleration $\boldsymbol{g}$ is aligned. The temperature field $T^*$ , averaged over the reference elementary volume, varies between $T^*_{\textit{min}}$ and $T^*_{\textit{max}}$ , and it is controlled by the advection–diffusion equation (Nield & Bejan Reference Nield and Bejan2017):

(2.1) \begin{align} [(1-\phi )(\rho c)_s + \phi (\rho c)_{\!f}] \frac {\partial T^*}{\partial t^*} + (\rho c)_{\!f} \boldsymbol{u}^*\boldsymbol{\cdot }\boldsymbol{\nabla} ^* T^* = \boldsymbol{\nabla} ^*\boldsymbol{\cdot }([(1-\phi )\lambda _s + \phi \lambda _{\!f}] \boldsymbol{\nabla} ^* T^*), \end{align}

where $s$ and $f$ refer to the solid and fluid phases, respectively, $c$ is the specific heat (considered at constant pressure for the fluid) and $\lambda$ the thermal conductivity (the superscript $^{*}$ indicates dimensional variables). In this case, we consider a system in local thermal equilibrium ( $T_{\!f} = T_s = T$ ) and without heat production. In addition, when the heat capacity of the solid is negligible ( $(\rho c)_s\to 0$ ), (2.1) reduces to (Hewitt Reference Hewitt2020)

(2.2) \begin{equation} \phi \frac {\partial T^*}{\partial t^*}+\boldsymbol{u}^{*}\boldsymbol{\nabla} ^{*} \boldsymbol{\cdot }T^* =\boldsymbol{\nabla} ^{*} \boldsymbol{\cdot }(\kappa \boldsymbol{\nabla} ^{*} T^*), \end{equation}

where $t^*$ is time, $\boldsymbol{u}^{*}=(u^{*},v^{*},w^{*})$ is the volume-averaged velocity field, $\phi$ is the porosity of the medium and $\kappa = [(1-\phi )\lambda _s + \phi \lambda _{\!f}]/(\rho c)_{\!f}$ is the thermal diffusivity of the medium, considered constant here. The fluid density, $\rho ^{*}$ , is assumed to depend linearly on temperature:

(2.3) \begin{equation} \rho ^{*}(T^*)=\rho ^{*}(T^*_{\textit{min}})-\Delta \rho ^{*}\frac {T^*-T^*_{\textit{min}}}{T^*_{\textit{max}}-T^*_{\textit{min}}} \text{ ,} \end{equation}

with $\Delta \rho ^{*}=\rho ^{*}(T^*_{\textit{min}})-\rho ^{*}(T^*_{\textit{max}})$ being the maximum density contrast within the domain. Assuming the validity of the Boussinesq approximation (Landman & Schotting Reference Landman and Schotting2007), the flow field is described by the continuity and the Darcy equations,

(2.4) \begin{equation} \boldsymbol{\nabla} ^{*}\boldsymbol{\cdot }\boldsymbol{u}^{*}=0,\quad \boldsymbol{u}^{*}=-\frac {K}{\mu }\left (\boldsymbol{\nabla} ^{*} P^{*}+\rho ^{*}g \boldsymbol{k}\right )\!, \end{equation}

respectively, with $\mu$ the (constant) fluid viscosity, $P^{*}$ the pressure and $\boldsymbol{k}$ the vertical unit vector. The walls are considered impermeable to the fluid, i.e. the following boundary conditions for the flow field apply:

(2.5) \begin{equation} w^*(z^*=0)=0,\quad w^*(z^*=H^{*})=0 . \end{equation}

Note that a slip at the walls is possible, due to the Darcy formulation adopted. Periodicity is considered across the vertical boundaries (directions $x,y$ ) for all flow variables.

Regarding the temperature field, at the upper and lower walls, both Dirichlet ( $T^*$ fixed) or Neumann ( $\partial T^*/\partial z^*$ fixed) boundary conditions are considered:

1. (i) Semi-infinite domain: the temperature is fixed at the top and no-flux condition is assumed at the bottom (see figure 1 a), namely

   (2.6) \begin{equation} T^*(z^*=H^*)=T^*_{\textit{min}}, \quad \frac {\partial T^*}{\partial z^*}\biggr |_{z^*=0} =0 . \end{equation}
   In addition, the domain is initially filled with still fluid and at temperature $T^*_{\textit{max}}$ . Thus, the velocity and temperature fields are initialised as $\boldsymbol{u^*}=0$ and
   (2.7) \begin{equation} T^*(x^*,y^*,z^*\lt H^*,t_0^*)=T^*_{\textit{min}} + (T^*_{\textit{max}}-T^*_{\textit{min}})\,\text{erf}\left [\sqrt {\frac {\phi }{\kappa }}\frac {(H^*-z^*)}{2\sqrt {t^*_0}}\right ]\!, \end{equation}
   respectively. The initial condition (2.7) is obtained evaluating at time $t_0^*$ the self-similar solution of (2.2) in the absence of convection ( $\boldsymbol{u}^*=0$ ) in a semi-infinite domain $(H^*\to \infty )$ (De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2017). The instant considered to initialise the simulation is $t_0^*=250 \phi \kappa (g \Delta \rho ^{*} K / \mu )^{-2}$ , corresponding to a dimensionless time of $\widehat {t}_0=250$ or $t_0=250/{Ra}$ , expressed in diffusive or convective units, respectively (see § 2.2 for further details). Finally, a random perturbation (white noise) is added. The amplitude of this noise controls the time at which the plumes form, i.e. the onset of convection occurs (De Paoli et al. Reference De Paoli, Zonta, Enzenberger, Coliban and Pirozzoli2025d ). This flow configuration, also labelled as ‘one-sided’, is representative of a domain cooled from above and insulated from below, and will be analysed in § 4.

2. (ii) Two-sided configuration: the temperature is imposed at both boundaries (see figure 1 b), such that the corresponding density profile is unstable:

   (2.8) \begin{equation} T^*(z^*=H^*)=T^*_{\textit{min}},\quad T^*(z^*=0)=T^*_{\textit{max}}. \end{equation}
   The velocity field is again initialised as $\boldsymbol{u^*}=0$ , while a linear profile is used for the temperature distribution,
   (2.9) \begin{equation} T^*(x^*,y^*,z^*,t^*=0) = T^*_{\textit{max}} - \frac {z^*}{H^*}(T^*_{\textit{max}}-T^*_{\textit{min}}) , \end{equation}
   corresponding to the initial condition labelled as IC1 in Hewitt et al. (Reference Hewitt, Neufeld and Lister2014). This configuration, labelled as ‘Rayleigh–Bénard’ or ‘two-sided’, is analysed in § 5.

2.2. Dimensionless equations

Depending on whether the domain height $H^*$ or an intrinsic diffusive length scale $\ell ^*$ (defined later) is used to make the flow quantities dimensionless, two different formulations can be derived (convective or diffusive units, both used in this work), described in §§ 2.2.1 and 2.2.2, respectively.

2.2.1. Convective units

Natural flow scales relevant to the convective system considered are the buoyancy velocity, $\mathscr{U}^{*}=g \Delta \rho ^{*} K / \mu$ , and the domain height, $H^{*}$ . Using the following set of dimensionless variables:

(2.10) \begin{equation} T=\frac {T^*-T^*_{\textit{min}}}{T^*_{\textit{max}}-T^*_{\textit{min}}},\quad x=\frac {x^{*}}{H^{*}},\quad \boldsymbol{u}=\frac {\boldsymbol{u}^{*}}{\mathscr{U}^{*}}, \end{equation}

(2.11) \begin{equation} t=\frac {t^*}{\phi H^{*}/\mathscr{U}^{*}},\quad p=\frac {p^{*}}{\Delta \rho ^{*}gH^{*}} , \end{equation}

and introducing the reduced pressure $p^{*}=P^*-\rho ^*(T^*_{\textit{min}})gz^*$ , we obtain the dimensionless form of the governing equations (2.2)–(2.4):

(2.12) \begin{equation} \frac {\partial T}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }T = \frac {1}{{Ra}} \boldsymbol{\nabla} ^2 T , \end{equation}

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

(2.14) \begin{equation} \boldsymbol{u}=-\left ( \boldsymbol{\nabla }\!p - T \boldsymbol{k}\right ) \!, \end{equation}

where

(2.15) \begin{equation} {Ra}=\frac {g \Delta \rho ^{*} \textit{KH}^{*} }{ \kappa \mu }=\frac {\mathscr{U}^{*} H^{*}}{\kappa }. \end{equation}

The flow is controlled by three dimensionless parameters: the Rayleigh number $Ra$ and the domain width in horizontal directions, $L_1=L^*_1/H^{*}$ and $L^*_2/H^{*}$ .

A distinction should be made depending on the behaviour of the solid phase concerning the transported scalar. When the solid phase is impermeable to the scalar field considered, e.g. in the case of mass transport problems, the definition (2.15) of the Rayleigh number changes as ${Ra}=\mathscr{U}^{\,*}H^*/(\kappa \phi )$ . However, the dimensionless formulation of the problem is the same for both permeable and impermeable behaviour, provided that when the permeable case is considered, the two phases are in local thermal equilibrium. Further details on this matter are provided by De Paoli (Reference De Paoli2023).

2.2.2. Diffusive units

An alternative non-dimensionalisation approach involves using $\ell ^*=\kappa /\mathscr{U}^{\,*}$ as the reference length scale. This strategy was introduced by Fu et al. (Reference Fu, Cueto-Felgueroso and Juanes2013) and is suitable for describing local flow dynamics that is not influenced by the largest scales of the flow. Defining the dimensionless variables as

(2.16) \begin{equation} T=\frac {T^*-T^*_{\textit{min}}}{T^*_{\textit{max}}-T^*_{\textit{min}}},\quad \widehat {x}=\frac {x^{*}}{\ell ^{*}},\quad \boldsymbol{u}=\frac {\boldsymbol{u}^{*}}{\mathscr{U}^{*}}, \end{equation}

(2.17) \begin{equation} \widehat {t}=\frac {t^*}{\phi \ell ^{*}/\mathscr{U}^{*}},\quad \widehat {p}=\frac {p^{*}}{\Delta \rho ^{*}g \ell ^{*}} , \end{equation}

we obtain the following dimensionless form of the governing equations (2.2)–(2.4):

(2.18) \begin{equation} \frac {\partial T}{\partial \widehat {t}} + \boldsymbol{u} \boldsymbol{\cdot }\widehat {\boldsymbol{\nabla }} T = \widehat {\boldsymbol{\nabla} }^2 T , \end{equation}

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

(2.20) \begin{equation} \boldsymbol{u}=-\left ( \widehat {\boldsymbol{\nabla }} \widehat {p} - T \boldsymbol{k}\right ) \!, \end{equation}

where $\widehat {\ \boldsymbol{\cdot }\ }$ indicates variables made dimensionless with respect to diffusive units. Note that in this case $Ra$ does not appear explicitly in the equations; instead, it represents the dimensionless height of the system, ${Ra}=\mathscr{U}^{\,*} H^{*}/ \kappa =H^*/\ell ^*$ .

2.3. Data analysed

The data analysed in this work consist of a collection of results presented in previous works (Pirozzoli et al. Reference Pirozzoli, De Paoli, Zonta and Soldati2021; De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022, Reference De Paoli, Zonta, Enzenberger, Coliban and Pirozzoli2025d ) and original simulations presented here for the first time (see Appendix A for a detailed description of the database). The relevant part of the data analysed is made freely available via De Paoli et al. (Reference De Paoli, Pirozzoli and Kondic2025a ).

The original simulations presented here (simulations A in table 1) are obtained as follows. The governing equations (2.12)–(2.14) are solved numerically with the aid of an (open-source) second-order finite-difference solver, AFiD-Darcy (De Paoli Reference De Paoli2025), which was previously employed to simulate convective porous media flows (De Paoli et al. Reference De Paoli, Yerragolam, Verzicco and Lohse2025b , Reference De Paoli, Yerragolam, Lohse and Verziccoc ). The code comprises an efficient parallel solver for simulating convective, incompressible and wall-bounded flows in porous media. This solver is based on the initial version of AFiD developed for turbulent flows (Van Der Poel et al. Reference Van Der Poel, Ostilla-Mónico, Donners and Verzicco2015). The algorithm utilises a pressure-correction scheme in conjunction with an efficient fast-Fourier-transform-based solver. Additional details of the algorithm are described by De Paoli et al. (Reference De Paoli, Yerragolam, Lohse and Verzicco2025c ).

3.1. Persistent homology: an example

Here we discuss a 2-D example. The basic geometric structures of interest are connected components and loops (holes), denoted by $\beta _0$ and $\beta _1$ , respectively. To illustrate the topological quantities describing such structures, we use the simple 2-D example shown in figure 2. The reader is referred to Kramár et al. (Reference Kramár, Goullet, Kondic and Mischaikow2013, Reference Kramár, Goullet, Kondic and Mischaikow2014b ) for a more in-depth discussion, along with simpler one-dimensional examples; additional toy examples in the context of granular media analysis can be found in Kramár et al. (Reference Kramár, Goullet, Kondic and Mischaikow2014a ) and Taghizadeh et al. (Reference Taghizadeh, Luding, Basak and Kondic2024) and in the context of porous media flow in Illingworth et al. (Reference Illingworth, Gu, Cummings and Kondic2026). Specifics related to computations of $\mathsf{PD}$ s are discussed later in § 3.2.

The example shown in figure 2 is a two-sided convection simulation from § 5, discussed later in the text (settings: ${Ra}=500$ , case C5 in table 2, field taken at the boundary layer), with the colour map showing the temperature. Figure 2 shows two threshold values, namely $0.5$ in figure 2(a) and $0.3$ in figure 2(b); the parts of the domain characterised by the temperature values lower than the specified threshold are not shown and appear light blue. By comparing figures 2(a) and 2(b), we observe a change in connectivity and structure as the threshold value changes.

Figure 2.

A 2-D example illustrating $\mathsf{PD}$ s for a function of two spatial coordinates, $x$ and $y$ (here we use the temperature results that will be discussed later in the text). The results for chosen thresholds of $T = 0.5$ (a) and $T = 0.3$ (b) – the light-blue areas are those where $T$ is less than the specified threshold value. Persistence diagrams corresponding to (c) $\beta _0$ (components) and to (d) $\beta _1$ (holes), and schematic examples of which are reported in (a); here $T = 0.5$ (black dashed line) and $T = 0.3$ (blue dashed line) show the threshold levels. The fields shown in (a,b) are two snapshots from simulation C5 (see table 2 for additional details). We refer to supplementary movie S1 for a visual interpretation of the threshold-dependent temperature field and of the corresponding $\mathsf{PD}$ s.

The change in connectivity, although visually apparent, needs to be quantified. Such quantification can be established by considering $\mathsf{PD}$ s. Let us consider first the PD corresponding to components, $\beta _0$ $\mathsf{PD}$ , shown in figure 2(c) (see also supplementary movie S1 available at https://doi.org/10.1017/jfm.2026.11804). As one descends from high values of temperature to lower ones, separated areas of high temperature appear (are born); each of these areas is called a component. Decreasing the threshold causes these components to merge; whenever two components merge, the ‘younger’ one (born at a lower temperature value) disappears (dies). The coordinates of each point (generator) in $\beta _0$ $\mathsf{PD}$ specify birth ( $b$ ) and death ( $d$ ) temperature for each component; the difference between birth and death values is called lifespan, $\mathscr L$ , and shows the range of temperatures over which a component ‘lived’ before merging with an older one. Therefore, the generators near the diagonal are associated with minor spatial variations in temperature. All key trends were verified to be insensitive to variations of the persistence threshold and to the number of temporal samples (Appendix B.2), confirming that the observed behaviour is not an artefact of post-processing choices. In contrast, those far from the diagonal are related to components that remained separated from others over an extensive range of temperatures (one could use a landscape analogy and say that there is a deep ‘valley’ surrounding the local temperature maximum). One generator that is treated slightly differently is the one that corresponds to the component that appeared first (at the largest value of temperature); this component never dies and therefore has a death coordinate equal to zero. Inspecting carefully $\beta _0$ $\mathsf{PD}$ shows that this particular generator has the birth coordinate close to unity in the example shown in figure 2. Large ‘clump’ of generators in $\beta _0$ $\mathsf{PD}$ , with coordinates close to $(0.8, 0.8)$ , shows that there are multiple regions of the temperatures larger than 0.8 which merged at similar values of temperature threshold; checking figure 2(a), we observe that these generators are due to large temperatures of the narrow regions connecting the temperature ‘nodes’ (inspecting supplementary movie S1 shows that this is indeed the case; the advantage of considering the $\mathsf{PD}$ instead is that this information is provided in a much simpler form).

Still considering $\beta _0$ $\mathsf{PD}$ , but shifting our focus to lower temperature values shows much less activity in the particular example shown in figure 2, suggesting there are few topological changes at these temperatures. This diagram also shows that there are no topological changes (considering components only) below the temperature value of $\approx 0.3$ ; the temperature field is fully connected. The reader familiar with the concept of Betti numbers (counting simply the number of components) will note that Betti number $\beta _0$ can be trivially obtained from the $\beta _0$ $\mathsf{PD}$ ; the opposite is, however, not true, since a $\mathsf{PD}$ contains significantly more detailed information about connectivity of the considered scalar field than the corresponding Betti number. Most importantly, a $\mathsf{PD}$ includes the information over all thresholds (temperature values), while Betti numbers are threshold-dependent. Furthermore, $\mathsf{PD}$ s are stable with respect to noise (small perturbations), while the same cannot be claimed for Betti numbers (Kramár et al. Reference Kramár, Goullet, Kondic and Mischaikow2014b ). We also note in passing that some works present the information contained in $\mathsf{PD}$ s using different visual representations, such as bar-codes (Carlsson Reference Carlsson2009), or sometimes consider sub-level thresholding instead of super-level as considered here; the connection between different graphical representations and thresholding approaches is straightforward.

In the 2-D example considered in figure 2, we also have $\beta _1$ $\mathsf{PD}$ , corresponding to holes (loops) in the temperature field, shown in figure 2(d). A loop appears (is born) when the temperature threshold is decreased sufficiently that a closed structure, enclosing the area of lower temperatures, forms. A loop dies when the threshold is decreased to the degree that the whole interior of the loop is filled up. As expected, both birth and death coordinates of loops are typically lower than those of components; this is observed by inspecting the two $\mathsf{PD}$ s shown in figure 2(c,d). Supplementary movie S1 helps in fully grasping the process of loop formation and disappearance.

We conclude the discussion of the presented example by noting that $\mathsf{PD}$ s provide significant data reduction, as they condense the information about the complicated two-variable function $T(x,y)$ into a point cloud ( $\mathsf{PD}$ ). This data reduction entails some information loss, primarily due to the loss of geometrical information ( $\mathsf{PD}$ includes connectivity information but not the size (spatial extent) of each component). If such information is of importance, then one wants to combine the information available by considering $\mathsf{PD}$ s by a complementary measure. Furthermore, although $\mathsf{PD}$ s provide reasonably complete information about the connectivity of the considered field, they are still point clouds, and one can come up with measures summarising a $\mathsf{PD}$ . While many elaborate approaches are being considered to achieve this goal, in this work we resort to a simple approach and account for three measures only (separately for $\beta _0$ and $\beta _1$ ): the number of generators $N_g$ , the average lifespan $\mathscr L$ (the average distance of the generators to the diagonal) and total persistence $\textit{TP} = {\mathscr L} N_g$ . One could think of $\textit{TP}$ as a simple measure describing (admittedly in a vague manner) the ‘intensity’ of a considered $\mathsf{PD}$ . We also note that to focus on the main features of the results, we exclude from consideration generators with ${\mathscr L} \lt 0.01$ , as the fluctuations of the temperature field on such a scale are not of interest. Such exclusion does not influence the statistics of the (relevant) long-lived generators. In §§ 4.2 and 4.3 we comment in more detail on the influence of the specific value (0.01) on the results.

3.2. Persistent homology: specifics

Over the past decade, several research groups have developed and released open-source software tools for computing $\mathsf{PD}$ s and related topological measures. All computations and results presented in this work were carried out using the GUDHI library (GUDHI 2025). In our calculations, we use the results obtained from the simulations directly, without any additional post-processing. Consistent with the simulation set-up, we assume periodic boundary conditions when computing topological measures (though using non-periodic boundaries yields only minor changes in the results). As illustrated by our toy example, figure 2, we use super-level filtration (meaning that the features that appear above the chosen threshold value are considered). For the considered data, the use of computing resources is very reasonable; e.g. for one data slice of simulation B7 (see table 1 in Appendix A), the computing time to build a cubical complex and compute persistence is $\approx 1.41$ s (on M2 Mac Studio with 128 GB of RAM). An example Python script for interacting with the GUDHI library is included in the supplementary material.

5.1. Convection regimes

We consider here a system heated from below and cooled from above, where a high (low) temperature is fixed at the bottom (top) wall, corresponding to the boundary conditions specified in (2.8). The domain is periodic in the horizontal directions, and the flow is initialised with a linear temperature distribution as in (2.9). After an initial transient phase, the flow attains a statistically steady state. The boundary conditions and an example of temperature distribution on the lateral boundaries of the domain are illustrated in figure 1(b).

In free fluids, i.e. in the absence of a porous medium, the two-sided flow is also called Rayleigh–Bénard flow, and it is particularly suitable to be analysed theoretically, due to its well-defined boundary conditions and its statistically steady nature: it represents a perfect candidate to develop new concepts on instabilities and dynamical systems (Lohse & Shishkina Reference Lohse and Shishkina2024). For the same reasons, the porous counterpart (also called Rayleigh–Bénard–Darcy or Rayleigh flow, where the free fluid is replaced by a fluid-saturated porous medium) has been widely investigated. In particular, extensive studies focused on the onset of the flow instabilities ( $0\lt {Ra}\lt O(10^2)$ ; see Horton & Rogers Jr Reference Horton and Rogers1945; Lapwood Reference Lapwood1948) and on the dynamics at moderate to high $Ra$ ( $O(10^2)\lt {Ra}\lt O(10^4)$ ; see Graham & Steen Reference Graham and Steen1994; Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004). In these systems, the global response parameter is the Nusselt number $\textit{Nu}$ defined as in (4.1) which, for ${Ra}\gt 10^3$ , scales as $\textit{Nu}\sim {Ra}$ plus sublinear corrections (Zhu et al. Reference Zhu, Fu and De Paoli2024). As a result, an increase of $Ra$ by a factor 10 would require substantially increasing the number of degrees of freedom required to resolve the flow, namely by a factor of $O(10^3)$ in three dimensions, corresponding to an even larger increase of the computational cost. In recent decades, however, numerical advancements have allowed one to explore in detail the flow dynamics in three dimensions and at large $Ra$ (Hewitt et al. Reference Hewitt, Neufeld and Lister2014; Pirozzoli et al. Reference Pirozzoli, De Paoli, Zonta and Soldati2021; De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022; Hu et al. Reference Hu, Xu and Yang2023), and a rich flow morphology has been observed. A review for ${Ra}\leqslant 740$ is provided by Hewitt et al. (Reference Hewitt, Neufeld and Lister2014), which we summarise and extend in the following.

Figure 13.

Temperature distribution over a portion of the domains considered, namely $0\leqslant x \leqslant 4$ , $0 \leqslant y \leqslant 4$ and $0\leqslant z \leqslant z_t$ (with $0.95\leqslant z_t \leqslant 0.999$ , depending on $Ra$ ). Distributions shown are relative to simulations C2 (a), C9 (b) and C17 (c) (see table 2 for additional details).

Figure 13 shows a few examples of the temperature distributions for three values of $Ra$ , illustrating the changes of the flow morphology. To highlight the complex near-wall flow pattern, the temperature field is shown on a horizontal slice located near the top boundary and located at $z=z_t$ ( $0.95\leqslant z_t \leqslant 0.999$ , depending on $Ra$ ), focusing on ${Ra} \geqslant 100$ . For smaller values (specifically, for ${Ra} \lt 4\pi ^2$ ) any instability is suppressed, and the flow is purely conductive. When $Ra$ is slightly increased, a 2-D unstable mode appears ( $4\pi ^2 \leqslant {Ra} \lesssim 4.5\pi ^2$ ), and dominates the flow also at larger $Ra$ , when multiple modes are present ( $4.5\pi ^2 \lesssim {Ra} \lesssim 97$ ). In these regimes, temperature structures spanning the entire system height, from one boundary layer to the other, populate the domain. The flow eventually develops a 3-D steady structure ( $97 \lesssim {Ra} \lesssim 300$ ; see figure 13 a), and multiple possible states exist. Finally, the flow is unsteady for ${Ra} \gtrsim 300$ : small flow instabilities develop and grow from the thermal boundary layers, move parallel to the horizontal walls and eventually merge into larger plumes, still covering the whole domain height, but now being unsteady (figure 13 b). For sufficiently large driving, namely ${Ra}\geqslant 1750$ , the system enters the so-called high- $Ra$ regime and the flow morphology does not exhibit any regular background structure (Hewitt et al. Reference Hewitt, Neufeld and Lister2014) (see e.g. figure 13 c). In this regime, the flow can be separated into three different parts: (i) a thin and time-dependent thermal boundary layer at the walls; (ii) an intermediate region populated by very dynamical sheet-like plumes, originated from the time-dependent wall layer; and (iii) the bulk, controlled by large columnar structures of hot (cold) rising (sinking) fluid, labelled megaplumes.

The dynamics at larger $Ra$ ( ${Ra}\lesssim 10^5$ ) has been recently explored in detail by Pirozzoli et al. (Reference Pirozzoli, De Paoli, Zonta and Soldati2021) and De Paoli et al. (Reference De Paoli, Pirozzoli, Zonta and Soldati2022). Those authors observed that when ${Ra}\gtrsim 10^4$ , the temperature distribution in the intermediate region is characterised not only by the presence of the small plumes originated from the boundary layer, but also by the existence of large structures, entraining many plumes and representing the footprint of the megaplumes populating the bulk of the flow (De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022). These structures, labelled as supercells, have been identified visually as the large loops gathering many smaller cells (see e.g. the dark loops on the horizontal cut in figure 13 c), and have been observed to be: (i) persistent in time and (ii) correlated in space with the megaplumes (De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022). Time persistence has been identified by averaging the temperature distribution in the near-wall region and observing that a pattern matching that of the supercells emerges. Spatial correlation with the megaplumes was demonstrated by filtering the near-wall temperature field to remove flow structures with wavelengths smaller than the dominant wavelength at the mid-plane.

Despite these efforts, a precise definition of the supercells and an accurate determination of the conditions required for their appearance remain elusive. The following questions remain unanswered: Can we provide a more detailed description of the supercells? Can we describe their presence/formation as a function of $Ra$ ? Do we have to rely on the bulk temperature field (i.e. the temperature field far from the top and bottom boundaries) to determine the existence and the morphology of the supercells? In this work, we aim at precisely answering these questions: we provide a robust criterion to identify supercells and describe their formation across a wide range of $Ra$ .

Figure 14.

Visualisation of the temperature fields on horizontal ( $x,y$ ) planes near the bottom wall (i) and at the mid-height ( $z=1/2$ ) (ii). The ‘C’ simulations and the Rayleigh numbers $Ra$ are indicated; for all other details we refer to table 2 (see Appendix A). The size of the domain shown is constant within each column, but varies from left to right for visualisation purposes. A scale bar (width $\widehat {W}$ expressed in diffusive units) is reported for reference. In some cases (e.g. C18), only a small portion of the domain is shown, while in others (e.g. C7), the field is replicated to match the desired size. Note the difference between the cell structure in the present results and in those of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) (see the results labelled as IC1 in Hewitt et al. (Reference Hewitt, Neufeld and Lister2014), corresponding to the same initial condition employed here), suggesting the strong effect of the domain size on the determination of the flow structure at low Rayleigh numbers $({Ra}\leqslant 500)$ . The data analysed relative to the ‘C’ simulations are available via De Paoli et al. (Reference De Paoli, Pirozzoli and Kondic2025a ).

Figure 14. Long description

Panel A: A heat map showing temperature fields near the bottom wall for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high. Panel B: A heat map showing temperature fields at mid-height for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high. Panel C: A heat map showing temperature fields near the bottom wall for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high. Panel D: A heat map showing temperature fields at mid-height for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high. Panel E: A heat map showing temperature fields near the bottom wall for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high. Panel F: A heat map showing temperature fields at mid-height for different Rayleigh numbers and simulation conditions. The heat map is divided into six rows, each representing a different Rayleigh number and simulation condition. The x-axis and y-axis represent spatial dimensions in diffusive units. The color scale ranges from dark red to bright yellow, indicating temperature variations from low to high.

5.2. Qualitative analysis of the flow morphology

We analyse the flow morphology for $5\times 10^1\leqslant {Ra}\leqslant 8\times 10^4$ , by considering the horizontal temperature distributions in the near-wall region (panels (i) in figure 14) and at the midplane (panels (ii) in figure 14). For ease of comparison, the size of the portion of domains shown in figure 14 is constant (in diffusive units) within each column (scale bar is reported as a reference). The simulations considered are either those by Pirozzoli et al. (Reference Pirozzoli, De Paoli, Zonta and Soldati2021) and De Paoli et al. (Reference De Paoli, Pirozzoli, Zonta and Soldati2022), or presented here for the first time, and the data are made available via De Paoli et al. (Reference De Paoli, Pirozzoli and Kondic2025a ).

At ${Ra}\leqslant 1750$ , the flow is also dependent on the initial condition, and it exhibits hysteresis effects (Otero et al. Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004). Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) explored the dynamics following two initial conditions: IC1, corresponding to a perturbed linear temperature profile, and IC2, obtained starting from a steady-state temperature field obtained from a simulation at a lower $Ra$ . Therefore, the initial condition adds here to the flow governing parameters ( ${Ra},\widehat {L}_1,\widehat {L}_2$ ). To reduce the parameter space, here we consider only one initial condition (corresponding to IC1 of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014)) characterised by a perturbed form of the stable state, (2.9), and we consider very large domains ( $\widehat {L}_1=\widehat {L}_2\geqslant 4\times 10^3$ ), such that the flow is independent of the values of $\widehat {L}_1$ and $\widehat {L}_2$ . As a result, the only remaining governing parameter is $Ra$ . Simulations performed in this work have been initialised with a linear temperature distribution, (2.9), which in dimensionless terms reads $T=1-z$ , and have been run for $t\geqslant 1000$ , to make sure the steady state (determined by keeping track of the time-averaged value of the Nusselt number) is achieved. Then, the temperature fields used to analyse the flow pattern are saved approximately every 10 or 20 convective time units. All the relevant simulation details are indicated in table 2 (see Appendix A).

Before quantifying the flow pattern using PH techniques, we notice that a visual inspection reveals some interesting features previously unobserved. In particular, at low Rayleigh numbers the flow organisation described in the literature (see figure 2b of Hewitt et al. (Reference Hewitt, Neufeld and Lister2014)) consists of cells forming a very regular pattern that for ${Ra}\leqslant 2\times 10^2$ is one-dimensional (rolls spanning across the entire domain width) and for $2\times 10^2\leqslant {Ra}\leqslant 3\times 10^2$ is 2-D and still very regular (equally squared cells). The domains considered in figure 14 are about 100 times larger than those in Hewitt et al. (Reference Hewitt, Neufeld and Lister2014), and the cells are organised in a more chaotic manner. For ${Ra}\leqslant 2\times 10^2$ , corresponding to figure 14(a–c), we still observe the formation of rolls, which in contrast to Hewitt et al. (Reference Hewitt, Neufeld and Lister2014) do not span across the entire domain in the horizontal directions. Also for ${Ra}= 3\times 10^2$ (figure 14 d) there are differences compared with previous works, where the cells were all equal: we do observe the transition from rolls to cells, but such cells have an irregular shape and are, in general, different from each other (polygons with 3–6 sides). These differences suggest that the previously observed regular flow pattern emerging at small $Ra$ is a genuine domain size effect, which is induced by the periodic boundary conditions applied in the horizontal directions (additional details are provided in Appendix B.1). These findings highlight the fact that at small $Ra$ the influence of the domain size on the flow morphology is considerable. Finally, we observe that at very large Rayleigh numbers ( ${Ra}\geqslant 2\times 10^4$ ; figure 14 p–r), where the supercells are clearly visible, their size increases with $Ra$ . This is not surprising, since the temperature structures in the bulk (megaplumes of size $\sim {Ra}^{1/2}$ in diffusive units; see Hewitt et al. Reference Hewitt, Neufeld and Lister2014; Pirozzoli et al. Reference Pirozzoli, De Paoli, Zonta and Soldati2021) are responsible for the formation of the supercells (De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022). However, this highlights once again the importance of domain size in capturing the large-scale structures in the flow.

Figure 15.

Loop $\mathsf{PD}$ s, $\beta _1\,{\mathsf{PD}}$ s, in two-sided convection simulation near the wall (simulations C in table 2). As the Rayleigh number ( $Ra$ ) increases, the morphology of the boundary develops significant changes, leading to variations in the corresponding $\beta _1\,{\mathsf{PD}}$ . Between ${Ra}=300$ and 500, we observe a morphological pattern change, with a clear transformation from a maze to a cellular network (see also figure 14). Proceeding from ${Ra}=1000$ , a division emerges, forming two distinct clusters of points in the $\mathsf{PD}$ s. This indicates an increasing number of loops within the boundary layer, appearing at high temperatures ( $T\approx 0.75$ ). For ${Ra} \geqslant 2500$ , the cluster at high temperatures becomes denser, suggesting an increased number of loops in the patterns shown in figure 14.

5.3. Quantification of pattern formation

We focus the presentation that follows on the identification of cells and supercells and, more generally, on the influence of $Ra$ on the pattern-formation process. We discuss loop-related measures only, as they provide a clear description of the patterns of interest.

Figure 15 shows the loop ( $\beta _1$ ) $\mathsf{PD}$ s for the simulations listed in table 2. The emerging behaviour, as $Ra$ increases, is an increase in the number of generators, and, in particular, the gradual appearance of a band of generators with the birth coordinates of $T\gtrsim 0.75$ starting at ${Ra} \gtrsim 1000$ . For larger $Ra$ , the number of points in the $\mathsf{PD}$ s increases substantially, and their structure becomes difficult to visualise. For this reason, we consider cumulative measures, which simplify the interpretation of the results, as discussed next.

Figure 16.

The p.d.f. of the birth values of $\beta _1$ topological generators, accumulated over 10 time outputs for each simulation (noise band 0.01 is removed). The birth coordinate corresponding to the maximum value of the p.d.f. is indicated by the dashed line and is reported in each panel. A transition from rolls (a–c) to cells/supercells (d–r) is clearly observed.

Figure 16. Long description

Panel A: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.21. Panel B: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.49. Panel C: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.40. Panel D: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.80. Panel E: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.77. Panel F: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.80. Panel G: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.74. Panel H: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.82. Panel I: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.83. Panel J: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.81. Panel K: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.76. Panel L: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.79. Panel M: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.80. Panel N: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.81. Panel O: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.82. Panel P: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.83. Panel Q: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.84. Panel R: A histogram showing the probability density function of birth values for topological generators. The horizontal axis represents the birth values (b) ranging from 0 to 1, and the vertical axis represents the probability density function (p.d.f.) ranging from 0 to 30. The histogram has a peak around 0.85.

Figure 17.

(a) The p.d.f. of $\beta _1$ topological generators as a function of the shifted temperature, with $\overline {T}$ corresponding to the space- and time-averaged temperature in the considered horizontal slice; (b) p.d.f of the lifespan, $\overline {\mathscr L}$ (noise band 0.01 is removed). Note that the lifespan is independent of the temperature shift. We observe approximate self-similarity of the results for large values of $Ra$ .

Figure 16 shows the normalised distribution (probability density function, p.d.f.) of the birth coordinate of the loops, using the data from $\mathsf{PD}$ s shown in figure 15, but averaged over 10 time outputs (saved approximately every 10 or 20 convective time units) in order to improve statistics. We recall that the birth coordinate represents the temperature at which the loops close and, therefore, at which the cells (or supercells) form. At low values of $Ra$ , we note a dramatic change in the p.d.f. as $Ra$ increases from 200 (figure 16 c) to 300 (figure 16 d), illustrating the topological change between rolls (figure 16 c) and cells (figure 16 d). This topological change can be visualised in figure 14: the rolls initially present in the system (figure 14 a–i) start to connect, forming cells (figures 14 b–i and 14 c–i). These cells fully develop and are clearly visible at larger $Ra$ (starting with figure 14 d–i). Further increase of $Ra$ leads to development of a clear and well-defined structure of the p.d.f.s shown in figure 16 with a dominant peak at $T\approx 0.8$ corresponding to supercells (see e.g. figures 14 o–i to 14 r–i), and broad distribution of the birth values for $0.5 \lesssim T \lesssim 0.7$ , corresponding to cells. We note that the portions of the flow in which $T \approx 0.75$ roughly correspond to the regions in which the horizontal temperature gradient is maximum (De Paoli et al. Reference De Paoli, Pirozzoli, Zonta and Soldati2022). We also comment that the results for p.d.f.s of connected components do not provide any useful information; therefore, considering loops is crucial for describing cells and supercells.

Further insight can be obtained by considering the loops’ birth coordinate as a function of $T-\overline {T}$ , where $\overline {T}$ is the temperature averaged in space (over the horizontal plane considered) and in time (over 10 time outputs as also implemented for the plots shown in figure 16). Figure 17 plots combined results, showing the $\beta _1$ birth coordinate distribution and the $\beta _1$ lifespan. Most importantly, the approximate collapse of the curves for ${Ra}\gtrsim 10^4$ for both birth values in figure 17(a) and lifespans in figure 17(b) illustrates the validity of the topological properties describing the temperature field across the wide range of parameters considered. More precisely, the collapse in figure 17(a) shows that the supercells appear at similar values of the shifted temperature $T - \bar {T}$ , while the collapse visible in figure 17(b) shows that the cells also die (merge with other cells) at similar values of the shifted temperature (since lifespan measures difference between the birth and death coordinates). For smaller values of $Ra$ , while the lifespans shown in figure 17 are too noisy to reach any useful insight, the birth coordinate distribution shows a clear shift in the peak of p.d.f.s from $T - \bar {T} \sim -0.6$ for ${Ra} = 50$ to $T - \bar {T} \sim -0.2$ for $100 \leqslant {Ra} \leqslant 200$ , and then to $T - \bar {T} \sim 0.2$ for larger values of $Ra$ . Additional details are provided in Appendix B.2, where we prove the robustness of our findings, with respect to the number of samples and the different aspect ratios considered.