Natalia Tomashenko

Salima Mdhaffar

Marc Tommasi

Yannick Estève

Jean-Franc ¸ois Bonastre

PRIVACY ATTACKS FOR AUTOMATIC SPEECH RECOGNITION ACOUSTIC MODELS IN A FEDERATED LEARNING FRAMEWORK

Keywords: Privacy, federated learning, acoustic models, attack models, speech recognition, speaker verification

This paper investigates methods to effectively retrieve speaker information from the personalized speaker adapted neural network acoustic models (AMs) in automatic speech recognition (ASR). This problem is especially important in the context of federated learning of ASR acoustic models where a global model is learnt on the server based on the updates received from multiple clients. We propose an approach to analyze information in neural network AMs based on a neural network footprint on the so-called Indicator dataset. Using this method, we develop two attack models that aim to infer speaker identity from the updated personalized models without access to the actual users' speech data. Experiments on the TED-LIUM 3 corpus demonstrate that the proposed approaches are very effective and can provide equal error rate (EER) of 1-2%.

INTRODUCTION

Federated learning (FL) for automatic speech recognition (ASR) has recently become an active area of research [START_REF] Cui | Federated acoustic modeling for automatic speech recognition[END_REF][START_REF] Dimitriadis | A federated approach in training acoustic models[END_REF][START_REF] Mdhaffar | Study on acoustic model personalization in a context of collaborative learning constrained by privacy preservation[END_REF][START_REF] Guliani | Training speech recognition models with federated learning: A quality/cost framework[END_REF][START_REF] Dimitriadis | Federated transfer learning with dynamic gradient aggregation[END_REF][START_REF] Yu | Federated learning in ASR: Not as easy as you think[END_REF]. To preserve the privacy of the users' data in the FL framework, the model is updated in a distributed fashion instead of communicating the data directly from clients to a server.

Privacy is one of the major challenges in FL [START_REF] Li | Federated learning: Challenges, methods, and future directions[END_REF][START_REF] Viraaji Mothukuri | A survey on security and privacy of federated learning[END_REF]. Sharing model updates, i.e. gradient information, instead of raw user data aims to protect user personal data that are processed locally on devices. However, these updates may still reveal some sensitive information to a server or to a third party [START_REF] Geiping | Inverting gradients-how easy is it to break privacy in federated learning?[END_REF][START_REF] Carlini | The secret sharer: Evaluating and testing unintended memorization in neural networks[END_REF]. According to recent research, FL has various privacy risks and may be vulnerable to different types of attacks, i.e. membership inference attacks [START_REF] Truex | Demystifying membership inference attacks in machine learning as a service[END_REF] or generative adversarial network (GAN) inference attacks [START_REF] Wang | Beyond inferring class representatives: User-level privacy leakage from federated learning[END_REF]. Techniques to enhance the privacy in a FL framework are mainly based on two categories [START_REF] Viraaji Mothukuri | A survey on security and privacy of federated learning[END_REF]: secure multiparty computation [START_REF] Bonawitz | Practical secure aggregation for federated learning on user-held data[END_REF] and differential privacy [START_REF] Dwork | Differential privacy[END_REF]. Encryption methods [START_REF] Manas | Privacy-preserving speech processing: cryptographic and string-matching frameworks show promise[END_REF][START_REF] Smaragdis | A framework for secure speech recognition[END_REF] such as fully homomorphic encryption [START_REF] Smaragdis | A framework for secure speech recognition[END_REF] and secure multiparty computation perform computation in the encrypted domain. These methods increase computational complexity. In a FL framework, this increase is not so significant compared to standard centralized training, since only the transmitted parameters need to be encrypted instead of large amounts of data, however with an increased number of participants, computational complexity becomes a critical issue. Differential privacy methods preserve privacy by adding noise to users' parameters [START_REF] Dwork | Differential privacy[END_REF][START_REF] Xie | Differentially private generative adversarial network[END_REF], however such solutions may degrade learning performance due to the uncertainty they introduce into the parameters.

Alternative methods to privacy protection for speech include deletion methods [START_REF] Cohen-Hadria | Voice anonymization in urban sound recordings[END_REF] that are meant for ambient sound analysis, and anonymization [START_REF] Tomashenko | Introducing the VoicePrivacy initiative[END_REF][START_REF] Tomashenko | The VoicePrivacy 2020 Challenge: Results and findings[END_REF] that aims to suppress personally identifiable information in the speech signal keeping unchanged all other attributes. These privacy preservation methods can be combined and integrated in a hybrid fashion into a FL framework.

Despite the recent interest in FL for ASR and other speechrelated tasks such as keyword spotting [START_REF] Leroy | Federated learning for keyword spotting[END_REF][START_REF] Hard | Training keyword spotting models on non-iid data with federated learning[END_REF], emotion recognition [START_REF] Latif | Federated learning for speech emotion recognition applications[END_REF], and speaker verification [START_REF] Granqvist | Improving on-device speaker verification using federated learning with privacy[END_REF], there is a lack of research on vulnerability of ASR acoustic models (AMs) to privacy attacks in a FL framework. In this work, we make a step towards this direction by analyzing speaker information that can be retrieved from the personalized AM locally updated on the user's data. To achieve this goal, we developed two privacy attack models that operate directly on the updated model parameters without access to the actual user's data. Parameters of neural network (NN) personalized AMs contain a wealth of information about the speakers [START_REF] Mdhaffar | Retrieving speaker information from personalized acoustic models for speech recognition[END_REF]. In this paper, we propose novel methods to efficiently and easily retrieve speaker information from the adapted AMs. The main idea of the proposed methods is to use an external Indicator dataset to analyze the footprint of AMs on this data. Another important contribution of this work is understanding how the speaker information is distributed in the adapted NN AMs.

This paper is structured as follows. Section 2 briefly introduces a considered FL framework for AM training. Section 3 describes the privacy preservation scenario and proposes two attack models. Experimental evaluation is presented in Section 4. We conclude in Section 5.

FEDERATED LEARNING FOR ASR ACOUSTIC MODELS

We consider a classical FL scenario where a global NN AM is trained on a server from the data stored locally on multiple remote devices [START_REF] Li | Federated learning: Challenges, methods, and future directions[END_REF]. The training of the global model is performed under the constraint that the training speech data are stored and processed locally on the user devices (clients), while only model updates are transmitted to the server from each client. The global model is learnt on the server based on the updates received from multiple clients.

The FL in a distributed network of clients is illustrated in Figure 1.

First, an initial global speech recognition AM Wg is distributed to the group of devices of N users (speakers). Then, the initial global model is run on every user si (i ∈ 1..N ) device and updated locally on the private user data. The updated models Ws i are then transmitted to the server where they are aggregated to obtain a new global model W * g . Typically, the personalized updated models are aggregated using federated averaging and its variations [START_REF] Bonawitz | Practical secure aggregation for federated learning on user-held data[END_REF][START_REF] Mcmahan | Communication-efficient learning of deep networks from decentralized data[END_REF]. Then, the updated global model W * g is shared with the clients. The process restarts and loops until convergence or after a fixed number of rounds. The utility and training efficiency of the FL AMs have been successfully studied in recent works [START_REF] Cui | Federated acoustic modeling for automatic speech recognition[END_REF][START_REF] Dimitriadis | A federated approach in training acoustic models[END_REF][START_REF] Mdhaffar | Study on acoustic model personalization in a context of collaborative learning constrained by privacy preservation[END_REF][START_REF] Guliani | Training speech recognition models with federated learning: A quality/cost framework[END_REF][START_REF] Dimitriadis | Federated transfer learning with dynamic gradient aggregation[END_REF][START_REF] Yu | Federated learning in ASR: Not as easy as you think[END_REF], and these topics are beyond the scope of the current paper. Alternatively, we focus on the privacy aspect of this framework.

Server

ATTACK MODELS

In this section, we describe the privacy preservation scenario and present two attack models.

Privacy preservation scenario

Privacy preservation is formulated as a game between users who share some data and attackers who access this data or data derived from it and aim to infer information about the users [START_REF] Tomashenko | Introducing the VoicePrivacy initiative[END_REF]. To preserve the user data, in FL, there is no speech data exchange between a server and clients, only model updates are transmitted between the clients and server (or between some clients). Attackers aim to attack users using information owned by the server. They can get access to some updated personalized models.

In this work, we assume that an attacker has access to the following data:

• An initial global NN AM Wg;

• A personalized model Ws of the target speaker s who is enrolled in the FL system. The corresponding personalized model was obtained from the global model Wg by fine-tuning Wg using speaker data. We consider this model as enrollment data for an attacker.

• Other personalized models of non-target and target speakers: Ws 1 ,...,Ws N . We will refer to these models as test trial data.

The attacker's objective is to conduct an automatic speaker verification (ASV) task by using the enrollment data model in the form of Ws and test trial data in the form of models Ws 1 ,...,Ws N .

Attack models

The motivation of the proposed approaches is based on the hypothesis that we can capture information about the identity of speaker s from the corresponding speaker-adapted model Ws and the global model Wg by comparing the outputs of these two neural AMs taken from hidden layers h on some speech data. We will refer to this speech data as Indicator data. Note, that the Indicator data is not related to any test or AM training data and can be chosen arbitrarily from any speakers.

Attack model A1

