# Predictions of Reynolds and Nusselt numbers in turbulent convection using
machine-learning models

Shashwat Bhattacharya<EMAIL_ADDRESS>Institut für
Thermo-und Fluiddynamik, Technische Universität Ilmenau, Postfach 100565,
D-98684 Ilmenau, Germany  Mahendra K. Verma  Department of Physics, Indian
Institute of Technology Kanpur, Kanpur 208016, India  Arnab Bhattacharya
Department of Computer Science and Engineering, Indian Institute of Technology
Kanpur, Kanpur 208016, India

###### Abstract

In this paper, we develop a multivariate regression model and a neural network
model to predict the Reynolds number (Re) and Nusselt number in turbulent
thermal convection. We compare their predictions with those of earlier models
of convection: Grossmann-Lohse [Phys. Rev. Lett. 86, 3316 (2001)], revised
Grossmann-Lohse [Phys. Fluids 33, 015113 (2021)], and Pandey-Verma [Phys. Rev.
E 94, 053106 (2016)] models. We observe that although the predictions of all
the models are quite close to each other, the machine learning models
developed in this work provide the best match with the experimental and
numerical results.

##  I Introduction

Thermal convection is encountered in natural phenomena, such as stars and
planets, as well as in various engineering applications. A simplified version
of thermal convection is Rayleigh-Bénard convection (RBC), which is a model of
a fluid enclosed between two horizontal plates with the bottom plate kept at a
higher temperature than the top one. RBC is governed primarily by two
nondimensional parameters: the Rayleigh number Ra, which is the ratio of the
buoyancy and the dissipative force, and the Prandtl number Pr, which is the
ratio of kinematic viscosity and thermal diffusivity of the fluid
Chandrasekhar (1981); Ahlers, Grossmann, and Lohse (2009); Chillà and
Schumacher (2012). Two important response parameters of RBC are the Nusselt
number (Nu) and the Reynolds number (Re), which are respective measures of
large-scale heat transport and velocity in turbulent RBC Ahlers, Grossmann,
and Lohse (2009); Chillà and Schumacher (2012). The Reynolds number is defined
as $\mathrm{Re}=Ud/\nu$, where $U$ is the large-scale velocity, $\nu$ is the
kinematic viscosity, and $d$ is the distance between the thermal plates. The
Nusselt number is given by $\mathrm{Nu}=1+\langle u_{z}T\rangle/(\kappa\Delta
d^{-1})$, where $u_{z}$ is the vertical component of velocity, $T$ is the
temperature field, $\kappa$ is the thermal diffusivity, $\Delta$ is the
temperature difference between the thermal plates, and $\langle~{}\rangle$
denotes the volume average.

Researchers have attempted to model Re and Nu on the governing parameters, Ra
and Pr Ahlers, Grossmann, and Lohse (2009); Chillà and Schumacher (2012);
Siggia (1994); Xia (2013); Verma (2018). Early theoretical, numerical, and
experimental studies of RBC reveal a power-law scaling of Nu and Re, i.e.,
$\mathrm{Nu}\sim\mathrm{Ra}^{\alpha}\mathrm{Pr}^{\beta}$ and
$\mathrm{Re}\sim\mathrm{Ra}^{\gamma}\mathrm{Pr}^{\delta}$. The exponents
$\alpha,\beta,\gamma$, and $\delta$ vary for different regimes of Pr and Ra.
For the scaling of Nu, the exponent $\alpha$ ranges from 1/4 for
$\mathrm{Pr}\ll 1$ to approximately 1/3 for $\mathrm{Pr}\gtrsim 1$, Malkus
(1954); Shraiman and Siggia (1990); Cioni, Ciliberto, and Sommeria (1997);
Scheel and Schumacher (2017); Castaing  _et al._ (1989); Chavanne  _et al._
(1997); Horn, Shishkina, and Wagner (2013); Wagner and Shishkina (2013);
Kaczorowski and Xia (2013); Niemela and Sreenivasan (2003); Funfschilling  _et
al._ (2005); Stevens, Verzicco, and Lohse (2010); Dong  _et al._ (2020);
Madanan and Goldstein (2020); Vial and Hernández (2017). However, Nu has a
weak dependence on Pr, with the exponent $\beta$ ranging from approximately
zero for $\mathrm{Pr}\gtrsim 1$ to $0.14$ for small $\mathrm{Pr}$ Xia, Lam,
and Zhou (2002); Verzicco and Camussi (1999). Regarding the scaling of Re, the
exponent $\gamma$ was observed to be approximately $2/5$ for $\mathrm{Pr}\ll
1$, $1/2$ for $\mathrm{Pr}\sim 1$, and $3/5$ for $\mathrm{Pr}\gg 1$ Castaing
_et al._ (1989); Chavanne  _et al._ (1997); Cioni, Ciliberto, and Sommeria
(1997); Niemela  _et al._ (2001); Lam  _et al._ (2002); Emran and Schumacher
(2008); Silano, Sreenivasan, and Verzicco (2010); Verma  _et al._ (2012);
Wagner and Shishkina (2013); Scheel and Schumacher (2017); Horn, Shishkina,
and Wagner (2013); Pandey and Verma (2016); Pandey  _et al._ (2016); and
$\delta$ has been observed to range from $-0.7$ for $\mathrm{Pr}\lesssim 1$ to
$-0.95$ for $\mathrm{Pr}\gg 1$. Xia, Lam, and Zhou (2002); Brown,
Funfschilling, and Ahlers (2007). In the limit of very large Ra (the ultimate
regime), Kraichnan Kraichnan (1962) argued that
$\mathrm{Nu}\sim\sqrt{\mathrm{RaPr}}$, $\mathrm{Re}\sim\sqrt{\mathrm{Ra/Pr}}$
for $\mathrm{Pr}\leq 0.15$, and $\mathrm{Nu}\sim\sqrt{\mathrm{RaPr^{-1/2}}}$,
$\mathrm{Re}\sim\sqrt{\mathrm{Ra/Pr^{3/2}}}$ for $0.15<\mathrm{Pr}\leq 1$,
with logarithmic corrections. However, the existence of the aforementioned
regime in RBC is still under debate Verma  _et al._ (2012); Lohse and Toschi
(2003); Pawar and Arakeri (2016); Chavanne  _et al._ (1997); Schmidt  _et al._
(2012); Roche  _et al._ (2001); Ahlers  _et al._ (2012); He  _et al._ (2012);
Niemela  _et al._ (2000); Urban  _et al._ (2012).

Since the scaling of Re and Nu depends crucially on the regime of Ra and Pr,
researchers felt a need for a unified model that encompasses all the regimes.
Grossmann and Lohse Grossmann and Lohse (2000, 2001) (henceforth referred to
as GL) derived the first set of formulae for Re and Nu. Recently,
Bhattacharya, Verma, and Samtaney (2021a) revised the GL model by
incorporating the observed suppression in the dissipation rates in turbulent
convection Emran and Schumacher (2008); Pandey and Verma (2016); Pandey  _et
al._ (2016); Bhattacharya  _et al._ (2018); Bhattacharya, Samtaney, and Verma
(2019). Pandey and Verma Pandey and Verma (2016); Pandey  _et al._ (2016) also
derived formulae for Nu and Re. All the above models have been largely
successful in predicting Nu and Re over a range of Ra and Pr.

In this paper, we employ machine-learning (ML) tools to construct models for
predicting Re and Nu. This approach, which is increasingly gaining popularity
in fluid mechanics Parish and Duraisamy (2016); Fonda  _et al._ (2019); Pandey
and Schumacher (2020); Pandey, Schumacher, and Sreenivasan (2020); Brunton,
Noack, and Koumoutsakos (2020), involves building and improving prediction
algorithms by “learning” from the existing data Goodfellow, Bengio, and
Courville (2016); Hastie, Tibshirani, and Friedman (2009); Burkov (2019).
Hence, in the present work, we do not delve much into the physics behind Re
and Nu relations; instead, we use the data from previous works to develop the
prediction models. We construct a multivariate regression (MR) model and a
neural network (NN) model for the predictions of Nu and Re. We compare the
predictions of the above ML models with those of the earlier models. We
observe that the predictions of all the models are close to each other, but
the machine learning models developed in this study provide better match with
the experimental and numerical results.

The outline of the paper is as follows. In Sec. II, we briefly describe GL and
Pandey-Verma models. In Sec. III, we discuss the MR and NN models developed in
the present study. We compare the predictions of the above models with the
existing ones in Sec. IV. We conclude in Sec. V.

##  II A summary of existing models

In this section, we summarize GL and Pandey-Verma models briefly.

###  II.1 Grossmann-Lohse (GL) model

