Laurent Soucasse

Bérengère Podvin

Philippe Rivière

Anouar Soufiani

Low-order models for predicting radiative transfer effects on Rayleigh-Bénard convection in a cubic cell at different Rayleigh numbers

Keywords:

This paper presents low-order models of Rayleigh-Bénard convection of a radiating gas in a cubic cell, in the Rayleigh number range Ra ∈ [10 6 -10 8 ]. Numerical simulations are carried out for an air/H 2 O/CO 2 mixture assumed to be radiating (coupled case) or transparent (uncoupled case). When coupling with radiation, it is shown that the kinetic energy of the flow and the thermal energy increase. At Ra = 10 6 , planar flow states are observed when radiation is taken into account, while diagonal flow states prevail in the uncoupled case. From Ra 3 × 10 7 , quasi-stable diagonal flows are observed in both coupled and uncoupled simulations, with occasional brief reorientations. The reorientation frequency seems to decrease with the Rayleigh number and seems to increase with radiation. A Proper Orthogonal Decomposition (POD) analysis reveals that eleven of the first twelve POD eigenfunctions are preserved over the Rayleigh number range, whatever the coupling conditions. However, POD eigenvalues are higher with radiation. POD-based low-order models are derived at different Rayleigh numbers, for both coupled and uncoupled cases. Radiative transfer effects are added in the model in an a priori fashion, from uncoupled simulation data. Coupled POD models predict the energy increase with radiation and the loss of stability of the diagonal rolls at Ra = 10 6 . Uncoupled and coupled models correctly reproduce reorientation frequencies over the Rayleigh number range. Finally, a generalized model is derived, solely based on uncoupled simulation data at Ra = 10 7 and energy scaling laws. This generalized model captures the change in dynamics associated with radiation effects and variations in Rayleigh number, except at Ra = 10 6 .

Introduction

The problem of Rayleigh-Bénard convection (RBC) consists of a fluid layer confined between two horizontal plates, heated from below and cooled from above. This problem has been widely studied to model buoyancy-driven flows encountered in geophysics or engineering, and to understand the fundamental properties of turbulence. In turbulent regime, the complexity comes from the chaotic nature of the flow and the wide range of time and length scales on which velocity and temperature fields vary.

The dynamics of turbulent RBC becomes more complicated when the working fluid is a radiating gas, which emits and absorbs thermal radiation. It is the case, for instance, in the atmosphere or in confined environments such as buildings, where the fluid (the air) contains radiating molecules in the infrared such as water vapour or carbon dioxide. Radiative transfer effects on turbulent convection have been first investigated by [START_REF] Spiegel | The smoothing of temperature fluctuations by radiative transfer[END_REF], who showed that radiation acts as a dissipation mechanism of temperature fluctuations. This damping effect prevails over conductive dissipation for large and intermediate turbulent scales and vanishes for small turbulent scales from a critical length scale depending on the radiative properties of the medium [START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF]. However, in the case of thermally-driven flow, radiation also affects the mean temperature gradient and buoyant motion. Several researchers have reported an increase of the mean kinetic energy of the flow due to an increase of the mean potential energy in various configurations involving a radiating gas: in a differentially heated cavity [START_REF] Kogawa | Influence of radiation effect on turbulent natural convection in cubic cavity at normal temperature atmospheric gas[END_REF], in a Rayleigh-Bénard cavity (Soucasse et al. 2014a) or for a confined plume generated by a linear heat source [START_REF] Wang | Numerical study of an unsteady confined thermal plume under the influence of gas radiation[END_REF]. The cost of radiative transfer computations in a turbulent medium restricts the numerical investigation of radiation effects due to the angular, spectral, and spatial dependence of the radiation field. Experimental works are also challenging to undertake, as non-intrusive techniques that do not perturb the radiation field are required. Therefore, the question arises whether a low-order model could capture radiation effects on RBC, given that radiative transfer rather affects the large and intermediate scales of the flow.

In the case of a non-radiating gas, the large-scale motion in RBC, also referred to as Large Scale Circulation (LSC), has aroused a great interest as it intermittently changes orientation. In cylindrical cells with aspect ratio of unity (diameter equal to height), azimuthal rotation of the LSC and reversal of the LSC (sudden change of direction in the same circulation plane) have been reported by 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 numerical studies [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]. In cubic cells, low-frequency reorientations of the LSC have been observed between four quasi-stable states, corresponding to the LSC lying in the two diagonal vertical planes of the cube associated with clockwise and anticlockwise motion [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]. For both cylindrical and rectangular cells, it has been shown that the LSC dynamics is very sensitive to the aspect ratio [START_REF] Xi | Azimuthal motion, reorientation, cessation, and reversal of the large-scale circulation in turbulent thermal convection: A comparative study in aspect ratio one and one-half geometries[END_REF][START_REF] Vasiliev | Reversals of large-scale circulation in turbulent convection in rectangular cavities[END_REF][START_REF] Ni | Reversals of the large-scale circulation in quasi-2D RayleighBénard convection[END_REF].

Several models have been proposed to explain the LSC dynamics in RBC for nonradiating gases. [START_REF] Brown | Large-scale circulation model for turbulent Rayleigh-Bénard convection[END_REF], 2008) derived a stochastic two-equation model to predict the time evolution of the strength and the azimuthal orientation of the LSC in cylindrical geometries. This model has been extended for cubical geometries, by the addition of a potential term which drives the azimutal angle towards the vertical edges of the cube and the LSC in the diagonal planes [START_REF] Bai | Ability of a low-dimensional model to predict geometrydependent 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]. Another phenomenological model for reversals in a square cell has been proposed by [START_REF] Araujo | Wind reversals in turbulent Rayleigh-Bénard convection[END_REF]. From this model, the authors have established a range of Prandtl and Rayleigh numbers where reversals occur and a scaling law for the reversal frequency. Other modelling approaches rely on modal decomposition of temperature and velocity fields based on Fourier modes [START_REF] Chandra | Dynamics and symmetries of flow reversals in turbulent convection[END_REF], Koopmann modes [START_REF] Giannakis | Koopman analysis of the long-term evolution in a turbulent convection cell[END_REF] or Proper Orthogonal Decomposition (POD) modes [START_REF] Bailon-Cuba | Aspect ratio dependece of heat transfer and large-scale flow in turbulent convection[END_REF][START_REF] Podvin | Proper orthogonal decomposition investigation of turbulent Rayleigh-Bénard convection in a rectangular cavity[END_REF].

POD is often used in fluid mechanics to extract large-scale coherent structures from numerical or experimental data. The POD modes form an orthogonal basis, which is optimal to represent the flow regarding its energy content. A low-order model for the flow can be easily derived from Galerkin projection of Navier-Stokes equation onto the orthogonal POD basis. In a previous study we have developed a POD model for capturing reorientations of the LSC in a cubic Rayleigh-Bénard cell at a Rayleigh number Ra = 10 7 [START_REF] Soucasse | Proper orthogonal decomposition analysis and modelling of large-scale flow reorientations in a cubic Rayleigh-Bénard cell[END_REF]). In the presence of radiating species, we have shown that radiation effects can be taken into account by a POD model, from a rigorous projection of the radiative source term of the energy balance onto the POD basis [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF]. However, a frequent pitfall of POD models is that their prediction capabilities are restricted in the neighbourhood of the flow parameters associated with the data. This paper aims to examine at which extent radiative transfer effects can be predicted from uncoupled simulation data across a wide range of Rayleigh numbers in a cubic cell. Direct numerical simulations (DNS), coupled or uncoupled with radiation, have been performed in the range Ra ∈ [10 6 -10 8 ], where reorientations are more likely to be observed [START_REF] Sugiyama | Flow reversals in thermally driven turbulence[END_REF]. A radiating air/H 2 O/CO 2 mixture at room temperature has been considered to make our study relevant for building applications. Numerical simulations are presented in section 2 and associated POD analysis is performed in section 3. Uncoupled POD models (non-radiating case) and predicted coupled POD models (radiating case) are derived in sections 4 and 5 at different Rayleigh numbers.

Finally, a general model across the Rayleigh range, solely based on uncoupled simulation data at Ra = 10 7 and energy scaling laws, is developed in section 6.

Direct numerical simulations

Problem setup

We consider the natural convection flow of an air/H 2 O/CO 2 mixture confined in a cubic Rayleigh-Bénard cell, with top and bottom isothermal walls at T cold and T hot and adiabatic vertical walls. The six walls are characterised by uniform grey emissivities ε, the horizontal isothermal walls being black (ε = 1), and the vertical adiabatic walls being perfectly diffuse reflecting (ε = 0). The parameter controlling the flow regime is the Rayleigh number defined by