The ASV task with the proposed attack model is performed in several steps as illustrated in Figure 2. Let denote a set of utterances in the Indicator dataset I as u1, . . . , uJ ∈ I; a sequence of vectors in utterance uj as {u

W h s i (uj) = {w h,t s i ,j } T j t=1 and W h g (uj) = {w h,t g,j } T j
t=1 , and per-frame differences between corresponding outputs:

∆ h s i (uj) = {∆ h,t s i ,j } T j t=1 , (1)
where ∆ h,t s i ,j = w h,t s i ,j -w h,t g,j , t ∈ 1..Tj. 2. For each personalized model, we compute mean and standard deviation vectors for ∆ h,t s i ,j over all speech frames in the Indicator dataset I:

µ h s i = I j=1 T j t=1 ∆ h,t s i ,j I j=1 Tj , (2)
σ h s i = I j=1 T j t=1 (∆ h,t s i ,j -µ h s i ) 2 I j=1 Tj 1 2 (3)
3. For a pair of personalized models Ws i and Ws k , we compute a similarity score ρ at hidden level h on the Indicator dataset based on the L2-normalised Euclidean distance between the corresponding vector pairs for means and standard deviations:

ρ(W h s i , W h s k ) = αµ ∥µ h s i -µ h s k ∥2 ∥µ h s i ∥2∥µ h s k ∥2 + ασ ∥σ h s i -σ h s k ∥2 ∥σ h s i ∥2∥σ h s k ∥2 , (4)
where αµ, ασ are fixed parameters in all experiments.

4. Given similarity scores for all matrix pairs, we can complete a speaker verification task based on these scores.

Pool of personalized models

Indicator speech data

𝑊 𝑠 1 ℎ (𝑢) 𝑊 𝑔 ℎ (𝑢) 𝑊 𝑠 𝑁 𝑊 𝑠 1 𝑊 𝑠 2 𝑊 𝑔 𝑢 𝑢 𝑊 𝑠 𝑁 ℎ (𝑢) Δ 𝑠 1 ℎ (𝑢) Δ 𝑠 𝑁 ℎ (𝑢) [μ 𝑠 1 ℎ , 𝜎 𝑠 1 ℎ ] [μ 𝑠 𝑁 ℎ , 𝜎 𝑠 𝑁 ℎ ]
… Fig. 2. Statistic computation for the attack model A1.

Attack model A2

For the second attack model, we train a NN model as shown in Fig- ure 3. This NN model uses personalized and global models and the speech Indicator dataset for training. It is trained to predict a speaker identity provided the corresponding personalized model. When the model is trained, we use it in evaluation time to extract speaker embeddings similarly to x-vectors and apply probabilistic linear discriminant analysis (PLDA) [START_REF] Snyder | X-vectors: Robust DNN embeddings for speaker recognition[END_REF][START_REF] Ioffe | Probabilistic linear discriminant analysis[END_REF]. As shown in Figure 3, the model consists of two parts (frozen and trained). The outputs of the frozen part are ∆ h s i sequences of vectors computed per utterance of the Indicator data as defined in Formula (1). For every personalized model Ws i , we compute ∆ h s i for all the utterances u of the Indicator corpus; then ∆ h s i (u) is used as input to the second (trained) part of the NN which comprises several time delay neural network (TDNN) layers [START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF] and one statistical pooling layer. The experiments were conducted on the speaker adaptation partition of the TED-LIUM 3 corpus [START_REF] Franc ¸ois Hernandez | TED-LIUM 3: twice as much data and corpus repartition for experiments on speaker adaptation[END_REF]. This publicly available data set contains TED talks that amount to 452 hours speech data in English from about 2K speakers, 16kHz. Similarly to [START_REF] Mdhaffar | Study on acoustic model personalization in a context of collaborative learning constrained by privacy preservation[END_REF], we selected from the TED-LIUM 3 training dataset three datasets: Train-G, Part-1, Part-2 with disjoint speaker subsets as shown in Table 1. The Indicator dataset was used to train an attack model. It is comprised of 320 utterances selected from all 32 speakers of test and development datasets of the TED-LIUM 3 corpus. The speakers in the Indicator dataset are disjoint from speakers in Train-G, Part-1, and Part-2.

For each speaker in the Indicator dataset we select 10 utterances. The size of the Indicator dataset is 32 minutes. The Train-G dataset was used to train an initial global AM Wg. Part-1 and Part-2 were used to obtain two sets of personalized models. 1

ASR acoustic models

The ASR AMs have a TDNN model architecture [START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF] and trained using the Kaldi speech recognition toolkit [START_REF] Povey | The Kaldi speech recognition toolkit[END_REF]. 40-dimensional Mel-frequency cepstral coefficients (MFCCs) without cepstral 1 Data partitions and scripts will be available online: https:// github.com/Natalia-T/privacy-attacks-asr 1. Data sets statistics truncation appended with 100-dimensional i-vectors were used as the input into the NNs. Each model has thirteen 512-dimensional hidden layers followed by a softmax layer where 3664 triphone states were used as targets 2 . The initial global model Wg was trained using the lattice-free maximum mutual information (LF-MMI) criterion with a 3-fold reduced frame rate as in [START_REF] Povey | Purely sequence-trained neural networks for ASR based on latticefree MMI[END_REF]. The two types of speech data augmentation strategies were applied for the training and adaptation data: speed perturbation (with factors 0.9, 1.0, 1.1) and volume perturbation, as in [START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF]. Each model has about 13.8M parameters. The initial global model Wg was trained on the Train-G. Personalized models Ws i were obtained by fine-tuning all the parameters of Wg on the speakers' data from Part-1 and Part-2 as described in [START_REF] Mdhaffar | Study on acoustic model personalization in a context of collaborative learning constrained by privacy preservation[END_REF]. For all personalized speaker models, we use approximately the same amount of speech data to perform fine-tuning (speaker adaptation) -about 4 minutes per model. For most of the speakers (564 in Part-1, 463 in Part-2) we obtained two different personalized models (per speaker) on disjoint adaptation subsets, for the rest speakers we have adaptation data only for one model.

Attack models

We investigate two approaches for attack models: A1 -a simple approach based on the comparative statistical analysis of the NN outputs and the associated similarity score between personalized models, and A2 -a NN based approach. For the test target trials, we use comparisons between different personalized models of the same speakers (564 in Part-2 and 1027 in the Part-1+Part-2), and for the non-target trials we randomly selected 10K pairs of models from different speakers in a corresponding dataset.

Attack model A1

The first attack model was applied as described in Section 3.2.1. The parameters αµ, ασ in Formula (4) equal to 1 and 10 respectively. This model was evaluated on two datasets of personalized models corresponding to Part-2 and combined Part-1+Part-2 datasets. The Indicator dataset is the same in all experiments.

Attack model A2

For training the attack model A2, we use 1300 personalized speaker models corresponding to 736 unique speakers from Part-1. When we applied the frozen part of the architecture shown in Figure 3 to the 32-minute Indicator dataset for each speaker model in Part-1, we obtained the training data with the amount corresponding to about 693h (32×1300). The trained part of the NN model, illustrated in Figure 3, has a similar topology to a conventional x-vector extractor [START_REF] Snyder | X-vectors: Robust DNN embeddings for speaker recognition[END_REF]. However, unlike the standard NN x-vector extractor, that is trained to predict speaker id-s by the input speech segment, our proposed model learns to predict a speaker identity from the W h s part of a speaker personalized model. We trained 2 attack models corresponding to the two values of parameter h ∈ {1, 5} -a hidden layer in the ASR neural AMs at which we compute the activations.

Values h were chosen based on the results for the attack model A1.

The output dimension of the frozen part is 512. The frozen part is followed by the trained part that consists of seven hidden TDNN layers and one statistical pooling layer introduced after the fifth TDNN layer. The output is a softmax layer with the targets corresponding to speakers in the pool of speaker personalized models (number of unique speakers in Part-1).

Results

The attack models were evaluated in terms of equal error rate (EER).

Denoting by Pfa(θ) and Pmiss(θ) the false alarm and miss rates at threshold θ, the EER corresponds to the threshold θEER at which the two detection error rates are equal, i.e., EER = Pfa(θEER) = Pmiss(θEER).

Results for the attack model A1 are shown in Figure 4 for Part-2 and combined Part-1 and Part-2 datasets. Speaker information can be captured for all values h with varying success: EER ranges from 0.86% (for the first hidden layer) up to 20.51% (for the top hidden layer) on Part-2. For the Part-1+Part-2 we observe similar results. To analyze the impact of each component in Formula (4) on the ASV performance, we separately compute similarity score ρ either using only means (ασ = 0) or only standard deviations (αµ = 0). Results on the Part-2 dataset are shown in Figure 5. Black bars correspond to ασ = 0 when only means were used to compute similarity scores ρ between personalized models. Blue bars represent results for αµ = 0 when only standard deviations were used to compute ρ. Orange bars correspond to the combined usage of means and standard deviations as in Figure 4 (αµ = 1, ασ = 10). The impact of each component in the sum changes for different hidden layers. When we use only standard deviations, we observe the lowest EER on the first layer. In case of using only means, the first layer is, on the contrary, one of the least informative for speaker verification. For all other layers, combination of means and standard deviations provided superior results over the cases when only one of these components were used. These surprising results for the first hidden layer could possibly be explained by the fact that the personalized models incorporate i-vectors in their inputs, thus the speaker information can be easily learnt at this level of the NN. We plan to investigate this phenomena in detail in our future research.

We choose two values h ∈ {1, 5} which demonstrate promising results for the model A1, and use the corresponding outputs to train two attack models with the configuration A2. The comparative results for the two attack models are presented in Table 2. For h = 5, the second attack model provides significant improvement in performance over the first one and reduces EER from 7% down to 2%. For h = 1, we could not obtain any improvement by training a NN based attack model: the results for A2 in this case are worse than for the simple approach A1. One explanation for this phenomenon could be the following. The first layers of the AMs provide highly informative features for speaker classification, however, training the proposed NN model on these features results in overfitting because training criterion of the NN is speaker accuracy, but not the target EER metric, and the number of targets is relatively small, hence, the NN overfits to classify the seen speakers in the training dataset. Table 2. EER, % evaluated on Part-2, h -indicator of a hidden layer

CONCLUSIONS

In this work, we focused on the privacy protection problem for ASR AMs trained in a FL framework. We explored to what extent ASR AMs are vulnerable to privacy attacks. We developed two attack models that aim to infer speaker identity from the locally updated personalized models without access to any speech data of the target speakers. One attack model is based on the proposed similarity score between personalized AMs computed on some external Indicator dataset, and another one is a NN model. We demonstrated on the TED-LIUM 3 corpus that both attack models are very effective and can provide EER of about 1% for the simple attack model A1 and 2% for the NN attack model A2. Another important contribution of this work is the finding that the first layer of personalized AMs contains a large amount of speaker information that is mainly contained in the standard deviation values computed on Indicator data. This interesting property of NN adapted AMs opens new perspectives also for ASV, and in future work, we plan to use it for developing an efficient ASV system.

Fig. 1 .

 1 Fig. 1. Federated learning in a distributed network of clients: 1) Download of the global model Wg by clients. 2) Speaker adaptation of Wg on the local devices using user private data. 3) Collection and aggregation of multiple personalized models Ws 1 ,...,Ws N on the sever. 4) Sharing the resulted model W * g with the clients.