Grossmann and Lohse Grossmann and Lohse (2000, 2001) provided the first
unifying model for thermal convection. Grossmann and Lohse Grossmann and Lohse
(2000, 2001) derived this model by splitting the total viscous and thermal
dissipation rates ($\tilde{D}_{u}$ and $\tilde{D}_{T}$ respectively) into bulk
and the boundary layer (BL) contributions as

$\displaystyle\tilde{D}_{u}$ $\displaystyle=$
$\displaystyle\tilde{D}_{u,\mathrm{bulk}}+\tilde{D}_{u,\mathrm{BL}},$ (1)
$\displaystyle\tilde{D}_{T}$ $\displaystyle=$
$\displaystyle\tilde{D}_{T,\mathrm{bulk}}+\tilde{D}_{T,\mathrm{BL}}.$ (2)

Using the Prandtl-Blasius theory of boundary layers, and assuming homogeneous
and isotropic turbulence-like properties of dissipation rates in the bulk,
Grossmann and Lohse Grossmann and Lohse (2000, 2001) expressed the bulk and
boundary layer contributions in terms of Re, Nu, and Pr. These contributions
were then substituted in the exact relations relating the total dissipation
rates with the Nusselt number Shraiman and Siggia (1990):

$\displaystyle\tilde{D}_{u}$ $\displaystyle=$
$\displaystyle\frac{\nu^{3}}{d^{4}}\mathrm{(Nu}-1)\frac{\mathrm{Ra}}{\mathrm{Pr}^{2}},$
(3) $\displaystyle\tilde{D}_{T}$ $\displaystyle=$
$\displaystyle\frac{\kappa\Delta^{2}}{d^{2}}\mathrm{Nu}.$ (4)

Following the above procedure, Grossmann and Lohse arrived at general
expressions for Nu and Re valid for all Pr and Ra. The formulae for Re and Nu
were later updated by Stevens  _et al._ (2013) using more up-to-date data of
RBC. The GL model provides good estimates for the Nusselt and Reynolds numbers
for a wide range of Ra and Pr.

Table 1: Summary of different models of RBC for predicting the Renolds number (Re) and the Nusselt number (Nu). Note: SS = Shraiman and Siggia Shraiman and Siggia (1990), $\epsilon_{u}$ = viscous dissipation rate, $\epsilon_{T}$ = thermal dissipation rate, ML = Machine learning. Model | Year | Type | Prediction algorithm | Comments
---|---|---|---|---
Grossmann and Lohse | 2001, | Analytical | SS formulae | $\epsilon_{u}\sim U^{3}/d$
(GL) Grossmann and Lohse (2000, 2001); Stevens  _et al._ (2013) | updated in 2013 | \+ data |  | $\epsilon_{T}\sim U\Delta^{2}/d$
Revised Grossmann- | 2021 | Analytical | SS formula with | $\epsilon_{u}\sim(U^{3}/d)f(\mathrm{Ra,Pr})$
Lohse (RGL) Bhattacharya, Verma, and Samtaney (2021a) |  | + | modified expression | $\epsilon_{T}\sim(U\Delta^{2}/d)f(\mathrm{Ra,Pr})$
|  | data | for dissipation rates |
Pandey and Verma | 2016 | Analytical | Model momentum | Prefactors for the
(PV) Pandey and Verma (2016); Pandey  _et al._ (2016) |  | + | equation | momentum equation
|  | data |  | terms
Multivariate regression | Present study | ML + data | Multivariate linear |
(MR) |  |  | regression |
Neural network (NN) | Present study | ML + data | Neural networks |

###  II.2 Revised Grossmann-Lohse (RGL) model

Although the GL model has been quite successful in explaining past
experimental and numerical results, it based on certain assumptions that are
not very accurate for RBC. First, the model assumes that the viscous and
thermal dissipation rates in the bulk scale as $U^{3}/d$ and $U\Delta^{2}/d$
(for $\mathrm{Pr}\leq 1$) respectively, as in passive scalar turbulence with
open boundaries Lesieur (2008); Verma (2019). Recent studies, however, show
that these rates get an additional correction of approximately
$\mathrm{Ra}^{-0.2}$ for $\mathrm{Pr}\sim 1$ Emran and Schumacher (2008);
Pandey and Verma (2016); Pandey  _et al._ (2016); Scheel and Schumacher
(2017); Bhattacharya  _et al._ (2018); Bhattacharya, Samtaney, and Verma
(2019). The above correction is because of the inhibition of nonlinear
interactions due to the presence of walls and buoyancy Verma (2018); Pandey
and Verma (2016); Pandey  _et al._ (2016); Bhattacharya  _et al._ (2019). In
addition, the above studies reveal that the thickness of the viscous boundary
layer in RBC deviates from $\mathrm{Re}^{-1/2}$, contrary to what is assumed
in GL model Scheel, Kim, and White (2012); Shi, Emran, and Schumacher (2012);
Bhattacharya  _et al._ (2018).

Recently, Bhattacharya, Verma, and Samtaney (2021a) incorporated the above
corrections in the dissipation rates and the boundary layer thickness in the
equations of the GL model. In the revised model (henceforth referred to as
Revised Grossmann-Lohse or RGL model), the Reynolds and Nusselt numbers are
obtained by solving a cubic polynomial equation consisting of four functions
[$f_{i}(\mathrm{Ra,Pr})$] that are prefactors for the dissipation rates in the
bulk and boundary layers. The functional forms of the prefactors in the RGL
model were determined using multivariate linear regression on 60 sets of
simulation data of RBC. Bhattacharya, Verma, and Samtaney (2021a) observed
that the predictions of the RGL model were marginally better than the GL
model, especially for extreme Prandtl numbers Bhattacharya, Verma, and
Samtaney (2021a).

###  II.3 Pandey-Verma (PV) model

Pandey and Verma Pandey and Verma (2016); Pandey  _et al._ (2016) used a
different approach from the aforementioned models. To compute the Péclet
number ($\mathrm{Pe}=\mathrm{RePr}$) as a function of Ra and Pr, Pandey and
Verma Pandey and Verma (2016); Pandey  _et al._ (2016) constructed a model
(henceforth referred to as PV model) by estimating the relative strengths of
the various terms of the momentum equation of RBC using data from their
numerical simulations. Based on these estimates, they arrived at a quadratic
formula for the Péclet number. The predictions of Re by the PV model were as
accurate as the those by the GL model Verma (2018). Pandey and Verma Pandey
and Verma (2016); Pandey  _et al._ (2016) also derived a general formula for
the Nusselt number using correlations. Note that the relation for Nu was
obtained using simulation data for $\mathrm{Pr}\geq 1$ and thus is not
accurate for $\mathrm{Pr}\ll 1$ regime as shown later in this paper.

The main aspects of GL, RGL, and PV models are summarized in Table 1. In the
next section, we describe the machine-learning models developed in the present
study.

##  III Machine-Learning models

Having discussed three important models of the past, we now describe two
machine-learning models to determine the Nusselt and Reynolds numbers for
given Rayleigh and Prandtl numbers. The first one, referred to as Multivariate
Regression model (MR), employs linear regression to the available data
directly. The second approach, referred to as Neural Network (NN) model,
employs artificial neural networks for the predictions.

To construct the ML models, we use 62 sets of simulation data as training
datasets. Sixty of these datasets were generated by Bhattacharya, Verma, and
Samtaney (2021a) using high-resolution direct numerical simulations of RBC.
The Rayleigh number for these datasets ranges from $5\times 10^{5}$ to
$5\times 10^{9}$, and the Prandtl number ranges from $0.02$ to $100$. All the
above simulations were conducted using the finite-difference code SARAS Samuel
_et al._ (2021); Verma  _et al._ (2020) on a cubical domain of unit dimension
with no-slip walls on all the sides. Isothermal boundary conditions were
imposed on the horizontal walls. The grid resolution varied between $257^{3}$
to $1025^{3}$, depending on the resolution requirements. For further details
on the simulations and their validations, refer to Bhattacharya, Verma, and
Samtaney (2021a, b). Additionally, we conduct two simulations of RBC for
$\mathrm{Pr}=0.02$ with $\mathrm{Ra}=10^{5}$ and $2\times 10^{5}$; these
simulations were performed on a $513^{3}$ grid.  Here, we remark that the
output parameters — Re and Nu — are global quantities that are based on
velocity and heat fluxes that are averaged over the entire domain. Further, Nu
and Re used for training the model are not computed for a specific timeframe;
rather they are computed for multiple timeframes (ranging from 14 to 263; see
Bhattacharya, Verma, and Samtaney (2021a) for details) and then averaged over
these timeframes. Thus, we take care of the spatial and temporal fluctuations
of quantities that are typically encountered in turbulent flows. The
parameters of all the data samples, along with their corresponding Re and Nu,
are tabulated in Table 2.

