1. Introduction
When a severe accident (SA) occurs in a nuclear power plant, the radioactive fuel and reactor metallic components melt and form a fluid called corium. The corium relocates from the core to the lower plenum of the reactor vessel, where non-miscible oxidic and metallic phases separate: the oxide phase contains the majority of the decay heat from radioactive elements and heats from below the less dense liquid metal phase floating at the surface (figure 1a). This top metal layer is usually thinner compared with the oxide layer and, thus, concentrates the power from the oxide to the vessel wall. This phenomenon is often referred to as the ‘focusing effect’ in the nuclear safety literature. When the external cooling of the reactor vessel is implemented as a SA management strategy (Theofanous et al. Reference Theofanous, Liu, Additon, Angelini, Kymäläinen and Salmassi1997; Carénini, Fichot & Seignour Reference Carénini, Fichot and Seignour2018), addressing the heat transfer through the top metal layer is fundamental to predict the failure of the vessel or to justify its integrity. This work focuses on this issue: a liquid metal layer considered as a cylindrical layer heated from below, cooled at the side with a constant temperature (assuming that the wall is being ablated and, thus, maintained at its melting temperature) and cooled at the top by radiative heat transfer. One difficulty of the problem lies in the top boundary condition, due to the interdependence of variables: the radiative heat flux from the upper surface depends on the surface temperature, which is affected by heat transfer within the metal layer, which in turn depends on the efficiency of the radiative heat flux. The approach proposed is to prescribe a uniform heat flux leaving from the top of the layer and focus on analysing the temperature profiles, fluid behaviour and heat transfers. This will allow correlating these outputs with the input control parameters and encompass all possible configurations within the reactor. Indeed, depending on the height of the metal layer and on the state of the reactor vessel structures above the pool, the radiative heat transfer can either play a major role for the power dissipation from the metal or be negligible. This approach has also the advantage of decoupling the study of the metal layer from the modelling of the radiative heat transfer. As illustrated in Le Guennic et al. (Reference Le Guennic, Skrzypek, Skrzypek, Bigot, Peybernes and Le Tellier2020), the consideration of the top radiative heat transfer introduces assumptions on its modelling directly in correlations established for the behaviour of the metal layer. With the present approach, coupling with more detailed radiative heat transfer models will be possible (Rein et al. Reference Rein, Carénini, Fichot, Le Bars and Favier2023). The impact of considering a uniform heat flux compared with a more realistic radiative exchange will be evaluated in future studies.
Figure 1. (a) Sketch of the corium phase stratification in the lower plenum of a nuclear reactor vessel during a SA. (b) Sketch of the modelled fluid metal layer in a vertical plane through the cylinder.
Given the specific boundary conditions, a mixture of different types of convection can be expected. Bottom heating and top cooling is reminiscent of Rayleigh–Bénard configurations with an imposed flux often investigated in the literature (e.g. Hurle et al. Reference Hurle, Jakeman, Pike and Sutton1967; Chapman & Proctor Reference Chapman and Proctor1980; Otero et al. Reference Otero, Wittenberg, Worthing and Doering2002; Verzicco & Sreenivasan Reference Verzicco and Sreenivasan2008; Johnston & Doering Reference Johnston and Doering2009; Fantuzzi Reference Fantuzzi2018). The lateral cooling is additionally expected to sustain natural convection, which has also been the subject of numerous studies (e.g. Batchelor Reference Batchelor1954; Churchill & Chu Reference Churchill and Chu1975; Bejan & Tien Reference Bejan and Tien1978; George & Capp Reference George and Capp1979;Wells & Worster Reference Wells and Worster2008; Ng et al. Reference Ng, Ooi, Lohse and Chung2015; Shishkina Reference Shishkina2016). In integral SA codes, like the ASTEC code developed by IRSN (Chatelard et al. Reference Chatelard2014), the entire process of a reactor core meltdown accident is simulated, from initiating events to the release of radioactive materials. Different modules address different aspects of the accident, the corium behaviour in the lower plenum of the vessel being one of them (Carénini, Fleurot & Fichot Reference Carénini, Fleurot and Fichot2014). In such a code, the focusing effect evaluation is based on a simplified approach proposed by Theofanous et al. (Reference Theofanous, Liu, Additon, Angelini, Kymäläinen and Salmassi1997), which combines correlations from both Rayleigh–Bénard and natural convection. It is assumed that the fluid in the bulk is thoroughly mixed and that the vertical heat transfer is symmetrical, meaning that the temperature difference between the bottom and the bulk is equal to the temperature difference between the bulk and the top. The validity of this approach was checked, in particular with the BALI-metal facility (Bonnet & Seiler Reference Bonnet and Seiler1999). Water was used to simulate the corium, and the top boundary condition was controlled by conduction through an epoxy plate and a temperature-regulated heat exchanger. This set-up was designed to closely resemble the conditions in a reactor with radiative heat transfer at the top. The BALI-metal tests have shown that, for a shallow layer thickness (aspect ratio above 10), the 0D model overestimates the side heat flux. However, CFD simulations of the metal layer (Shams et al. Reference Shams2020) showed that fluid properties (water versus steel) have a significant effect on the global behaviour, especially the Prandtl number. With steel, the lateral heat flux is up to 
$50\,\%$ higher compared with water under similar conditions. This questions the validity of using water as a simulant for molten steel, and consequently, the previously derived model for integral SA codes.
More generally, the competition between forced convection involving dominant vertical heat fluxes and horizontal or natural convection involving dominant horizontal heat fluxes is at the core of many geophysical situations. The competition between Rayleigh–Bénard and horizontal modes of convection is important for the dynamics within subglacial lakes in Greenland and Antarctica (Couston, Nandaha & Favier Reference Couston, Nandaha and Favier2022; Livingstone et al. Reference Livingstone2022). Planetary oceans are another example, since they receive latitudinally dependent solar radiations while being heated from the bottom by the geothermal flux (Wang et al. Reference Wang, Huang, Zhou and Xia2016). Finally, heterogeneous heat fluxes along the core-mantle boundary in the Earth's core, which are due to large-scale convective patterns within the solid mantle, can sustain large-scale azimuthal flows (Sumita & Olson Reference Sumita and Olson1999; Mound & DaviesReference Mound and Davies2017).
This paper presents three-dimensional (3-D) direct numerical simulations (DNS) of a liquid layer with a fixed Prandtl number of 
$0.1$ motivated by nuclear safety issues involving liquid metals (Carénini et al. Reference Carénini, Fichot and Seignour2018), with the aim of (i) getting a better understanding of the fluid behaviour and heat transfers for different characteristics found in reactors, (ii) establishing scaling laws, and (iii) deriving a one-dimensional (1-D) model suitable for use in integral SA codes.
This paper is divided into four sections. Section 2 introduces the governing equations and the numerical simulation tool used. Section 3 focuses on identifying the heat transfer mechanisms for the asymptotic regimes (dominant side or top heat flux) by analysing scaling laws for mean temperature variables. We conclude that a minima, a 1-D radial temperature description of the turbulent regime is necessary for nuclear safety evaluations. Hence, we delve into the fluid flow structure to determine this temperature profile. Finally, § 4 discusses the implications of our results for nuclear safety evaluations and outlines future developments of our work.
2. Mathematical and numerical formulation
2.1. Governing equations
We consider the flow of an incompressible fluid with buoyancy effects being included using the Boussinesq approximation. The fluid is confined within a cylinder of thickness 
$H$ and radius 
$R$ (see figure 1b) and gravity is pointing downward 
$\boldsymbol {g}=-g\boldsymbol {e}_z$. It is heated from below with a uniform heat flux per unit area 
$\phi _{in}$ and cooled from above by a uniform outgoing flux 
$\phi _{out}$. Note that we are interested in the cases where 
$\phi _{in}\neq \phi _{out}$ so that the residual heat flux is necessarily escaping the domain through the side boundary. The dimensional temperature on the side boundary is fixed at 
$\theta _0$. We model the bottom interface between the oxide layer and the liquid metal layer by a no-slip rigid boundary (a rigid crust forms at the oxide surface due to cooling, Carénini et al. Reference Carénini, Fichot and Seignour2018); we model the upper free surface of the liquid layer by a rigid stress-free boundary, neglecting free surface deformations. The side boundary is a rigid no-slip boundary. Lengths are rescaled using the height of the cylinder 
$H$ while time is rescaled using the vertical diffusive time scale 
$H^2/\kappa$, with 
$\kappa$ the constant thermal diffusivity. The dimensionless temperature 
$T$ is defined relatively to the imposed side temperature and rescaled using the imposed bottom flux 
$\phi _{in}$,
\begin{equation} T=\frac{k}{\phi_{in}H}\left(\theta-\theta_0\right) , \end{equation}
where 
$k$ is the thermal conductivity. The dimensionless conservation equations of momentum, mass and energy are then
$$\begin{gather} \frac{1}{Pr}\left(\frac{\partial\boldsymbol{u}}{\partial t}+ \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right) ={-}\boldsymbol{\nabla}P + Ra_{\phi} T\boldsymbol{e}_z+{\nabla}^2\boldsymbol{u} , \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0 , \end{gather}$$
$$\begin{gather}\frac{\partial T}{\partial t}+ \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla} T=\nabla^2 T . \end{gather}$$
Here 
$\boldsymbol {u}$, 
$P$ and 
$T$ are the dimensionless velocity, pressure and temperature of the fluid, respectively. The problem is characterised by four dimensionless parameters: the aspect ratio 
$\varGamma$, the flux ratio 
$R_F$, the Rayleigh–Roberts number 
$Ra_{\phi }$ that is based on the heat flux imposed at the bottom 
$\phi _{in}$ and the Prandtl number 
$Pr$ fixed to 
$0.1$ throughout the paper. They are defined by
\begin{equation} \varGamma = \frac{R}{H},\quad R_F=\frac{\phi_{out}}{\phi_{in}},\quad Ra_{\phi}=\frac{\beta g \phi_{in} H^4}{k \nu \kappa},\quad Pr=\frac{\nu}{\kappa}=0.1 , \end{equation}
where 
$\beta$ is the thermal expansion coefficient and 
$\nu$ is the kinematic viscosity, both assumed to be constant. The dimensionless boundary conditions can be written as
\begin{equation} \left.\begin{gathered} \boldsymbol{u}(z=0)=\boldsymbol{0}\quad\mbox{and}\quad\frac{\partial T}{\partial z}(z=0)={-}1 ,\\ \frac{\partial u}{\partial z}(z=1)=\frac{\partial v}{\partial z}(z=1)=w(z=1)=0\quad\mbox{and}\quad \frac{\partial T}{\partial z}(z=1)={-}R_F ,\\ \boldsymbol{u}(r=\varGamma)=\boldsymbol{0}\quad\mbox{and}\quad T(r=\varGamma)=0 , \end{gathered}\right\} \end{equation}
where 
$\boldsymbol {u}=(u,v,w)$ are the velocity components in cylindrical coordinates 
$\boldsymbol {e}_r,\boldsymbol {e}_{\varphi },\boldsymbol {e}_z$, respectively.
2.2. Numerical approach
2.2.1. Nek5000
The governing equations (2.2)–(2.4) with boundary conditions (2.6) are solved numerically using Nek5000 (Fischer Reference Fischer1997; Deville, Fischer & Mund Reference Deville, Fischer and Mund2002), which has been used extensively in thermal convection studies (e.g. Scheel, Emran & Schumacher Reference Scheel, Emran and Schumacher2013; Léard et al. Reference Léard, Favier, Le Gal and Le Bars2020; Terrien, Favier & Knobloch Reference Terrien, Favier and Knobloch2023). The entire cylindrical geometry is discretised using up to 
$\mathcal {E}=36\,608$ hexahedral elements that have been refined close to all boundaries to properly resolve viscous and thermal boundary layers. The velocity is discretised within each element using Lagrange polynomial interpolants based on tensor-product arrays of Gauss–Lobatto–Legendre quadrature points. The polynomial order 
$N$ on each element varies between 
$6$ and 
$10$ in this study. We use the 
$3/2$ rule for dealiasing with extended dealiased polynomial order 
$3N/2$ to compute nonlinear products. A third-order time stepping using a mixed explicit–implicit backward difference approach is used. A summary of the simulations physical and numerical parameters is provided in table 1 in the Appendix.
2.2.2. Numerical protocol and statistics
We initialise all simulations with a fluid at rest and a uniform temperature field 
$T=0$ everywhere. Infinitesimal temperature perturbations of amplitude 
$10^{-3}$ are introduced. Thermal convection grows during a transient that typically lasts for approximately five vertical diffusive times, and which is longer as the aspect ratio 
$\varGamma$ increases. Once the system has reached a statistically stationary state, various spatio-temporal averages are computed. Note that we have tested that different initial conditions (for example, starting with the equilibrium diffusive temperature distribution) eventually lead to the same statistically stationary state.
We first define the temporal and volume average operator 
$\langle.\rangle$ over the whole fluid domain volume 
$V$ and over time 
$\tau$ as
\begin{equation} \left\langle T\right\rangle=\frac{1}{\tau V}\int_{t_0}^{t_0+\tau}\int_{V}T\,\mathrm{d}V\,\mathrm{d}t . \end{equation}
The typical 
$\tau$ value ranges between 
$2$ and 
$0.2$ diffusive times for the lowest and the largest Rayleigh–Roberts numbers, respectively. Additionally, adding specific variables as a subscript means that an average along those specific directions is made. We always consider temporal averages during the statistically steady state so that the time variable is never explicitly written. For instance, 
$\langle.\rangle _{\varphi }$ indicates an average in time and along the azimuthal direction only.
2.2.3. Filtered simulations
While most of the results discussed below are obtained using DNS, some extreme cases were only accessible via filtered simulations following the approach described in Fischer & Mullen (Reference Fischer and Mullen2001). To distinguish between DNS and filtered simulations, a viscous dissipation criterion has been used. The mean viscous dissipation rate 
$\epsilon$, is defined by
\begin{equation}
\epsilon=2Pr\left\langle \boldsymbol{\mathsf{S}}:\boldsymbol{\mathsf{S}}\right\rangle\quad\mbox{with}\
\boldsymbol{\mathsf{S}}=\tfrac{1}{2}\left(\boldsymbol{\nabla}
\boldsymbol{u} + \boldsymbol{\nabla} \boldsymbol{u}^{\rm
T}\right). \end{equation}
A simulation with polynomial order 
$N$ is considered to be a DNS when the time and volume-averaged viscous dissipation 
$\epsilon$ is varying by less than 
$5\,\%$ when compared with the same simulation but using 
$N+2$ polynomial order. A simulation failing to satisfy this criteria is labelled as filtered and numerical stability is ensured by using a 
$1\,\%$ filter on the last two polynomials (Fischer & Mullen Reference Fischer and Mullen2001). Note that in the system studied, the Prandtl number is 
$0.1$, meaning that the filter does not impact the more diffusive temperature field but only the viscous dissipation at small scales.
Alternatively, we also used the criterion discussed in Scheel et al. (Reference Scheel, Emran and Schumacher2013) that compares the isotropic Kolmogorov dissipative scale with the numerical grid size. For all the DNS simulations presented in this study, the numerical grid size (
$L$) is below the Kolmogorov dissipative scale (
$\eta _K$) (see table 1 in Appendix).
3. Results
3.1. Qualitative overview
Let us start with a qualitative description of the different flow regimes. First, irrespective of the control parameters, no motion-less steady state exists in this system. Maintaining a constant temperature at the side generates a radial temperature gradient, which cannot be balanced by the hydrostatic pressure gradient. This leads to natural convection in the form of a downward recirculation along the side boundary and then towards the centre along the bottom boundary. At low 
$Ra_{\phi }$, this flow is axisymmetric, while at larger 
$Ra_{\phi }$, an instability develops and breaks the symmetry, leading to a drifting thermal branches pattern associated to the most extreme heat fluxes and temperatures found in the system. While these observations call for a linear stability analysis to identify the symmetry-breaking mechanism, we leave this aspect to future works. The focus of the present work is to identify scaling laws in the turbulent regime at large Rayleigh–Roberts numbers 
$Ra_{\phi }$, irrespective of the underlying linear instability mechanism.
We first present 3-D visualisations showing both the temperature field at the surface and streamlines coloured with the velocity amplitude in figure 2. We focus on the representative case 
$Ra_{\phi }=10^7$ and compare two flux ratios, 
$R_F=0.1$ and 
$0.9$, and three aspect ratios, 
$\varGamma =4$, 
$8$ and 
$16$. At low flux ratio (see figure 2a–c), a large-scale temperature pattern is clearly visible, with a number of azimuthal branches increasing with 
$\varGamma$. This thermal pattern is associated with intense radially outward flows. As we will see, these regimes are dominated by convective motions reminiscent of horizontal convection (Hughes & Griffiths Reference Hughes and Griffiths2008). At high flux ratio (see figure 2d–f), the system is more azimutally symmetrical and more intense fluctuations cover most of the domain. We also observe a clear temperature gradient between the core of the cylinder and the isothermal side boundary. As we will see, these regimes are dominated by thermal structures reminiscent of Rayleigh–Bénard convection (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000).
Figure 2. Each of the six subplots shows the 3-D view of (i) on the left half, the temperature field at the top surface and (ii) on the right half, some velocity streamlines. Plots (a–c) show results for 
$R_F=0.1$ and (d–f) for 
$R_F=0.9$, with 
$\varGamma =4$ for (a,d), 
$\varGamma =8$ for (b,e) and 
$\varGamma =16$ for (c, f). In all cases, 
$Pr=0.1$ and 
$Ra_{\phi }=10^7$. Note that the temperature colour scale is the same for all plots.
3.2. Low flux ratio regime
This section is devoted to the low flux ratio regime for which we fix 
$R_F=0.1$ as a representative value. We first discuss the mean temperature scaling as a function of the two other input parameters 
$Ra_{\phi }$ and 
$\varGamma$ before providing a theoretical explanation based on simple dimensional arguments.
3.2.1. Scaling for 
$R_F=0.1$
In this section we fix 
$R_F=0.1$, meaning that 
$90\,\%$ of the power transferred through the lower boundary goes out through the side, and we seek scaling laws for the mean temperature systematically varying 
$Ra_{\phi }$ and 
$\varGamma$. For each simulation, the statistically stationary state is reached, which typically takes 
$5$ diffusive time scales, and we compute the mean temperature of the system noted 
$\langle T\rangle$ using the space–time average operator defined in (2.7).
The mean temperature evolution with 
$Ra_{\phi }$ and 
$\varGamma$ is reported in figure 3. When the value of 
$Ra_{\phi }$ increases, the mean non-dimensional temperature of the system decreases as expected (see figure 3a). This is because a higher Rayleigh–Roberts number results in more efficient heat transfers within the system, leading to an average temperature getting closer to the side temperature that is zero in our dimensionless units (2.1). Note however that the dimensional temperature obviously increases when we increase the heat flux. When 
$Ra_{\phi }$ exceeds 
$10^5$, a power-law behaviour emerges with an exponent that appears to be unaffected by 
$\varGamma$. Upon estimating the power exponent (best fit using the least square method), it has been found that 
$\langle T\rangle \sim Ra_{\phi }^{-0.20\pm 0.03}$. To determine the mean exponent, we compute the average of the exponents associated with 
$\varGamma =4,5,8,16$ considering 
$Ra_{\phi } \geqslant 10^5$ data only, while the variability is quantified by the largest difference between these exponents. This scaling law is derived by considering DNS and filtered simulation data. The exclusion of the filtered data from the analysis only leads to a slight alteration in the scaling law, resulting in 
$\langle T\rangle \sim Ra_{\phi }^{-0.21\pm 0.03}$.
Figure 3. Log-log plot of the mean temperature evolution with Rayleigh number for different aspect ratios in (a) and with the aspect ratio for different Rayleigh numbers in (b). The width of the thick light lines is equal to 
$3$ times the standard deviation of the mean temperature time series at the statistically stationary state. The dash-dotted lines show the best power-law fit for each 
$\varGamma$ in (a) and 
$Ra_{\phi }$ in (b). Empty/full symbols respectively indicate filtered/DNS simulations. Here 
$Pr=0.1$ and 
$R_F=0.1$.
Additionally, when 
$\varGamma$ increases at constant 
$Ra_{\phi }$, the mean temperature of the system increases. Indeed, when 
$\varGamma$ increases for a fixed 
$Ra_{\phi }$, the ratio between the heating bottom surface and the cooling lateral surface also increases, leading to a larger global energy input into the system. Figure 3(b) shows that when 
$Ra_{\phi }$ is greater than 
$10^5$, a power-law behaviour in 
$\varGamma$ also emerges. Considering both DNS and filtered simulation data, and employing a consistent methodology for exponent estimation as applied to the 
$Ra_{\phi }$ dependence, we obtain 
$\langle T\rangle \sim \varGamma ^{0.83 \pm 0.11}$ whereas we find that 
$\langle T\rangle \sim \varGamma ^{0.89 \pm 0.08}$ without taking filtered data into account. Once again these results indicate a minor influence of the filtered simulations on the overall scaling behaviour.
The two scalings can be combined leading to the following final power law:
\begin{equation} \left\langle T\right\rangle \sim Ra_{\phi}^{{-}1/5}\varGamma^{4/5} . \end{equation}Here the particular exponent values will be theoretically justified below (see § 3.2.2) and fall well into the fitted range of exponents found from our simulations.
Figure 4 shows the mean temperature, compensated by scaling (3.1) as a function of 
$Ra_{\phi }$. All rescaled data converge to the same constant close to unity. For 
$\varGamma =4$, the scaling law seems to apply when 
$Ra_{\phi }>10^5$, whilst it is necessary to wait until 
$Ra_{\phi }=10^8$ when 
$\varGamma =16$. In the following section, we use dimensional analysis to gain insight into the physics underlying this scaling.
Figure 4. Compensated mean temperature following (3.1) as a function of 
$Ra_{\phi }$. Different symbols/colours correspond to different aspect ratios. The coloured areas indicate 
$3$ times standard deviations of the mean temperature time series at the statistically stationary state multiplied by 
$Ra_{\phi }^{1/5}\varGamma ^{-4/5}$. Empty/full symbols respectively indicate filtered/DNS simulations. Here 
$Pr=0.1$ and 
$R_F=0.1$.
3.2.2. Theoretical analysis for $R_F=0.1$
We focus on the sidewall considering the dimensionless steady axisymmetrical governing equations using cylindrical coordinates
$$\begin{gather} u\frac{\partial u}{\partial r}+ w\frac{\partial u}{\partial z}={-}Pr\frac{\partial P}{\partial r} +Pr\nabla^2 u , \end{gather}$$
$$\begin{gather}u\frac{\partial w}{\partial r} + w\frac{\partial w}{\partial z}={-}Pr\frac{\partial P}{\partial z} + Ra_{\phi}PrT+ Pr\nabla^2 w , \end{gather}$$
$$\begin{gather}\frac{1}{r}\frac{\partial ru}{\partial r}+\frac{\partial w}{\partial z}=0 , \end{gather}$$
$$\begin{gather}u\frac{\partial T}{\partial r} + w\frac{\partial T}{\partial z} = \nabla^2 T . \end{gather}$$
Due to the low Prandtl regime, the viscous boundary layer is nested within the thermal one. A diagram of this configuration is shown in figure 5. The behaviour of the vertical (
$w_s$) and radial (
$u_s$) velocities at the edge of the viscous boundary layer is derived from the conservation of mass (3.4) and energy (3.5). We assume that the typical scale of variation in the radial direction is the dimensionless thickness of the side boundary layer (
$\partial /\partial _r \sim 1/\delta$). This thickness may either represent the thermal boundary layer (
$\delta _{th}$) or the viscous boundary layer (
$\delta _v$), depending on whether the radial variation under consideration pertains to temperature or velocity. Moreover, we assume that the scale of variation in height is the dimensionless height equal to 
$1$ (
$\partial /\partial _z \sim 1$). Then mass conservation (3.4) leads to
\begin{equation} \frac{u_s}{\delta_v}\sim w_s. \end{equation}
Figure 5. Sketch of the side boundary layers when 
$Pr<1$.
On the left-hand side of (3.5), both advection terms, 
$u\partial _r T$ and 
$w\partial _z T$, scale as 
$w_s T$, due to mass conservation (3.6). Because radial variations are much larger than height variations near the lateral boundaries, we assume that the dominant term in the Laplacian operator scales with the dimensionless thickness of the side boundary layer squared (
$\nabla ^2 \sim 1/\delta ^2$). Therefore, (3.5) leads to
\begin{equation} w_s\sim 1/\delta_{th}^2. \end{equation}
Let us now approximate the temperature variation across the thermal boundary layer (
$\delta T_s$) by the average temperature of the system 
$\langle T\rangle$ (recall that the dimensionless temperature on the sidewall is zero). Within the thermal boundary layer, we expect a force balance between the inertia term and the buoyancy term, so that the vertical momentum balance (3.3) reduces to
\begin{equation} u \frac{\partial w}{\partial r} + w \frac{\partial w}{\partial z} \simeq RaPrT. \end{equation}
On the left-hand side of (3.8), each advection term, 
$u\partial _r w$ and 
$w\partial _z w$, scales like 
$w_s^2$, due to mass conservation (3.6). Then (3.8) reduces to 
$w_s^2 \sim Ra_{\phi }Pr\langle T\rangle$ or, equivalently, using (3.7),
\begin{equation} \delta_{th}\sim\left(Ra_{\phi}Pr\left\langle T\right\rangle\right)^{{-}1/4} . \end{equation}Finally, to link the heat flux applied on the bottom surface to the thermal characteristics of the side, a global flux balance is required. Integrating (3.5) over the volume at steady state leads to
\begin{equation} {\rm \pi}\varGamma^2\left(1-R_F\right)=2{\rm \pi}\varGamma\frac{\left\langle T\right\rangle}{\delta_{th}} , \end{equation}
where the left-hand side corresponds to the power mismatch between the lower and upper boundaries, which is balanced by the conducting flux across the thermal boundary layer on the right-hand side. Thus, the averaged temperature 
$\langle T\rangle$ is proportional to the aspect ratio and to the thickness of the thermal boundary layer
\begin{equation} \left\langle T\right\rangle\sim\delta_{th}(1-R_F)\varGamma . \end{equation}Combining (3.9) and (3.11), the mean temperature can therefore be expressed in terms of the control parameters through the relationship
\begin{equation} \left\langle T\right\rangle\sim Ra_{\phi}^{{-}1/5}\varGamma^{4/5}\left(1-R_F\right)^{4/5} Pr^{{-}1/5} . \end{equation} These simple dimensional arguments allow us to recover the scaling (3.1) obtained via numerical simulations. In addition, this reveals the dependencies on 
$R_F$ and 
$Pr$, which we did not observe since both these parameters have been fixed for now. Note that assuming that the average temperature is representative of the temperature difference across the radial thermal boundary layer is presumably only valid when the Rayleigh–Roberts number is large enough to mix efficiently the bulk of the convective system. In order to compare our results with existing literature, it is standard to use the Nusselt number (the ratio of convective to diffusive heat flux). However, this measurement is only meaningful when the heat transfer can be unambiguously defined, that is, when the average isothermal surfaces are parallel. Or in other words, when temperature varies only along one dimension in the system. In our system there is no specific direction for heat transfer except in the asymptotic regime of high or low flux ratio. In these scenarios, it is conceivable to determine the Nusselt value, which indicates the main heat flux direction (vertical for 
$R_F \rightarrow 1$ and horizontal for 
$R_F=0$).
In the low flux ratio regime, the heat flux mainly goes in the horizontal direction; therefore, we define the Nusselt number by the relationship
\begin{equation} \phi_{in} \equiv \frac{k}{R}\left(\left\langle\theta\right\rangle-\theta_0\right)Nu . \end{equation}With the choice made for the non-dimensional temperature, the Nusselt number reads as the inverse of the mean temperature, hence
\begin{equation} Nu\sim Ra_{\phi}^{1/5}\varGamma^{1/5}\left(1-R_F\right)^{{-}4/5} Pr^{1/5}. \end{equation}
Although we find a 
$1/5$ exponent for the Rayleigh dependency, reminiscent of the horizontal convection scaling of the Rossby (Reference Rossby1965) regime, it is important to note that the force balance in the boundary layer is different. In Rossby (Reference Rossby1965) buoyancy and viscosity at the bottom both play a role, while in our case, the balance is between inertia and buoyancy on the vertical side boundary. The equivalent Rayleigh exponent based on a temperature scale (instead of a flux scale as done here) is 
$1/4$ and corresponds to the scaling of vertical convection identified by Batchelor (Reference Batchelor1954). Furthermore, in a two-dimensional rectangular system with identical thermal boundary conditions to the present case (except at the top boundary where 
$R_F$ was set to 
$0$), Ganzarolli & Milanez (Reference Ganzarolli and Milanez1995) identified a Nusselt scaling law. We obtain identical exponents for the Rayleigh–Roberts number and the aspect ratio as those found by Ganzarolli & Milanez (Reference Ganzarolli and Milanez1995). Both exhibited 
$1/5$ exponents. It is worth mentioning that the validity of the Rayleigh scaling was recently shown to be limited to laminar boundary layers by Shishkina (Reference Shishkina2016), so that a different scaling is expected at even large Rayleigh–Roberts numbers (we considered 
$Ra_{\phi }\leqslant 10^9$). Previous works (George & Capp Reference George and Capp1979) also pointed out this possibility. Finally, our results show the same 
$1/5$ Prandtl exponent dependence predicted by Shishkina (Reference Shishkina2016) and recently emphasized byZwirner et al. (Reference Zwirner, Emran, Schindler, Singh, Eckert, Vogt and Shishkina2022). However, we recall here that we have not checked this Prandtl scaling with our numerical simulations, which were all performed at 
$Pr=0.1$.
To validate the underlying physical approximations required to derive our scaling law, we now examine the secondary variables, specifically the thickness of the thermal boundary layer and the vertical velocity. Thanks to the relations (3.7) and (3.11), the scaling laws for those quantities can be written as
\begin{equation} \left.\begin{gathered} \delta_{th}\sim Ra_{\phi}^{{-}1/5}\varGamma^{{-}1/5}\left(1-R_F\right)^{{-}1/5} Pr^{{-}1/5} ,\\ w_s\sim Ra_{\phi}^{2/5}\varGamma^{2/5}\left(1-R_F\right)^{2/5} Pr^{2/5} . \end{gathered}\right\} \end{equation} The thickness of the thermal boundary layer is estimated by determining the radius at which 
$90\,\%$ of the maximum temperature near the edge is reached, while the vertical velocity is based on seeking the minimum vertical velocity near the edge. More precisely, we first compute the vertical velocity at mid-height of the domain to avoid the influence of the top and bottom boundaries, and we average in time and in the azimuthal direction. We then search for the first peak of negative vertical velocity, starting from the boundary and moving toward the centre of the domain. We conducted these measurements for aspect ratios of 
$4\leqslant \varGamma \leqslant 16$ and for 
$10^5\leqslant Ra_{\phi }\leqslant 10^8$. Results shown in figure 6 are in excellent agreement with the scaling laws (3.15), hence, further validating our approach.
Figure 6. Log-log plot of (a) the size of the thermal boundary layer and (b) the absolute value of the vertical peak velocity near the edge (defined as the first minimum of the negative vertical velocity at mid-height, starting from the boundary and moving toward the centre of the domain) versus the determined scaling laws (3.15) for different 
$Ra_{\phi }$ and aspect ratios. The dash-dotted lines correspond to the scaling laws and for all the simulations. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameters are 
$Pr=0.1$ and 
$R_F=0.1$.
3.3. High flux ratio regime
We now consider the other limiting case of a flux ratio close to unity. We follow the same approach as in the previous section and we start with scalings obtained from numerical simulations followed by a theoretical explanation.
3.3.1. Scaling for $R_F=0.9$
We now fix 
$R_F=0.9$, meaning that 
$90\,\%$ of the heating power goes out through the top surface. Similarly to the precedent section, we seek power laws systematically varying 
$Ra_{\phi }$ and 
$\varGamma$. However, in this second regime closer to the classical Rayleigh–Bénard configuration, heat transfers are mainly along the vertical direction. We therefore focus on the mean temperature difference between the top and bottom surfaces, denoted 
$\Delta T_v$ and computed as
\begin{equation} \Delta T_v=\left\langle T(t,r,\varphi,z=0)-T(t,r,\varphi,z=1)\right\rangle_{r\varphi} , \end{equation}
where the average corresponds to a temporal and surfacic average along radius and azimuthal angle. As we will see below, this averaged quantity is more relevant than the mean temperature of the system that was more representative of the radial temperature differences between the bulk and the side boundary when 
$R_F=0.1$. The evolution of the mean temperature difference with 
$Ra_{\phi }$ and 
$\varGamma$ is plotted in figure 7. We observe that, as 
$Ra_{\phi }$ increases, 
$\Delta T_v$ decreases, suggesting that the system becomes more homogeneous vertically. The decrease in the top–bottom temperature difference appears to be independent of the aspect ratio, which suggests a local mechanism. Note also that this vertical temperature difference is the same at different locations, as will be further discussed below in § 3.5. When 
$Ra_{\phi }$ is larger than 
$10^5$, a power-law scaling emerges with an exponent 
$Ra_{\phi }^{-0.20\pm 0.02}$ independently of 
$\varGamma$. As we will show below, the closest relevant scaling is given by
\begin{equation} \Delta T_v\sim Ra_{\phi}^{{-}1/5} . \end{equation}
Figure 7. Log-log plot of the bottom–top temperature difference for different 
$\varGamma$ values. The error bars are computed by taking 3 times the standard deviation of the bottom–top temperature difference time series at the statistically stationary state. The dash-dotted line corresponds to the best-fit scaling law 
$Ra_{\phi }^{-1/5}$, the blue diamond (
$\Diamond$) plot shows the horizontally homogeneous Chapman & Proctor (Reference Chapman and Proctor1980) configuration (
$R_F=1$ in a doubly periodic Cartesian box) and the continuous black line corresponds to the Grossmann & Lohse (Reference Grossmann and Lohse2000) 
$I_l$ regime written in terms of flux Rayleigh number 
$Ra_{\phi }$. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameters are 
$Pr=0.1$ and 
$R_F=0.9$.
3.3.2. Theoretical analysis for $R_F=0.9$
In this configuration close to Rayleigh–Bénard convection, one might initially assume that the system is controlled by the heat flux across a thin thermal boundary layer (Malkus Reference Malkus1954). Let us define the Rayleigh number based on the temperature difference 
$Ra_{\Delta T_v}=(\beta g \Delta T_v H^3)/(\nu \kappa )$, which is an output parameter in our case since the temperature difference is not known a priori. Note that in the asymptotic limit where 
$R_F$ tends to one, 
$Ra_{\Delta T_v}$ and 
$Ra_{\phi }$ are getting proportional to each other with the Nusselt number (as defined in Rayleigh–Bénard convection) being the proportionality coefficient (Cioni, Ciliberto & Sommeria Reference Cioni, Ciliberto and Sommeria1997). The classical scaling 
$\Delta T_v\sim Ra_{\Delta T_v}^{-1/3}$ would equivalently give 
$Ra_\phi ^{-1/4}$, which is not compatible with our result (3.17). Our regime is closer to the regime 
$I_l$ predicted by Grossmann & Lohse (Reference Grossmann and Lohse2000) in which an energetic approach is used to estimate the viscous and thermal dissipation rates for the 
$Pr<1$ case. In the 
$I_l$ regime, dissipation rates are dominated by their boundary layer contributions, therefore, the 
$I_l$ regime is obtained at relatively low 
$Ra_{\Delta T_v}$, i.e. when the turbulence is sufficiently underdeveloped, so that the dissipation is mostly concentrated within the boundary layers.
Considering a homogeneous bulk (reached at sufficiently high Rayleigh–Roberts numbers) and thanks to the 
$I_l$ regime of Grossmann & Lohse (Reference Grossmann and Lohse2000), the Nusselt scaling 
$Nu\sim \delta _b^{-1}\sim Ra_{\Delta T_v}^{1/4}$, where 
$\delta _b$ is the thickness of the bottom thermal boundary layer, leads to 
$\Delta T_v \sim Ra_{\phi }^{-1/5}$ consistent with (3.17). It is noteworthy that analogy with the 
$I_l$ regime only makes sense if one assumes an equivalence between the upper and lower thermal boundary layers, which is not necessarily the case here since the velocity boundary conditions are mixed (free slip at the top, no slip at the bottom). In addition, all the results of the study of Grossmann & Lohse (Reference Grossmann and Lohse2000) are obtained with an imposed temperature difference contrary to the imposed heat flux here.
When 
$R_F$ tends to 
$1$, our configuration is close to the one studied by Chapman & Proctor (Reference Chapman and Proctor1980), where a constant flux is imposed at the bottom surface and goes out by the top surface (
$R_F=1$) in a horizontally infinite domain (
$\varGamma = \infty$). This limiting case can be explored with a Cartesian box of size 
$2\times 2\times 1$ periodic in both horizontal directions, with the same imposed flux at the top and bottom and with no-slip and free-slip conditions, respectively, at the bottom and top. Figure 7 shows 
$\Delta T_v$ measured after the statistically stationary state is obtained. We also plot the temperature difference 
$\Delta T_v$ predicted by the 
$I_l$ regime scaling, keeping the prefactor determined for rigid boundaries and imposed temperatures (Grossmann & Lohse Reference Grossmann and Lohse2000). We find good agreement between the local Cartesian set-up and the 
$I_l$ regime, with the same Rayleigh exponent being found. A slight difference in the prefactor is nevertheless noted. Modifying the thermal and velocity boundary conditions did not have a significant impact, with only a minor change in the prefactor. Our 
$R_F=0.9$ case is also found to be in good agreement with the 
$I_l$ regime, but with a small difference. Here 
$10\,\%$ of the flux exits through the side, leading to a disruption in the vertical temperature gradient and, therefore, a difference in the temperature difference from bottom to top caused by the baroclinic flow.
Another way to verify the relevance of the 
$I_l$ regime is to check the bottom thermal boundary layer behaviour. The scaling law predicts that 
$\delta _{b} \sim Ra_{\phi }^{-1/5}$. A measure of the bottom thermal boundary layers has been done based on the temperature variance computed as
\begin{equation} \sigma_T=\left\langle\left(T(t,r=\varGamma/2,\varphi,z)-\left\langle T(t,r=\varGamma/2,\varphi,z)\right\rangle_{\varphi}\right)^2\right\rangle_{\varphi}^{1/2} . \end{equation}
We focus here on the temperature at 
$r=\varGamma /2$ to avoid the side and centre areas that are more significantly affected by the overall baroclinic circulation driven at the sidewall. In figure 8(a), 
$\sigma _T$ as a function of the altitude 
$z$ is plotted for different 
$Ra_{\phi }$ and for 
$\varGamma =16$. Near the top and bottom boundaries, the temperature variance is rapidly varying, which indicates the existence of boundary layers. In the bulk 
$0.2< z<0.8$, the system is more homogeneous especially when 
$Ra_{\phi }$ is high. Vertical dotted lines indicate the thickness of the bottom boundary layer 
$\delta _b$. It is estimated as the position 
$z$ for which the temperature variance has decreased by 
$90\,\%$ compared with its value at the boundary. We used the temperature variance rather than the vertical temperature profile because the recirculation flow (induced by the cold side) involves a heat transport mechanism similar to the horizontal convection (Mullarney, Griffiths & Hughes Reference Mullarney, Griffiths and Hughes2004) altering the thermal boundary layer structure. In figure 8(b) the corresponding boundary layer thickness is plotted as a function of 
$Ra_{\phi }$ for several 
$\varGamma$ at 
$R_F=0.9$. It is independent of 
$\varGamma$ and scales like 
$\delta _b \sim Ra_{\phi }^{-0.20\pm 0.05}$ (determined by the least square method). The 
$I_l$ regime law (Grossmann & Lohse Reference Grossmann and Lohse2000) is also plotted and is consistent with our measurements at 
$R_F=0.9$ in terms of exponent, with again an offset on the prefactor presumably due to residual effects of the large-scale circulation driven at the sidewall.
Figure 8. (a) Plot of the temperature variance with the altitude 
$z$ for different 
$Ra_{\phi }$ at 
$\varGamma =16$, 
$Pr=0.1$ and 
$R_F=0.9$. The irregular profile for 
$Ra_{\phi }=10^8$ is due to a lack of statistical samples for this very costly computation. The vertical dotted lines indicate the boundary layer altitudes based on the 
$90\,\%$ decrease of the temperature variance near the bottom. (b) Log-log plot of the bottom boundary layer with 
$Ra_{\phi }$ for different 
$\varGamma$. The continuous black line represents the Grossmann & Lohse (Reference Grossmann and Lohse2000) scaling corresponding to the 
$I_l$ regime. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameters are 
$Pr=0.1$ and 
$R_F=0.9$.
3.4. Intermediate regimes
In §§ 3.2 and 3.3 we have identified two different temperature averages that seem to characterise the system behaviour in the low/high flux ratio regimes. We now investigate the case 
$R_F=0.5$ where the heat flux equally goes out through the top and the side. As before, we seek power laws systematically varying 
$Ra_{\phi }$ and 
$\varGamma$. In the precedent limiting regimes, heat transfers were dominantly radial or vertical, but the outgoing flux is now evenly distributed between the top and side boundaries, so that a mixed state is expected. Thus, we now carry out an analysis using both the mean temperature 
$\langle T\rangle$ and the top/bottom temperature difference 
$\Delta T_v$.
Figure 9 shows plots of both the mean temperature 
$\langle T\rangle$ and the bottom–top temperature difference 
$\Delta T_v$ compensated by the low/high flux ratio scaling law, respectively, for a fixed value of 
$\varGamma =8$ and for the three flux ratios 
$R_F\in [0.1,0.5,0.9]$. The best-fit power laws for 
$\langle T\rangle$ and 
$\Delta T_v$ at the intermediate flux ratio of 
$R_F=0.5$ do not correspond to the scaling laws observed in previous regimes, exhibiting dependencies of 
$Ra_{\phi }^{-0.23\pm 0.02}$ and 
$Ra_{\phi }^{-0.13\pm 0.05}$, respectively. To clarify these observations, let us define the averaging operator
\begin{equation} \left\langle T\right\rangle(r_1< r< r_2)=\frac{1}{\tau {\rm \pi}(r_2^2-r_1^2)}\int_{t_0}^{t_0+\tau}\int_0^1\int_0^{2{\rm \pi}}\int_{r_1}^{r_2}T\,r\,\textrm{d}r\,\textrm{d}\varphi \,\textrm{d}z\,\textrm{d}t, \end{equation}
representing a temporal and volume average restricted to a particular radius range 
$r_1< r< r_2$. A similar definition without vertical integration is used for the top–bottom temperature difference. In figure 10 we show such local averages for various increasing values of the limiting radii 
$r_1$ and 
$r_2$. We fix 
$R_F=0.1$, 
$\varGamma =8$ and 
$10^3< Ra_{\phi }<10^9$. The mean temperature exhibits different power laws with varying exponents in different radial regions. Close to the side boundary (
$7.2< r<\varGamma =8$ typically), the local average 
$\langle T\rangle$ follows a 
$Ra_{\phi }^{-1/5}$ scaling, corresponding the low flux ratio regime as expected from the proximity with the sidewall circulation. In the bulk however (
$0< r<0.8$ typically), an exponent of 
$-0.27$ is obtained, reminiscent of the scaling observed for the high flux ratio regime. This implies that even at a flux ratio 
$R_F=0.1$, the high flux ratio mechanism persists in some form close to the centre. As shown in figure 10(b), 
$\Delta T_v(r_1< r< r_2)$ also exhibits a regime transition. For 
$r>4$, 
$\Delta T_v$ becomes negative at low 
$Ra_{\phi }$ and increasingly negative towards the edge, highlighting the sidewall impact on the dynamics. In the centre region however (
$0< r<0.8$), the power law for 
$\Delta T_v$ is consistent with the 
$I_l$ Grossmann & Lohse (Reference Grossmann and Lohse2000) regime, following a 
$Ra_{\phi }^{-1/5}$ dependency, as expected from the quasi-homogeneous behaviour of the bulk convection far from the sidewall.
Figure 9. (a) Compensated mean temperature, (b) compensated bottom–top temperature difference as a function of 
$Ra_{\phi }$ for 
$R_F=0.1\ (\circ )$, 
$R_F=0.5\ (\square )$, 
$R_F=0.9\ (\Diamond )$. The error bars are computed by taking 3 times the standard deviation of the time series at the statistically stationary state multiplied by 
$Ra_{\phi }^{1/5}\varGamma ^{-4/5}$ and 
$Ra_{\phi }^{1/5}$, respectively. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameters are 
$Pr=0.1$ and 
$\varGamma =8$.
Figure 10. (a) Mean temperature 
$\langle T\rangle$ and (b) bottom–top temperature difference 
$\Delta T_v$, defined on successive rings with 
$Ra_{\phi }$. The step radius is 
$0.8$, starting from 
$0$ to 
$0.8$ and ending from 
$7.2$ to 8. For (a), the power-law exponent (slope) is measured for each ring. The inset in (b) shows a log-log plot of 
$\Delta T_v$ and associated power-law exponents for the most inner rings. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameters are 
$Pr=0.1$, 
$\varGamma =8$ and 
$R_F=0.1$.
The system is therefore inherently inhomogeneous with a permanent interaction between two limiting regimes: one in the bulk defined by the 
$I_l$ Grossmann & Lohse (Reference Grossmann and Lohse2000) regime, and the other one near the edge defined by vertical natural convection. In the asymptotic regimes, the limiting mechanism (side and bottom–top boundary layers, respectively, for 
$R_F=0.1$ and 
$0.9$) controls the overall heat transfers and thermal structure. The variation of the flux ratio parameter involves a continuous transition with a gradual shift from the dominance of one regime to another. But for any flux ratio, signatures of both regimes can be seen in the radial profiles, as will be studied below. The regime interaction artificially arises from describing an inhomogeneous system using global variables (
$\langle T\rangle, \Delta T_v$). In fact, these two mechanisms operate simultaneously but in different locations. The high flux ratio regime is primarily located near the centre of the domain, while the low flux ratio regime is predominantly located at the side. Our study hence shows that the global mean temperature values 
$\langle T\rangle$ and 
$\Delta T_v$ are not adequate to fully characterise the system at any flux ratio due to the presence of two distinct mechanisms, and the radial inhomogeneities in the system statistics, as discussed in the next section.
3.5. Radial inhomogeneities
In this section we study the radial structure of the two temperature variables examined previously. In figure 11 we present the spatio-temporal averages 
$\langle T\rangle _{\varphi z}$ and 
$\langle \Delta T_v\rangle _{\varphi }$ as a function of radius, for different flux ratio values with 
$Ra_{\phi }=10^8$ and 
$\varGamma =8$. For 
$R_F=0.9$, 
$\Delta T_v$ is constant on a large part of the domain except near the sidewall where an inhomogeneity is clearly visible. The lower the 
$R_F$ the larger this inhomogeneous domain, as illustrated in figure 11(a). Indeed, as 
$R_F$ decreases, the sidewall circulation gets stronger and perturbs the bulk forced convection that tends to vertically homogenise the temperature. The radial evolution of the mean temperature (figure 11b) also reveals two distinct regions. One region is located near the lateral sidewall and displays a nearly constant temperature (excluding the thin boundary layer developing along the sidewall) in good agreement with the volume-averaged temperature at low flux ratio, shown as a thin horizontal line for 
$R_F=0.1$. The second, inner region shows a linear increase in temperature towards the centre. An increase in the prescribed heat flux at the top affects the radial temperature profile. Indeed, when 
$R_F=0.9$, the uniform temperature zone is getting pushed towards the sidewall, while the same temperature gradient is observed in the central region. Regarding the intermediate case 
$R_F=0.5$, the two distinct areas previously identified can still be seen, but the plateau zone is smaller compared with the case 
$R_F=0.1$. The size of this region is directly tied to the magnitude of the lateral heat flux and will be further studied in the next section by analysing the flow structure.
Figure 11. Radial profile of the bottom–top temperature difference (a) and of the mean temperature (b) for 
$R_F=0.1$ (
$\square$), 
$R_F=0.5$ (
$\Diamond$) and 
$R_F=0.9$ (
$\circ$). The horizontal lines indicate the radial-mean value for 
$R_F=0.9$ in (a) and 
$R_F=0.1$ in (b), while the dash lines indicate the 
$r_p$ value for 
$R_F=0.1$. The input parameters are 
$Ra_{\phi }=10^7$, 
$\varGamma =8$ and 
$Pr=0.1$.
In view of these two zones clearly identified in the radial temperature profile, we suggest to approximate it (outside the outer thermal boundary layer) by considering the three parameters model
\begin{equation} \left\langle T\right\rangle_{\varphi z}(r)=\left\{\begin{array}{@{}ll@{}} G_{T}r_p\left(1-\dfrac{r}{r_p}\right)+T_p & \mbox{for}\ r\leqslant r_p, \\[10pt] T_p & \mbox{for}\ r\geqslant r_p, \end{array} \right. \end{equation}
with the constant temperature near the sidewall 
$T_p$, the slope of the temperature profile in the linear region 
$G_T$ and the transition radius 
$r_p$. In the next sections we focus on the determination of each of these parameters and discuss their physical origin.
3.5.1. Radial temperature gradient
First, let us compute the slope of the radial temperature profile 
$G_T$ as 
$|\langle {\mathrm {d}\!\langle T\rangle _{\varphi z}}/{\mathrm {d}r}\rangle _{r} |$ and the radial velocity modulus computed as 
$\sqrt {\langle \langle u\rangle _{\varphi z}^2\rangle _r}$ restricted to 
$r\leqslant 5/8\varGamma$ and 
$r\leqslant 0.95\varGamma$, respectively, for 
$R_F=0.1$ and 
$0.9$, 
$10^5< Ra_{\phi }<10^8$ and considering 
$\varGamma =4,8$ and 
$16$. All data are plotted in figure 12. Both the mean temperature gradient and the radial velocity do not significantly depend on the flux ratio, nor on the aspect ratio, and they seem to be anti-correlated with each other. Indeed, they both follow a power-law behaviour with 
$Ra_{\phi }$, with the same approximate exponent 
$1/3$ but positive for the radial velocity and negative for the temperature gradient. This suggests that, to understand the physical origin of the thermal gradient, a closer look at the flow structure is necessary.
Figure 12. Log-log plot of the radial temperature gradient in (a) and of the RMS radial velocity in (b) as a function of 
$Ra_{\phi }$, with various values of 
$\varGamma$ shown with symbols and of 
$R_F$ shown with colours. The dash-dotted lines indicate the best-fit power law. The coloured areas are computed by taking the standard deviation of the radial temperature gradient and the square of the radial velocity in the bulk for (a,b), respectively. Empty/full symbols respectively indicate filtered/DNS simulations. The input parameter is 
$Pr=0.1$.
To do so, we consider the representative case 
$Ra_{\phi }=10^7$ and 
$\varGamma =8$. We compute in a 
$r,z$ plane and, for both regimes, 
$R_F=0.1$ and 
$R_F=0.9$, the norm of the velocity field denoted 
${U_n}$, shown in the top panels of figure 13, as well as the velocity fluctuations field denoted 
${u_n}^\prime$, shown at the bottom. Those fields are computed as follows:
\begin{align}
{U_n}(r,z)=\big[\langle u\rangle_{\varphi}^2+\langle
w\rangle_{\varphi}^2\big]^{1/2},\quad {u_n}^\prime(r,z)=
\big[\langle(u-\langle
u\rangle_{\varphi})^2\rangle_{\varphi}+ \langle(w-\langle
w\rangle_{\varphi})^2\rangle_{\varphi}\big]^{1/2}
.
\end{align}
In both regimes a global circulation surrounds the whole domain, with intense mean velocities close to the boundaries. In the bulk, fluctuations dominate the flow, at least for 
$0< r<5$ at 
$R_F=0.1$ and up to the side boundary layer for 
$R_F=0.9$. These regions correspond to the domain with a significant radial gradient of the mean temperature, shown in figure 11. In order to understand the emergence of the radial temperature gradient, let us look at the heat transport equation. We average it in time, in azimuth, but also over a particular thickness 
$h$. This thickness corresponds to the altitude at which the vertical profile of the radial velocity 
$\langle u\rangle _{\varphi r}$ changes sign. We denote this averaging operation as 
$\langle \bullet \rangle _{\varphi h}$. To clarify, we do not average over the whole depth because the terms related to radial advection would vanish owing to continuity.
Figure 13. Map in the 
$r,z$ plane of the mean (a,b) and fluctuating (c,d) velocity norm for (a,c) 
$R_F = 0.1$ and (b,d) 
$R_F = 0.9$. The dashed lines indicate the 
$r_p$ value. The input parameters are 
$Pr = 0.1$, 
$\varGamma = 8$ and 
$Ra_{\phi }=10^7$.
Our averaging procedure allows us to distinguish between the mean flow carrying hot fluid from the centre to the edge of the domain in the upper region 
$z>h$ and cold fluid advected from the edge to the centre in the bottom region 
$z< h$.
By decomposing the temperature and velocity fields between mean and fluctuating components, the heat equation can be written as
\begin{equation} \frac{\partial T'}{\partial t} +(\boldsymbol{\mathcal{U}}+\boldsymbol{u'})\boldsymbol{\cdot}\boldsymbol{\nabla}(\mathcal{T}+T')= {\nabla}^2 (\mathcal{T}+T') , \end{equation}
where 
$\boldsymbol {\mathcal {U}}=({U},{V},{W})^\textrm {T}=\langle \boldsymbol {u}\rangle _{\varphi h}$ and 
$\mathcal {T}=\langle T\rangle _{\varphi h}$ respectively represent the velocity and temperature mean component and 
$\boldsymbol {u'}=(u',v',w')^\textrm {T}$, 
$T'$ the fluctuating ones.
Thus, the heat equation reads
\begin{equation}
{U}\frac{\mathrm{d}\mathcal{T}}{\mathrm{d}r}+\left\langle
u'\frac{\partial T'}{\partial r}\right\rangle_{\varphi h}
+\left\langle({W}+w')\frac{\partial T'}{\partial
z}\right\rangle_{\varphi h}=\frac{1}{h} +
\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}
r}\left(r\frac{\mathrm{d}\mathcal{T}}{\mathrm{d} r}\right)+
\frac{1}{2{\rm \pi} h}\int_{0}^{2{\rm \pi}}\frac{\partial
T'}{\partial z}_{|z=h} \,\mathrm{d\varphi} .
\end{equation}Neglecting the mean vertical transport (third term on the left-hand side), diffusive terms (second and third terms on the right-hand side) and the transport induced by fluctuations (second and fourth terms on the left-hand side), (3.23) reduces to
\begin{equation} \frac{\mathrm{d}\mathcal{T}}{\mathrm{d}r} \simeq \frac{1}{\mathrm{U}h} . \end{equation} The radial temperature gradient appears to be inversely proportional to the radial velocity. Note that our approach here was guided and validated by empirical observations, and we did not perform any formal ordering between the various terms of (3.23). In figure 14 we plot the radial temperature profiles and the radial velocity modulus profiles for both regimes 
$R_F=0.1$ and 
$R_F=0.9$, with 
$\varGamma =16$ and 
$Ra_{\phi }$ between 
$10^5$ and 
$10^8$. We also plot the theoretical radial temperature gradient estimated by (3.24). To do so, 
${U}$ was estimated as the mean value of the radial velocity modulus (
$\langle |\langle u\rangle _{\varphi }|\rangle _z$) when 
$2< r<10$ and 
$2< r<12$, respectively, for 
$R_F=0.1$ and 
$R_F=0.9$, i.e. when the velocity is relatively constant. It can be seen in figure 14 that when 
$Ra_{\phi }$ increases, the values of the radial velocity modulus profile increases while the radial temperature gradient decreases. Furthermore, we can see that our estimate (3.24) matches very well to the radial temperature profiles observed in simulations.
Figure 14. Radial profile of the mean temperature (a,b) and of the radial velocity modulus (c,d) for different 
$Ra_{\phi }$ and for (a,c) 
$R_F=0.1$ and (b,d) 
$R_F=0.9$. The other input parameters are 
$Pr=0.1$ and 
$\varGamma =16$. The dash-dotted lines correspond to the theoretical slope (3.24) for the top row and to the mean value of the radial velocity modulus for the bottom row.
This model derives from the assumption that the mean flow is responsible for the radial temperature gradient. One could have argued that on the contrary, the radial temperature gradient creates the flow by a baroclinic torque: the temperature gradient would then be proportional (and not inversely proportional as we observe) to 
$U$. Thus, as will be further detailed in the next sections, we conclude that the side cold temperature generates the mean flow that then builds up the radial temperature gradient. Incidentally, the radial gradient scaling with 
$Ra_{\phi }$ shown in figure 12, 
$G_T\sim Ra_{\phi }^{-1/3}$, implies that it decreases faster with 
$Ra_{\phi }$ than the mean temperature of the system (
$\langle T\rangle \sim Ra_{\phi }^{-1/5}$, see (3.1)). This hints that the radial thermal gradient is a secondary aspect of the thermal structure. At first order, the plateau scaling is the dominant factor, which explains why the scaling of the low flux ratio regime works well at high 
$Ra_{\phi }$ values for all 
$R_F$, at least in the vicinity of the sidewall.
3.5.2. Temperature plateau value
The physical reason behind the uniform radial temperature profile near the sidewall, particularly visible at low 
$R_F$, is not trivial. Indeed, figure 15 shows that the temperature map averaged over azimuth and time, but not over depth, is very heterogeneous along the vertical direction. A cold zone is localized at the bottom close to the edge, in close connection with the localized, strong velocity zone seen in figure 13(a,c). Indeed, close to the sidewall, there is a strong downward flow that then turns into a cold jet with a characteristic radial extent, seemingly corresponding to the plateau area on the temperature profile. The larger 
$Ra_{\phi }$, the more extended this area. Actually, as we shall see, when 
$Ra_{\phi }$ increases, the cold jet that drives the global recirculation gets more inertia and propagates more into the domain. Averaging over depth, the 
$T_p$ plateau value is well predicted by (3.11) resulting from the analysis of the heat transfer through the side boundary layer made in the first regime. We now focus on understanding the cold jet dynamics in order to estimate the 
$r_p$ parameter.
Figure 15. Map in the 
$r,z$ plane of the temperature field averaged in time and azimuthal direction. The input parameters are 
$Ra_{\phi }=10^7$, 
$\varGamma =8$, 
$R_F=0.1$ and 
$Pr=0.1$.
3.5.3. Cold jet penetration length
We assume that the cold jet is a turbulent, self-similar structure, where a given amount of momentum initially injected at the side is propagating in the radial direction and heated from below by the lower plate. Following the seminal work of Morton, Taylor & Turner (Reference Morton, Taylor and Turner1956), it is well known that the thickness of jets and thermal plumes increases linearly along their propagation direction due to the turbulent entrainment of the ambient fluid (see also, e.g. List Reference List1982; Turner Reference Turner1986). Our structure can be seen as a mixture between a jet and a plume, and so it is natural to fit its thickness by a linear growth as
\begin{equation} b = \alpha (\varGamma-r) + r_0, \end{equation}
where 
$r_0$ is the radius close to the sidewall and 
$\alpha$ a constant known as the entrainment coefficient. In figure 16 the thickness of the cold jet is plotted as a function of the radius for the case where 
$Ra_{\phi }=10^8$, 
$\varGamma =8$, and for 
$R_F=0$ and 
$0.7$. We estimate the cold jet thickness at a given 
$r$ by seeking the depth where the ratio between the radial velocity (averaged in time and azimuthal direction) and the maximal value of this velocity located in the cold jet is equal to 
$0.5$. We see the linear behaviour corresponding to the cold jet spreading starting near the edge and stopping where 
$r\approx 5/8 \varGamma$ for 
$R_F=0$ and 
$r\approx 0.8\varGamma$ for 
$R_F=0.7$. The experimental entrainment coefficient ranges between 
$0.091$ and 
$0.121$. These values are consistent in terms of order of magnitude with values reported in the jets/plumes literature for other configurations (e.g. Carazzo, Kaminski & Tait Reference Carazzo, Kaminski and Tait2006; Fischer et al. Reference Fischer, List, Koh, Imberger and Brooks2007; Van Reeuwijk & Craske Reference Van Reeuwijk and Craske2015). Note that it does not make sense to be more quantitative here, since each configuration gives a different value for 
$\alpha$.
Figure 16. Plot of the cold plume thickness 
$b$ as a function of the radius. The blue points are the data, while the continuous line in orange is the best fit. The insets show the radial velocity normalized by the maximal velocity inside the jet, as a function of the height normalized by the jet thickness. Several radii indicated by the colour bar are shown going from near the edge (red) to the bulk (blue). Here 
$Ra_{\phi }=10^8$, 
$\varGamma =8$, 
$Pr=0.1$ and 
$R_F=0$ and 
$0.7$, respectively, for (a,b).
Furthermore, in the encapsulated graphs of figure 16, we plot the normalized radial velocity with the height normalized by the thickness of the jet. The radial velocity is normalized by the maximum values observed within the jet denoted 
$v_{max}$, which are determined using the expression
\begin{equation} v_{max}(r)=\mathrm{max}\left[|\left\langle u(t,r,\varphi,0\leqslant z\leqslant b) \right\rangle_{\varphi}|\right]. \end{equation}
This analysis is made for various radii ranging from the edge to the bulk. For both flux ratios in figure 16, all velocity profiles collapse on a unique profile that indicates a self-similar structure. In addition, when considering a fixed value of 
$Ra_{\phi }$ and 
$\varGamma$, a lower flux ratio results in a larger radial extent of the cold jet (not shown). The propagation of a jet on a heated plate is a well-studied topic (Schneider & Wasel Reference Schneider and Wasel1985; Steinrück Reference Steinrück1995; Higuera Reference Higuera1997; Fernandez-Feria, del Pino & Fernández-Gutiérrez Reference Fernandez-Feria, del Pino and Fernández-Gutiérrez2014; Fernandez-Feria & Castillo-Carrasco Reference Fernandez-Feria and Castillo-Carrasco2016), often examined in the context of a jet with uniform inlet velocity on a uniformly heated (imposed temperature) plate. The extension of the jet is determined by the competition between inertia and buoyancy forces. The velocity of the jet decreases due to entrainment of surrounding fluid, reducing its inertia, while lateral heating generates a vertical buoyancy force. The transition point is typically defined when the inertia is balanced by the buoyancy forces, often using a local Froude number to quantify this transition (Daniels & Gargaro Reference Daniels and Gargaro1993; Fernandez-Feria & Castillo-Carrasco Reference Fernandez-Feria and Castillo-Carrasco2016),
\begin{equation} \mathcal{F}r(r)=\left (\frac{v_{max}^2(r)}{\beta g \Delta T_{jet}(r)b^3(r)}\right)^{1/2} . \end{equation}
Here, 
$\Delta T_{jet}$ represents the temperature difference between the bottom and the outside of the cold plume 
$(\langle T\rangle _{\varphi }(r,z=0) -\langle T\rangle _{\varphi }(r,z=b) )$ at a given 
$r$. In figure 17 the local Froude number is plotted as a function of the radius for several 
$R_F$ going from 
$0$ to 
$0.9$. Two areas can be distinguished, near the edge where the cold jet inertia dominates (
$\mathcal {F}r>1$) and in the bulk where buoyancy finally dominates (
$\mathcal {F}r<1$). In order to reach the equilibrium between buoyancy and inertia (
$\mathcal {F}r\approx 1$), the cold jet travels a greater distance as the flux ratio is lower. When the Froude number falls below a certain threshold (function of 
$Pr$, 
$\varGamma$ and 
$Ra_{\phi }$), the cold jet is unable to be maintained, resulting in the radially inward motions being hindered by the adverse pressure gradient caused by buoyancy (Daniels & Gargaro Reference Daniels and Gargaro1993; Higuera Reference Higuera1997). In figure 18 we plot the radial temperature profile (
$\langle T\rangle _{\varphi z}(r)$) rescaled by the plateau temperature (
$T_p$) as a function of the radius rescaled by the radius where 
$\mathcal {F}r=1$, for various 
$R_F$. We observe that as 
$r/r(\mathcal {F}r=1)$ approaches 1, the rescaled temperature profiles converge towards the plateau temperature. This convergence is particularly clear when the flux ratio is low. As 
$R_F$ increases, the cold jet region becomes less meaningful, and in the asymptotic case where 
$R_F=1$, it effectively disappears since all the heat flux is evacuated at the top. The transition radius in (3.20) is thus well determined by the parameter 
$r_p = r(\mathcal {F}r=1)$. We might notice on figure 18 that accounting for this threshold is satisfying at first order only: while the critical Froude number is close to 1, its exact threshold value might also slightly depend on the aspect ratio and the Rayleigh number, which is left to future works.
Figure 17. Local Froude number as a function of radius for 
$Ra_{\phi }=10^8$, 
$\varGamma =8$, 
$Pr=0.1$ and 
$R_F$ going from 
$0$ to 
$0.9$. The dotted horizontal line indicates 
$\mathcal {F}r=1$ and the encapsulated plot represents a zoom on the area where 
$4.5< r<7$.
Figure 18. Radial temperature profile rescaled by the plateau temperature as a function of the radius rescaled by the radius where 
$\mathcal {F}r=1$, for 
$R_F$ going from 
$0$ to 
$0.9$. The input parameters are 
$Ra_{\phi }=10^8$, 
$\varGamma =8$ and 
$Pr=0.1$.
4. Conclusions and future works
In this paper a systematic numerical study was made of a system where a thin cylindrical layer of fluid (
$Pr=0.1$) is heated from below, one part of the heating power being extracted from the top surface, the other part being extracted from the side. This system is inherently heterogeneous in the radial direction: its spatial thermal structure results from the superposition of two asymptotic regimes corresponding, respectively, to forced and natural convection. Natural convection is of a Rayleigh–Bénard type and is modified by the presence of a convective structure, similar to a turbulent jet, along the bottom surface, originating from the sidewall. Combining various scaling laws, we have quantified these two regimes as well as their radial extent to propose a generic model of the radial temperature profile. This 1-D model is more relevant than existing models (developed mostly for nuclear safety analyses) that were based only on the average temperature of the system and neglected the radial variations of temperature (see Rein et al. (Reference Rein, Carénini, Fichot, Le Bars and Favier2023) for more details).
One of the extensions of this work is to relax the constraint of a uniform heat flux at the top and make it dependent on the radial position. This will make the analysis more consistent, since we have shown that the system cannot be considered as homogeneous in temperature, at least in the radial direction. This issue requires careful consideration and systematic analysis, following, for instance, the recent work of Clarté et al. (Reference Clarté, Schaeffer, Labrosse and Vidal2021).
The conclusions of this study allow us to better predict the average quantities of the system. However, relying solely on mean values is insufficient for nuclear safety analysis. Because of the turbulent nature of the system, important fluctuations are observed (in the calculations) at the sidewall (see e.g. the observed thermal branches in figure 2). Such events cannot be predicted by considering only the mean values. Therefore, it is also necessary to consider the statistics of fluctuations and the potential for extreme values of heat flux. In this view and for the extreme regimes relevant for the nuclear safety application (i.e. 
$Ra_\phi$ up to 
$10^{10}$), the direct simulation tool, even filtered, is limited, in particular to collect enough data for statistical convergence. Hence, our numerical study could be adequately complemented by an experimental study. All the above-mentioned points will be the focus of future works.
Acknowledgements
Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high-performance computing resources. This work was granted access to the HPC resources of IDRIS under the allocations A0120407543 and A0140407543 made by GENCI. Commissariat à l'Énergie Atomique et aux Énergies Alternatives (CEA) and Electricité de France (EDF) are gratefully acknowledged for their financial support. This work was granted access to the Very Large Computing Center (TGCC) of the CEA under the allocation 32004148 made by IRSN to the Centre de Calcul Recherche et Technologie (CCRT).
Funding
The authors declare the following interests regarding the funding and support received for this work: L'Institut de Radioprotection et de Sureté Nucléaire (IRSN), Commissariat à l’Énergie Atomique et aux Énergies Alternatives (CEA) and Electricité de France (EDF) provided financial support for this project. The aforementioned organizations had no influence on the design, data collection, analysis, interpretation of results or the decision to publish. The content of this work remains the sole responsibility of the authors.
Declaration of interests
The authors report no conflict of interest.
Appendix. Summary of the simulation parameters
Table 1. Simulations summary (DNS or filtered) according to the physical and numerical parameters.