( 1 )

 1 j , . . . , u (T j ) j }; a set of personalized models as Ws 1 , . . . , Ws N ∈ W; and an identifier of a hidden layer in the global or personalized AM as h. 1. ∀ Ws i ∈ W, ∀ uj ∈ I we compute activation values from the layer h for model pairs:

Fig. 3 .

 3 Fig. 3. Training a speaker embedding extractor for the attack model A2.4. EXPERIMENTS

Fig. 4 .

 4 Fig. 4. EER, % for the attack model A1 depending on the hidden layer h (in Wg and Ws i ) which was used to compute outputs, evaluated on Part-2 and on the combined Part-1+Part-2 dataset.

Fig. 5 .

 5 Fig. 5. EER, % for the attack model A1 depending on the hidden layer h, evaluated on Part-2 dataset. µ+σ -both means and standard deviations were used to compute similarity score ρ; µ -only means; and σ -only standard deviations were used.

Following the notation from[START_REF] Peddinti | A time delay neural network architecture for efficient modeling of long temporal contexts[END_REF], the model configuration can be described as follows: {-1,0,1} × 6 layers; {-3,0,3} × 7 layers.

This work was funded by VoicePersonae and DEEP-PRIVACY (ANR-18-CE23-0018) projects. This work was performed using HPC resources from GENCI-IDRIS (Grant 2021-AD011013331).


M Delort-Laval

L Soucasse
email: laurent.soucasse@centralesupelec.fr

Ph Rivière

A Soufiani

Rayleigh-Bénard convection in a cubic cell under the effects of gas radiation up to Ra=10 9

Keywords: turbulent convection, radiative transfer, cubic Rayleigh-Bénard cell

This paper investigates radiative transfer effects on Rayleigh-Bénard convection in a cubic cell over a large range of Rayleigh numbers, from Ra = 10 3 (below the onset of convection) to Ra = 10 9 in the turbulent regime. Coupled direct numerical simulations are carried out for a radiating air/H 2 O/CO 2 mixture at room temperature, using a Chebyshev spectral method for the flow and a ray-tracing method for the radiation field. For the highest Rayleigh numbers, a subgrid model is used to account for the radiation of the smallest, non-optically thin, turbulent scales. Symmetry and time-averaging (for unsteady solutions) are applied to compare coupled and uncoupled results, regardless of the multiple flow configurations that may be obtained. At low Rayleigh number, the potential energy decreases, and the onset of convection is delayed when radiation is taken into account. However, once convection settles, the potential energy increases with radiation, leading to a higher convective flux in the core and a higher kinetic energy. Specific contributions of radiative transfer to the potential energy balance and to the thermal energy balance are highlighted. It is also shown that the ratio of radiative and convective source terms in the energy balance roughly scales as Ra -1/2 and that radiative transfer effects weaken at high Rayleigh numbers. Finally, radiative transfer effects on turbulence budgets of mechanical and thermal fluctuations are analysed in the range 10 7 ≤ Ra ≤ 10 9 . The magnitude of each term of these budgets is stronger when radiation is taken into account. However, radiative dissipation has little influence on the temperature fluctuation budget.

Introduction

Rayleigh-Bénard (RB) convection, in which a fluid layer is heated from below and cooled from above, remains a challenging problem for both fundamental studies and for natural or engineering applications. These applications include atmospheric physics, convection in the Earth mantle, thermal management of buildings, and many industrial processes. A layer of transparent fluid heated from below is known to become unstable when the Rayleigh number reaches a critical value of about Ra cr =1708 for an infinite geometry in the directions perpendicular to gravity, and for rigid walls [START_REF] Chandrasekhar | Hydrodynamic and Hydromagnetic Stability[END_REF]. When the fluid is confined in a closed cavity, the critical Rayleigh number increases due to friction on the lateral walls. Fluid confinement also leads to a great variety of possible flow structures and corresponding heat transfer rates just above the transition, due to the symmetries of the system. Several studies have revealed a variety of steady and unsteady flows at moderate Rayleigh numbers for a cubic cavity heated from below. Depending on the boundary conditions on the lateral walls, up to seven or nine (for insulated lateral walls [START_REF] Pallares | Flow transitions in laminar Rayleigh-Bénard convection in a cubical cavity at moderate Rayleigh numbers[END_REF][START_REF] Puigjaner | Stability analysis of the flow in a cubical cavity heated from below[END_REF][START_REF] Bousset | Three-dimensional convection regimes in a cubical cavity[END_REF]) or even more (for conducting lateral walls [START_REF] Puigjaner | Bifurcation analysis of steady Rayleigh-Bénard convection in a cubical cavity with conducting sidewalls[END_REF]) structures have been observed and analyzed using numerical simulations or continuation technique identification. The existence of multiple solutions at Ra = 10 5 has also been observed experimentally depending on the initial conditions [START_REF] Leong | Experimental Nusselt numbers for a cubical-cavity benchmark problem in natural convection[END_REF]. Different sets of the Nusselt number have indeed been identified for the same Rayleigh number. Also, benchmark numerical simulations in the range 10 5 -10 8 of the Rayleigh number have shown an important disparity between the results from several numerical codes [START_REF] Pepper | Summary of benchmark numerical studies for 3-D natural convection in an air-filled enclosure[END_REF].

Subsequently, several studies have shown that, for confined enclosures and intermediate Ra numbers, typically in the range 10 6 -10 9 , the weakly turbulent RB flows are characterized by a Large-Scale Circulation (LSC) with low-frequency intermittent phenomena such as reversals and reorientations. Azimuthal rotation of the LSC and sudden reversals in cylindrical enclosures have been observed in experimental studies [START_REF] Sreenivasan | Mean wind and its reversal in thermal convection[END_REF][START_REF] Brown | Reorientation of the large-scale circulation in turbulent Rayleigh-Bénard convection[END_REF] and by numerical simulations [START_REF] Benzi | Numerical simulations of flow reversal in Rayleigh-Bénard convection[END_REF][START_REF] Mishra | Dynamics of reorientations and reversals of large-scale flow in Rayleigh-Bénard convection[END_REF]. Recently, large-scale flow reversals in circular enclosures have been demonstrated both experimentally [START_REF] Wang | Mechanism of large-scale flow reversals in turbulent thermal convection[END_REF] and numerically [START_REF] Xu | Tristable flow states and reversal of the large-scale circulation in two-dimensional circular convection cells[END_REF]. In cubic enclosures, LSC is generally organized in the diagonal vertical planes and reorientations from one plane to another have been observed [START_REF] Vasiliev | High Rayleigh number convection in a cubic cell with adiabatic sidewalls[END_REF][START_REF] Foroozani | Reorientations of the large-scale flow in turbulent convection in a cube[END_REF].

Several models have been developed to predict LSC low-frequency dynamics, either based on stochastic phenomenological approaches [START_REF] Araujo | Wind reversals in turbulent Rayleigh-Bénard convection[END_REF][START_REF] Brown | Large-scale circulation model for turbulent Rayleigh-Bénard convection[END_REF][START_REF] Brown | Azimuthal asymmetries of the large-scale circulation in turbulent Rayleigh-Bénard convection[END_REF][START_REF] Bai | Ability of a low-dimensional model to predict geometry-dependent dynamics of large-scale coherent structures in turbulence[END_REF][START_REF] Ji | Low-dimensional model of the large-scale circulation of turbulent Rayleigh-Bénard convection in a cubic container[END_REF] or modal decomposition [START_REF] Bailon-Cuba | Aspect ratio dependence of heat transfer and large-scale flow in turbulent convection[END_REF][START_REF] Chandra | Dynamics and symmetries of flow reversals in turbulent convection[END_REF][START_REF] Giannakis | Koopman analysis of the long-term evolution in a turbulent convection cell[END_REF][START_REF] Soucasse | Proper orthogonal decomposition analysis and modelling of large-scale flow reorientations in a cubic Rayleigh-Bénard cell[END_REF][START_REF] Xu | Correlation of internal flow structure with heat transfer efficiency in turbulent Rayleigh-Bénard convection[END_REF]. The development of such predictive models is a very attractive topic since Direct Numerical Simulations (DNS) must be run for very long integration times to capture these low-frequency phenomena.