Figure 1: (color online) For $\mathrm{Ra}=10^{6}$ (black) and $\mathrm{Ra}=10^{7}$ (red): (a) plots of Re versus Pr, and (b) plots of Nu versus Pr. Also shown are the piecewise linear fits for Re and Nu with the breakpoints (exhibited by vertical dashed lines) computed using the “pwlf” module Jekel and Venter (2019) of python. Figure 2: (color online) Plots of the matching functions $H_{i}(\mathrm{Pr})$ vs. Pr. $H_{1}$, $H_{2}$, and $H_{3}$ become unity in the regimes given by $\mathrm{Pr}<0.45$, $0.45<\mathrm{Pr}<10$, and $\mathrm{Pr}>10$ respectively. They attain the value of $1/2$ at the regime boundaries and become negligible outside their respective regimes. Table 2: Details of the training dataset for the Multivariate Regression (MR) and the Neural Network (NN) models. The input parameters are the Prandtl number (Pr) and the Rayleigh number (Ra), and the output parameters are the Nusselt number (Nu) and the Reynolds number (Re). S. No | Pr | Ra | Nu | Re | S. No | Pr | Ra | Nu | Re
---|---|---|---|---|---|---|---|---|---
1 | 0.02 | $1\times 10^{5}$ | 3.01 | 1180 | 32 | 6.8 | $1\times 10^{6}$ | 7.90 | 24.9
2 | 0.02 | $2\times 10^{5}$ | 3.52 | 1640 | 33 | 6.8 | $2\times 10^{6}$ | 9.46 | 35.6
3 | 0.02 | $5\times 10^{5}$ | 4.48 | 2440 | 34 | 6.8 | $5\times 10^{6}$ | 12.9 | 59.7
4 | 0.02 | $1\times 10^{6}$ | 5.78 | 3200 | 35 | 6.8 | $1\times 10^{7}$ | 15.9 | 89.2
5 | 0.02 | $2\times 10^{6}$ | 6.91 | 4290 | 36 | 6.8 | $2\times 10^{7}$ | 19.5 | 128
6 | 0.02 | $5\times 10^{6}$ | 8.85 | 6650 | 37 | 6.8 | $5\times 10^{7}$ | 26.1 | 217
7 | 0.02 | $1\times 10^{7}$ | 10.3 | 9420 | 38 | 6.8 | $1\times 10^{8}$ | 31.6 | 314
8 | 0.1 | $5\times 10^{5}$ | 6.11 | 749 | 39 | 6.8 | $2\times 10^{8}$ | 38.5 | 452
9 | 0.1 | $1\times 10^{6}$ | 7.34 | 1030 | 40 | 6.8 | $5\times 10^{8}$ | 50.5 | 729
10 | 0.1 | $2\times 10^{6}$ | 8.85 | 1380 | 41 | 6.8 | $1\times 10^{9}$ | 65.7 | 1070
11 | 0.1 | $5\times 10^{6}$ | 11.3 | 2090 | 42 | 6.8 | $2\times 10^{9}$ | 77.0 | 1520
12 | 0.1 | $1\times 10^{7}$ | 13.9 | 2870 | 43 | 6.8 | $5\times 10^{9}$ | 101 | 2400
13 | 0.1 | $2\times 10^{7}$ | 16.4 | 3870 | 44 | 50 | $1\times 10^{6}$ | 8.17 | 3.53
14 | 0.1 | $5\times 10^{7}$ | 20.8 | 6020 | 45 | 50 | $2\times 10^{6}$ | 9.66 | 5.19
15 | 0.1 | $1\times 10^{8}$ | 26.7 | 8140 | 46 | 50 | $5\times 10^{6}$ | 13.8 | 9.38
16 | 0.5 | $1\times 10^{6}$ | 8.38 | 285 | 47 | 50 | $1\times 10^{7}$ | 16.7 | 14.0
17 | 0.5 | $3\times 10^{6}$ | 11.4 | 482 | 48 | 50 | $2\times 10^{7}$ | 20.2 | 21.1
18 | 0.5 | $1\times 10^{7}$ | 15.9 | 874 | 49 | 50 | $5\times 10^{7}$ | 26.4 | 35.2
19 | 0.5 | $3\times 10^{7}$ | 21.6 | 1480 | 50 | 50 | $1\times 10^{8}$ | 31.8 | 50.8
20 | 0.5 | $1\times 10^{8}$ | 30.6 | 2610 | 51 | 50 | $2\times 10^{8}$ | 38.7 | 76.4
21 | 1 | $1\times 10^{6}$ | 8.18 | 147 | 52 | 50 | $5\times 10^{8}$ | 51.8 | 137
22 | 1 | $2\times 10^{6}$ | 10.1 | 213 | 53 | 50 | $1\times 10^{9}$ | 61.5 | 202
23 | 1 | $5\times 10^{6}$ | 13.3 | 340 | 54 | 100 | $1\times 10^{6}$ | 7.94 | 1.80
24 | 1 | $1\times 10^{7}$ | 16.3 | 491 | 55 | 100 | $2\times 10^{6}$ | 10.4 | 2.78
25 | 1 | $2\times 10^{7}$ | 19.8 | 702 | 56 | 100 | $5\times 10^{6}$ | 13.9 | 4.90
26 | 1 | $5\times 10^{7}$ | 26.0 | 1100 | 57 | 100 | $1\times 10^{7}$ | 16.8 | 7.02
27 | 1 | $1\times 10^{8}$ | 31.4 | 1530 | 58 | 100 | $2\times 10^{7}$ | 20.1 | 9.91
28 | 1 | $2\times 10^{8}$ | 38.6 | 2170 | 59 | 100 | $5\times 10^{7}$ | 26.1 | 17.1
29 | 1 | $5\times 10^{8}$ | 49.2 | 3330 | 60 | 100 | $1\times 10^{8}$ | 31.8 | 26.0
30 | 1 | $1\times 10^{9}$ | 61.2 | 4700 | 61 | 100 | $2\times 10^{8}$ | 39.1 | 37.5
31 | 1 | $2\times 10^{9}$ | 76.8 | 6580 | 62 | 100 | $5\times 10^{8}$ | 49.7 | 71.4
Figure 3: (color online) A schematic diagram of the neural network models for
the predictions of (a) Re and (b) Nu. Here, the input, hidden, and output
layers are depicted by blue, orange and red nodes respectively. The models for
Re and Nu have 6 and 16 hidden nodes respectively.

###  III.1 Multivariate regression (MR) model

We construct multivariate regression (MR) model using machine-learning
software WEKA Frank  _et al._ (2009) and obtain the functional forms of Re and
Nu. We look for a power-law relation for Re and Nu of the form
$A\mathrm{Ra}^{\alpha}\mathrm{Pr}^{\beta}$, take logarithms of this
expression, and employ linear regression to obtain $A$, $\alpha$, and $\beta$.
However, as discussed in Sec. I, the Pr and Ra dependences of Re and Nu are
different at small, moderate, and large Prandtl number regimes. Hence, the
aforementioned power-law will not be accurate of all the regimes of Pr.

To overcome the above problem, we split our parameter space into three regimes
– small Pr, moderate Pr, and large Pr – such that in each regime, the
quantities $A$, $\alpha$, and $\beta$ do not vary significantly. We obtain the
boundaries (or breakpoints) between the regimes as follows. For
$\mathrm{Ra}=10^{6}$ and $10^{7}$, we plot the logarithms of Nu and Re versus
the logarithm of Pr and fit a piecewise linear curve with three segments [see
Figs. 1(a) and (b)]. We employ the “PiecewiseLinFit” function from the “pwlf”
module of python Jekel and Venter (2019). The above module determines the
breakpoints of the piecewise linear curve using an optimization process that
solves a least squares fit several times for different breakpoint locations.
The two breakpoints of Pr versus Nu and Pr versus Re curves determined by the
“pwlf” module for $\mathrm{Ra}=10^{6}$ and $10^{7}$ are tabulated in Table 3.
Finally, we average the breakpoints over the two Rayleigh numbers and using
this average, we obtain the parameter regimes for Re as

Small Pr $\displaystyle:$ $\displaystyle\mathrm{Pr\leq 0.45},$ Moderate Pr
$\displaystyle:$ $\displaystyle 0.45\leq\mathrm{Pr}\leq 10,$ Large Pr
$\displaystyle:$ $\displaystyle\mathrm{Pr}\geq 10,$

and the parameter regimes for Nu as

Small Pr $\displaystyle:$ $\displaystyle\mathrm{Pr\leq 0.25},$ Moderate Pr
$\displaystyle:$ $\displaystyle 0.25\leq\mathrm{Pr}\leq 4.6,$ Large Pr
$\displaystyle:$ $\displaystyle\mathrm{Pr}\geq 4.6.$