Ra = gβ∆T L 3 ν f a , (2.1)
where g is the gravitational acceleration, β = 1/T 0 is the thermal expansion coefficient (T 0 being the mean temperature), L is the size of the cavity, ν f is the kinematic viscosity, a is the thermal diffusivity and ∆T = T hot -T cold is the temperature difference between hot and cold walls. We investigate the Rayleigh range Ra ∈ [10 6 -10 8 ] in which reorientations are likely to be observed [START_REF] Sugiyama | Flow reversals in thermally driven turbulence[END_REF]. We vary the Rayleigh number by changing the temperature difference ∆T and other parameters are fixed. 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, of molar composition X H2O = 0.02 and X CO2 = 0.001, and the cavity size is set to L = 1 m, so that the temperature difference vary in the range ∆T ∈ [1.1 × 10 -2 -1.1] K. 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 = ν f /a = 0.707).

Mass, momentum and energy balance are made dimensionless using the cavity size L, the reference time L 2 /(a √ Ra) and the reduced temperature θ = (T -T 0 )/∆T and write under Boussinesq approximation (2.4) where u = (u, v, w) is the dimensionless velocity vector and p the dimensionless motion pressure. The velocity is zero on all walls. The temperature is set to 0.5 and -0.5 respectively on the bottom and top walls. The conductive flux -∇θ • n is zero on the four lateral walls, as these walls have a zero radiative emissivity.

∇ • u = 0, (2.2) ∂u ∂t + u • ∇u = -∇p + Pr θ e z + Pr √ Ra ∇ 2 u, (2.3) ∂θ ∂t + u • ∇θ = 1 √ Ra ∇ 2 θ + P rad ,
The dimensionless radiative power P rad in equation (2.4) accounts for absorption and emission of thermal radiation by the medium and is defined by

P rad (r) = L 2 λ∆T ν κ ν 4π I ν (Ω, r)dΩ -4πI • ν (T (r)) dν, (2.5)
where I ν (Ω, r) is the actual radiative intensity at wavenumber ν, direction Ω and position r. I • ν (T (r)) is the Planck equilibrium intensity at temperature T and κ ν is the absorption coefficient of the medium. In accordance with the Boussinesq approximation, the absorption coefficient is assumed to be spatially uniform. The Planck equilibrium intensity is given by

I • ν (T (r)) = 2hν 3 c 2 0 1 exp hν k B T (r) -1 , (2.6)
where h is the Planck constant, k B is the Boltzmann constant and c 0 the speed of light in vacuum. The radiative intensity field is obtained by solving the radiative transfer equation

Ω •∇ I ν (Ω, r) = κ ν (I • ν (T (r)) -I ν (Ω, r)) .
(2.7)

The associated boundary condition at boundary points r w for gray diffuse reflecting walls writes

I ν (Ω, r w ) = εI • ν (T (r w )) + 1 -ε π Ω •n<0 I ν (Ω , r w )|Ω •n |dΩ , (2.8)
for directions Ω such that Ω•n > 0, n being the wall normal directed towards the inside of the domain.

It is worth noting that the flow equations are written and solved in dimensionless form, while the radiative transfer equations are treated in dimensional form since we consider an actual molecular radiating gas. A key parameter for radiation is the optical thickness τ ν = κ ν L that varies over several orders of magnitude in our model. Considering a grey fluid (wavelength-independent absorption) would facilitate a parametric study of radiation effects but would fail to represent the behaviour of actual radiating gases.

Numerical methods

The numerical methods used for solving the coupled system of equations (2.2)-(2.4) and equation (2.7) are presented and validated in details by [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] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF]. We briefly mention here the main features of the coupled algorithm.

Navier-Stokes equations are solved using a Chebyshev collocation method [START_REF] Xin | An extended Chebyshev pseudo-spectral benchmark for the 8:1 differentially heated cavity[END_REF]. Domain decomposition along the vertical direction is carried out by

Ra

Convection mesh Radiation mesh δt × 10 3 ∆t 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 (a) 2/1.5 (b) 5,000 10 8 161x161x(8x21) 40x40x40 (a) 1 5,000

Table 1. Simulation parameters: convection mesh, radiation mesh, convection time step δt and integration time interval ∆t in the asymptotic regime. 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) Radiation subgrid model is used. (b) δt is 2×10 -3 in the uncoupled simulation and 1.5 × 10 -3 in the coupled simulation.

the Schur complement method to make the computations parallel [START_REF] Xin | 3D spectral parallel multi-domain computing for natural convection flows[END_REF]. Time integration is performed through a second-order semi-implicit scheme. The velocity divergence free condition is enforced using a projection method. The radiative transfer equation is solved using a ray-tracing algorithm, made parallel by distributing the rays among the different processors. The 4π angular domain is uniformly discretized using 900 rays from volume cell centers and 450 rays from boundary cell centers. 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]) is used to take into account the spectral variations of the absorption coefficient of the air/H 2 O/CO 2 mixture: it 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 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 by [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 [START_REF] Soucasse | Numerical study of coupled molecular gas radiation and natural convection in a differentially heated cubical cavity[END_REF]. DNS have been performed, considering the air/H 2 O/CO 2 mixture as radiating (coupled case) or transparent (uncoupled case, X H2O = 0, X CO2 = 0 and thus P rad = 0). Simulation parameters are given in Table 1. Five different Rayleigh numbers have been considered: Ra = {10 6 ; 3 × 10 6 ; 10 7 ; 3 × 10 7 ; 10 8 }. The convection mesh is built from Chebyshev collocation points. We have checked that the number of points in the boundary layers is sufficient regarding the criterion proposed by [START_REF] Shishkina | Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution[END_REF]. 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 and Ra = 10 8 , the radiation mesh is coarsened by a factor of four in each direction of space compared to the convection mesh and we use a radiation subgrid model (Soucasse et al. 2014b) 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. It has been used for the simulation of coupled natural convection and radiation in a differentially heated cavity at Ra = 3 × 10 9 [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]. 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 (the flow time step is imposed by numerical stability constraints and does not correspond to significant variations of the temperature field). 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 coupled case).

Radiative transfer effects

Radiative transfer effects on RBC at Ra = 10 7 have been discussed in a previous work [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF]. When the fluid emits and absorbs radiation, heat transfer is no longer restricted to the boundary layer region. An energy exchange between convection and radiation in the core of the cavity leads to a significant increase of the convective flux compared to the uncoupled case. The LSC is strengthened and both mean and turbulent kinetic energies increase. Temperature fluctuations also increase but to a lesser extent because of radiative damping.

The same effects are observed at other Rayleigh numbers. The total kinetic energy e k = 0.5 u • u dr, the total thermal energy e θ = 0.5 θθ dr, the conductive flux at the walls q cond = ∂θ/∂z dxdy and the convective flux in the core (z = 0.5) q conv = √ Ra wθ dxdy are displayed in figure 1 as a function of the Rayleigh number ( • denotes the time average). It can be noticed that the kinetic energy and the convective flux significantly increase in the presence of radiation while the thermal energy and the conductive flux are not much affected. However, radiation effects diminish with the Rayleigh number. This can be explained by the following scaling analysis. If the temperature dependance of the Planck equilibrium intensity is linearised around the mean temperature T 0 , an order of magnitude for the dimensional radiative power is 16κ P σT 3 0 ∆T , where κ P = ν κ ν I • ν (T 0 )dν × (π/σT 4 0 ) is the Planck mean absorption coefficient of the mixture, σ is the Stefan-Boltzmann constant. Therefore, the radiative source term in the energy balance (2.4) roughly scales as

P rad √ Ra ∼ O 1 √ Ra 16κ P σT 3 0 L 2 λ , (2.9)
and decreases in Ra -1/2 while the order of magnitude of the convective term u • ∇θ remains constant.