The brief and non-exhaustive picture of the state of the art described above is related to RB convection in non-radiating fluids (transparent or opaque). However, it is well known that the small amounts of water vapour and carbon dioxide in ambient air are responsible of significant radiative transfer that could affect the mean and fluctuating temperature fields and the resulting thermal convection. This radiative transfer is at the origin of the greenhouse effect in Earth atmosphere. The effects of radiation on RB convection, with application to atmospheric physics, have been investigated in the pioneering studies [START_REF] Spiegel | The smoothing of temperature fluctuations by radiative transfer[END_REF][START_REF] Spiegel | The convective instability of a radiating fluid layer[END_REF][START_REF] Gille | Convection in a radiating gas[END_REF]. Using linear stability analyses, radiative transfer was shown to delay the onset of instability due to two mechanisms, namely, the homogenization of the temperature field in the core of the fluid layer (effect on the base state) and the damping of temperature disturbances [START_REF] Bdéoui | The onset of Rayleigh-Bénard instability in molecular radiating gases[END_REF]. Numerical simulations [START_REF] Bdéoui | A numerical study of Rayleigh-Bénard convection in radiating gases[END_REF][START_REF] Lan | Stability analysis for three-dimensional Rayleigh-Bénard convection with radiatively participating medium using spectral methods[END_REF] and experimental studies [START_REF] Gille | Convection in a radiating gas[END_REF] have also confirmed this behavior. However, the above mentioned studies were limited to the first transition to the convection regime or addressed the flow structure and heat transfer just above this transition.

For higher Rayleigh numbers, although several studies have addressed the effects of radiation on natural convection in differentially heated cavities [START_REF] Soucasse | Transitional regimes of natural convection in a differentially heated cavity under the effects of wall and molecular gas radiation[END_REF][START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF][START_REF] Kogawa | Influence of radiation effect on turbulent natural convection in cubic at normal temperature atmospheric gas[END_REF][START_REF] Kogawa | Experimental evaluation of thermal radiation effects on natural convection with a Rayleigh number of 10 8 -10 9 by using an interferometer[END_REF], there are very few studies related to radiation effects on RB convection in chaotic and turbulent regime. First attempts have been performed in Refs. [START_REF] Sakurai | Radiative heat transfer analysis in a turbulent natural convection obtained from direct numerical simulation[END_REF][START_REF] Sakurai | Radiation effects on mixed turbulent natural and forced convection in a horizontal channel using direct numerical simulation[END_REF] for Rayleigh numbers of about 10 6 -10 7 . However, no full coupling was considered in Ref. [START_REF] Sakurai | Radiative heat transfer analysis in a turbulent natural convection obtained from direct numerical simulation[END_REF] and, in Ref. [START_REF] Sakurai | Radiation effects on mixed turbulent natural and forced convection in a horizontal channel using direct numerical simulation[END_REF], devoted to mixed convection, the gray gas assumption was used for a non-confined thin layer of radiating fluid. A preliminary investigation [START_REF] Soucasse | Effects of molecular gas radiation on Rayleigh-Bénard convection in a 3D cubical cavity[END_REF] of radiation effects on RB convection in a cubic cell was limited to Ra = 10 7 and, due to limitations in computational resources, the integration time was not sufficient to capture LSC reorientations. More recently, low order models based on Proper Orthogonal Decomposition (POD) have been developed to analyse these LSC reorientations with and without radiation [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF][START_REF] Soucasse | Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different rayleigh numbers[END_REF] in the range 10 6 Ra 10 8 . Similar first eigenmodes were found with and without radiation but with higher energies (eigenvalues) when radiation was taken into account. It was found also that the frequency of reorientation events tends to increase with radiative transfer. However, the effects of radiation on heat transfer and flow dynamics were not fully analyzed in these studies.

The aim of this paper is to analyse molecular radiation effects on heat transfer and flow dynamics in a large range of the Rayleigh number, from 10 3 to 10 9 , in a cubic cavity heated from below and cooled from the top at constant temperatures, with insulated lateral walls. Mixture of air with small amounts of H 2 O and CO 2 as radiating gases, typical of ambient humid air conditions, is considered. Given the variety of flow patterns (and associated reorientations in unsteady regime) that may be observed in the large range of Rayleigh numbers studied, we will focus the analysis on symmetry-averaged and time-averaged (for unsteady solutions) quantities. The governing equations and the employed numerical methods for DNS are described in Sec. 2. The symmetries of the considered problem are also discussed in this section. The effects of radiation on mean flow fields and on mean conductive, convective and radiative fluxes are discussed in Sec. 3. Then we analyse the budgets of kinetic energy, potential energy, and thermal energy (square of the temperature differences) in Sec. 4. Finally, in Sec. 5, second order statistics are processed in the range 10 7 Ra 10 9 to analyse the budgets of turbulent kinetic energy and of the variance of temperature fluctuations.

Direct numerical simulations

Problem set-up and governing equations

The studied configuration is displayed in Fig. 1. We consider a cubic cavity of size L, heated from below and cooled from above. Top and bottom walls are isothermal at temperature T cold and T hot , and are black (emissivity ε = 1). The four side walls are adiabatic and perfectly diffuse reflecting (emissivity ε = 0). The cavity is filled with a radiating air/H 2 O/CO 2 mixture at a mean temperature T 0 .

In the case of a non radiating fluid the problem would be fully controlled by two parameters: the Prandtl number Pr = ν f /a and the Rayleigh number

Ra = gβ∆T L 3 ν f a , ( 1
)
where g is the gravitational acceleration, β = 1/T 0 is the thermal expansion coefficient, ∆T = T hot -T cold is the temperature difference, ν f is the kinematic viscosity and a is the thermal diffusivity. In the presence of gas radiation, it is not possible to define a limited set of dimensionless parameters due to the complex structure of molecular absorption spectra. In order to make our study relevant for building applications, we consider an air/H 2 O/CO 2 mixture at a mean temperature of T 0 = 300 K, at atmospheric pressure, of molar composition X H 2 O = 0.02 and X CO 2 = 0.001. The cavity size is set to L = 1 m and we vary the Rayleigh number by changing the temperature difference ∆T , other parameters being fixed. Thermophysical properties are assumed to be uniform (low temperature differences), not affected by the small amount of water vapour and carbon dioxide, and constant at all Rayleigh numbers (thermal conductivity λ = 0.0263 W m -1 K -1 , thermal diffusivity a = 2.25 × 10 -5 m 2 s -1 , Prandtl number Pr = 0.707). Mass, momentum and energy balance write under Boussinesq approximation

∂u i ∂x i = 0, (2)
∂u i ∂t + u j ∂u i ∂x j = - ∂p ∂x i + Prθδ i3 + Pr √ Ra ∂ 2 u i ∂x j ∂x j , (3)
∂θ ∂t + u i ∂θ ∂x i = 1 √ Ra ∂ 2 θ ∂x i ∂x i + P rad , (4)
where x i is the i th dimensionless Cartesian coordinate, u i is the i th dimensionless velocity vector component, p is the dimensionless motion pressure, θ = (T -T 0 )/∆T is the dimensionless temperature and P rad is the dimensionless radiative power. Eqs. ( 2)-( 4) are made dimensionless using the reference time L 2 /(a √ Ra) and the reference length L. A non-slip boundary condition is imposed on the velocity at the six walls and the thermal boundary conditions write θ = 0.5 on x 3 = 0, θ = -0.5 on

x 3 = 1, ∇θ•n = 0 on x 1 = 0, x 1 = 1, x 2 = 0, x 2 = 1.    (5)
Note that there is no radiative flux on the reflecting adiabatic sidewalls. The Boussinesq approximation is fully justified in the present study since, at the highest considered value of the Rayleigh number Ra = 10 9 and with the above numerical values, the relative temperature difference is ∆T /T 0 = 0.036. This would imply insignificant non-Boussinesq effects according to Ref. [START_REF] Sameen | Specific roles of fluid properties in non-Boussinesq thermal convection at the Rayleigh number of 2 × 10 8[END_REF] for instance.

The dimensionless radiative power is defined by

P rad (x) = L 2 λ∆T ν κ ν 4π I ν (x, Ω) dΩ -4πI • ν (T (x)) dν, (6)
where I ν (x, Ω) is the actual radiative intensity at wavenumber ν, position x, and direction Ω and I • ν (T (x)) is the Planck equilibrium intensity (blackbody intensity) at temperature T . The absorption coefficient κ ν is assumed to be uniform in accordance with the Boussinesq approximation. The high resolution spectrum of the absorption coefficient of the considered mixture is shown in Fig. 1. It contains thousands of spectral lines which makes expensive the computation of the integral over the wavenumbers in Eq. [START_REF] Leong | Experimental Nusselt numbers for a cubical-cavity benchmark problem in natural convection[END_REF]. Therefore, we use the Absorption Distribution Function (ADF) model [START_REF] Pierrot | A fictitious-gas-based absorption distribution function global model for radiative transfer in hot gases[END_REF] which consists in substituting the integration over the wavenumber with an integration over the values of the absorption coefficient, for which a coarse logarithmic discretisation is sufficient. In the present study, the values of the absorption coefficient of Fig. 1 have been logarithmically discretised in 16 classes and the accuracy of the model has been shown to be better than 1 % [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF]. Model parameters and computational details for the considered mixture are given in Refs. [START_REF] Soucasse | Effets des transferts radiatifs sur les écoulements de convection naturelle dans une cavité différentiellement chauffée en régimes transitionnel et faiblement turbulent[END_REF][START_REF] Soucasse | Transitional regimes of natural convection in a differentially heated cavity under the effects of wall and molecular gas radiation[END_REF].

Using the ADF model, the dimensionless radiative power writes

P rad (x) = L 2 λ∆T k κ k 4π I k (x, Ω) dΩ -4w k σT 4 (x) , (7)
where κ k and w k are respectively the absorption coefficient and the weight associated with the k th ADF class and σ is the Stefan-Boltzmann constant.

The radiative intensity I k (x, Ω) is obtained by solving the radiative transfer equation for a non scattering medium, for each class k

Ω • ∇I k (x, Ω) = κ k w k σT 4 (x) π -I k (x, Ω) . (8)
The associated boundary condition at wall positions x w for propagation directions Ω such that Ω•n > 0, n being the wall normal directed towards the interior of the domain, writes

I k (x w , Ω) = εw k σT 4 (x w ) π + 1 -ε π Ω •n<0 I k (x w , Ω )|Ω •n |dΩ . (9)

Numerical methods