Table 3: For $\mathrm{Ra}=10^{6}$ and $10^{7}$: boundaries separating the small, moderate and large Pr regimes. Here, $\mathrm{Pr}_{l,\mathrm{Re}}$ = boundary between the small and moderate Pr regimes for the predictions of Re; $\mathrm{Pr}_{u,\mathrm{Re}}$ = boundary between the moderate and large Pr regimes for the predictions of Re; $\mathrm{Pr}_{l,\mathrm{Nu}}$ = boundary between the small and moderate Pr regimes for the predictions of Nu; $\mathrm{Pr}_{u,\mathrm{Nu}}$ = boundary between the moderate and large Pr regimes for the predictions of Nu. Ra | $\mathrm{Pr}_{l,\mathrm{Re}}$ | $\mathrm{Pr}_{u,\mathrm{Re}}$ | $\mathrm{Pr}_{l,\mathrm{Nu}}$ | $\mathrm{Pr}_{u,\mathrm{Nu}}$
---|---|---|---|---
$10^{6}$ | 0.25 | 6.8 | 0.28 | 2.5
$10^{7}$ | 0.65 | 14 | 0.22 | 6.8
Average | 0.45 | 10 | 0.25 | 4.6

To ensure continuity across the parameter regimes, we employ six matching
functions $H_{i}(\mathrm{Pr})$; $i=1$ to $6$. The functions $H_{1}$, $H_{2}$,
$H_{3}$, $H_{4}$, $H_{5}$, and $H_{6}$ are such that they become unity inside
the regimes defined by $\mathrm{Pr}<0.45$, $0.45<\mathrm{Pr}<10$,
$\mathrm{Pr}>10$, $\mathrm{Pr}<0.25$, $0.25<\mathrm{Pr}<4.6$, and
$\mathrm{Pr}>4.6$, respectively, and become negligible outside their regimes.
The value of these functions is 1/2 at the boundaries of their respective
regimes. The following expressions of $H_{i}$ satisfy the aforementioned
conditions:

$\displaystyle H_{1}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{1}(0.45-\mathrm{Pr})}},$ (5) $\displaystyle
H_{2}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{1}(\mathrm{Pr}-0.45)}}-\frac{1}{1+e^{-k_{2}(\mathrm{Pr}-10)}},$
(6) $\displaystyle H_{3}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{2}(\mathrm{Pr}-10)}},$ (7) $\displaystyle
H_{4}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{1}(0.25-\mathrm{Pr})}},$ (8) $\displaystyle
H_{5}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{1}(\mathrm{Pr}-0.25)}}-\frac{1}{1+e^{-k_{2}(\mathrm{Pr}-4.6)}},$
(9) $\displaystyle H_{6}(\mathrm{Pr})$ $\displaystyle=$
$\displaystyle\frac{1}{1+e^{-k_{2}(\mathrm{Pr}-4.6)}},$ (10)

where $k_{1}=50$ and $k_{2}=7$. The above values of $k_{1}$ and $k_{2}$
provide the optimal steepness for the $H_{i}$ curves at the boundary points.
See Fig. 2 for an illustration of the behavior of the matching functions
$H_{1}$, $H_{2}$, and $H_{3}$. We then determine the logarithms of the
prefactor $A$ and the exponents $\alpha$ and $\beta$ for each regime using the
linear regression function of WEKA, and then combine these regimes using the
matching functions $H_{i}$.

We make use of cross-validation technique to to get an insight on how the MR
model will generalise to an independent dataset and thus to estimate how
accurately it will perform in practice. We employ ten-fold cross-validation,
which involves random partitioning of the complete dataset into 10 subdatasets
of equal size. Of these 10 subdatasets, a single subdataset is retained as the
validation data for testing the model, and the remaining 9 subdatasets are
used as training data. This process is then repeated 10 times, with each of
the 10 subdatasets used exactly once as the validation data. The mean of the
absolute percentage error, given by

$\mathcal{D}=\left|\frac{\mbox{Predicted value}-\mbox{Actual
value}}{\mbox{Actual value}}\right|\times 100,$ (11)

is computed for all the samples of each fold and then averaged over all the
folds. The above average is observed to be 4% for the predictions of both Re
and Nu, thus ensuring that the model is robust. After having done the cross-
validation, we invoke the learning algorithm for a final (11th) time on the
entire dataset and obtain the following expressions for Re and Nu.

$\displaystyle\mathrm{Re}$ $\displaystyle=$ $\displaystyle
0.40H_{1}\mathrm{Ra}^{0.45}\mathrm{Pr}^{-0.72}+0.11H_{2}\mathrm{Ra}^{0.52}\mathrm{Pr}^{-0.86}$
(12) $\displaystyle+0.055H_{3}\mathrm{Ra}^{0.58}\mathrm{Pr}^{-0.99},$
$\displaystyle\mathrm{Nu}$ $\displaystyle=$ $\displaystyle
0.24H_{4}\mathrm{Ra}^{0.27}\mathrm{Pr}^{0.16}+0.15H_{5}\mathrm{Ra}^{0.29}\mathrm{Pr}^{0.018}$
(13) $\displaystyle+0.13H_{6}\mathrm{Ra}^{0.30}\mathrm{Pr}^{0.0084}.$

Thus, interested researchers can use the above equations to obtain the Nusselt
and Reynolds numbers for an arbitrary set of Rayleigh and Prandtl numbers.

The above training algorithm was run on a personal computer consisting of an
8-gigahertz Core i5 processor with 8 gigabytes of RAM. The time taken for
executing the training algorithm is 0.3 second for both Re and Nu. It is
important to note that the present model is an application of machine-learning
model to the raw data itself, in contrast to RGL model where ML model was used
for parameter estimation of the prefactors of the dissipation equations [Eqs.
(1) and (2)].

###  III.2 Neural Network (NN) model

Lastly, we discuss the neural network model. Neural network is known to be a
more robust prediction algorithm than regression. This is due to the intrinsic
nonlinearity of neural networks and also by the virtue of universal
approximation theorems Hornik, Stinchcombe, and White (1989) associated with
these networks. A schematic diagram of the neural networks for predicting Re
and Nu are exhibited in Fig. 3(a,b). We employ the Keras module Gulli and Pal
(2017) of Python to construct the neural networks. The input features for both
the networks are the logarithms of Rayleigh and Prandtl numbers. The networks
consist of one hidden layer, with 6 nodes for Re predictions and 16 nodes for
Nu predictions. We arrived at the above configuration after trial and error as
it provided us the best results.  It must be noted from our discussion in Sec.
I that the range of the scaling exponent for Nu ($0.25\leq\alpha\leq 0.33$) is
shorter than that for Re ($0.40\leq\gamma\leq 0.60$); hence, more care is
needed for predicting Nu. Thus, more neurons are required in the neural
network for Nu.

The hidden layers of both the networks are activated by the Exponential Linear
Unit (ELU) function Goodfellow, Bengio, and Courville (2016). The networks are
set to minimize the mean-squared error between the actual values in the
training set and the predictions. We employ the Adam algorithm Gulli and Pal
(2017) to optimize our models. The networks are trained using the mini-batch
gradient descent algorithm with a batch size of 8, that is, eight training
examples are used to estimate the error gradient. The gradient descent
algorithm is run for 1000 epochs. The training algorithms were run on the same
personal computer on which the MR algorithm was executed. The computational
time taken to train the model is 11 seconds for both Re and Nu; thus the
neural networks take more time than the multivariate regression for training.
However, the time required for training is still not very large and the
necessary computations can easily be carried out on a laptop. The relevant
segment and functions of the code describing the model (for Re) are given
below.

def baseline_model():
model = Sequential()
’ ’ ’ Specify the number of input parameters (input_dim),
which in this case is 2. Add the hidden layer with
6 nodes, activated using ELU function.’ ’ ’
model.add(Dense(6, input_dim=2,
kernel_initializer = ’normal’, activation = ’elu’))
# Add the output layer containing one output node.
model.add(Dense(1, kernel_initializer =
’normal’))
# Specify the loss function and the optimizer
model.compile(loss=’mean_squared_error’,
optimizer=’adam’)
return model
# Specify the batch size and the number of epochs
estimator = KerasRegressor(build_fn = baseline_model, epochs = 1000,
batch_size = 8, verbose = 0)
#Functions for training and predicting:
#Train the model (X = input, Y = output)
estimator.fit(X,Y)
#Predict the values for all X
predictions = estimator.predict(X)