The unsteady dynamics of Rayleigh-Bénard convection in a cubic container is characterised by low-frequency reorientations of the Large Scale Circulation (LSC) (Bai et al. 2016; [START_REF] Foroozani | Reorientations of the large-scale flow in turbulent convection in a cube[END_REF], that can be monitored using the time evolution of the x and y components of the angular momentum with respect to the cavity centre r 0 L = (r -r 0 ) × u dr.

(2.10)

Figure 2 shows the time evolution of components L x and L y at the different Rayleigh numbers for the coupled and uncoupled cases. Note that time integration is shorter for the two highest Rayleigh numbers. The uncoupled case is characterised at each Rayleigh number by quasi-stable diagonal states, with abrupt reorientations between these states. A diagonal state means that the LSC lies in one of the two diagonal planes x = y or x = 1 -y, with clockwise or anticlockwise motion (four diagonal states are available), and is characterized by a non-zero equilibrium value for both L x and L y components. Sudden reorientations between two diagonal states occur when either L x or L y changes sign, which corresponds to a rotation of π/2 of the LSC around the vertical axis. A detailed description of the reorientation process in the cubic cell and associated flow patterns is provided by [START_REF] Foroozani | Reorientations of the large-scale flow in turbulent convection in a cube[END_REF] and [START_REF] Vasiliev | Transient flows and reorientations of large-scale convection in a cubic cell[END_REF]. Although the overall dynamics is similar over the Rayleigh number range, it can be noticed that the stability of the diagonal states increases with the Rayleigh number: the flow spends less time around zero angular momentum and reorientations are less frequent when the Rayleigh number increases. At Ra = 10 6 , the dynamics is more chaotic with quick passing around zero of either L x or L y components. At Ra = 10 8 , the oscillation amplitude of L x and L y around the equilibrium value is weaker and only one reorientation event is observed during the sequence (only two of the four diagonal states are visited in this case).

In the coupled case, a significant change in dynamics compared to the uncoupled case is noticeable at Ra = 10 6 , where quasi-stable planar sates are observed. A planar state means that the LSC lies either in x planes or y planes, with clockwise or anticlockwise motion (four planar states are available), and is characterized by a zero equilibrium value for one of the two components L x or L y . Figure 3 shows a snapshot of temperature and velocity fields at Ra = 10 6 . In the uncoupled case (left part of figure 3) a diagonal state is observed as the LSC lies in the diagonal plane x = 1-y with L x > 0, L y > 0. The fluid flows up along the right vertical edge (x; y) = (0; 1) and flows down along the left vertical edge (x; y) = (1; 0). This main diagonal roll is slightly tilted and small counter-rotating rolls are noticeable in the top right corner (x; y; z) = (0; 1; 1) and in the bottom left corner ((x; y; z) = (1; 0; 0). In the coupled case (right part of figure 3) a planar state is observed as the LSC lies in y planes with L x 0, L y < 0. The fluid flows up along the front plane x = 1 and flows down along the rear plane x = 0. Counter-rotating structures are noticeable near the top horizontal edge (x; z) = (1; 1) and the bottom horizontal edge (x; z) = (0; 0). For Ra 3×10 7 , quasi-stable diagonal states are observed with radiation. However, for a given Rayleigh number, the dynamics is more chaotic and reorientations seem to be more frequent when the flow is coupled with radiation. It is worth noticing here that, although stable planar states are not observed in the uncoupled case, they have been reported at low Rayleigh numbers [START_REF] Puigjaner | Stability analysis of the flow in a cubical cavity heated from below[END_REF]).

In order to quantify radiative transfer effects on the temporal dynamics we have computed two characteristic frequency scales: the circulation frequency f c (or circulation time τ c = 1/f c ) and the reorientation frequency f r (or reorientation time τ r = 1/f r ). The circulation frequency is a high frequency associated with the rotation frequency of the LSC roll. Frequencies are reported in table 2. The circulation frequency is estimated using a reference ellipsoid path length in the diagonal plane and a reference velocity. In the uncoupled case, the circulation frequency is nearly constant with the Rayleigh number as convective time units are used. In the coupled case, the increase of the kinetic energy leads to an increase of the circulation frequency. This increase compared to the uncoupled case diminishes with the Rayleigh number. The reorientation frequency is estimated by tracking the zeros of the filtered time evolution of L x and L y to avoid small-scale noise. Data in table 2 confirm the observations of figure 2: the reorientation frequency seems to decrease with the Rayleigh number and, at a given Rayleigh number, seems to be higher when radiation is taken into account.

It should be noted here that the uncertainty associated with reorientation frequencies is significant, owing to the few reorientations. If we assume a Poisson distribution of reorientation events (as it has been observed for reversals in cylindrical cells [START_REF] Sreenivasan | Mean wind and its reversal in thermal convection[END_REF])), the uncertainty on the number of reorientations N r would be ± √ N r , and the uncertainty on the reorientation frequency f r would be ± (f r /∆t), ∆t being the integration time. Uncertainties on reorientation frequencies are given in table 2. Relative uncertainties range from 20 % (Ra = 3 × 10 6 , uncoupled case) to 70 % (Ra = 3 × 10 7 , uncoupled case). These uncertainties moderate the conclusions on radiative transfer effects on the reorientation frequency, especially at Ra = 10 7 where the increase of the reorientation frequency with radiation is not statistically significant. However, the decrease of the reorientation frequency with the Rayleigh number in the uncoupled case seems to be statistically significant in the range 3 × 10 6 Ra 10 8 and has been reported in other works [START_REF] Sugiyama | Flow reversals in thermally driven turbulence[END_REF].

Proper Orthogonal Decomposition

Methodology

The POD in fluid mechanics aims at finding an optimal basis of spatial eigenfunctions φ n (r) to represent an unsteady flow variable vector U (r, t) of size M on a lowdimensional subspace. The POD eigenfunctions or POD modes φ n (r) are solution of the following eigenvalue problem [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF])

M k=1 U m (r, t)U k (r , t) φ k n (r )dr = λ n φ m n (r), (3.1)
where • denotes the time average. They form an orthonormal basis allowing the decomposition

U (r, t) = ∞ n=1 a n (t)φ n (r), (3.2)
where the projection coefficient a n (t) are statistically uncorrelated and their energy is equal to the eigenvalue such that a n (t)a m (t) = δ nm λ n , where δ nm is the Kronecker symbol. The eigenvalue associated with a POD mode is thus a measure of its energy content and the objective is to restrict the decomposition (3.2) to a few modes with the largest eigenvalues so that the associated low-order subspace captures most of the energy of the field U (r, t). To take into account the coupling between velocity and temperature fields in thermal convection, we define the flow variable vector as U = {u, γθ} (M = 4), where γ is a rescaling factor

γ = u(r, t) • u(r, t) dr θ 2 (r, t) dr , (3.3)
such that velocity and temperature fields have the same energy [START_REF] Podvin | Low-order models for the flow in a differentially heated cavity[END_REF]. The POD modes thus combine velocity and temperature: φ n = {φ u n , γφ θ n }. Equation (3.1) is solved in practice using the method of snapshots [START_REF] Sirovich | Turbulence and the dynamic of coherent structures. Part I: Coherent structures[END_REF]. A snapshot set of 1,000 samples U (r, t i ) is extracted from each simulation at discrete times t i with a fixed sampling period of 10 dimensionless time units for 10 6 Ra 10 7 and a fixed sampling period of 5 dimensionless time units for 3 × 10 7 Ra 10 8 . This sampling period might seem large compared to the transition time between two states and the few reorientations, but according to [START_REF] Podvin | Precursor for wind reversal in a square Rayleigh-Bénard cell[END_REF], precursor events for reorientations are associated with large-scale interactions during relatively large time scales (of the order of the reorientation period). However, we have seen in the time evolution of the angular momentum components in figure 2 that the four possible quasi-stable flow states (planar or diagonal) were not necessarily visited or not equally represented during the time sequence, although there is no physical evidence suggesting that these states are not equiprobable. This is an artefact owing to the relatively short simulation time compared with the time scale separating two reorientations, especially at high Rayleigh number. In order to enforce an equal statistical weight for each flow state and to improve the convergence of the POD method, we have built enlarged snapshot sets, obtained by the action of the symmetry group of the problem on the original snapshot sets. The use of statistical symmetries of the flow in POD is discussed at length by [START_REF] Moin | Characteristic-eddy decomposition of turbulence in a channel[END_REF] and [START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF].

In the uncoupled case, the problem satisfies four independent reflection symmetries S x , S y , S z and S d with respect to the planes x = 0.5, y = 0.5, z = 0.5 and x = y [START_REF] Puigjaner | Bifurcation analysis of steady Rayleigh-Bénard convection in a cubical cavity with conducting sidewalls[END_REF]). In the coupled case, radiative transfer should break the S z symmetry (radiative emission being proportional to T 4 ), but owing to the weak temperature gradients, nonlinear effects are negligible so that we can consider that the S z symmetry is still satisfied. These four elementary symmetries generate a symmetry group of sixteen elements. This allows us to multiply the number of snapshot by a factor of sixteen. In conclusion, for each simulation (coupled, uncoupled, 10 6 Ra 10 8 ), we work on an enlarged snapshot set made of 16,000 samples (1,000 original snapshots multiply by sixteen after applying all the symmetries).

Energy spectra

The POD spectra obtained in the uncoupled case at different Rayleigh numbers are shown in the left part of figure 4. Each spectrum is normalized such that n λ n = 1 but since the total POD energy is nearly constant with the Rayleigh number in the uncoupled case, it is relevant to compare the spectra between them. The mode ordering roughly corresponds to a ranking of the eigenfunctions in terms of a characteristic spatial scale, and we can therefore associate the low-order modes with the largest spatial scales and the high-order modes with the smallest spatial scales. Examination of figure 4 suggests that the POD spectra can be split into three parts: (i) the large scales for n 10 which will correspond to the modes retain in the models in section 4, (ii) the intermediate scales for 10 n N , with N such that the N first POD modes capture 95 % of the total energy (N increases with the Rayleigh number, it goes from 250 at Ra = 10 6 to about 5,000 at Ra = 10 8 ); (iii) the small scales for n N . At high Rayleigh number Ra 10 6 3 × 10 6 10 7 3 × 10 7 10 8 uncoupled 0.95 1.10 1.30 1.50 1.75 coupled 1.17 1.24 1.42 1.62 1.85

Table 3. γ factor for uncoupled and coupled cases at different Rayleigh numbers.

(Ra = 3 × 10 7 , Ra = 10 8 ), there is a different behaviour of large scales and intermediate scales, with a very fast decay of the low-order modes followed by a very slow decay in the rest of the spectrum owing to the turbulent nature of the flow. On the contrary, at low Rayleigh number (Ra = 10 6 , Ra = 3 × 10 6 ), the energy contained in the intermediate scales is proportionately more important and a fast decay of the energy of the small scales is observed. Interestingly, the spectrum decay is roughly log-linear in the intermediate range and this will be used in section 6 to model the variations of the POD spectrum with the Rayleigh number.