Direct Numerical Simulations (DNS) have been performed, considering the air/H 2 O/CO 2 mixture as radiating (radiation case) or transparent (noradiation case, X H 2 O = 0, X CO 2 = 0 and thus P rad = 0) for 12 Rayleigh number values in the range 10 3 ≤ Ra ≤ 10 9 . The coupled algorithm for solving flow equations ( 2)-( 4) and radiative transfer equation ( 8) has been extensively used and described in previous works, for both differentially heated cavities [START_REF] Soucasse | Transitional regimes of natural convection in a differentially heated cavity under the effects of wall and molecular gas radiation[END_REF][START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF] and Rayleigh-Bénard cavities [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF][START_REF] Soucasse | Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different rayleigh numbers[END_REF]. We will briefly recall here the main features of this algorithm. Simulation parameters are given in Tab. 1.

Flow equations are solved using a collocation spectral method [START_REF] Xin | An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF][START_REF] Xin | 3D spectral parallel multi-domain computing for natural convection flows[END_REF] based on a Chebyshev polynomial expansion for the three directions of space. Time integration is performed through a second-order semi-implicit scheme. The pressure-flow coupling is ensured by a projection method that forces the velocity divergence free condition. Domain decomposition is applied along the vertical direction to perform parallel computations. The convection mesh is built from Chebyshev collocation points and is given in Tab. 1 for each Rayleigh number. We have checked that the number of points in the boundary layers is sufficient regarding the criterion proposed in Ref. [START_REF] Shishkina | Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution[END_REF]. Although this criterion was established for non-radiating fluids, calculation results show a posteriori that radiation does not affect so much the thickness of boundary layers. We therefore adopt the same meshing for flow filed calculations with and without radiation. The radiative transfer equation is solved using a ray-tracing algorithm [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF][START_REF] Soucasse | Transitional regimes of natural convection in a differentially heated cavity under the effects of wall and molecular gas radiation[END_REF], made parallel by distributing the rays among the different processors. The 4π steradian angular domain is uniformly discretized using between 900 rays to 3600 rays from volume cell centers, depending on the spatial mesh size. Up to Ra = 10 7 , the radiation mesh is built from the convection mesh, coarsened by a factor of two in each direction of space. For Ra ≥ 3 × 10 7 , the radiation mesh is built from the convection mesh, coarsened by a factor of four to eight in each direction of space and we use a subgrid model [START_REF] Soucasse | Subgrid-scale model for radiative transfer in turbulent participating media[END_REF][START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF] to account for the radiation of small spatial scales. This subgrid model has been validated in various configurations and its accuracy is about a few percent on radiative power and wall fluxes.

Finally, an explicit coupling is carried out between flow and radiation calculations and the radiation source term is updated every 10 convection time steps δt (every 5 convection time steps at Ra = 10 9 ). Indeed, the flow time step is imposed by numerical stability constraints and does not correspond to significant variations of the temperature field. Note that the computed radiative source term is based on instantaneous temperature fields and will be discussed itself in terms of mean and fluctuating fields (see Sec. 5). For unsteady solutions, time integration is carried out over a period ∆t once the asymptotic regime (statistically steady) is reached. It should be mentioned here that radiation calculations are much more computationally expensive than convection calculations (the CPU time is about 30 times greater in the radiation case).

Validation of the numerical tools

Validations of the pseudo-spectral solver were provided first in Ref. [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF], where steady solutions of the natural convection problem for differentially heated cavities at different Rayleigh numbers varying from Ra = 10 5 to 10 7 , without radiation, were calculated and successfully compared with the 3D benchmark solutions of Ref. [START_REF] Tric | A first incursion into the 3D structure of natural convection of air in a differentially heated cubic cavity, from accurate numerical solutions[END_REF]. For weakly turbulent natural convection without radiation, the predictions from the present solver were compared with the results provided by Refs. [START_REF] Trias | Direct numerical simulation of a two-and three-dimensional turbulent natural convection flows in a differentially heated cavity of aspect ratio 4[END_REF][START_REF] Trias | Direct numerical simulation of a differentially heated cavity of aspect ratio 4 with Rayleigh numbers up to 10 11 . Part I: Numerical methods and time-averaged flow[END_REF][START_REF] Trias | Direct numerical simulation of a differentially heated cavity of aspect ratio 4 with Rayleigh numbers up to 10 11 . Part II: Heat transfer and flow dynamics[END_REF] in the case of a parallelepipedic differentially heated cavity of aspect ratio L z /L x = 4, L y /L x = 1, filled with air (Pr = 0.71), at a Rayleigh number based on the height of the cavity Ra = 2×10 9 . Excellent agreement was found for both mean fields and second order statistics [START_REF] Soucasse | Transitional regimes of natural convection in a differentially heated cavity under the effects of wall and molecular gas radiation[END_REF]. Moreover, the grid convergence was discussed in Ref. [START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF] for differentially heated cubic cavities at Rayleigh numbers up to Ra= 3 × 10 9 .