The relevant code segment for describing the model for Nu is similar to that
for Re, except that the number of nodes in the hidden layer is specified to be
16.  The user can pass a two-dimensional array of arbitrary Rayleigh and
Prandtl numbers as an argument to the “estimator.predict()” function; this
function returns the corresponding array for predictions. Refer to the
supplemental materials for the complete codes (“suppl1.txt” for predicting Re
and “suppl2.txt” for predicting Nu) and the input dataset (“suppl3.csv”).

We employ 10-fold cross-validation for our above models (similar to that done
for the MR model described earlier) to get an estimate of how the models
generalize to an independent dataset. We observe the mean absolute error from
the cross-validation to be 4% for the predictions of Re and 3% for the
predictions of Nu. The small values of the errors indicate that our neural-
network model has a robust predictive capability. It must be noted, however,
that cross-validation may be faulty in determining the prediction reliability
of a machine-learning model in certain cases Bishop (2006) and a more
appropriate tool is to compare the absolute percentage error between the final
model (developed using entire training dataset) and the actual values of the
test datasets. This will be covered in detail in Sec. IV.

To ascertain the absence of overfitting, we add weighted regularizing terms
(of $\mathrm{L}^{2}$ norm Goodfellow, Bengio, and Courville (2016)) to our
neural network models and observe how the models fair against the test data.
In the present work, the test data consists of the simulation results of
Scheel and Schumacher (2017) ($\mathrm{Pr}=0.005$ and $0.02$), Wagner and
Shishkina (2013) ($\mathrm{Pr}=0.7$), Emran and Schumacher (2008)
($\mathrm{Pr}=0.7$), Kaczorowski and Xia (2013) ($\mathrm{Pr}=4.38$), and
Horn, Shishkina, and Wagner (2013) ($\mathrm{Pr}=2547.9$); and the
experimental results of Cioni, Ciliberto, and Sommeria (1997)
($\mathrm{Pr}=0.02$) and Niemela  _et al._ (2001) ($\mathrm{Pr}=0.7$).

Table 4: Effect of adding regularizing hyperparameters to the neural networks on the performance of the models on test data. Here, $D_{\mathrm{Re}}$ and $D_{\mathrm{Nu}}$ respectively represent the percentage deviation of the predictions from the actual values of the test data. Regularizer hyperparameter | $D_{\mathrm{Re}}$ | $D_{\mathrm{Nu}}$
---|---|---
None | 19.2% | 9.76%
0.001 (Hidden layer) | 20.4% | 14.0%
0.001 (Hidden layer) | 20.5% | 14.4%
0.001 (Output layer) |  |
0.01 (Hidden layer) | 21.7% | 25.8%
0.01 (Hidden layer) | 23.5% | 26.1%
0.01 (Output layer) |  |

In Table 4, we list the regularizing hyperparameters added in the hidden and
the output layers of the neural networks and their effects on the accuracy of
the networks. The table shows that the addition of regularizing terms does not
improve the forecast; rather the forecast worsens as the values of the
hyperparameters increase. Thus, in the final neural network models, no
regularizing counterterms are used.

Table 5: Different trial configurations of neural networks and their corresponding performances. Here, HL = hidden layer, $D_{\mathrm{Re}}$ = percentage deviation of the predictions from the actual Re, $D_{\mathrm{Nu}}$ = percentage deviation of the predictions from the actual Nu. Configuration | $D_{\mathrm{Re}}$ | $D_{\mathrm{Re}}$ | Configuration | $D_{\mathrm{Nu}}$ | $D_{\mathrm{Nu}}$
---|---|---|---|---|---
(Re) | (training) | (testing) | (Nu) | (training) | (testing)
1 HL, |  |  | 1 HL, |  |
6 nodes per HL | 2.98% | 19.2% | 16 nodes per HL | 2.31% | 9.76%
2 HL, |  |  | 2 HL, |  |
3 nodes per HL | 2.98% | 36.0% | 8 nodes per HL | 3.04% | 17.9%
3 HL, |  |  | 4 HL |  |
2 nodes per HL | 5.93% | 81.6% | 4 nodes per HL | 4.64% | 17.1%

Finally, we remark that we also tried to use less width and more depth to our
neural networks. As shown in Table 5, we examined the performance of the
neural network using two and three hidden layers for predicting Re and two and
four hidden layers for predicting Nu. In all the above trials, we kept the
same total number of hidden nodes (that is, 6 for Re and 16 for Nu).
Interestingly, we observed that the performance of the networks decreases with
the increase in depth. This is likely because the deeper neural networks have
a tendency to overfit and hence their predictions are less accurate.

The important aspects of the MR and NN models are summarized in Table 1. In
the next section, we discuss and compare the predictions of these models with
respect to the GL, RGL, and PV models.

##  IV Performance of the ML models

Table 6: Absolute percentage error in the predictions of Re by the GL model, the RGL model, the PV model, the multivariate regression (MR) model, and the neural network (NN) model for different sets of simulation and experimental data. The asterisks (*) refer to the numerical data of Bhattacharya, Verma, and Samtaney (2021a), used for training the MR and NN models. Pr | GL | RGL | PV | MR | NN
---|---|---|---|---|---
0.005 | 48% | 11% | 11% | 13% | 4.5%
$0.02^{*}$ | 52% | 4.7% | 31% | 3.5% | 8.0%
0.02 | 53% | 11% | 19% | 8.8% | 9.0%
$0.1^{*}$ | 30% | 1.9% | 11% | 3.2% | 3.4%
$0.5^{*}$ | 14% | 1.3% | 3.3% | 4.5% | 3.1%
0.7 | 24% | 6.7% | 23% | 13% | 9.2%
$1^{*}$ | 20% | 2.8% | 8.6% | 5.1% | 1.4%
$6.8^{*}$ | 27% | 3.4% | 8.4% | 5.2% | 2.8%
$50^{*}$ | 84% | 6.0% | 12% | 4.5% | 3.4%
$100^{*}$ | 150% | 3.8% | 15% | 4.6% | 3.6%
2547.9 | 560% | 85% | 77% | 36% | 37%
Overall | 82% | 12% | 20% | 9.3% | 7.7%
Table 7: Absolute percentage error in the predictions of Nu by the GL model, the RGL model, the PV model, the multivariate regression (MR) model, and the neural network (NN) model for different sets of simulation and experimental data. The asterisks (*) refer to the numerical data of Bhattacharya, Verma, and Samtaney (2021a), used for training the MR and NN models. Pr | GL | RGL | PV | MR | NN
---|---|---|---|---|---
0.005 | 17% | 9.6% | 93% | 3.3% | 2.7%
$0.02^{*}$ | 34% | 8.9% | 150% | 2.2% | 5.1%
0.02 | 15% | 11% | 84% | 6.1% | 12%
$0.1^{*}$ | 5.0% | 3.1% | 22% | 6.2% | 2.5%
$0.5^{*}$ | 5.4% | 1.4% | 6.0% | 1.5% | 1.6%
0.7 | 11% | 9.1% | 9.7% | 6.5% | 3.0%
$1^{*}$ | 5.8% | 3.6% | 9.2% | 1.1% | 0.85%
4.38 | 6.3% | 5.7% | 7.6% | 1.0% | 0.95%
$6.8^{*}$ | 6.4% | 7.2% | 3.8% | 5.1% | 3.3%
$50^{*}$ | 7.2% | 3.2% | 7.3% | 4.9% | 1.9%
$100^{*}$ | 3.9% | 2.7% | 6.4% | 4.9% | 2.1%
2547.9 | 17% | 2.3% | 7.6% | 4.7% | 14%
Overall | 11% | 5.7% | 34% | 4.0% | 4.2%