In the right part of figure 4 is shown the ratio between the coupled POD spectrum and the uncoupled POD spectrum at each Rayleigh number. This ratio is always greater than one (except for one mode at Ra = 10 6 ) which confirms that the coupled POD spectrum captures the energy increase associated with radiation effects. However, this energy increase depends on the Rayleigh number and on the part of the spectrum. At low Rayleigh number, the energy increase is proportionately higher in the large-scale range, while at high Rayleigh number, the energy increase is proportionately higher in the intermediate scale range. In the small-scale range, the spectrum ratio is nearly constant and remains greater than 1. As expected from the discussion in section 2.3, the ratio between coupled and uncoupled total POD energies decreases with the Rayleigh number as radiation effects weaken: it is equal to 1.59 at Ra = 10 6 and equal to 1.19 at Ra = 10 8 .

A last parameter characterising the POD spectra is the factor γ, the values of which are reported in table 3. This factor γ (corresponding to the ratio between mechanical energy and thermal energy, see equation (3.3)) increases with the Rayleigh number and is always greater in the coupled case at a given Rayleigh number, which is consistent with the observations of figure 1.

Spatial eigenfunctions

3.3.1. Ra = 10 6 , uncoupled case Figure 5 shows the first twelve POD eigenfunctions at Ra = 10 6 in the uncoupled case, which correspond to nine spatial structures (three of the modes are doubly degenerated). The figure displays the contribution of each structure to the mean convective heat flux in the vertical direction which can be written as √ Ra n λ n φ θ n φ w n owing to the POD decomposition. Five spatial structures have been already highlighted in previous works [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] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF] and correspond to the first seven POD modes at Ra = 10 7 . They are labelled M , L x/y , D, BL x/y and C. We briefly mention their properties and physical meaning:

(i) The M mode corresponds to the mean flow: it is made of two counter rotating torus and it is thermally stratified.

(ii) The L x and L y modes form a pair of degenerate modes (only the L x mode is shown in figure 5, the L y mode is its image by a rotation of π/2 around the vertical axis). They correspond to a single large-scale roll lying in either x planes or y planes. When combined, the L x and L y modes from a single large-scale diagonal roll.

(iii) The D mode is an 8-roll mode that transports fluid from one corner to the other and strengthens the circulation along the diagonal.

(iv) The BL x and BL y modes (pair of degenerate modes) are constituted of two longitudinal co-rotating structures around the x axis or the y axis. They connect the core of the cell with the horizontal boundary layers.

(v) The C mode is a corner-roll mode which favours planar flow states and was found to promote reorientations between diagonal flow states.

In addition, four other spatial structures are noticeable in figure 5. They are labelled D * , C * , BL x/y * and M * because they share common features with modes D, C, BL x/y and M .

(i) The D * mode is a 4-rolls mode that shares some similarities with the D mode, as it transports fluid from one corner to the other. However, the rolls of the D * mode extend from bottom to top, while the rolls of the D mode are confined in either the upper half or the lower half of the cell.

(ii) The C * mode is another corner-roll mode. Unlike the C-mode these corner rolls extend from top to bottom and are not confined in a half cell.

(iii) The BL x * and BL y * modes (pair of degenerate modes) are made of two longitudinal counter-rotating structures around the x or the y axis. A 2D-version of these POD modes has been highlighted by [START_REF] Podvin | Precursor for wind reversal in a square Rayleigh-Bénard cell[END_REF] in the square cell (symmetrybreaking mode labelled S).

(iv) The M * mode is a single toroidal structure linking the top and bottom walls. The fluid flows up along the adiabatic vertical walls and flows down in the center of the cell.

Ra = 10 8 , uncoupled case

Figure 6 shows the first twelve POD eigenfunctions at Ra = 10 8 in the uncoupled case. Despite the large difference in Rayleigh number, most of the modes observed at Ra = 10 6 are retrieved. One can identify again modes M , L x/y , D, BL x/y , C, as well as modes D * , C * and BL x/y * . However, the mode ordering is not exactly the same. In addition, streamlines and iso-surfaces of convective heat flux in figure 6 are closer to the walls than those in figure 5, as boundary layers are thinner when the Rayleigh number increases. The M * mode is no longer observed, but a new mode, the L z mode, appears. This mode is constituted of a single roll lying in z plane. This is a pure mechanical mode ( φ u 1) and one can note that the convective heat flux associated with the L z mode is much weaker than that of the other modes.

Other coupled and uncoupled cases, similarities

For the other cases studied (other Rayleigh numbers, coupled with radiation), the first twelve POD modes belong to the set of the thirteen POD modes described above: M , L x/y , D, BL x/y , C, D * , C * , BL x/y * , M * and L z . The list of the modes, ordered according to their eigenvalue, is given in table 4. A remarkable feature is that the eigenfunctions are globally preserved when radiation is taken into account, although the associated eigenvalues are higher than in the uncoupled case. In the range 10 7

Ra 10 8 , there are very few differences in the mode ordering, whatever the coupling conditions. It can be noted than the L z mode becomes more important when the Rayleigh number increases and when radiation is considered. In the range 10 6 Ra 3 × 10 6 , the mode ordering is less stable. Indeed, the eigenvalues of the modes are closer to each other, given the 10 6

3 × 10 6 10 7 3 × 10 7 10 8

uncoupled coupled uncoupled coupled uncoupled coupled uncoupled coupled uncoupled coupled slower decay of the POD spectrum at low Rayleigh number. The ordering observed in the coupled case at Ra = 10 6 is the most different from the others: the D mode is no longer in fourth position, the BL x/y modes are also downshift far from the first modes and both M * and L z modes are present. This is not surprising as this case is the only one associated with quasi-stable planar flow states. The question arising is how close are two eigenfunctions of the same nature but belonging to different POD bases (for instance the D mode in the uncoupled case at Ra = 10 6 and the D mode in the coupled case at Ra = 3 × 10 7 ). To answer it, we have taken the POD basis of the uncoupled case at Ra = 10 7 as a reference and we have computed, for a given POD basis B, the differences compared to this reference using the following error measure

M M M M M M M M M M L x/y L x/y L x/y L x/y L x/y L x/y L x/y L x/y L x/y L x/y D C * D D D D D D D D D * M * D * C * BL x/
e B = 12 n=1 1 12 1 - φ B n (r)φ ref n (r) dr , (3.4)
where φ ref n are the eigenfunctions of the reference POD basis and φ B n are the eigenfunctions of the POD basis B, reordered according to the mode ranking of the reference case. For the L z mode, which is not contained in the reference basis, we have taken the L z mode obtained at Ra = 10 7 in the coupled case as the reference. Values of indicator e B are reported in table 4. The differences compared to the reference basis are rather small and are about 1 % or less, except for the coupled case at Ra = 10 6 . This result will be key for the development of a POD model across the Rayleigh range (see section 6) and for predicting radiative transfer effects.

A last comment on the 13 identified modes can be made regarding their contribution to the global angular momentum and their symmetry properties, given in table 5. The symmetry analysis sheds light on mode interactions in the low-order models and the angular momentum is used to monitor the flow reorientations. The modes contributing to the x and y components of the angular momentum are of course the LSC modes L x and L y but also, in a lesser extent the boundary layer modes BL x and BL y . Furthermore the L x and the BL x modes have the same symmetry properties, as well as the L y and the BL y modes. It denotes strong interactions between these two pairs of modes, corresponding to

mode name Lx Ly Lz Symmetries M - - -Sx Sy Sz S d Lx X - -Sx ASy ASz - Ly -X -ASx Sy ASz - D - - -ASx ASy Sz S d BLx X - -Sx ASy ASz - BLy -X -ASx Sy ASz - C - - -Sx Sy Sz AS d D * - - -ASx ASy ASz S d BLx * - - -Sx ASy Sz - BLy * - - -ASx Sy Sz - C * - - -Sx Sy ASz AS d M * - - -Sx Sy ASz S d Lz - -X ASx ASy Sz AS d
Table 5. Name, angular momentum and symmetry properties of the 13 identified modes. X indicates a non zero angular momentum. S x/y/z/d and AS x/y/z/d denote respectively a symmetry and an anti-symmetry with respect to the planes x = 0.5, y = 0.5, z = 0.5 and x = y.

the connection between the LSC and the boundary layers. Interestingly, all the degenerate modes break the S d symmetry and each pair possesses opposite symmetry properties in S x and S y . In addition, one can note that D * , C * , BL x/y * and M * modes have the same symmetry properties in S x , S y and S d than their companion mode D, C, BL x/y and M , but have opposite symmetry properties in S z .

Uncoupled models

In this section we derive low-order models for the uncoupled case from uncoupled DNS data. We consider two different truncations: a 6-D truncation, with modes L x/y , D, BL x/y and C, and a 11-D truncation, with in addition modes D * , C * , BL x/y * , M * or L z modes. The mode M , corresponding to the mean flow, is always taken as constant.