For Rayleigh-Bénard convection in cubic cavities and without radiation, we compare in Tab. 2 the Nusselt number (dimensionless conductive flux av- The flow strengths are characterized by the two Reynolds numbers Re rms = u i u i v Ra/Pr and Re rms u 3 = u 3 u 3 v Ra/Pr, based either on the total velocity magnitude or the vertical velocity magnitude (• denotes the statistical average and • v the average over the entire volume of the cavity). The comparisons given in Tab. 3 show an excellent agreement with less than 1% differences in all cases.

As for radiative transfer computations, the ray-tracing algorithm was successfully compared to Monte Carlo results in the case of prescribed analytical temperature fields in a cubical enclosure [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF]. The accuracy of the refined angular discretization and of the ADF model were studied in Ref. [START_REF] Soucasse | Effets des transferts radiatifs sur les écoulements de convection naturelle dans une cavité différentiellement chauffée en régimes transitionnel et faiblement turbulent[END_REF] and Ref. [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF], respectively, and the radiation subgrid-scale model was validated, considering a snapshot of turbulent natural convection in a differentially heated cavity at Ra = 3 × 10 9 [START_REF] Soucasse | Subgrid-scale model for radiative transfer in turbulent participating media[END_REF].

Symmetries

Without radiative source term, Eqs. ( 2)-(4) satisfy four independent reflection symmetries S x 1 , S x 2 , S x 3 and S d with respect to the planes x 1 = 0.5, x 2 = 0.5, x 3 = 0.5 and x 1 = x 2 [START_REF] Puigjaner | Bifurcation analysis of steady Rayleigh-Bénard convection in a cubical cavity with conducting sidewalls[END_REF]. These symmetries act on the velocity and temperature fields as follows

S x 1 : (x 1 , x 2 , x 3 ) → (1 -x 1 , x 2 , x 3 ) (u 1 , u 2 , u 3 , θ) → (-u 1 , u 2 , u 3 , θ) (10)
S x 2 : (x 1 , x 2 , x 3 ) → (x 1 , 1 -x 2 , x 3 ) (u 1 , u 2 , u 3 , θ) → (u 1 , -u 2 , u 3 , θ) (11)
S x 3 : (x 1 , x 2 , x 3 ) → (x 1 , x 2 , 1 -x 3 ) (u 1 , u 2 , u 3 , θ) → (u 1 , u 2 , -u 3 , -θ) (12)
S d : (x 1 , x 2 , x 3 ) → (x 2 , x 1 , x 3 ) (u 1 , u 2 , u 3 , θ) → (u 2 , u 1 , u 3 , θ) (13)
These four elementary symmetries generate a symmetry group of sixteen elements. Radiation emission being proportional to T 4 , radiative transfer should break the S x 3 symmetry as the mean temperature gradient is directed along the x 3 axis. However, owing to the weak temperature gradients (∆T 10 K for the highest Rayleigh number), non-linear effects are negligible (namely 1 -

4T 3 0 ∆T T 4
hot -T 4 cold 3 × 10 -4 for the highest Rayleigh number) so that we can consider that the S x 3 symmetry is still satisfied in the radiation case. Symmetries ( 10)-( 13) act on the radiative power field in the same way as they act on the temperature field.

These many symmetries are responsible for multiple flow patterns in the studied Rayleigh range. In the no-radiation case, symmetry breaking solutions are observed from the onset of convection at a critical Rayleigh number Ra c 3400 [START_REF] Catton | The effect of insulating vertical walls on the onset of motion in a fluid heated from below[END_REF][START_REF] Puigjaner | Stability analysis of the flow in a cubical cavity heated from below[END_REF]: one observes a single roll flow pattern around either the x 1 axis or the x 2 axis (which breaks either the S x 2 or S x 1 symmetry, and also the S d symmetry). When the Rayleigh number increases up to the onset of unsteadiness, multiple stable or unstable flow patterns have been identified, as highlighted in Sec. 1. A detailed description of these flow patterns and associated bifurcation diagram, from the onset of convection and up to Ra = 1.5 × 10 5 , for Pr = 0.71 in a cubic cell, is provided in Ref. [START_REF] Puigjaner | Stability analysis of the flow in a cubical cavity heated from below[END_REF] for the no-radiation case. At higher Rayleigh numbers Ra ≥ 10 6 in unsteady regime, large-scale single roll flow patterns are observed, aligned along the diagonal planes x 1 = x 2 or x 1 = 1 -x 2 and breaking both S x 1 and S x 2 symmetries. However, all symmetries could be recovered for sufficiently long integration time, because of low-frequency reorientations of this large-scale roll in the horizontal plane, that are rotations of π/2 around the vertical axis [START_REF] Bai | Ability of a low-dimensional model to predict geometry-dependent dynamics of large-scale coherent structures in turbulence[END_REF][START_REF] Foroozani | Reorientations of the large-scale flow in turbulent convection in a cube[END_REF]. Reorientations of the large-scale roll in the horizontal plane have also been reported in the case of a radiating gas in the range 10 6 ≤ Ra ≤ 10 8 [START_REF] Soucasse | Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different rayleigh numbers[END_REF], radiation causing an increase of the reorientation frequency.

Given the variety of flow patterns that may be observed, it might seem difficult to compare the results of the radiation case and the no-radiation case. Therefore, we apply the symmetries to any flow variable φ after the simulation and we define the statistical average as follows

φ = 1 ∆t ∆t 1 16 16 n=1 S n (φ) dt, (14)
S n (φ) being the action of the n th symmetry (among 16 elements) of the symmetry group on the variable φ. The time average in Eq. ( 14) is only applied in the asymptotic regime of unsteady solutions.

Flow fields and heat transfer

In this section, we analyse the vertical distribution of key quantities averaged over horizontal planes. Statistical average (time average and symmetry average) as defined in Eq. ( 14) is applied to each solution.

Low Rayleigh numbers: onset of convection

Figure 2 shows the temperature, the kinetic energy and the radiative power (defined as the difference between the absorbed and emitted powers) obtained at the different Rayleigh numbers, in the radiation case and in the no-radiation case. At the lowest Rayleigh number (10 3 ), there is no kinetic energy for both the radiation and no-radiation case, which means that the convection has not started. The configuration is purely conductive in the noradiation case, which means that the temperature profile is linear. However, the presence of radiative energy transfer in the radiation case diminishes the temperature gradient in the center of the cavity: gas-gas radiative exchanges tend to homogenise the temperature field in the cavity. When the Rayleigh number increases (between Ra = 3 × 10 3 and Ra = 10 4 ), the convection settles in the no-radiation case, thus mixing the fluid and flattening the temperature profile. In the radiation case, though, the convection does not settle: the radiation has a delaying effect on the onset of the flow in the cavity. The gas is rather emitting near the hot wall and absorbing near the cold wall, and these radiative exchanges make the temperature field more stable.

But when the convection starts in the radiation case, between Ra = 10 4 and Ra = 3 × 10 4 , the kinetic energy in the cavity is higher than in the noradiation case. This is mainly due to an increase of potential energy in the radiation case. This will be discussed in details in Sec. 4. As the convection starts, the center of the cavity becomes more isothermal, but the temperature profile is steeper in the radiation case. Let us consider the lower half of the cavity: near the wall, the fluid is hotter than in the nearly isothermal center of the cavity. Therefore, the gas is emitting near the wall and then is absorbing up to the cavity center, creating a hump in the radiative power distribution, which is absent without convection. Because of the S x 3 symmetry, it is the exact opposite in the upper half of the cavity.

Higher Rayleigh numbers: unsteady convection

The onset of unsteadiness in the flow field occurs at different Rayleigh numbers, depending on the presence of radiation or not (between Ra = 10 5 and Ra = 3 × 10 5 for the no-radiation case, and between Ra = 3 × 10 5 and Ra = 10 6 for the radiation case). Radiation seems to delay the onset of unsteadiness, though a detailed investigation of the stability of all possible flow patterns would be required to confirm this point. As the Rayleigh number increases, the convection is more efficient at mixing the fluid. Therefore, the cavity is quasi-isothermal and thermal boundary layers form near the upper and lower walls. The kinetic energy is more and more homogeneous in the center, with peaks near the walls, where the temperature gradient is stronger. As before, the presence of gas radiation increases the kinetic energy, even if the relative effect of radiation diminishes as the Rayleigh number increases [START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF]. In the energy budget, the radiative power roughly scales as

1 √ Ra P rad = O 1 √ Ra κ P σT 3 0 L 2 λ ( 15
)
where κ P = κ ν I 0 ν (T 0 )dν × π/(σT 4 0 ) is the Planck mean absorption coefficient. Thus, the radiative power remains of the same order of magnitude but is penalized by the 1/ √ Ra factor, whereas the order of magnitude of the convective term u • ∇θ remains the same regardless of the Rayleigh number.

In order to better show the thermal boundary layer in the range 10 7 Ra 10 9 , Fig. 3 displays the temperature and radiative power distributions near the lower wall. Only the radiation case is shown as radiation has little effect on the boundary layer thickness for this Rayleigh number range. It can be seen that the thickness of the thermal boundary layer roughly scales as Ra -1/3 , in agreement with literature results in turbulent regime [START_REF] Malkus | The heat transport and spectrum of thermal turbulence[END_REF][START_REF] Niemela | Turbulent convection at very high Rayleigh numbers[END_REF]. The radiative power profile is similar to the one described at Ra = 10 5 : in the lower half of the cavity, the gas is rather absorbing except in a thin layer near the lower hot wall where the gas is emitting (the opposite happens nearby the upper wall). As the Rayleigh number increases, the convection mixing is stronger, the temperature gradient at the wall is sharper and the radiative power gradient at the wall increases. The thickness of the emitting layer and the position of the absorption peak also roughly scale as Ra -1/3 .

Heat fluxes

Figure 4 shows the three components of the total heat flux q tot q tot = -∂ θ s ∂z where q cond , q conv and q rad denote respectively the conductive, the convective and the radiative flux. • s denotes the spatial average over the horizontal plane and • denotes the statistical average as defined in Eq. ( 14). The total flux q tot is constant with height as the side walls are adiabatic. Values of wall fluxes and mid-height fluxes are reported in Tabs. 4 and 5 for the no-radiation case and the radiation case, respectively. In the no-radiation case, there is a balance between the conductive and convective fluxes. The radiative flux, reported in Tab. 4, represents the energy exchange between the isothermal black upper and lower walls and is constant with height and across the range of Rayleigh numbers. Figure 4 shows that, at high Rayleigh numbers (Ra 10 6 ), variations of the conductive and convective fluxes are restricted to the boundary layer regions. In the center, convection prevails and conduction vanishes.

q cond + √ Ra θu 3 s qconv + L λ∆T k 4π I k Ω • e x 3 dΩ s q rad , (16)
When gas radiation is taken into account, however, it modifies the balance between the fluxes. At low Rayleigh numbers, before the onset of convection in the cavity, the conduction-radiation coupling enables a greater conductive flux near the wall, as seen is Fig. 4. With the onset of convection, the presence of a radiating gas allows a greater convective flux in the center of the cavity: contrary to the no-radiation case, the variations of the convective wall mid-height Ra q cond q rad q cond q conv q rad 10 3

1.00 125.05 1.00 -125.05 3 × 10 3 1.00 125. 05 wall mid-height Ra q cond q rad q cond q conv q rad 10 3

2. [START_REF] Bejan | Heat Transfer Handbook[END_REF] flux are not restricted to the boundary layers due to a convection-radiation coupling in the quasi-isothermal part of the cavity. Though, the influence of radiative transfer decreases as the Rayleigh number increases, as seen in Eqs. ( 16): the convective flux scales as √ Ra, but not the radiative nor the conductive fluxes, whose relative effects become less important at higher Rayleigh numbers (q rad = O(κ P σT 3 0 L 2 /λ), the cavity size L being constant). Comparing Tabs. 4 and 5, it can also be noticed that the radiative flux at the wall decreases in the radiation case (screening effect due to gas radiation) and that gas radiation has little effect on the conductive flux at the wall for Ra 10 6 .

Energy budgets

We investigate in this section the effect of gas radiation on the balance of three quantities: the kinetic energy e k = 1 2 u i u i , the potential energy e p = -Prθ(x 3 -0.5) and the thermal energy e θ = 1 2 θ 2 , where • denotes the statistical average as defined in Eq. ( 14). For the latter quantity e θ , the word energy is used in a statistical sense and not in a physical sense and we use italics to emphasize it. The volume average of these three quantities at different Rayleigh numbers is shown in Fig. 5.

As stated in Sec. 3.1, gas radiation in the cavity delays the onset of convection: in the no-radiation case, the convection starts between Ra = 3 × 10 3 and Ra = 10 4 , whereas with gas radiation it starts between Ra = 10 4 and Ra = 3 × 10 4 . Kinetic energy increases with the Rayleigh number and reaches a peak at Ra = 3 × 10 5 . In the no-radiation case, the kinetic energy remains quite constant across the unsteady regime. The velocity was made dimensionless using the reference velocity u ref = a √ Ra L , which corresponds to a balance between the inertial and buoyancy forces. The quasi-plateau of the kinetic energy in the unsteady regime confirms that the reference velocity is well suited in this study. Gas radiation clearly increases the kinetic energy from Ra 10 5 : the kinetic energy peak is almost 50% higher in the radiation case than in the no-radiation case. Like the fluxes, the relative effect of radiation is less important as the Rayleigh number increases and the kinetic energy then decreases and becomes closer to the no-radiation case.

Potential energy and thermal energy are related to differences between the local temperature of the fluid and the mean temperature in the cavity. Both decrease with the Rayleigh number, given that the boundary layer zones (with the highest temperature gradients) narrow. As stated in Sec. the onset of convection, gas radiation decreases the temperature gradient in the core, and thus the potential energy and the thermal energy. However, after the onset, both are higher in the radiation case because of a higher temperature gradient in the core. The kinetic energy balance can be written in the form

∂e k ∂t + ∂u j e k ∂x j = - ∂u i p ∂x i + Prθu 3 + Pr √ Ra ∂ ∂x j u i ∂u i ∂x j - Pr √ Ra ∂u i ∂x j ∂u i ∂x j , (17)
After applying volume average and statistical average, this equation writes

Pr θu 3 v τu v = Pr √ Ra ∂u i ∂x j ∂u i ∂x j v εu v , (18)
where τ u v is the total production of kinetic energy (work of buoyant forces) and ε u v denotes the total mechanical dissipation. Radiation does not affect directly this balance. The potential energy balance writes

∂e p ∂t + ∂u j e p ∂x j = -Prθu 3 - Pr √ Ra ∂ ∂x i (x 3 -0.5) ∂θ ∂x i + Pr √ Ra ∂θ ∂x 3 - Pr √ Ra (x 3 -0.5)P rad . (19)
After applying volume average and statistical average, this equation writes

Pr √ Ra (Nu -1) ω θ v - Pr √ Ra (x 3 -0.5)P rad v ω rad v = Pr θu 3 v τu v , ( 20
)
where Nu is the Nusselt number (conductive flux at the wall). Terms ω θ v and ω rad v correspond to sources / sinks of potential energy associated with conduction and radiation transport of temperature variations. These sources / sinks are balanced with kinetic energy production ( τ u v , see Eq. ( 18)).

The three terms of the balance of Eq. ( 20) are displayed in Fig. 6. At the lowest Rayleigh numbers, τ u v = 0 because there is no fluid movement, and thus Nu = 1 in the no-radiation case, which is consistent with the profile of the conductive flux in Fig. 4. However, the conductive flux is sharper near the walls in the radiation case, even in the absence of convection (see 3.3). Thereby, the Nusselt number is greater than one and ω θ v is positive at low Rayleigh numbers. Interestingly the sign of the term ω rad v changes with the Rayleigh number: it is negative at low Rayleigh numbers and positive from Ra = 10 5 , when the kinetic energy is greater in the radiation case than in the no-radiation case. Therefore, the sign of this term determines whether gas radiation slows down convection or promotes convection. The sign change of ω rad v is correlated with the modification of the shape of the radiative power in Fig. 2. In the absence of convection, the lower half of the cavity (x 3 < 0.5) is exclusively emitting (P rad < 0) and the upper half (x 3 > 0.5) is exclusively absorbing (P rad > 0), which corresponds to a destruction of potential energy ( ω rad v < 0). With the onset of convection, the lower half becomes absorbing, except for a thin layer near the isothermal wall that is still emitting, and the opposite happens in the upper half: this corresponds to a production of potential energy ( ω rad v > 0).

The thermal energy balance writes

∂e θ ∂t + ∂u j e θ ∂x j = 1 √ Ra ∂ ∂x i θ ∂θ ∂x i - 1 √ Ra ∂θ ∂x i ∂θ ∂x i + 1 √ Ra θP rad . (21)
After applying volume average and statistical average, this equation writes Nu √ Ra

τ θ v = 1 √ Ra ∂θ ∂x i ∂θ ∂x i v ε θ v - 1 √ Ra θP rad v ε rad v . ( 22
)
The production of thermal energy ( τ θ v ) is dissipated through two mechanisms: conduction ( ε θ v ) and radiation ( ε rad v ). The three terms of the balance of Eq. ( 22) are displayed in Fig. 7.

At the lowest Rayleigh numbers, the production of thermal energy is much higher in the radiation case than in the no-radiation one, due to significant differences in conductive fluxes. The dissipation distributes quite evenly between conduction and radiation until the onset of convection, where the importance of radiation reduces compared to conduction. At higher Rayleigh numbers, the radiative dissipation is negligible.

Turbulence budgets

We analyse in this section turbulence budgets of velocity and temperature fluctuating fields in the range 10 7 Ra 10 9 . Mean and fluctuating fields are defined according to the statistical average • given by Eq. ( 14). Details on the computation of second order and third order statistics with this statistical mean are given in Appendix A.

The turbulent kinetic energy balance e k (u ) = 1 2 u i u i writes

∂e k (u ) ∂t + u j ∂e k (u ) ∂x j = - ∂ ∂x j u j p + 1 2 u i u i u j - Pr √ Ra ∂e k (u ) ∂x j δ u +Pru 3 θ τ u -u i u j ∂u i ∂x j ζ u - Pr √ Ra ∂u i ∂x j ∂u i ∂x j -ε u . (23)
In the above equation, τ u is the production term of turbulent kinetic energy by the work of buoyancy forces, ζ u , is an energy exchange term between mean and fluctuating flow, ε u is the molecular dissipation term (always positive), while δ u appears as a diffusion term -which vanishes when integrating over the whole spatial domain for non slip boundary conditions. Table 6 provides these volume integrated production ( τ u v ), transfer ( ζ u v ), and dissipation ( ε u v ) terms, in the radiation and no-radiation cases for the three Rayleigh numbers 10 7 , 10 (by a factor ∼ 10) contrary to what was observed in differentially heated configurations [START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF]. Radiation always increases this production term (and also the dissipation one necessarily), the magnitude of the increase being reduced when the Rayleigh number increases. This increase is consistent with the convection enhancement in the radiation case discussed in the previous sections.

Figure 8 shows the vertical distribution of each term in the right hand side of Eq. ( 23), averaged over the horizontal plane. Contrary to volume integrated terms, all terms can significantly contribute to the local balance of turbulent kinetic energy. In addition, all sources of mechanical fluctuations are significant in the core of the cavity which means the dynamics of turbulent fluctuations is not restricted to the boundary layers near the hot and cold walls. Note that the transfer term like the diffusion one have a varying sign through all the domain, which explains that they can contribute significantly to the local balance and in a much lesser extent to the global balance. Interestingly, the transfer term is negative in the boundary layer, while the mean flow shear at the walls is usually a source of turbulent fluctuations in most flows.

The balance of the variance of temperature fluctuations e θ (θ ) = writes

∂e θ (θ ) ∂t + u j ∂e θ (θ ) ∂x j = - ∂ ∂x j 1 2 u j θ θ - 1 √ Ra ∂e θ (θ ) ∂x j δ θ -u j θ ∂θ ∂x j ζ θ - 1 √ Ra ∂θ ∂x j ∂θ ∂x j -ε θ - 1 √ Ra θ P rad -ε rad (24)
Here again appear a diffusion term δ θ and a transfer term ζ θ between the mean temperature field and the fluctuating temperature field. ε θ corresponds to the conductive dissipation term. In the case of a radiating mixture, a supplementary radiative dissipation term ε rad appears. After integration over the whole spatial domain, this balance of temperature fluctuation variance reduces to an equilibrium between the transfer source term ζ θ v , and the conductive ε θ v and radiative (if relevant) ε rad v dissipation terms -the contribution of the diffusion term vanishes since temperature fluctuations at isothermal horizontal walls and temperature gradients at adiabatic perfectly reflecting walls vanish. Table 7 provides these volume integrated terms in the radiation and no-radiation cases for the three Rayleigh numbers 10 7 , 10 8 , and 10 9 . Here again, radiation is shown to increase the transfer term, but to a lesser extent compared to the increase observed for velocity fluctuations.

The radiative dissipation term contributes to 10% of the transfer term at Ra=10 7 but reduces to 5 and 1% at Ra =10 8 and 10 9 respectively. This behaviour was already observed in differentially heated configurations [START_REF] Soucasse | Natural convection in a differentially heated cubical cavity under the effects of wall and molecular gas radiation at Rayleigh numbers up to 3× 10 9[END_REF] and was explained by the ability of radiative transfer to dissipate large scale thermal structures while conduction becomes more efficient for smaller scale structures that occur at higher Rayleigh numbers [START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF].

Figure 9 shows the vertical distribution of each term in the right hand side of Eq. ( 24), averaged over the horizontal plane. There is a slight increase of each term when radiation is taken into account, as observed above on volume integrated terms: this is particularly noticeable at Ra = 10 7 . Let note again that the diffusion term, though vanishing in volume integrated balances, may have a significant local contribution in the redistribution of temperature fluctuations, in particular in the boundary layers near horizontal hot and cold walls. Contrary to mechanical fluctuations, the transfer term is positive in the boundary layers and the dynamics of fluctuations in the core Ra 10 is very weak.

Conclusion

Direct numerical simulations of Rayleigh-Bénard convection in a cubic cell for a radiating air-H 2 O-CO 2 mixture have been performed over a large range of Rayleigh numbers (Ra ∈ [10 3 -10 9 ]), from the onset of convection to the turbulent regime.

At low Rayleigh numbers, conduction-radiation energy exchanges make the temperature field more stable. This coupling delays the onset of convection. Once convection settles, convection-radiation exchanges in the cavity core have been highlighted. When gas radiation is taken into account energy transfer is no longer restricted to the boundary layers: convection is enhanced and so is the total kinetic energy of the flow. Gas radiation also increases the conductive flux and decreases the radiative flux at the two opaque isothermal walls, though differences with the no-radiation case are below 5 % from Ra 10 6 . Because the Rayleigh number is varied through the temperature difference ∆T , radiative transfer effects weaken when the Rayleigh number increases.

Specific contributions of radiative transfer to the potential energy balance and the thermal energy balance have been highlighted. It has been shown that gas radiation acts as a sink of potential energy at low Rayleigh numbers and as a source of potential energy from Ra 10 5 . In the thermal energy balance, radiation effects correspond to a dissipation of temperature differences. This radiative dissipation is comparable to the conductive dissipation up to Ra = 10 5 , but decreases and becomes negligible for higher Rayleigh numbers.

At high Rayleigh numbers, in the range 10 7 Ra 10 9 , the analysis of turbulence statistics has shown that all sources of turbulent fluctuations are stronger in the radiation case, especially velocity fluctuations. This turbulence intensification is associated with the convection enhancement observed when gas radiation is accounted for. Sources of temperature fluctuations are less affected by gas radiation: the convection enhancement is mitigated by radiative dissipation of temperature disturbances, although this term remains weak compared to other sources. Again, radiation effects on turbulence statistics weaken when the Rayleigh number increases.

To the best of our knowledge, the present coupled simulations are the first to address Rayleigh-Bénard convection of a radiating gas at such high Rayleigh numbers and such wide Rayleigh number range. Results will be useful for the development and validation of reduced order models or large-eddy simulation models, accounting for radiation effects. However, the simulations were restricted to a single gas mixture and a single cavity size. Sensitivity of radiation effects on these parameters will be the topic of future research.

Figure 1 :

 1 Figure 1: (a) Cubic Rayleigh-Bénard cavity filled with a radiating air/H 2 O/CO 2 mixture. Top and bottom walls are isothermal and black while side vertical walls are adiabatic and perfectly diffuse reflecting. (b) Absorption coefficient spectrum of the considered air/H 2 O/CO 2 mixture (T 0 = 300 K, X H2O = 0.02, X CO2 = 0.001, atmospheric pressure) and Planck function I • ν (T 0 ).

Figure 2 :

 2 Figure 2: Flow fields at different Rayleigh numbers. Temperature, kinetic energy and radiative power, averaged over the horizontal plane and statistically averaged. Radiation case (black lines) and no-radiation case (red lines).

Figure 3 :

 3 Figure 3: Temperature and radiative power close to the hot wall, averaged over the horizontal plane and statistically averaged. Rayleigh number range 10 7 Ra 10 9 . Radiation case. On the right plots, the abscissa is scaled by Ra 1/3 .

Figure 4 :

 4 Figure 4: Heat fluxes at different Rayleigh numbers. Conductive flux, convective flux and radiative flux, averaged over the horizontal plane and statistically averaged. Radiation case (black lines) and no-radiation case (red lines).

Figure 5 :

 5 Figure 5: Total kinetic energy e k v , potential energy e p v and thermal energy e θ v as a function of the Rayleigh number. Radiation case (black symbols) and no-radiation case (red symbols).

  radiation) ω θ (no-radiation) ω rad (radiation)