We test the performance of the multivariate regression and neural network
models with the GL, RGL, and PV models. We compare their predictions for the
numerical results of Bhattacharya, Verma, and Samtaney
(2021a)($\mathrm{Pr}=0.02$, 0.1, 0.5, 1, 6.8, 50, and 100), Scheel and
Schumacher (2017) ($\mathrm{Pr}=0.005$ and $0.02$), along with the test
dataset Wagner and Shishkina (2013) ($\mathrm{Pr}=0.7$), Emran and Schumacher
(2008) ($\mathrm{Pr}=0.7$), Kaczorowski and Xia (2013) ($\mathrm{Pr}=4.38$),
and Horn, Shishkina, and Wagner (2013) ($\mathrm{Pr}=2547.9$); and the
experimental results of Cioni, Ciliberto, and Sommeria (1997)
($\mathrm{Pr}=0.02$), and Niemela  _et al._ (2001) ($\mathrm{Pr}=0.7$). The
simulations of Wagner and Shishkina (2013), Kaczorowski and Xia (2013), and
Bhattacharya, Verma, and Samtaney (2021a) involve a cubical cell, whereas the
others employ a cylindrical cell. The aspect ratio of the RBC cell is unity in
all the above works. We compute the absolute percentage error between the
estimated and actual values as per Eq. (11). For every Pr, we average the
aforementioned errors over all the Rayleigh numbers and list this average
error in Table 6 (for the predictions of Re) and Table 7 (for the predictions
of Nu). Note that for $\mathrm{Pr}=0.02$, we provide two entries in the above
tables: one corresponding to the data of Bhattacharya, Verma, and Samtaney
(2021a) (which were used for training), and the other corresponding to those
of Scheel and Schumacher (2017) and Cioni, Ciliberto, and Sommeria (1997). In
Fig. 4(a,b), we plot the predictions of the normalized Re by the different
models discussed in the paper along with the numerically computed values of
Refs. Bhattacharya, Verma, and Samtaney (2021a); Scheel and Schumacher (2017);
Wagner and Shishkina (2013); Emran and Schumacher (2008); Horn, Shishkina, and
Wagner (2013). See Fig. 4(c, d) for the corresponding plots for Nu. Note that
in the above figures, we normalize Nu with $\mathrm{Gr}^{0.3}$, where
$\mathrm{Gr}=\mathrm{Ra}/\mathrm{Pr}$ is the Grashoff number. This is because
Nu has only a weak dependence on Pr, due to which the Nu versus Ra curves for
different Pr’s would have tended to overlap on each other. Further, we do not
show the plots for $\mathrm{Pr}=6.8$ (only for Nu) and $\mathrm{Pr}=0.5$ in
the above figures because their datapoints lie very close to those of
$\mathrm{Pr}=4.38$ and $0.7$ respectively.

Figure 4: (color online) Comparison between Grossmann and Lohse (GL, blue
lines), Pandey and Verma (PV, orange lines), Revised Grossmann and Lohse (RGL,
black lines), multivariate regression (MR, green lines), and the neural
network (NN, red lines) models for the predictions of Re for (a)
$\mathrm{Pr}<1$, and (b) $\mathrm{Pr}\geq 1$; and of Nu for (c)
$\mathrm{Pr}<1$, and (d) $\mathrm{Pr}\geq 1$. The predictions are compared
with the past numerical and experimental results for $\mathrm{Pr}=0.005$
(triangles), $\mathrm{Pr}=0.02$ (squares), $\mathrm{Pr}=0.1$ (pentagons),
$\mathrm{Pr}=0.7$ (circles), $\mathrm{Pr}=1$ (left-pointing triangles),
$\mathrm{Pr}=4.38$ (right-pointing triangles), $\mathrm{Pr}=6.8$ (hexagons),
$\mathrm{Pr}=50$ (inverted triangles), and $\mathrm{Pr}=2547.9$ (diamonds).
The filled markers correspond to the training dataset for the MR and NN models
that includes the numerically computed values of Bhattacharya, Verma, and
Samtaney (2021b).

###  IV.1 Predictions of the Reynolds number

First, we analyse in detail the performance of the models in their predictions
of the Reynolds number. From Table 6 and Figs. 4(a,b), we observe that the GL
model’s predictions for Re are close to the actual values for moderate Prandtl
numbers. Note, however, that Re is sensitive to the way it is computed (that
is, based on root mean square velocity, maximum plume velocity, etc.). The RGL
model further improves the predictions of Re, especially for extreme Prandtl
numbers. The PV model also marginally improves the predictions of the Reynolds
number compared to the GL model but is not as accurate as the RGL model.

Now, we will analyze the performances of the MR and NN models developed in the
present study. Both the above models exhibit very good accuracy in their
predictions of Re and provide significantly better predictions than the GL and
the PV models. As exhibited in Table 6, the MR and NN models provide better
predictions than even the RGL model for large Prandtl numbers. For small Pr’s,
however, the MR model fares marginally worse than the NN and RGL models, with
the accuracy of the latter two models being comparable. It must be noted,
however, that the MR model performs better than the RGL model if we consider
all the Prandtl numbers into account (as shown in Table 6).

###  IV.2 Predictions of the Nusselt number

Having analyzed the predictions of Re, we will now examine the preformance of
the models in predicting Nu. Figures. 4(c,d) and Table 7 show that the GL
model predicts Nu with a good accuracy. The RGL model marginally improves the
predictions of Nu for moderate Pr but significantly improves the predictions
for extreme Pr. The PV model’s predictions of Nu are close to those by the GL
model for moderate and large Pr. However, for small Prandtl numbers, Nu
predicted by the PV model deviates significantly from the actual values. This
is because Pandey and Verma Pandey and Verma (2016); Pandey  _et al._ (2016)
derived their expression for Nu using the data from simulations with
$\mathrm{Pr}\geq 1$. As a result, the final expression for the Nusselt number
did not contain any Pr dependence, which is applicable only for moderate and
large Prandtl numbers. Hence, for small Prandtl numbers, the errors between
the predictions and the actual Nu values are very high for the PV model. Thus,
as far as the previously developed models are concerned, the RGL model
exhibits the best performance.

The MR and NN models developed in this study exhibit very good accuracy in
their predictions of Nu and are observed to mostly perform better than the GL
and PV models. The predictions of the MR and NN models are even better, albeit
marginally, than the RGL model. For small and moderate Pr’s
($\mathrm{Pr}\lesssim 1$), both MR and NN models perform better than the RGL
model, with the MR model marginally outperforming the NN model. For larger
Pr’s, with the exception of $\mathrm{Pr}<2547.9$, the NN model performs better
than the MR and RGL models. The performances of the latter two models are
comparable in the aforementioned regime. For $\mathrm{Pr}=2547.9$, the ML
models, especially the NN model, do not perform as well as they do for the
remaining Pr’s; this is possibly because of lack of data points for very large
Pr’s.

###  IV.3 Significance of the ML models

In Sec. IV.1 and Sec. IV.2, we have shown that the ML models developed in this
study exhibit better accuracy in predictions compared to GL and PV models,
Further, although the ML models offer better predictions than the RGL model as
well, the improvement is only marginal. In spite of the above, the ML models
are more promising than the RGL model as they easier to construct provided we
have sufficient training dataset. Fortunately, there is a large amount of
datasets available in the literature (see, for example, Ahlers, Grossmann, and
Lohse (2009), Chillà and Schumacher (2012)). Also, experiments and numerical
simulations of the past provide a wide range of Ra and Pr for training the
neural network. For example, Ahlers  _et al._ (2012), He  _et al._ (2012), and
Urban  _et al._ (2012) performed experiments for Ra as large as $10^{15}$,
while Niemela  _et al._ (2000) reached Ra up to $10^{17}$. In addition, Zhu
_et al._ (2018) and Iyer  _et al._ (2020) recently carried out numerical
simulations for Ra up to $10^{14}$ and $10^{15}$ respectively. Thus, the
machine learning models developed in this study has a good scope for
improvement with inclusion of these datasets.

On the other hand, the training dataset for the RGL model is more difficult to
construct because it requires the prefactors for the viscous and thermal
dissipation rates in the bulk and boundary layers as output parameters. The
data on these prefactors is quite limited; in fact, to the best of our
knowledge, only Bhattacharya, Verma, and Samtaney (2021a) have computed the
above quantities using numerical data of RBC. Further, there is no available
data on $f_{i}$’s for extreme Rayleigh numbers ($\gg 10^{10}$). Note that most
of the work on such large Rayleigh numbers are experimental because
simulations of RBC in that regime require large computational resources.
However, computation of prefactors for the dissipation rates requires
multipoint measurements, which is difficult in experiments. Due to these
constraints of the RGL model, we believe that NN and MR models hold a lot of
promise for a robust model for forecasting Re and Nu.  We conclude in the next
section.

##  V Summary and conclusions

In summary, we construct two ML models–multivariate regression (MR) model and
neural network (NN) model –to predict Re and Nu for given Ra and Pr. We
compared the models’ predictions with existing models. We observe that
although all the above models predict Re and Nu with good accuracy, the ML
models exhibit the best performace. In general, we observed that Re is more
sensitive to modeling parameters compared to Nu; this is because Re depends
more critically on Pr compared to Nu. We also discussed that the machine
learning models developed in this study can be enhanced further by using data
for a wider range of parameters.

The positive results from the machine learning models reinforce the importance
of data-driven methods for modelling turbulence, a practice that is gaining
popularity Pandey and Schumacher (2020); Pandey, Schumacher, and Sreenivasan
(2020); Fonda  _et al._ (2019); Parish and Duraisamy (2016). An accurate
estimation of the large-scale velocity and heat transport helps scientists and
engineers to better understand thermally driven flows encountered in nature
and in engineering applications. For example, our results will be useful in
understanding atmospheric flows and thus can aid in more accurate climate
modelling and weather predictions. On the other hand, our findings will help
engineers to better design buildings, heat sinks, solar collectors, and
chimneys. Further, since the improvement in the predictions by our models is
more pronounced for extreme Prandtl numbers, we believe that our results will
be useful in modelling natural flows and designing applications where such
extreme Prandtl numbers are involved. Examples of above include the mantle
flow (very large Prandtl number) and flows in liquid metal batteries (small
Prandtl number). We hope that in future, more machine-learning tools will be
constructed for thermal convection that will help model natural and
engineering flows.