Construction

POD low-order models are derived from Galerkin projection of Navier-Stokes equations (2.3-2.4) onto a POD basis set. Using the decomposition (3.2), this yields a system of ordinary differential equations for the mode amplitudes a n (t) of the form

da n (t) dt = (L B nm + L D nm )a m (t) + Q nmp a m (t)a p (t) + T n (t), (4.1)
where L B nm and L D nm are linear coefficients associated with buoyancy and diffusion and Q nmp are quadratic coefficients associated with advection. They are directly determined from spatial modes φ n , see appendix A. T n (t) is a closure term, which models the effect of the truncation and corresponds to an equivalent dissipation term. Following [START_REF] Podvin | A large-scale investigation of wind reversal in a square Rayleigh-Bénard cell[END_REF], the closure term is expressed as a combination of linear and cubic terms, and by the addition of noise such that

T n (t) = L A nm   1 + 1 k p 2 |a p (t)| 2   a m (t) + σ n (t), (4.2)
where L A nm are adjustable linear parameters, k is the time average of the energy of the fluctuating modes in the truncation and n (t) is a random noise perturbation of amplitude one.

The adjustable coefficients L A nm were determined for equilibria a n (t) = a eq n corresponding to diagonal states from the DNS (da eq n /dt = 0). We retain the following equilibrium: a eq M =

√

λ M ; a eq L x/y = λ L x/y ; a eq D = √ λ D ; a eq BL x/y = -η BL λ BL x/y and a eq n = 0 for the other modes. Note that boundary layer modes BL x/y are negatively correlated with LSC modes L x/y and their equilibrium value in absolute value |a eq BL x/y | is smaller than the standard deviation of their fluctuations (η BL < 1). The coefficient η BL is about 0.15-0.2 except at Ra = 10 8 where it is much weaker (η BL = 0.05). A large value of η represents a strong connection between the boundary layers and the LSC in the diagonal state. The coefficients L A nm can only be determined unambiguously for modes M , L x,y , D and BL x,y , which are non-zero at equilibrium. For modes which are zero at equilibrium, the sum of linear contributions was adjusted to balance quadratic interactions with the mean mode, in order to avoid excessive energy. We checked that the model behavior did not change when small changes were made in the adjustable parameters.

The noise level was estimated as follows

σ = C 2000 8 λ n 7 2 λ n , (4.3)
where C = 0.016 is a constant, determined at Ra = 10 7 from a scaling analysis [START_REF] Soucasse | Proper orthogonal decomposition analysis and modelling of large-scale flow reorientations in a cubic Rayleigh-Bénard cell[END_REF]. The ratio in equation ( 4.3) gives a rough estimate of the energy transferred from the low-order truncation to the next higher-order scales. This ratio reaches a maximum at Ra = 3 × 10 6 and is lowest for the highest Rayleigh numbers. The noise level is indicated in appendix A.

Time integration

The dynamics of modes L x/y ,D, C and BL x/y have already been discussed in previous work for the case Ra = 10 7 [START_REF] Soucasse | Proper orthogonal decomposition analysis and modelling of large-scale flow reorientations in a cubic Rayleigh-Bénard cell[END_REF]. We checked that the general features of the phase portraits held over the range of Rayleigh numbers considered. In particular, it was established that (i) Quasi-steady states constituting the large-scale circulation are characterized by roll L x/y and boundary layer modes BL x/y connecting the boundary layer and the core region.

(ii) The mode D indicates along which diagonal the large-scale circulation takes place and tends to stabilize the large-scale circulation.

(iii) The corner mode C is essential to reproduce reorientations in the model. The frequency of reorientations also depends on the coupling between the boundary layer modes and the roll modes, as well as on the noise level. A higher intensity of the boundary layer mode makes the large-scale circulation less likely to become unstable.

Figure 7 shows histograms of modes L x and L y at different Rayleigh numbers obtained with the 11-D model and with the DNS. They have been computed from time series of mode amplitudes a Lx and a Ly and represent the probability density function of the flow in the phase space L x -L y . It can be seen that the model captures the change in dynamics with the Rayleigh number: the planar states are more frequently visited at low Rayleigh number and the stability of the diagonal states increases at high Rayleigh number. This is further confirmed by the values of reorientation frequencies. Frequencies of zero crossings for modes L x , L y and D, for the uncoupled DNS and model are reported in table 6. Frequency f D corresponds to the reorientation frequency f r defined in section 2.3, 2). At the higher Rayleigh numbers, we found that it is not always possible to determine reorientation frequencies in the DNS as the available simulation time is relatively small compared to the characteristic time between reorientations. However it was possible to obtain values for the model since the integration time could be easily extended.

Time spectra of roll mode L x and higher-order mode D * are represented in figure 8 and compared with the DNS at Ra = 3×10 6 . A good agreement is observed for frequencies up to 10 -2 . However the model is not able to predict the local increase about the circulation frequency f c ∼ 2×10 -2 and the energetic content at higher frequencies is underpredicted, which is not surprising given the small size of the truncation.

Effects of the additional modes

We study in this section the effects of the additional modes included in the 11-D truncation: modes D * , C * , BL x/y * , M * (for Ra 10 7 ) and L z (for Ra 3×10 7 ). Modes L z and M * respectively correspond to a toroidal and poloidal circulation in the cube. These differences correspond to different interactions with other modes: the evolution of mode L x,y is determined by the interaction of M * with mode BL x,y * , as well as mode L z with mode BL x,y . Moreover, the dynamics of M * depends on interactions between modes L and BL, modes BL and BL * with the same orientation (x or y), as well as interactions of mode C and C * , and D and D * . In contrast, the evolution of L z depends on cross-interactions of modes C and D, as well as C * and D * , and modes L and BL corresponding to different directions (x and y). Phase portraits in figure 9 shows that the correlations observed in the DNS are reproduced by the model. More quantitatively, the correlation coefficient between modes C and product DL z at Ra = 10 8 is equal to 0.66 in both the model and the DNS. At Ra = 10 7 , the correlation coefficient between modes M * and product CC * is equal to 0.5 in the model ad 0.3 in the DNS.

To get further insights on the effects of the additional modes, we have compared the dynamics of the model using either the 11-D truncation or the 6-D truncation. Generally speaking, the results obtained with the 6-D truncation are very similar to those obtained with the 11-D truncation, which means the additional modes do not alter the general dynamics of the reorientations. Reorientation frequencies given in table 6 are retrieved with the 6-D model within the same accuracy. This small influence of the additional modes is further confirmed by a linear stability analysis performed around diagonal equilibria for both truncations at each Rayleigh number. The resulting eigenvalue spectra are given in figure 10. First, it can be noted that the shape of the spectrum is globally preserved with the Rayleigh number. Comparison of the 6-D and the 11-D truncations shows that the 6-D reproduces the essential features of the 11-D truncation, that is: (i) a pair of very stable eigenvalues associated with an eigenvector with a strong component

Predicted coupled models

We derive in this section low-order models that include radiation effects. These models are referred to as predicted coupled models because they only rely on uncoupled DNS data and correspond to an a priori attempt to model radiation effects.

Construction

In the framework of the Boussinesq approximation (weak temperature differences), the temperature dependence of the radiation field can be assumed to be linear. We are therefore able to define a modal-radiative power P rad n (r), corresponding to the radiative response to a temperature eigenfunction φ θ n (r) and such that P rad (r, t) = n a n (t)P rad n (r). Using this decomposition in the energy balance and applying Galerkin projection, a linear radiation term is obtained and the coupled POD model takes the general form

da n (t) dt = (L B nm + L D nm + L R nm )a m (t) + Q nmp a m (t)a p (t) + T n (t), (5.1)
where L R nm are linear radiation coefficients. Definition and computation details for coefficients L R nm and for the modal radiative power P rad n (r) are given in appendix B and can be also found in [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF]. We will assume that the form of the closure law (equation (4.2)) remains valid in the presence of radiation.

Our goal is to build an a priori model of radiative transfer effects based on uncoupled DNS data and uncoupled POD eigenfunctions. Therefore, model coefficients are determined as follows:

(i) Coefficients L B nm , L D nm and Q nmp are taken from the uncoupled model. (ii) Coefficients L R nm are added and computed from uncoupled temperature eigenfunctions φ θ n and uncoupled factor γ. (iii) Adjustable coefficients L A nm are taken from the uncoupled model. This means we assume that the energy transfer from the large scales to the small scales is unchanged.

However, equilibrium values of mode amplitude a eq,rad n are required in the coupled model to estimate the contribution of the mean mode M , to estimate the fluctuating energy k in equation (4.2) and to estimate the noise level. They are not known a priori and need to be predicted. To do so, we use equilibrium relations (da eq,rad n /dt = 0) for diagonal states, with the additional radiation terms, the adjustable coefficients being known. Predicted equilibria ∼ a eq,rad n