Figure 6 :

 6 Figure 6: Potential energy balance: τ u v , ω θ v and ω rad v a function of the Rayleigh number. Radiation case (black symbols) and no-radiation case (red symbols).

  radiation) ε θ (no-radiation) ε rad (radiation)

Figure 7 :

 7 Figure 7: Thermal energy balance: τ θ v , ε θ v and ε rad v as a function of the Rayleigh number. Radiation case (black symbols) and no-radiation case (red symbols).

Table 1 :

 1

	Ra	Convection mesh Radiation mesh	δt × 10 3	∆t
	10 3	41x41x(4x11)	20x20x20	10	N/A
	3 × 10 3	41x41x(4x11)	20x20x20	10	N/A
	10 4	41x41x(4x11)	20x20x20	10	N/A
	3 × 10 4	41x41x(4x11)	20x20x20	10	N/A
	10 5	65x65x(4x17)	32x32x32	10	N/A
	3 × 10 5	65x65x(4x17)	32x32x32	10	10,000 | N/A (a)
	10 6	81x81x(4x21)	40x40x40	5		10,000
	3 × 10 6	81x81x(4x21)	40x40x40	2.5	10,000
	10 7	81x81x(4x21)	40x40x40	2.5	10,000
	3 × 10 7 121x121x(6x21)	30x30x30 (b)	2 | 1.5 (a)	5,000
	10 8	161x161x(8x21)	40x40x40 (b)	1		5,000
	10 9	321x321x(16x41)	40x40x40 (b)	0.5 | 0.25 (a)	100
	Ra				10 6	10 7	10 8
	Kaczorowski et al. (2013) [53], DNS, Pr = 0.7 8.32 16.30 31.30
	Giannakis et al. (2018) [23], DNS, Pr = 0.7	-	16.57	-
	Foroozani et al. (2017) [15], LES, Pr = 0.7	8.10	-	31.60
	Xu et al. (2019) [54], LB, Pr = 0.7		8.33 16.22	-
	Our results, DNS, Pr = 0.707		8.30 16.24 31.06