## Acknowledgements

We thank J. Schumacher, K. R. Sreenivasan, A. Pandey, M. Sharma, R. Samuel,
and S. Alam for useful discussions. The present work was mostly conducted at
Indian Institute of Technology Kanpur, India, and was funded by the research
grant number SPO/STC/PHY/2018037 from Indian Space Research Organization,
India. The simulations of convection were performed on Shaheen II of King
Abdullah University of Science and Technology, Saudi Arabia, under the project
k1416. Shashwat Bhattacharya is currently funded by a postdoctoral fellowship
of Alexander von Humboldt Foundation (Germany).

## Conflicts of interest

The authors have no conflicts to disclose.

## Supplementary Material

See supplementary material for the python codes for predicting Re (suppl1.txt)
and Nu (suppl2.txt) using neural network, as well as the training dataset
(suppl3.csv).

## References

## References

  * Chandrasekhar (1981) S. Chandrasekhar,  _Hydrodynamic and Hydromagnetic Stability_ (Dover publications, Oxford, 1981).
  * Ahlers, Grossmann, and Lohse (2009) G. Ahlers, S. Grossmann, and D. Lohse, “Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection,” Rev. Mod. Phys. 81, 503–537 (2009).
  * Chillà and Schumacher (2012) F. Chillà and J. Schumacher, “New perspectives in turbulent Rayleigh-Bénard convection,” Eur. Phys. J. E 35, 58 (2012).
  * Siggia (1994) E. D. Siggia, “High Rayleigh number convection,” Annu. Rev. Fluid Mech. 26, 137–168 (1994).
  * Xia (2013) K.-Q. Xia, “Current trends and future directions in turbulent thermal convection,” Theor. App. Mech. Lett. 3, 052001 (2013).
  * Verma (2018) M. K. Verma,  _Physics of Buoyant Flows: From Instabilities to Turbulence_ (World Scientific, Singapore, 2018).
  * Malkus (1954) W. V. R. Malkus, “The Heat Transport and Spectrum of Thermal Turbulence,” Proceedings of the Royal Society of London. Series A 225, 196–212 (1954).
  * Shraiman and Siggia (1990) B. I. Shraiman and E. D. Siggia, “Heat transport in high-Rayleigh-number convection,” Phys. Rev. A 42, 3650–3653 (1990).
  * Cioni, Ciliberto, and Sommeria (1997) S. Cioni, S. Ciliberto, and J. Sommeria, “Strongly turbulent Rayleigh–Bénard convection in mercury: comparison with results at moderate Prandtl number,” J. Fluid Mech. 335, 111–140 (1997).
  * Scheel and Schumacher (2017) J. D. Scheel and J. Schumacher, “Predicting transition ranges to fully turbulent viscous boundary layers in low Prandtl number convection flows,” Phys. Rev. Fluids 2, 123501 (2017).
  * Castaing  _et al._ (1989) B. Castaing, G. Gunaratne, Kadanoff, L. P., A. Libchaber, and F. Heslot, “Scaling of hard thermal turbulence in Rayleigh-Bénard convection,” J. Fluid Mech. 204, 1–30 (1989).
  * Chavanne  _et al._ (1997) X. Chavanne, F. Chillà, B. Castaing, B. Hebral, B. Chabaud, and J. Chaussy, “Observation of the ultimate regime in Rayleigh-Bénard convection,” Phys. Rev. Lett. 79, 3648–3651 (1997).
  * Horn, Shishkina, and Wagner (2013) S. Horn, O. Shishkina, and C. Wagner, “On non-Oberbeck–Boussinesq effects in three-dimensional Rayleigh–Bénard convection in glycerol,” J. Fluid Mech. 724, 175–202 (2013).
  * Wagner and Shishkina (2013) S. Wagner and O. Shishkina, “Aspect-ratio dependency of Rayleigh-Bénard convection in box-shaped containers,” Phys. Fluids 25, 085110 (2013).
  * Kaczorowski and Xia (2013) M. Kaczorowski and K.-Q. Xia, “Turbulent flow in the bulk of Rayleigh–Bénard convection: small-scale properties in a cubic cell,” J. Fluid Mech. 722, 596–617 (2013).
  * Niemela and Sreenivasan (2003) J. J. Niemela and K. R. Sreenivasan, “Confined turbulent convection,” J. Fluid Mech. 481, 355–384 (2003).
  * Funfschilling  _et al._ (2005) D. Funfschilling, E. Brown, A. Nikolaenko, and G. Ahlers, “Heat transport in turbulent Rayleigh-Bénard convection in cylindrical samples with aspect ratio one and larger,” J. Fluid Mech. 536, 145–154 (2005).
  * Stevens, Verzicco, and Lohse (2010) R. J. A. M. Stevens, R. Verzicco, and D. Lohse, “Radial boundary layer structure and Nusselt number in Rayleigh–Bénard convection,” J. Fluid Mech. 643, 495–507 (2010).
  * Dong  _et al._ (2020) D.-L. Dong, B.-F. Wang, Y.-H. Dong, Y.-X. Huang, N. Jiang, Y.-L. Liu, Z.-M. Lu, X. Qiu, Z.-Q. Tang, and Q. Zhou, “Influence of spatial arrangements of roughness elements on turbulent rayleigh-bénard convection,” Phys. Fluids 32, 045114 (2020).
  * Madanan and Goldstein (2020) U. Madanan and R. J. Goldstein, “High-Rayleigh-number thermal convection of compressed gases in inclined rectangular enclosures,” Phys. Fluids 32, 017103 (2020).
  * Vial and Hernández (2017) M. Vial and R. H. Hernández, “Feedback control and heat transfer measurements in a Rayleigh-Bénard convection cell,” Phys. Fluids 29, 074103 (2017).
  * Xia, Lam, and Zhou (2002) K.-Q. Xia, S. Lam, and S.-Q. Zhou, “Heat-flux measurement in high-Prandtl-number turbulent Rayleigh-Bénard convection,” Phys. Rev. Lett. 88, 064501 (2002).
  * Verzicco and Camussi (1999) R. Verzicco and R. Camussi, “Prandtl number effects in convective turbulence,” J. Fluid Mech. 383, 55–73 (1999).
  * Niemela  _et al._ (2001) J. J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly, “The wind in confined thermal convection,” J. Fluid Mech. 449, 169–178 (2001).
  * Lam  _et al._ (2002) S. Lam, X.-D. Shang, S.-Q. Zhou, and K.-Q. Xia, “Prandtl number dependence of the viscous boundary layer and the Reynolds numbers in Rayleigh-Bénard convection,” Phys. Rev. E 65, 066306 (2002).
  * Emran and Schumacher (2008) M. S. Emran and J. Schumacher, “Fine-scale statistics of temperature and its derivatives in convective turbulence,” J. Fluid Mech. 611, 13–34 (2008).
  * Silano, Sreenivasan, and Verzicco (2010) G. Silano, K. R. Sreenivasan, and R. Verzicco, “Numerical simulations of Rayleigh–Bénard convection for Prandtl numbers between 10-1 and 104 and Rayleigh numbers between 105 and 109,” J. Fluid Mech. 662, 409–446 (2010).
  * Verma  _et al._ (2012) M. K. Verma, P. K. Mishra, A. Pandey, and S. Paul, “Scalings of field correlations and heat transport in turbulent convection,” Phys. Rev. E 85, 016310 (2012).
  * Pandey and Verma (2016) A. Pandey and M. K. Verma, “Scaling of large-scale quantities in Rayleigh-Bénard convection,” Phys. Fluids 28, 095105 (2016).
  * Pandey  _et al._ (2016) A. Pandey, A. Kumar, A. G. Chatterjee, and M. K. Verma, “Dynamics of large-scale quantities in Rayleigh-Bénard convection,” Phys. Rev. E 94, 053106 (2016).
  * Brown, Funfschilling, and Ahlers (2007) E. Brown, D. Funfschilling, and G. Ahlers, “Anomalous Reynolds-number scaling in turbulent Rayleigh–Bénard convection,” J. Stat. Mech. Theor. Exp. 2007, P10005 (2007).
  * Kraichnan (1962) R. H. Kraichnan, “Turbulent thermal convection at arbitrary prandtl number,” Phys. Fluids 5, 1374–1389 (1962).
  * Lohse and Toschi (2003) D. Lohse and F. Toschi, “Ultimate state of thermal convection,” Phys. Rev. Lett. 90, 034502 (2003).
  * Pawar and Arakeri (2016) S. S. Pawar and J. H. Arakeri, “Two regimes of flux scaling in axially homogeneous turbulent convection in vertical tube,” Phys. Rev. Fluids 1, 042401(R) (2016).
  * Schmidt  _et al._ (2012) L. E. Schmidt, E. Calzavarini, D. Lohse, and F. Toschi, “Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell,” J. Fluid Mech. 691, 52–68 (2012).
  * Roche  _et al._ (2001) P.-E. Roche, B. Castaing, B. Chabaud, and B. Hebral, “Observation of the 1/2 power law in Rayleigh-Bénard convection,” Phys. Rev. E 63, 045303(R) (2001).
  * Ahlers  _et al._ (2012) G. Ahlers, X. He, D. Funfschilling, and E. Bodenschatz, “Heat transport by turbulent Rayleigh-Bénard convection for $Pr\approx 0.8$ and $3\times 10^{12}\lesssim Ra\lesssim 10^{15}$: aspect ratio T = 0.50,” New J. Phys. 12, 103012 (2012).
  * He  _et al._ (2012) X. He, D. Funfschilling, H. Nobach, E. Bodenschatz, and G. Ahlers, “Transition to the Ultimate State of Turbulent Rayleigh-Bénard Convection,” Phys. Rev. Lett. 108, 024502 (2012).
  * Niemela  _et al._ (2000) J. J. Niemela, L. Skrbek, K. R. Sreenivasan, and R. J. Donnelly, “Turbulent convection at very high Rayleigh numbers,” Nature 404, 837–840 (2000).
  * Urban  _et al._ (2012) P. Urban, P. Hanzelka, T. Kralik, V. Musilová, A. Srnka, and L. Skrbek, “Effect of Boundary Layers Asymmetry on Heat Transfer Efficiency in Turbulent Rayleigh-Bernard Convection at Very High Rayleigh Numbers,” Phys. Rev. Lett. 109, 154301 (2012).
  * Grossmann and Lohse (2000) S. Grossmann and D. Lohse, “Scaling in thermal convection: a unifying theory,” J. Fluid Mech. 407, 27–56 (2000).
  * Grossmann and Lohse (2001) S. Grossmann and D. Lohse, “Thermal convection for large Prandtl numbers,” Phys. Rev. Lett. 86, 3316–3319 (2001).
  * Bhattacharya, Verma, and Samtaney (2021a) S. Bhattacharya, M. K. Verma, and R. Samtaney, “Revisiting Reynolds and Nusselt numbers in turbulent thermal convection,” Phys. Fluids 33, 015113 (2021a).
  * Bhattacharya  _et al._ (2018) S. Bhattacharya, A. Pandey, A. Kumar, and M. K. Verma, “Complexity of viscous dissipation in turbulent thermal convection,” Phys. Fluids 30, 031702 (2018).
  * Bhattacharya, Samtaney, and Verma (2019) S. Bhattacharya, R. Samtaney, and M. K. Verma, “Scaling and spatial intermittency of thermal dissipation in turbulent convection,” Phys. Fluids 31, 075104 (2019).
  * Parish and Duraisamy (2016) E. J. Parish and K. Duraisamy, “A paradigm for data-driven predictive modeling using field inversion and machine learning,” J. Comput. Phys. 305, 758–774 (2016).
  * Fonda  _et al._ (2019) E. Fonda, A. Pandey, J. Schumacher, and K. R. Sreenivasan, “Deep learning in turbulent convection networks,” Proc. Natl. Acad. Sci. U.S.A. 116, 8667–8672 (2019).
  * Pandey and Schumacher (2020) S. Pandey and J. Schumacher, “Reservoir computing model of two-dimensional turbulent convection,” Phys. Rev. Fluids 5, 113506 (2020).
  * Pandey, Schumacher, and Sreenivasan (2020) S. Pandey, J. Schumacher, and K. R. Sreenivasan, “A perspective on machine learning in turbulent flows,” J. of Turbulence 21, 567–584 (2020).
  * Brunton, Noack, and Koumoutsakos (2020) S. L. Brunton, B. R. Noack, and P. Koumoutsakos, “Machine Learning for Fluid Mechanics,” Annu. Rev. Fluid Mech. 52, 477–508 (2020).
  * Goodfellow, Bengio, and Courville (2016) I. Goodfellow, Y. Bengio, and A. Courville,  _Deep Learning_ (The MIT Press, Cambridge, Massachusetts, 2016).
  * Hastie, Tibshirani, and Friedman (2009) T. Hastie, R. Tibshirani, and J. Friedman,  _The Elements of Statistical Learning_ (Springer, New York, 2009).
  * Burkov (2019) A. Burkov,  _The Hundred-Page Machine Learning Book_ (Andriy Burkov, Québec, 2019).
  * Stevens  _et al._ (2013) R. J. A. M. Stevens, E. P. van der Poel, S. Grossmann, and D. Lohse, “The unifying theory of scaling in thermal convection: the updated prefactors,” J. Fluid Mech. 730, 295–308 (2013).
  * Lesieur (2008) M. Lesieur,  _Turbulence in Fluids_ (Springer-Verlag, Dordrecht, 2008).
  * Verma (2019) M. K. Verma,  _Energy trasnfers in Fluid Flows: Multiscale and Spectral Perspectives_ (Cambridge University Press, Cambridge, 2019).
  * Bhattacharya  _et al._ (2019) S. Bhattacharya, S. Sadhukhan, A. Guha, and M. K. Verma, “Similarities between the structure functions of thermal convection and hydrodynamic turbulence,” Phys. Fluids 31, 115107 (2019).
  * Scheel, Kim, and White (2012) J. D. Scheel, E. Kim, and K. R. White, “Thermal and viscous boundary layers in turbulent Rayleigh–Bénard convection,” J. Fluid Mech. 711, 281–305 (2012).
  * Shi, Emran, and Schumacher (2012) N. Shi, M. S. Emran, and J. Schumacher, “Boundary layer structure in turbulent Rayleigh–Bénard convection,” J. Fluid Mech. 706, 5–33 (2012).
  * Samuel  _et al._ (2021) R. J. Samuel, S. Bhattacharya, A. Asad, S. Chatterjee, M. K. Verma, R. Samtaney, and S. F. Anwer, “SARAS: A general-purpose PDE solver for fluid dynamics,” J. Open Source Softw. 6, 2095 (2021).
  * Verma  _et al._ (2020) M. K. Verma, R. J. Samuel, S. Chatterjee, S. Bhattacharya, and A. Asad, “Challenges in fluid flow simulations using exascale computing,” S.N. Comput. Sci. 1, 178 (2020).
  * Bhattacharya, Verma, and Samtaney (2021b) S. Bhattacharya, M. K. Verma, and R. Samtaney, “Prandtl number dependence of the small-scale properties in turbulent Rayleigh-Bénard convection,” Phys. Rev. Fluids 6, 063501 (2021b).
  * Jekel and Venter (2019) C. F. Jekel and G. Venter, “Pwlf: a python library for fitting 1d continuous piecewise linear functions,” URL: https://github. com/cjekel/piecewise_linear_fit_py (2019).
  * Frank  _et al._ (2009) E. Frank, M. Hall, G. Holmes, R. Kirkby, B. Pfahringer, I. H. Witten, and L. Trigg, “Weka-a machine learning workbench for data mining,” in  _Data mining and knowledge discovery handbook_ (Springer, 2009) pp. 1269–1277.
  * Hornik, Stinchcombe, and White (1989) K. Hornik, M. Stinchcombe, and H. White, “Multilayer feedforward networks are universal approximators,” Neural Netw. 2, 359–366 (1989).
  * Gulli and Pal (2017) A. Gulli and S. Pal,  _Deep Learning with Keras_ (Packt Publishing, Birmingham, 2017).
  * Bishop (2006) C. M. Bishop,  _Pattern Recognition and Machine Learning_ (Springer, Singapore, 2006).
  * Zhu  _et al._ (2018) X. Zhu, V. Mathai, R. J. A. M. Stevens, R. Verzicco, and D. Lohse, “Transition to the Ultimate Regime in Two-Dimensional Rayleigh-Bénard Convection,” Phys. Rev. Lett. 120, 144502 (2018).
  * Iyer  _et al._ (2020) K. P. Iyer, J. D. Scheel, J. Schumacher, and K. R. Sreenivasan, “Classical 1/3 scaling of convection holds up to Ra $=10^{15}$,” Proc. Natl. Acad. Sci. U.S.A. 117, 7594–7598 (2020).