for modes M , L x/y , D, and BL x/y are given in table 7 (a eq,rad n = 0 for the other modes) and compared with actual coupled equilibria a eq,rad n and uncoupled equilibria a eq n extracted from the DNS. We can see that the effect of radiative coupling is to increase the energy of the equilibria, by lesser amounts as the Rayleigh number increases. Finally, the noise level in the coupled model σ rad is determined from the energy increase of the large-scale modes compared to the uncoupled case

σ rad σ = n ∼ a eq,rad n 2 n (a eq n ) 2 (5.2)
Noise level values are reported in appendix A. As expected, the noise level with radiation σ rad is always greater than the noise level without radiation due to the energy increase.

Mode energies

The predicted coupled model was then integrated in time, using an 11-D truncation, leading to predicted amplitudes a rad i . The predicted energy (or predicted eigenvalue) of the modes with radiative coupling ∼ λ rad was then determined using:

• the ratio η n determined above for the non-zero equilibrium values η n = ( ∼ a eq,rad n /a eq n ) 2 (with n corresponding to modes M, L x,y , D, BL x,y ) • the ratio ηn predicted from the model ηn = (a rad n ) 2 / a 2 n for modes which are zero at equilibria (with n corresponding to modes D * , C, C * , BL *

x,y , L z or M z ). Figure 11 compares for each mode and each Rayleigh number the predicted energy ∼ λ rad n with the eigenvalues λ n and λ rad n found in the uncoupled and coupled DNS. In addition, the values of the total energy contained in the first twelve modes are reported in table 8. Except for the lowest Rayleigh number, corresponding to a large radiative increase which is overpredicted by the model, the total energy increase is correctly predicted by the model (within about 10 % maximum). Overall, the model tends to overpredict the amplitude of the BL x/y modes (which are associated with the largest relative increase) and underpredict that of the L x/y modes. However the agreement is generally good, although the ordering of the modes observed in the coupled DNS is not always captured.

Figure 12 shows histograms of modes L x and L y at different Rayleigh numbers obtained with the coupled predicted model and the coupled DNS. At Ra = 10 6 , the radiative coupling makes diagonal states unstable, and the large-scale circulation in the simulation tends to become parallel to the cavity sides. This change in dynamics is qualitatively predicted by the model: it can be seen in the histograms of figure 12 that the system spends considerable time near the roll states (L x = 0 or L y = 0). This will be further 9. Average frequencies fn = 1/τn in the coupled DNS and in the coupled model where τn is the average time between zeros of an (restricted to times larger than 5τc ). The value N/A is indicated when it was not possible to obtain a value (less than two switches observed at higher Rayleigh numbers, many rapid switches at Ra = 3 × 10 6 ). Integration time for the model is 40,000 time units.

confirmed by the linear stability analysis in section 5.4. At Ra = 3 × 10 6 , the coupled system seems to spend more time near the roll states although the diagonal states are still dominant. This is also correctly captured by the predicted model. Changes due to radiative coupling are less important at Ra 10 7 in both the model and the DNS. At Ra 3 × 10 7 , very few reorientations are observed and the amplitude of the oscillations about the equilibrium states decrease.

Reorientation frequencies

Table 9 compares frequencies of zero crossings for modes L x , L y and D, for the coupled DNS and the predicted coupled model. Generally speaking, a good agreement is observed between the predicted model and the coupled DNS (discrepancies of 10 % to 30 %, of the order of the uncertainty on DNS frequencies), except at Ra = 10 8 where the reorientation frequency is underestimated. We note that results are not given for the coupled case at Ra = 10 6 since the diagonal states are no longer observed. We note that at Ra = 3 × 10 6 large and fast in the D mode do not make it possible to determine a frequency in both the coupled simulation and the model. As in section 4.2, values for the model reported in the table were obtained for an integration time of 40,000 time units.

Model stability at Ra = 10 6

A remarkable feature of the predicted coupled model is its ability to predict the loss of stability of the diagonal rolls for the benefit of the planar rolls at Ra = 10 6 . In fact, linear stability analysis of the diagonal equilibria for both 6-D and 11-D truncations (see figure 13) shows that the least stable pair of eigenvalues associated with mode C becomes markedly unstable with radiative coupling. Integration of the models without noise shows that in the presence of radiation the now unstable fixed point destabilises towards another fixed point in the L x = 0 or L y = 0 space (also in figure 13) corresponding to roll modes. We emphasize that both 6-D and 11-D truncations provide very similar results, indicating that the essential dynamics are captured by the 6-D truncation.

At other Rayleigh numbers, we found that for both 6-D and 11-D truncations the pair of eigenvectors associated with the corner mode C became less stable when radiation was taken into account, but their eigenvalue remained negative.

Generalised models over the Rayleigh number range

In this section, we develop a generalised model, for both uncoupled and coupled cases, that infers the dynamics over the Rayleigh number range from uncoupled DNS data at Ra = 10 7 and energy scaling laws. We restrict the models to a 6-D truncation with modes L x/y , D, BL x/y and C.

Uncoupled case

The basic assumption to build the generalised model is that the variations of the eigenfunctions can be neglected over the Rayleigh number range. The key ingredient is thus to estimate correctly the distribution of the POD spectrum in order to determine the equilibria, the adjustable coefficients and the noise level. The form of the uncoupled generalised model is the same as that of the uncoupled model (equation (4.1) together with closure law (4.2)). Model parameters are determined as follows:

(i) Linear and quadratic coefficients L B nm , L D nm and Q nmp are directly determined from a set of POD eigenfunctions φ n = {φ u n , γφ θ n }. Here, φ u n and φ θ n are assumed to be constant with the Rayleigh number and taken at Ra = 10 7 . The coefficient γ 2 is assumed to vary with the Rayleigh number according to the following law γ 2 ∼ 0.032Ra 1/4 , (6.1) which is consistent with the scaling of the turbulent fluctuations in the bulk θ * ∼ Ra -0.14 [START_REF] Castaing | Scaling of hard thermal turbulence in Rayleigh-Bénard convection[END_REF]).

(ii) The adjustable coefficients L A nm are determined from equilibrium relations for a diagonal state. Equilibrium values a eq n strongly depend on the Rayleigh number and are estimated using the following ansatz for the energy λ n of the large-scale modes n ∈ {M, L x/y , D, BL x,y , C}:

λ n = h n log(Ra) + g n , (6.2)
which is a good approximation at all Rayleigh numbers except the lowest one (Ra = 10 6 ), as it can be seen in figure 14. The lowest Rayleigh number was therefore not taken into account when fitting the coefficients h n and g n . Note that the ratio between the equilibrium value of modes BL x/y and their amplitude was taken to be constant and equal to η BL = 0.2. (iii) As in section 4, the noise level is determined using equation (4.3), from the energy of the large scales and the energy of the intermediate scales. We have highlighted in section 3.2 that the spectrum decay of the intermediate scales is log-linear and we adopt the following modelling for their energy Ra) with log(C(Ra)) = 1.9 -0.5 log(Ra) α(Ra) = -2.0 + 0.07 log(Ra) (6.3)

λ n = C(Ra)n α(
The fit was determined from the variations of the spectrum λ n with 8 n N/2, with N such that the N first POD modes capture 95 % of the total energy. We have checked that the asymptotic decay was correctly captured. It can also be noted in figure 14 that the intermediate-scale fit provides a correct estimation of the energy of mode C, with a discrepancy between the large-scale and the intermediate-scale fit of 25 % maximum at the lowest Rayleigh number. Noise levels are reported in appendix A. The uncoupled generalised model was integrated in time for the different Rayleigh numbers. The measured reorientation frequencies are given in table 10. Overall, the model correctly captures the amplitude of the modes and the associated time scales, except at Ra = 10 6 . At this Rayleigh number, the energy of the large-scale modes is underestimated by the large-scale fit, which leads to an overestimation of the noise level and an overestimation of the reorientation frequency. At Ra = 10 8 , no reorientation event has been detected by the model despite a large integration time. 15 and compared with values given by the coupled DNS and the coupled predicted model. In addition, the energy increase for each largescale mode M, L x,y , D, BL x,y is given in figure 16. In terms of total large-scale energy, a very good agreement is observed for the generalized model (discrepancies about 5 %), even at Ra = 10 6 . In this case, as figure 16 shows, it is due to the compensation of estimation errors (M is overpredicted by 25 %, while modes L are underpredicted by 30 %). At other Rayleigh numbers, the evolution of the energy of the different modes is relatively well predicted (between 5 % to 20 %), except for an underprediction of modes L. Overall a good agreement is observed between the different models.

The measured reorientation frequencies obtained with the coupled generalised model are given in table 11 and compared with the coupled DNS. Generally speaking, the model fairly reproduces the coupled DNS results (discrepancies of 10 % to 30 %) and predicts the increase of the reorientation frequency with radiation, which tends to subside with the Rayleigh number. However, it should be pointed that the change in dynamics observed at Ra = 10 6 (loss of stability of the diagonal rolls) is not predicted with the generalized coupled model, which is consistent with the lack of validity of the modelling assumptions at this Rayleigh number.

Conclusion

DNS of Rayleigh-Bénard convection in a cubic cell for a radiating air/H 2 O/CO 2 mixture has been performed in the Rayleigh number range Ra ∈ [10 6 -10 8 ]. When radiative transfer is ignored (uncoupled case), the flow dynamics is always characterised by four quasi-stable diagonal states, with occasional brief reorientations between them. The stability of the diagonal states increases and the reorientation frequency decreases with the Rayleigh number in the range 3 × 10 6 Ra 10 8 . When radiative transfer is taken into account (coupled case), the kinetic energy of the flow and the temperature fluctuations increase. At Ra = 10 6 , quasi-stable planar states are observed while for 0.9 1.1 0.7 0.9 1.7 1.5 3 × 10 7 0.4 0.6 0.4 0.5 1.3 1.0 10 8

1.1 0.9 N/A 1.0 1.1 1.4

Table 11. Reorientation frequencies obtained with the coupled generalised model. Integration time for the model is 40,000 time units.

Ra 3 × 10 6 diagonal states are retrieved. Radiation seems to increase the reorientation frequency.

A POD analysis of DNS results has been conducted and shows that the first twelve POD eigenfunctions are globally preserved over the whole Rayleigh number range, whatever the coupling conditions. However, the POD eigenvalues are higher in the coupled case at each Rayleigh number. The first seven modes identified in previous works at Ra = 10 7 [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] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF] are retrieved (modes M , L x/y , D, BL x/y and C), though not necessarily in the same order depending on the Rayleigh number and the coupling conditions. Additional modes have been highlighted: (i) modes C * , D * , BL * x/y that share common features with the previous set ; (ii) modes M * and L z associated with toroidal and poloidal circulations and present either at low or high Rayleigh number, respectively.