Simulation parameters: Convection mesh, radiation mesh, convection time step δt and integration time interval ∆t for unsteady solutions. For the convection mesh, numbers in parenthesis indicate the number of spatial domains times the number of collocation points in the vertical in each domain.

(a)

No-radiation case | radiation case.

(b)

Radiation subgrid model is used.

Table 2 :

 2 Comparison of the averaged wall conductive Nusselt number with literature. eraged over the hot and cold walls) computed in the present study with the few values found in the literature. A good agreement is obtained with differences smaller than 1% with the DNS results of Kaczorowski et al. (2013) [53], 2% with the DNS results of Giannakis et al. (2018) [23], 2.5% with the Large-Eddy Simulation (LES) results of Foroozani et al. (2017) [15], and 0.4% with the Lattice-Boltzmann (LB) results of Xu et al. (2019)[START_REF] Xu | Lattice Boltzmann simulations of threedimensional thermal convective flows at high Rayleigh number[END_REF]. Note that the 1% difference for the values of the Prandtl number used by the different authors are not expected to affect the Nusselt numbers by more than 1% according to the classical correlations that may be found in the literature[START_REF] Bejan | Heat Transfer Handbook[END_REF].In addition to the Nusselt number, we compare in Tab. 3 the flow strengths obtained in the present study with the LB computations of[START_REF] Xu | Lattice Boltzmann simulations of threedimensional thermal convective flows at high Rayleigh number[END_REF] [START_REF] Xu | Lattice Boltzmann simulations of threedimensional thermal convective flows at high Rayleigh number[END_REF].

	Re rms	Re rms u 3

Table 3 :

 3 Comparison of the flow strength with results of Xu et al. (2019) [54].

Table 4 :

 4 Wall heat fluxes and mid-height heat fluxes at different Rayleigh numbers. Conductive flux, convective flux and radiative flux, averaged over the horizontal plane and statistically averaged. No-radiation case.

		1.00	-	125.05
	10 4	1.94 125.05 0.12 1.82 125.05
	3 × 10 4 2.67 125.05 -0.13 2.80 125.05
	10 5	4.24 125.05 -0.11 4.35 125.05
	3 × 10 5 5.92 125.05 -0.06 5.99 125.05
	10 6	8.30 125.05 0.18 8.12 125.05
	3 × 10 6 11.48 125.05 0.15 11.34 125.05
	10 7	16.24 125.05 0.08 16.14 125.05
	3 × 10 7 22.06 125.05 0.04 22.02 125.05
	10 8	31.06 125.05 0.02 31.04 125.05
	10 9	61.38 125.05 0.002 61.87 125.05

Table 5 :

 5 Wall heat fluxes and mid-height heat fluxes at different Rayleigh numbers. Conductive flux, convective flux and radiative flux, averaged over the horizontal plane and statistically averaged. Radiation case.

Table 6 :

 6 Production τ u v , transfer ζ u v , and dissipation ε u v terms of the turbulent kinetic energy balance, averaged on the whole domain, in the radiation (rad) and noradiation (no-rad) cases for the three Rayleigh numbers Ra= 10 7 , 10 8 , and 10 9 . Term variations due to radiation are provided in parentheses integrated production term τ u v highly dominates the transfer term ζ u v

  Turbulent kinetic energy budget in the range 10 7 ≤ Ra ≤ 10 9 . Production τ u (green lines), transfer ζ u (blue lines), dissipation u (red lines) and diffusion δ u (black lines) of turbulent kinetic energy for radiation and no-radiation cases.

				no-radiation					radiation
		0.01						0.01
	7
	Ra=10	0						0
		-0.01						-0.01
	Ra=10 8	-0.01 0 0.01						-0.01 0 0.01							ζ u' (transfer) δ u' (diffusion) -ε u' (dissipation) τ u' (production)
	0.005						0.005
	9
	Ra=10	0						0
	-0.005						-0.005
		0	0.2	0.4	x 3	0.6	0.8	1	0	0.2	0.4	x 3	0.6	0.8	1
	Figure 8:
															1 2 θ θ

Table 7 :

 7 Transfer ζ θ v , conductive ε θ v and radiative ε rad v dissipation terms of the balance equation for the variance of temperature fluctuations, averaged on the whole domain, in the radiative (rad) and non-radiative (no-rad) cases for the three Rayleigh numbers Ra= 10 7 , 10 8 , and 10 9 .

			7		10 8		10 9
	terms ×10 3 no-rad	rad	no-rad	rad	no-rad	rad
	ζ θ v	1.73	1.99 (15%)	1.16	1.23 (6%)	0.774 0.780 (1%)
	ε θ v	1.73	1.81 (5%)	1.16	1.23 (6%)	0.797 0.797 (--)
	ε rad v	-	0.20	-	0.06	-	0.008

Acknowledgements

This work benefited from the financial support of the "Agence de l'environnement et de la maîtrise de l'énergie" (ADEME, France). This work was granted access to the HPC resources of IDRIS under the allocation 2020-A0062B00209 attributed by GENCI (Grand Equipement National de Calcul Intensif). This work was also performed using HPC resources from the Mésocentre computing center of CentraleSupélec and École Normale Supérieure Paris-Saclay supported by CNRS and Région Île-de-France (http://mesocentre.centralesupelec.fr/).

Appendix A. Turbulence statistics using symmetry average

This section details the computation of second order and third order turbulence statistics in turbulence budgets of Sec. 5 with the statistical average defined in Eq. ( 14), which implies an ensemble average of symmetrised time evolutions of each flow variable of interest.

The first step is to compute the symmetry average (ensemble average over symmetrised states) of physical quantities of the form φ, φψ or φφψ.

To do this, we need to express the action of each elementary symmetry S x 1 , S x 2 , S x 3 and S d (Eqs. ( 10)-( 13)) on each flow variable and each product of flow variables of interest. This is given in Tab. A.8. Note that for third

Table A.8: Action of elementary symmetries defined in Eqs. ( 10)-( 13)) on each flow variable or product of flow variable of interest.

order statistics, only products of the type φφu 3 are required to compute the diffusion terms δ u and δ θ in Eqs. ( 23)-( 24), averaged over horizontal planes.

The symmetry average can be applied before the time average or after the time average, indifferently. Once averaged quantities φ, φψ or φφψ are known, the second step consists in computing second order statistics φ ψ or third order statistics φ φ ψ with the following relationships