POD-based low-order models have been derived for the uncoupled case in an a posteriori fashion, from uncoupled DNS data at each Rayleigh number. Truncations using either the first twelve modes or the set of key seven modes have been considered. It is shown that the set of key seven modes is sufficient to model the flow dynamics and reproduces correctly reorientation frequencies. For the coupled case, low-order models including radiation have been derived in an a priori fashion, solely based on uncoupled DNS data. Radiative transfer effects are linear in the model owing to the weak temperature differences. A remarkable feature of this coupled predicted model is its ability to predict the loss of stability of the diagonal states at Ra = 10 6 . The coupled predicted model also foresees the energy increase and reorientation frequency increase associated with radiative transfer effects. Finally, a generalised model for both coupled and uncoupled cases has been proposed in order to capture variations in Rayleigh number from uncoupled DNS data at Ra = 10 7 and scaling laws of the uncoupled POD eigenspectrum. This generalised approach yields satisfactory results but is not able to predict the change in dynamics at Ra = 10 6 , even while it captures changes in energies. This shows the limit of POD model predictions when extrapolating in the range of parameters far from the original data.

In summary, this study shows that radiative transfer effects on the large-scale dynamics of Rayleigh-Bénard convection can be predicted satisfactorily from uncoupled simulation data over a wide range of Rayleigh number. The coupled predicted POD model derived in this paper is a valuable tool to investigate radiative transfer effects, given the high computational cost of radiation calculations. In addition, the evolution of the dynamics with the Rayleigh number can be foreseen, at least in part, from the dynamics at Ra = 10 7 and using scaling laws of the POD spectrum. Nevertheless, this study is restricted to the cubical geometry, given that the cavity shape alters the POD eigenfunctions and the LSC dynamics. Extension of these models to account for aspect ratio dependencies will be the topic of future works. 1.5 1.3 1.1 0.9 σ rad (coupled predicted model) 1.9 1.5 1.2 1.0 0.8 σ rad (coupled generalized model) 2.6 1.9 1.4 1.1 0.9

Table 13. Noise level σ ×10 3 for the different models and Rayleigh numbers.

a eq M = √ λ M ; a eq L x/y = λ L x/y ; a eq D = √ λ D ; a eq BL x/y = -η BL λ BL x/y and a eq n = 0 for the other modes. Coefficients L A nm are thus computed from equation (A 6) with L R nm = 0. In the coupled case, equilibrium values are not known. However, we assume that radiation does not affect the energy transfer from the large scales to the small scales and take the adjustable coefficients from the uncoupled models. Equilibrium values with radiation can be thus predicted from equilibrium relations (A 6), the adjustable coefficients being known. More details on the equilibrium value computation in the coupled case are given in [START_REF] Soucasse | Reduced-order modelling of radiative transfer effects on Rayleigh-Bénard convection in a cubic cell[END_REF].

Finally the noise level σ is determined either from equation (4.3) using the POD spectrum (taken from the DNS or estimated from the fit in the generalised uncoupled model) or from equation (5.2) using equilibrium values in the coupled models. Noise level values are reported in table 13.

Figure 1 .

 1 Figure1. Total kinetic energy e k , total thermal energy e θ , conductive flux at the wall q cond and convective flux in the core qconv as a function of the Rayleigh number for coupled (red squares) and uncoupled (black triangles) cases.

Figure 2 .

 2 Figure 2. Time evolution of x and y components of angular momentum Lx (blue lines) and Ly (red line).

Figure 3 .

 3 Figure 3. Instantaneous flow field at Ra = 10 6 and t = 5, 000 for uncoupled (left) and coupled (right) cases. Streamlines and isotherms θ = {0; ±0.05; ±0.1}. In the uncoupled case (left) a diagonal state is observed (Lx > 0, Ly > 0) while in the coupled case (right) a planar state is observed (Lx 0, Ly < 0).

Figure 4 .

 4 Figure 4. POD eigenspectrum obtained from uncoupled simulations, normalized such that n λn = 1 (left) and POD eigenspectrum ratio between coupled and uncoupled results (right).

Figure 5 .

 5 Figure 5. First twelve modes at Ra = 10 6 (uncoupled case). Streamlines and iso-surfaces of the contribution to the mean convective heat flux φ θ φ w = 0.25, coloured by mode temperature.For degenerate modes, only the x-oriented one is displayed. Colour map for mode temperature ranges from -0.5 (blue) to 0.5 (red).

Figure 6 .

 6 Figure6. First twelve modes at Ra = 10 8 (uncoupled case). Streamlines and iso-surfaces of the contribution to the mean convective heat flux φ θ φ w = 0.25 (0.025 for the Lz mode), coloured by mode temperature. For degenerate modes, only the x-oriented one is displayed. Colour map for mode temperature ranges from -0.5 (blue) to 0.5 (red).

Figure 7 .

 7 Figure 7. Histograms of modes Lx and Ly: DNS (left) and model (right). Uncoupled case.Integration time for the model equal to ∆t (DNS).

Figure 8 .

 8 Figure 8. DNS and model spectra at Ra = 3 × 10 6 for mode Lx (left) and mode D * (right).

Figure 9 .

 9 Figure 9. Top: Phase portrait of M * and CC * in the DNS (left) and in the model (right) at Ra = 10 7 ; Bottom: Phase portrait of C and in the DNS (left) and in the model (right) at Ra = 10 8 .

Figure 10 .

 10 Figure 10. Linear stability analysis spectrum around the (Lx, Ly) equilibria. Comparison between the 6-D and 11-D truncations of the uncoupled model.

Figure 11 .

 11 Figure 11. Predicted energies compared with POD eigenvalues in the uncoupled and coupled DNS. Energies are ranked according to the mode ordering in the uncoupled simulation.

Figure 12 .

 12 Figure 12. Histograms of modes Lx and Ly: DNS (left) and model (right). Coupled case.Integration time for the model equal to ∆t (DNS).

Figure 13 .

 13 Figure 13. Left: Linear stability analysis spectrum around the (Lx, Ly) diagonal equilibria for the predicted coupled model at Ra = 10 6 . Right: Trajectory in the (Lx, Ly) space of the coupled predicted model without noise from unstable diagonal equilibria: 6-D truncation (black dashed lines) and 11-D truncation (red dotted lines). Symbols indicate stable fixed point for both coupled and uncoupled models.

Figure 14 .

 14 Figure14. Evolution of the energy λn of the large-scale modes given by the DNS and the large-scale fit (6.2).

Figure 16 .

 16 Figure 16. Relative increase of the amplitude of the large-scale modes due to radiation with the Rayleigh number -comparison of the DNS, the coupled predicted model and the coupled generalised model.

Table 2 .

 2 Circulation frequency fc and reorientation frequency fr for the uncoupled and coupled cases at different Rayleigh numbers. The value N/A is indicated when it was not possible to obtain a value (less than two switches observed at Ra = 10 8 -uncoupled case, many rapid switches at Ra = 3 × 10 6 -coupled case) or when planar flow states are observed at Ra = 10 6coupled case.

		fc × 10 2	fr × 10 3
	Ra	uncoupled coupled uncoupled coupled
	10 6	2.12	2.86	1.3 ± 0.3	N/A
	3 × 10 6	2.06	2.47	2.2 ± 0.4	N/A
	10 7	2.05	2.33	1.4 ± 0.4 1.7 ± 0.4
	3 × 10 7	2.03	2.27	0.4 ± 0.3 1.3 ± 0.5
	10 8	2.03	2.19	N/A	1.1 ± 0.5

Table 4 .

 4 List of the first twelve POD modes for each case studied, ordered from top to bottom. The last line of the table gives the error measure e B in percentage between the POD basis and the reference POD basis (Ra = 10 7 , uncoupled) as defined in equation 3.4.

	y	BL x/y	BL x/y	BL x/y	BL x/y	Lz

Table 6 .

 6 Average frequencies fn = 1/τn in the uncoupled DNS and in the uncoupled model where τn is the average time between zeros of an (restricted to times larger than 5τc ). The value N/A is indicated when it was not possible to obtain a value (less than two switches observed at higher Rayleigh numbers). Integration time for the model is 40,000 time units.as the D mode changes sign at each reorientation. Values for the model reported in the table were obtained for an integration time of 40,000 time units. Uncertainties on model frequencies are about 4 to 8 times smaller than uncertainties on DNS frequencies, according to the respective integration times. It can be noted in table 6 that the model captures the decrease of reorientation frequencies with the Rayleigh number and fairly reproduces the frequencies observed in the DNS. Discrepancies between DNS and model results are of the order of the uncertainty on DNS frequencies (see table

			uncoupled model
	Ra	fL x × 10 3	fL y × 10 3	fD × 10 3
		DNS model DNS model DNS model
	10 6	0.8	1.1	0.9	1.1	1.3	1.7
	3 × 10 6 1.4	1.4	1.0	1.1	2.2	1.7
	10 7	0.8	0.8	0.7	1.0	1.4	1.3
	3 × 10 7 N/A 0.4 N/A 0.3	0.4	0.8
	10 8	N/A N/A N/A 0.05 N/A 0.14

Table 8 .

 8 Predicted energy (∼) of the first twelve modes compared with eigenvalues in the uncoupled simulation λn and in the coupled simulation λ rad n . Relative ratios compared to the uncoupled case are indicated in parentheses.Table

	Ra	10 6	3 × 10 6	10 7		3 × 10 7	10 8
	12 n=1 λn 0.0320 12 n 0.0548 (1.71) 0.0425 (1.51) 0.0401 (1.44) 0.0353 (1.28) 0.0323 (1.21) 0.0281 0.0277 0.0275 0.0267 n=1 λ rad ∼ rad 12 n=1 λ n 0.0986 (3.08) 0.0405 (1.44) 0.0380 (1.37) 0.0388 (1.41) 0.0320 (1.20)
			coupled predicted model
		Ra	fL x × 10 3	fL y × 10 3	fD × 10 3
			DNS model DNS model DNS model
		3 × 10 6 1.6	1.5	1.2	1.4 N/A N/A
		10 7	0.9	1.2	0.7	1.0	1.7	1.7
		3 × 10 7 0.4	0.5	0.4	0.5	1.3	0.9
		10 8	1.1 N/A N/A 0.1	1.1	0.2

Table 10 .

 10 Reorientation frequencies obtained with the uncoupled generalized model. Integration time for the model is 40,000 time units.Following the procedure described in section 5.1, radiation effects are incorporated in the generalised model. Model parameters for the coupled generalised model are thus determined as follows:(i) Coefficients L B nm , L D nm and Q nmp are taken from the uncoupled generalised model. (ii) Coefficients L R nm are added and computed from uncoupled temperature eigenfunctions φ θ n at Ra = 10 7 and γ factor modelled by equation (6.1). (iii) Adjustable coefficients L A nm are taken from the uncoupled generalised model. This allows us to estimate equilibrium values

		uncoupled generalised model
	Ra	fL x × 10 3	fL y × 10 3	fD × 10 3
		DNS model DNS model DNS model
	10 6	0.8	1.7	0.9		1.6	1.3 N/A
	3 × 10 6 1.4	1.4	1.0		1.2	2.2	2.0
	10 7	0.8	0.6	0.7		0.6	1.4	1.1
	3 × 10 7 N/A 0.2 N/A 0.2	0.4	0.4
	10 8	N/A N/A N/A N/A N/A N/A
		6.2. Coupled case
			∼ a	eq,rad n	and energies	∼ λ	rad n	for the modes in the
	truncation (except mode C which is zero at equilibrium).
	(iv) Noise level is determined using equation (5.2).
	Differences between all the models (uncoupled model, predicted coupled model, uncou-
	pled generalised model and coupled generalised model) are highlighted in table 12 of
	appendix A.
	The total energy of the large-scale modes	7 n=1	∼ rad λ n	predicted by the coupled gener-
	alised model is represented in figure

  Total energy of the large-scale modes 7 n=1 as a function of the Rayleigh number -comparison of the DNS, the coupled predicted model, and the coupled generalised model.

				0.1
						uncoupled DNS
				0.08		predicted coupled model generalized coupled model
			n T n	0.06		coupled DNS
				0.04
				10 6 0.02			10 7	10 8
							Ra
	1 1.2 1.4 1.6 1.8 2.2 Figure 15. 10 6 M rad /M 2	10 7 coupled DNS predicted coupled model generalized coupled model	10 8	L rad /L	10 6 1 1.1 1.2 1.3 1.4 1.6 1.5	10 7 coupled DNS predicted coupled model generalized coupled model	10 8
			Ra					Ra
		1.6
			coupled DNS					coupled DNS
		1.5	predicted coupled model generalized coupled model			2	predicted coupled model generalized coupled model
	/D	1.4				/BL
	D rad	1.2 1.3				BL rad	1.5
		1.1					1
		10 6 1	10 7		10 8		10 6	10 7	10 8
			Ra					Ra

Table 12 .

 12 Model parameter determination.Ra 10 6 3 × 10 6 10 7 3 × 10 7 10 8

		uncoupled model	predicted coupled		uncoupled	coupled generalised
			model		generalised model	model
	L B nm , L D nm , L R nm , Qnmp	computed from uncoupled DNS data at each	computed from uncoupled DNS data at each	computed from uncoupled DNS data at Ra = 10 7 -	computed from uncoupled DNS data at Ra = 10 7 -
		Rayleigh number	Rayleigh number	eq. (6.1) for γ	eq. (6.1) for γ
		(L R nm = 0)			(L R nm = 0)
	L A nm	determined from	taken from the		determined from	taken from the
		equilibrium	uncoupled model	equilibrium	uncoupled
		relations (A 6)			relations (A 6)	generalised model
	a eq	determined from	determined from	determined from	determined from
		POD spectrum	equilibrium		POD spectrum	equilibrium
			relations (A 6)		large-scale fit,	relations (A 6)
					eq. (6.2)
	σ	computed from	computed from		computed from	computed from
		eq. (4.3) using	eq. (5.2) using		eq. (4.3) using	eq. (5.2) using
		DNS POD	equilibria a eq		fitted POD	equilibria a eq
		spectrum			spectrum
		σ (uncoupled model)	1.2	1.4	1.2	1.0	0.8
		σ (uncoupled generalized model) 1.5

† Email address for correspondence: laurent.soucasse@centralesupelec.fr

Acknowledgements

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 Suprieure Paris-Saclay supported by CNRS and Rgion Île-de-France (http://mesocentre.centralesupelec.fr/).

Declaration of Interests

The authors report no conflict of interest.

Appendix A. Model coefficients

The low-order models for the time evolution of the POD mode amplitude a n (t), derived in sections 4-6, all take the general form

with L R nm = 0 when radiation effects are ignored. The way to determine the model coefficients is summarised in table 1 for the uncoupled model (section 4), the predicted coupled model (section 5) and the uncoupled and coupled generalised models (section 6).

Coefficients L B nm , L D nm , L R nm and Q nmp are respectively associated with buoyancy, diffusion, radiation and advection and are defined as

where P rad m is a modal-radiative power associated with mode m, whose definition is given in appendix B. These coefficients can be fully computed from a set of POD eigenfunctions {φ u , γφ θ }.

The adjustable coefficients L A nm can be determined from equilibria a n (t) = a eq n corresponding to a diagonal state of the flow (da eq n /dt = 0). Equilibrium relations write (L B nm + L D nm + L R nm + 2L A nm )a eq m + Q nmp a eq m a eq p = 0. (A 6)

In the uncoupled case, equilibrium values a eq n are taken from the POD spectrum (uncoupled model) or from the large-scale fit of equation (6.2) (uncoupled generalised model):

Appendix B. Computation of radiative terms

This appendix details the computation of the linear coefficients L R nm (equation (A 4)) associated with radiation effects in the low-order models.

Owing to the weak temperature gradients, the Planck function I • ν (T (r, t)) can be linearised around the mean temperature T 0 and the POD decomposition of the temperature field can be introduced

Because of the linearity of the radiative transfer equation (2.7), the radiative intensity field I ν (Ω, r, t) can be decomposed similarly

where ψ θ ν,n (Ω, r) is a modal-intensity field, associated with the POD temperature eigenfunction φ θ n (r) and satisfying the following transport equation

and boundary condition for It can be noted that the total radiative power P rad (r, t) is the sum of all modal-radiative powers, weighted by the associated POD coefficient P rad (r, t) = n a n (t)P rad n (r). (B 6)