Amine Benchrifa

Andrew Buchan

Romain Ceolato

Steven Dargaville

Maxime Delort- Laval

Bruno Dias

Frédéric Doumenc

Béatrice Guerrier

Frédérique Laurent-Nègre

Lorenn Le Turnier

Thierry Magin

Sophie Mergui

Peder Olesen

Christopher Pain

Lucas Paulien

Baraa Qaddah

Jean-Baptiste Salmon

James B Scoggins

Clara Velte

Laurent Zimmer

Pascal Boulet

Cyril Caliot

Pedro Coelho

Franck Enguehard

Chrisophe Ménézo

Jörg Schu- Macher. J'

Bérengère Podvin

Philippe Rivière

Anouar Soufiani

mes travaux de recherches ces dix dernières années.

Titre : Modélisation des transferts radiatifs : propriétés radiatives, simulation numérique et interactions couplées

Mots clés : transfert radiatif, modélisation des propriétés radiatives, simulations numériques couplées, convection turbulente, modélisation d'ordre réduit Résumé : Ce manuscrit d'habilitation synthétise des travaux portant sur la modélisation du transfert radiatif à l'aide de la simulation numérique. La première partie porte sur la modélisation des propriétés radiatives des milieux gazeux ou particulaires. Un modèle statistique à bande étroite hybride est établi pour prédire le rayonnement hors équilibre impactant les véhicules spatiaux lors des entrées atmosphériques. Par ailleurs, une méthodologie pour la modélisation de la rétrodiffusion d'agrégats fractals de suie pour les applications LIDAR est présentée.

La deuxième partie étudie deux stratégies pour traiter la dépendance angulaire du champ de rayonnement dans les simulations de transfert radiatif : un algorithme adaptatif-cible utilisant une base angulaire d'ondelettes de Haar et un modèle d'ordre réduit angulaire basé sur une décomposition en modes propres orthogonaux.

La troisième partie est consacrée à l'étude des interactions couplées entre le rayonnement thermique et la convection turbulente. Des simulations numériques directes couplées sont entreprises et analysées pour des mélanges air/H2O/CO2 dans des cellules cubiques de Rayleigh-Bénard. Un modèle de bas ordre est développé : il capture la dynamique à grande échelle de l'écoulement et les effets associés au rayonnement du gaz. Enfin, les interactions couplées convection-rayonnement dans des grottes peu profondes, cavités peu ventilées et humides sous la surface du sol, sont étudiées. Title: Radiative transfer modelling: radiative property modelling, numerical simulation and coupled interactions Keywords: radiative transfer, radiative property modelling, coupled numerical simulations, turbulent convection, reduced-order modelling Abstract : This professorial thesis summarises research contributions to the field of radiative transfer modelling using numerical simulation. The first part focuses on modelling the radiative properties of gaseous or particulate media. A Hybrid Statistical Narrow Band model is derived to predict nonequilibrium radiation in atmospheric entries of space vehicles. Further, a methodology for the backscattering modelling of soot fractal aggregates in LIDAR applications is presented.

The second part investigates two strategies to address the angular dependence of the radiation field in radiative transfer simulations: a goal-based adaptivity algorithm using Haar-wavelet angular discretisation and an angular reduced-order model using proper orthogonal decomposition.

The third part is dedicated to the study of coupled interactions between thermal radiation and turbulent convection. Coupled direct numerical simulations are performed and analysed for air/H2O/CO2 mixtures in cubic Rayleigh-Bénard cells. A low-order model is then derived: it captures the large-scale dynamics of the convective flow and associated radiation effects. Finally, coupled convection-radiation interactions in shallow caves, that are weakly ventilated and humid cavities below the ground surface, are investigated.

Radiative transfer modelling:

radiative properties, numerical simulation and coupled interactions Radiative transfer is an energy transfer mode associated with electromagnetic waves (or photons). The heat we feel in front of the sun or in front of a fire are everyday examples of thermal radiation, which, unlike the other energy transfer modes, conduction and convection, can travel over a long distance. Understanding radiative transfer is of primary importance for climate sciences as radiation interactions between the sun, the atmosphere and the Earth mostly determine the temperature of our planet. Radiative transfer plays a key role in various applications in the field of energy production or energy efficiency: solar panels, solar furnaces, combustion systems, high-temperature processes or energy transfer in buildings. Also, radiation is a major tool for physicists to probe the matter, observe the Earth or get information from space. Many medical diagnostics are based on radiation transport.

The physical model governing radiative transfer propagation is the Radiative Transfer Equation (RTE), which is a Boltzmann transport equation on the radiative intensity field. The radiative intensity I ν (r, Ω) is the energy flux (per unit surface) transported by electromagnetic waves (or photons) at position r, along propagation direction Ω and at wavenumber ν. For propagating medium of optical index equal to one, the RTE writes

Ω •∇ I ν (Ω, r) = η ν (r) emission + σ ν (r) 4π p ν (r, Ω, Ω ′ )I ν (Ω ′ , r)dΩ ′ scattering source -κ ν (r)I ν (Ω, r) absorption -σ ν (r)I ν (Ω, r) scattering extinction , (1.1)
where η ν , κ ν , σ ν and p ν are respectively the emission coefficient, the absorption coefficient, the scattering coefficient and the phase function, characterising the radiative properties of the matter. When the medium is at Local Thermal Equilibrium (LTE) the emission coefficient writes η ν = κ ν I • ν (T ), where I • ν is the Planck function at temperature T . A boundary condition must specify the incoming intensity at the contour of the domain of interest to complete Eq. (1.1). In the case of an opaque diffuse reflecting boundary of emissivity ε ν at LTE it writes

I ν (Ω, r b ) = ε ν (r b )I • ν (T (r b )) + 1 -ε ν (r b ) π Ω ′ •n<0 I ν (Ω ′ , r b )|Ω ′ •n |dΩ ′ , (1.2)
for boundary points r b and for propagation directions Ω • n > 0, n being the boundary normal directed towards the inside of the spatial domain. Prescribed incoming intensity or pure specular reflecting boundaries will also be considered in the following.

Challenges of numerical simulation

The numerical simulation of radiative transfer aims at increasing the scientific knowledge of this phenomenon and mastering natural or engineering processes where radiation occurs. But solving the RTE is a numerical challenge because of the intricate dependence of the radiative intensity on wavelength, propagation direction, space, and potentially time. We briefly discuss here the main obstacles.

• Spectral dependence The absorption spectrum of common radiating molecules in the infrared such as water vapour or carbon dioxide are made of millions of lines associated with energy transitions between internal molecular energy levels [START_REF] Goody | Atmospheric radiation: theoritical basis[END_REF]. In ionised plasma, absorption spectra are further complicated by energy transitions between electronic energy levels, bound-free and free-free radiative processes. In propagating media involving dropplets or solid particles, absorption and scattering properties are a complex function of the size, the shape and the optical index of each particle.

INTRODUCTION

• Angular dependence The angular variation of the radiation field strongly depends on the optical thickness ζ ν (L) = (κ ν + σ ν )L, where L is a characteristic optical path length L. At low optical thicknesses (ζ ν << 1) radiation propagates directionally and requires accurate high-order discrete ordinates methods or ray tracing algorithms to capture it. At large optical thicknesses (ζ ν >> 1), radiation propagates diffusely and approximate methods such as the Rosseland approximation or the spherical harmonics expansion can be appropriate [START_REF] Modest | Radiative heat transfer[END_REF][START_REF] Howell | Thermal Radiation Heat Transfer[END_REF]. For the propagating media described above, the optical thickness varies of several orders of magnitude, so both regimes must be captured.

• Space dependence The spatial discretisation of the radiative intensity depends on the optical thickness of the mesh (each spatial cell should be optically thin) and on the length scale over which the matter properties (temperature, pressure, composition) vary. In turbulent flows, all length scales up to the Kolmogorov scale must be resolved in theory, which is usually not feasible because of computational time and memory constraints. Simplified approaches based on time averaging or space filtering are used to release these discretisation constraints but require the modeling of closure terms referred to as turbulence-radiation interactions [START_REF] Coelho | Numerical simulation of the interaction between turbulence and radiation in reactive flows[END_REF].

• Time dependence Radiation propagates at the speed of light, which is often much larger than any reference velocity characterising the state of the matter, so that radiation propagation can be considered instantaneous. However, the radiative properties of the matter vary in time and the intensity field varies accordingly. In a turbulent medium, time discretisation constraints meet those of the space discretisation, highlighted above.

Moreover, the numerical difficulties specific to radiation add to those of other interacting physical processes, such as the heat transport in turbulent flows, solid materials and particulate media, or the chemistry in reacting media. Examples of multiphysics problems involving thermal radiation are the design of fire safety equipments [START_REF] Parent | Radiation attenuation and opacity in smoke and water sprays[END_REF], the development of air receivers in concentrated solar power applications [START_REF] Wu | Numerical simulation of convective heat transfer between air flow and ceramic foams to optimise volumetric solar air receiver performances[END_REF], the optimisation of combined photovoltaic-thermal collectors integrated into buildings [START_REF] Muresan | Numerical simulation of a vertical solar collector integrated in a building frame: Radiation and turbulent natural convection coupling[END_REF] or the study of the heat island effect in urban environments [START_REF] Caliot | Model of spectral and directional radiative transfer in complex urban canopies with participating atmospheres[END_REF]. Certainly, developments in High-Performance Computing (HPC) over the past 20 years have contributed to overcoming some of the limitations set out. Still, the simulation of multiphysics problems involving highly turbulent flows and/or large-scale spatial domains is not accessible to reference numerical methods minimising modelling uncertainties. Also, the environmental impact of HPC questions its use by the scientific community. Physicists and engineers need to assess the appropriate level of detail for the simulation with respect to the actual needs and issues that must be responded to. Developing accurate models at low computational costs remains a key issue for research in the field of radiative transfer.

Quantities of interest Two scenarios arise. First, in the field of remote sensing or optical diagnostics, the quantity of interest is the radiative intensity itself as one wants to model the signal received by a detector. Detectors measure incoming radiation within a small solid angle (or field of view) and within several spectral ranges, as many remote sensing techniques take advantage of the spectral variation of the received signal.

Second, in the field of heat transfer, the quantities of interest are the total contributions of radiation to the energy balance of the matter. One wants to model the radiative power P rad (power per unit volume) defined as

P rad (r) = ν 4π
κ ν (r)I ν (r, Ω)dΩdν -4π ν η ν (r)dν, (1.3) and the radiative flux q rad (power per unit surface) at boundaries which writes

q rad (r b ) = π ν ε ν (r b )I • ν (T (r b ))dν - ν Ω•n<0 ε ν (r b )I ν (r b , Ω)|Ω • n|dΩdν (1.4) 1 INTRODUCTION
for an opaque wall at LTE. It is worth noting that scattering does not come into play in Eqs. (1.3)- (1.4) although the medium may scatter radiation. The accuracy and efficiency of numerical methods and models will be assessed with respect to the prediction of these two quantities.

Scope of the manuscript This professorial thesis summarises my contribution to the field of radiative transfer modelling and coupled interactions from research work performed in the period 2014-2023. It addresses part of the questions raised above related to the characterisation of the radiative properties of reacting or particulate media, the efficient and accurate simulation of radiation propagation and the understanding of coupled interactions between radiation and other heat and mass transport phenomena. The presented work was performed in various application contexts such as: the heat shield design of space vehicles, the LIDAR monitoring of soot particles, the study of pollution dispersion in urban environments, the improvement of energy efficiency and thermal comfort in buildings, the climate prediction within shallow caves.

It comprises studies carried out as a postdoctoral research assistant at the Von Karman Institute (2014)(2015), as a research fellow at Imperial College London (2015-2017) and as assistant professor at CentraleSupélec (2017-2023). Research work associated with my PhD thesis [Th1], dedicated to coupled convection-radiation interactions in differentially heated cavities, will not be discussed. In terms of computational facilities, this work used the ARCHER UK National Supercomputing Service, the HPC resources of IDRIS attributed by GENCI (Grand Equipement National de Calcul Intensif), and the HPC resources from the Mésocentre computing centre of CentraleSupélec and École Normale Supérieure Paris-Saclay.

The manuscript is organised in three parts:

• Section 2 gathers research work on radiative property modelling of plasma flow and clouds of soot fractal aggregates.

• Section 3 discusses numerical methods designed to address the angular dependence of the radiation field efficiently.

• Section 4 is dedicated to the understanding and modelling of coupled interactions between thermal radiation and turbulent natural convection flows.

The conclusion summarises the main findings and draws perspectives for orienting future research in the field. A presentation of my career history, scientific responsibilities and publications I have co-authored is given in Appendix A.

Radiative property modelling

Absorption and emission of radiation by a molecule come from energy transitions between the internal energy levels corresponding to electronic, vibrational and rotational states of the molecule. This results in line spectra that can be computed from spectroscopic database gathering line parameters (position, intensity, width) for each energy transition of the considered chemical species. Line spectra of typical radiating molecules in gaseous mixtures are made of million of peaky lines, which need a very fine spectral discretisation to describe.

Accounting for radiative properties with such high spectral resolution is referred to as the Line By Line (LBL) approach, but the computational cost of this approach is prohibitive for most radiative transfer simulations, especially in 3D spatial domains. Approximate models of radiative properties aim at coarsening the spectral discretisation at the scale of a narrow band (spectral interval where the Planck function can be assumed to be constant) or at the scale of the whole spectral range of interest. These approaches are referred to as band models and global (of full spectrum) models, respectively. Among band models, two main strategies arise in the literature: the Statistical Narrow Band (SNB) [START_REF] Malkmus | Random lorentz band model with exponential-tailed S -1 line intensity distribution function[END_REF] model, which provides an expression of averaged transmissivities from statistical assumptions on line properties, and the correlated-k (CK) [START_REF] Goody | Atmospheric radiation: theoritical basis[END_REF] model, which is based on a reordering of the spectrum in terms of absorption coefficient values to facilitate spectral integrations. Global models are usually an extension of the CK model over the whole spectral range of interest. Generally speaking, global models are more efficient than band models in reducing the spectral dimension but are also less accurate [START_REF] Taine | Gas IR radiative properties: from spectroscopic data to approximate models[END_REF]. Inaccuracies come from the dependence of radiative properties on temperature and pressure, which is not monotonous with frequency.

In particulate media, dispersed particles in gases like liquid droplets in clouds or solid particles in smoke emit and absorb radiation. They also significantly scatter radiation, contrary to molecules whose scattering effects are usually negligible in the field of heat transfer. The interaction between an incident radiation and an isolated particle can be computed from the electromagnetic theory [START_REF] Kahnert | Numerical methods in electromagnetic scattering theory[END_REF]. Depending on the shape of the particle, one may rely on the Mie therory [START_REF] Bohren | Absorption and scattering of light by small particles[END_REF] for spherical particles, the Superposition T-Matrix (STM) [START_REF] Mackowski | A multiple sphere t-matrix fortran code for use on parallel computer clusters[END_REF] method for aggregates made of non-overlapping touching spheres or the Discrete Dipole Approximation (DDA) [191,[START_REF] Enguehard | Mie Theory and the Discrete Dipole Approximation[END_REF] for particles of arbitrary shape to determine the absorption and scattering cross sections. But the challenge in modelling radiative properties of particulate media is to determine the collective effect of an ensemble of particles given the uncertainties in the optical index, size and shape of each particle, especially in the case of soot fractal aggregates.

This section contributes to the radiative property modelling of gases, plasma and particulate media. Section 2.1 briefly introduces the global model ADF for radiative properties of nearly uniform gases. A band model strategy is presented in Sec. 2.2 to model the radiation of hypersonic plasma flows in disequilibrium conditions. Finally, Sec. 2.3 discusses the backscattering modelling of soot fractal aggregate ensembles with application to the simulation of LIDAR signals.

ADF modelling of nearly uniform media

This section presents the global ADF model for radiative properties of gaseous mixtures in the case of weak temperature gradients and/or weak composition gradients such that radiative properties can be assumed to be uniform in space. The medium is assumed to be non scattering (σ ν = 0) and at local thermal equilibrium, which means that the emission coefficient η ν is related to the absorption coefficient κ ν such that η ν = κ ν I • ν (T ). There is no specific novelty in the following developments but the ADF model will be extensively used in Secs. 3 and 4 and is presented for the sake of clarity.

The Absorption Distribution Function (ADF) [START_REF] Pierrot | Accuracy of narrow-band and global models for radiative transfer in H 2 O, CO 2 , and H 2 O-CO 2 mixtures at high temperature[END_REF] model consists in substituting the integration of the radiation fields over the wavenumber with an integration over the values k of the absorption coefficient, for which a coarse discretisation is sufficient. To this aim, the following absorption distribution function [START_REF] Denison | An absorption-line blackbody distribution function for efficient calculation of total gas radiative transfer[END_REF] is introduced

F(k) = π σ SB T 4 0 ν, κν (T 0 )≤k I • ν (T 0 )dν. (2.1)
This function reorders the absorption spectrum according to the values k using the Planck function at a reference temperature T 0 as a weight to account for its variations accross the whole spectral interval. This function is discretised in intervals [k - i ; k + i ] of averaged value k i . The weights of this distribution associated with each interval i are defined as

w i = F(k + i )-F(k - i ) ( i w i = 1
). The RTE (1.1) and associated boundary condition (1.2) then become

Ω • ∇I i (r, Ω) = k i w i σ SB T 4 (r) π -I i (r, Ω) , (2.2)
I i (r b ) = ε(r b )w i σ SB T 4 (r b ) π - 1 -ε(r b ) π Ω ′ •n<0 I i (r b , Ω ′ )|Ω ′ •n |dΩ ′ , (2.3)
where σ SB is the Stefan-Boltzmann constant. The total intensity integrated over the wavenumber is simply retrieved by summing the contribution of each k-class: I = I ν dν = i I i . The radiative power and the radiative flux at boundary points are retrieved according to

P rad (r) = i 4π k i I i (r, Ω)dΩ - i 4w i k i σ SB T 4 (r), (2.4)
q rad (r b ) = ϵ(r b )σ SB T 4 (r b ) -ϵ(r b ) i Ω•n<0 w i I i (r b , Ω)|Ω • n|dΩ. (2.5)
It should be noted that the model is restricted to gray emissivities of the opaque walls at the boundaries of the spatial domain. Figure 2.1 shows the absorption spectrum of a typical air/H 2 O/CO 2 mixture at room temperature and the associated function F. This high resolution spectrum has been computed with the HITRAN 2016 database [START_REF] Gordon | The HITRAN2016 molecular spectroscopic database[END_REF]. The distribution function is monotonous and smooth and is easily captured with a coarse logarithmic discretisation in k. Here, 16 k-classes are used, evenly spaced in log-scale except the first class associated with the lowest k values (optically thin limit). As soon as the absorption coefficient is uniform, the ADF model is exact and accuracy only depends on the discretisation in k. Accuracy in computing gas emissivities with the discretised function F in Fig. 2.1 compared to the LBL approach is better than 1 % for a wide range of uniform optical path lengths at T 0 . However, in the case of heterogeneous media when temperature, pressure and composition gradients become significant, additional assumptions are required to apply the ADF model [START_REF] Denison | The spectral-line-based weighted-sum-of-gray-gases model in nonisothermal nonhomogeneous media[END_REF][START_REF] Pierrot | A fictitious-gas-based absorption distribution function global model for radiative transfer in hot gases[END_REF][START_REF] Modest | The full-spectrum correlated-k distribution for thermal radiation from molecular gas-particulate mixtures[END_REF].

HSNB modelling of reacting plasma flows

References [A6], [P4, P5, P6], [T1] Radiative property modelling in hypersonic entry flows is essential for designing spacecraft heat shields and the success of planetary missions. Spacecrafts undergo severe convective and radiative heating during atmospheric entry and radiative heating usually dominates early in the entry when the velocity is high and the density of the atmosphere is low. Incident radiation comes from the shock layer ahead of the vehicle, where the matter is potentially in disequilibrium conditions with dissociation of molecules and ionisation.

The LBL description of the emission and absorption properties of such plasma flow depends on populations of internal energy levels of atoms and molecules and on fundamental spectroscopic data gathered in spectral databases such as NEQAIR [START_REF] Whiting | NEQAIR96, nonequilibrium and equilibrium radiative transport and spectra program: User manual[END_REF], SPRADIAN [START_REF] Fujita | SPRADIAN, structural package for radiation analysis: theory and application[END_REF], MONSTER [START_REF] Kuznetsova | Absorption cross sections of diatomic molecules for problems of radiative heat transfer in low-temperature plasma[END_REF], SPECAIR [START_REF] Laux | Radiation and nonequilibrium collisinal-radiative models[END_REF]. In the present study, we use the HTGR (High Temperature Gas Radiation) database [START_REF] Chauveau | Contributions of diatomic molecular electronic systems to heated air radiation[END_REF][START_REF] Babou | Spectroscopic data for the prediction of radiative transfer in CO 2 -N 2 plasmas[END_REF][START_REF] Soufiani | High temperature gas radiation (HTGR) database and models[END_REF], developed for O 2 -N 2 and CO 2 -N 2 plasma applications. This database gathers up-to-date atomic spectroscopic data from various sources (such as NIST [START_REF] Kramida | NIST Atomic Spectra Database (version 5.2)[END_REF] and TOPbase [START_REF] Cunto | TOPbase at the CDS[END_REF]) together with ab-initio calculations of diatomic molecular spectra and atomic line shapes. It includes bound-bound atomic and molecular transitions, bound-free transitions resulting from various mechanisms, and free-free transitions. The covered spectral range is [1000 -200,000 cm -1 ], and the targeted maximum temperature is 30,000 K. However, although the LBL method is very accurate, the large number of radiative transitions that have to be taken into account makes it very computationally expensive and impractical for coupled simulations in complex geometries.

The Smeared-Rotational-Band (SRB) model or the multiband opacity-binning are common ways to simplify the calculation of molecular radiation but their accuracy is restricted to small optical thicknesses [START_REF] Johnston | Spectrum modeling for air shock-layer radiation at lunar-return conditions[END_REF][START_REF] Johnston | Extension of multiband opacity-binning to molecular, non-boltzmann shock layer radiation[END_REF]. More sophisticated approaches for radiative property modelling include the k-distribution methods over the whole spectrum or spectral narrow bands. They have been used for instance to model visible, utraviolet and vacuum ultraviolet radiation of astrophysical [START_REF] Mihalas | Foundations of radiation hydrodynamics[END_REF] or thermal [145] plasmas. In the framework of hypersonic nonequilibrium flows, the Full Spectrum Correlated-k approach (FSCK) has been efficiently implemented to model radiation of air mixtures [START_REF] Bansal | Multi-scale k-distribution model for gas mixtures in hypersonic nonequilibrium flows[END_REF] and carbonaceous atmospheres [START_REF] Bansal | Modeling of radiative heat transfer in carbonaceous atmospheres using k-distribution models[END_REF]. Successful comparisons with LBL results demonstrated the accuracy of this approach. But implementing the method becomes tedious when considering a large number of overlapping, non-weak molecular electronic systems. Moreover, the spectral information is completely lost when using such full-spectrum approaches.

This section investigates a novel approach to modelling radiative properties of hypersonic plasma flows following the work of Lamet et al. [START_REF] Lamet | Narrow-band model for nonequilibrium air plasma radiation[END_REF] and based on the Statistical Narrow Band (SNB) model. Band parameters have been computed using the HTGR database and tabulated against translational-rotational and vibrational temperatures. The model can easily include new radiating species and electronic systems, using the uncorrelation assumption inside narrow bands. Earth entry conditions corresponding to the FIRE 2 experiment, as well as Titan entry conditions corresponding to the Huygens probe, are considered to assess the accuracy and the efficiency of the model.

HSNB model for radiative transfer

The spectral radiative intensity at an arbitrary point s of an optical path starting at point s = 0 is given by the solution of the RTE for a non-scattering medium, such that

I ν (s) = I ν (0)τ ν (0, s) + s 0 η ν (s ′ )τ ν (s ′ , s)ds ′ , (2.6)
where τ ν (s ′ , s) = exps s ′ κ ν (s ′′ )ds ′′ is the spectral transmissivity between points s ′ and s. We search for an expression of the averaged intensity I ν (s)

∆ν over a spectral narrow band ∆ν. First, the radiative mechanisms are grouped into different contributions: e.g. a molecular electronic system, a set of atomic lines, or a continuum process. These contributions are assumed to be statistically uncorrelated, which allows us to write

I ν (s) ∆ν = I ν (0) ∆ν k τ k ν (0, s) ∆ν + k s 0 η k (s ′ )τ k ν (s ′ , s) ∆ν k ′ ̸ =k τ k ′ ν (s ′ , s) ∆ν ds ′ , (2.7)
where the index k refers to a radiative contribution. Note that we also assume the mean intensity at the starting point of the path (s = 0) to be uncorrelated with the total transmissivity. The evaluation of the term η k ν (s ′ )τ k ν (s ′ , s) ∆ν , in Eq. (2.7), depends on the type of radiative contribution k. The procedures used for each contribution are detailed below. When all contribution types are included, the resulting method is named the Hybrid Statistical Narrow Band (HSNB) model.

Narrow band model for optically thick molecular systems For an optically thick molecular system, the emission coefficient and transmissivity are strongly correlated. However, the ratio η ν /κ ν can be considered uncorrelated to the transmissivity τ ν , such that

η k ν (s ′ )τ k ν (s ′ , s) ∆ν = η k ν (s ′ ) κ k ν (s ′ ) ∂τ k ν (s ′ , s) ∂s ′ ∆ν ≃ η k ν (s ′ ) κ k ν (s ′ ) ∆ν ∂τ k ν (s ′ , s) ∂s ′ ∆ν . (2.8)
This assumption is valid at thermal equilibrium where the ratio η ν /κ ν is equal to the Planck function which is nearly constant within a narrow band. For thermal nonequilibrium conditions, Lamet et al. [START_REF] Lamet | Narrow-band model for nonequilibrium air plasma radiation[END_REF] checked that this assumption remains satisfactory for atmospheric entry flow applications. Discretising the optical path into homogeneous cells of size ∆s i = s i+1 -s i , the contribution of optically thick molecular systems to the mean intensity is thus written as

I thick ν (s j ) ∆ν = k∈T j-1 i=0 τ k ν (s i+1 , s j ) ∆ν -τ k ν (s i , s j ) ∆ν η k ν κ k ν i ∆ν k ′ ∈S k ′ ̸ =k τ k ′ ν (s * i , s j ) ∆ν , (2.9)
where T is the set of optically thick molecular systems, S is the set of all the systems and s * i is a mean equivalent point between s i and s i+1 introduced to simplify the spatial integration.

From statistical assumptions concerning the intensity and the position of the lines within a narrow band ∆ν, the SNB model [START_REF] Malkmus | Random lorentz band model with exponential-tailed S -1 line intensity distribution function[END_REF] provides an expression for the mean transmissivity of a homogeneous column of the form

τ ν (l) ∆ν = 1 ∆ν ∆ν exp (-κ ν l) dν = exp - W δ , (2.10)
where δ is the mean spacing between the line positions within ∆ν and W is the mean black equivalent width of these lines. The mean black equivalent width can be set as a function of three band parameters: the mean absorption coefficient κ ν ∆ν of the absorbing species and two overlapping parameters β D ∆ν and β L ∆ν related to Doppler and Lorentz broadening. For addressing non homogeneous optical paths, the Lindquist-Simmons approximation [START_REF] Young | Nonisothermal band model theory[END_REF] has been used.

Box model for optically thin molecular systems and continua If a molecular system k is optically thin for all wavenumbers ν, the mean transmissivity over a spectral narrow band can be simply expressed by τ k ν (s ′ , s) ∆ν = exps s ′ κ k ν (s ′′ ) ∆ν ds ′′ . In addition, the correlation between the emission coefficient and the transmissivity is weak such that one can write η ν τ ν ∆ν ≃ η ν ∆ν τ ν ∆ν . These two simplifications also hold for continua because of their weak spectral dynamics.

Therefore, the contribution to the mean intensity of the set B of optically thin molecular systems and continua, is then written as

I box ν (s j ) ∆ν = k∈B j-1 i=0 η k ν i ∆ν k ′ ∈S τ k ν (s * i , s j ) ∆ν ∆s i . (2.11)
The required band parameters for the box model are thus the spectral average of the emission coefficient η ν ∆ν and the absorption coefficient κ ν ∆ν . For bound-free processes, both contributions of true absorption κ abs ν ∆ν and induced emission κ ie,eq ν ∆ν to the absorption coefficient are tabulated separetly to account for potential chemical disequilibrium. The spectral average of the absorption coefficient is retrived following

κ k ν ∆ν = κ k,abs ν ∆ν -κ k,ie,eq ν ∆ν χ neq , (2.12)
where χ neq is a disequilibrium coefficient and is a function of the ratio between equilibrium and actual composition of chemical species involved the bound-free process.

Line-by-line treatment of atomic lines For atmospheric entry applications, many atomic radiative transitions are optically thick. Attempts have been made to derive a SNB model for atoms but the results were not sufficiently accurate due to the weak spectral density of atomic lines [START_REF] Lamet | Narrow-band model for nonequilibrium air plasma radiation[END_REF]. Therefore, the contribution of atomic lines to the mean intensity is treated in a LBL manner,

I at ν (s j ) ∆ν = j-1 i=0 η at ν κ at ν i (τ at ν (s i+1 , s j ) -τ at ν (s i , s j )) ∆ν × k ′ ∈S k ′ ̸ =at τ k ′ ν (s * i , s j ) ∆ν , (2.13)
where the first average is computed exactly over the narrow band. High resolution atomic radiative properties are computed according to

η ν = ul n u A ul 4π f se ul (ν -ν ul ), (2.14)
κ ν = ul n l B lu f a ul (ν -ν ul ) -n u B ul f ie ul ((ν -ν ul )) hν, (2.15)
where A ul , B ul and B lu are the Einstein coefficients related to spontaneous emission, induced emission and absorption of the transition u → l, n u and n l are the number densities of the upper and lower levels and ν ul is the wavenumber of the transition, f se ul , f a ul and f ie ul are the line profiles associated to spontaneous emission, absorption and induced emission, respectively. The line shapes are related to one another to retrieve equilibrium at T ve . A Voigt profile is considered for the spontaneous emission line shape f se ul . Einstein coefficients are taken from the NIST database and collisional broadening data are taken from Ref. [START_REF] Chauveau | Contributions of diatomic molecular electronic systems to heated air radiation[END_REF].

The LBL treatment of atomic lines is not too penalizing because of the small number of atomic lines (of the order of few thousand) compared to the number of molecular lines (of the order of several million). Furthermore, the spectral grid dedicated to atomic radiation can be much smaller as compared to the LBL spectral grid including all radiative contributions. For this work, an adaptive spectral grid which combines the 11 point stencil per line proposed in Ref. [START_REF] Da | An adaptive line-by-line statistical model for fast and accurate spectral simulations in low-pressure plasmas[END_REF] and a refinement procedure between two lines in order to accurately capture the far wing regions has been implemented.

Conclusion

The molecular systems and continua included in the database are listed in Fig. 2.2. Band parameters for thick and thin molecular systems have been tabulated according to two temperatures T (translational-rotational) and T ve (vibrational-electronic). It should be noted that the parameters can be converted for treating arbitrary electronic level populations as each electronic molecular system is treated independently. The criterion retained to decide whether a molecular system is thick or thin is based on the maximum value of the optical depth for a plasma at thermodynamic equilibrium with T =8000 K, p = 2 atm and l = 10 cm. If the maximum value of κ ν l is greater than 0.1, the molecular system is considered as thick. Band parameters for continua have been tabulated according to one temperature (T ve for free-free processes, photoionization and photodetachment, and T for O 2 photodissociation). We have considered 199 spectral bands of constant size ∆ν = 1000 cm -1 in the range [1000 -200,000 cm -1 ].

Application to atmospheric entries

Stagnation-line flow modelling A 1D plasma flow model along the stagnation-line of a spherical body has been derived following the methodology of Ref. [START_REF] Klomfass | Calculation of stagnation streamline quantities in hypersonic blunt body flows[END_REF] to simulate realistic entry flowfields. The thermal state of the plasma is described according to a two-temperature model in which the translation of heavy species and rotation of molecules are assumed to follow a Boltzmann distribution at the temperature T , and the translation of electrons, vibration of molecules, and electronic excitation of heavy species are assumed to follow a Boltzmann distribution at the temperature T ve . Thermodynamic properties, transport properties, chemical reaction rates and energy transfer rates are taken from the Mutation ++ Library [START_REF] Scoggins | Muta-tion++: Multicomponent thermodynamic and transport properties for ionized gases in c++[END_REF].

Coupled flow-radiation calculations and uncoupled radiation calculations are carried out. For solving the RTE and computing the radiative source terms in the flow energy balance, we make use of the 1D tangent slab approximation. The radiative properties of the medium are assumed to vary only along the stagnation-line direction and are assumed to be constant in planes that extend to infinity, perpendicular to this direction. It is a reasonable approximation as we will consider shock layers of which the size is small compared to the radius of the vehicle.

R 0 (m) T w (K) T ∞ (K) u ∞ (km/s) ρ ∞ (kg/m 3 ) y ∞ (

Comparison between LBL and HSNB results

In order to assess the accuracy and the efficiency of the HSNB model, a comparison with the LBL method is carried out. In LBL calculations, radiative properties of the plasma are computed from the spectroscopic HTGR database [START_REF] Soufiani | High temperature gas radiation (HTGR) database and models[END_REF] on a high resolution spectral grid of 4.4 × 10 6 points in order to capture correctly all the atomic and molecular lines. The HSNB model is also compared with the SRB model, which is often used as a simple model to treat molecular radiation but may lead to a strong overestimation of radiative fluxes. For this purpose, we implemented a model similar to the SRB that will be called hereafter HSNB-Weak. It consists in computing the mean transmissivity of thick molecular systems according to

τ k ν (s ′ , s) ∆ν = exp - s s ′ κ k ν (s ′′ )
∆ν ds ′′ instead of Eq. (2.10).

Two hypersonic entry conditions have been studied: an Earth entry condition corresponding to the trajectory point t = 1634 s of the FIRE 2 experiment and one condition of Titan entry corresponding to the trajectory point t = 191 s of the Huygens probe entry. The FIRE 2 experiment has been the subject of many numerical studies because of the availability of flight data [START_REF] Cauchon | Radiative heating results from the FIRE II flight experiment at a reentry velocity of 11.4 km/s[END_REF] and the early trajectory point t = 1634 s is associated with strong thermal nonequilibrium effects. The Titan test case has been chosen to show strong molecular radiation effects coming from the CN-Violet system [START_REF] Caillault | Radiative heating predictions for Huygens entry[END_REF]. Boundary conditions for the numerical simulations are given in Tab. 2.1. The wall of the vehicle is assumed to be non ablative, non catalytic and isothermal at T = T ve = T w and a no slip condition for the velocity is prescribed. In the free stream, temperature, velocity and mass densities are imposed. For radiation, we assume that the boundaries of the computational domain are black walls at T w and T ∞ . Finally, chemical species and radiative systems taken into account are listed in Tab. 2.2.

Calculations have been performed with the LBL, HSNB-Weak and HSNB models from the same flowfield corresponding to the coupled result obtained with the HSNB model. Table 2.3 gives the incident radiative flux at the wall, together with the total computational time for one radiation calculation, for the different combinations of models and test cases. Compared to the reference LBL solutions, the HSNB model provides an accurate prediction of the incident radiative flux, with an error between 3 and 5 % and a speed up factor around 80 for the computational time. Most of the computational gain comes from the calculation of LBL molecular spectra which is very expensive due to the large number of molecular lines. The HSNB-Weak model provides reasonably accurate results for Earth entry with a difference of 3.5 %. However, for the Titan entry case, the incident radiative flux is over-predicted by 26 %. Concerning the computational times, the HSNB-Weak model is much faster than the HSNB model. The Lindquist-Simmons approximation used for computing mean transmissivities over non homogeneous paths for thick molecular systems is responsible for the lower computational efficiency of the HSNB model.

The spectral and cumulated incident radiative fluxes at the wall are displayed in Fig. 2.3 for the Earth test case. It can be seen that the complex structure of the LBL spectral flux is correctly captured by both HSNB and HSNB-Weak models, with a good agreement on the total cumulated flux. The incoming radiation mostly arises from molecular and atomic transitions in the vacuum ultraviolet. The accuracy of the HSNB model should also be assessed regarding the total radiative power along the stagnation line which is the radiative source term in the total energy balance. Figure 2.3 shows this distribution together with the difference with LBL calculations normalized by the absolute maximum value of the total radiative power. The differences do not exceed 5 % for both HSNB and HSNB-Weak models. The highest discrepancies are located near the shock position, where the radiation emission is at a maximum.

The spectral and cumulated incident radiative fluxes at the wall and the total radiative power along the stagnation line are shown for the Titan test case on Figure 2.4. While the HSNB model reproduces with a good accuracy the LBL calculation, both the spectral flux and the total radiative power are strongly over-predicted by the HSNB-Weak model. This failure comes from an incorrect treatment of the CN-Violet molecular system in the spectral range [25,000 -29,000 cm -1 ].

Comparison with literature results

The accuracy of the HSNB model has been further assessed in comparison with the LBL radiation simulation of Bansal and Modest [START_REF] Bansal | Modeling of radiative heat transfer in carbonaceous atmospheres using k-distribution models[END_REF] of the trajectory point t = 189 s of the Huygens probe entry.

In order to reproduce their results we have computed the radiative source term along the stagnation line with both the HSNB model and the LBL approach from uncoupled flowfield taken from Johnston [START_REF] Johnston | Nonequilibrium Shock-Layer Radiative Heating for Earth and Titan Entry[END_REF]. For the CN electronic states, both Boltzmann populations and nonequilibrium populations based on the QSS model of Bose et al. [START_REF] Bose | Modeling and experimental assessment of CN radiation behind a strong shock wave[END_REF] have been taken into account. The QSS model uses simple correlations to estimate non-Boltzmann population of electronic levels from total number densities. The comparison is shown in Fig. 2.5 for the CN red and violet systems. First of all, we can see that for a given population assumption, our LBL and HSNB results give similar results. When the QSS model is used, we obtain a good agreement with the LBL results of Bansal and Modest [START_REF] Bansal | Modeling of radiative heat transfer in carbonaceous atmospheres using k-distribution models[END_REF] for both LBL and HSNB. Note that nonequilibrium effects are not negligible in the considered simulation and lead to about 16 % difference at the peak value of the radiative source term.

Entries of ablative materials In addition to an efficient modelling of thick molecular systems, the interest in the HSNB framework lies in the ability to easily add new species or radiative contributions. The formulation of the model in terms of transmissivity takes advantage of the lack of correlation between the different radiative contributions.

Carbon-phenolic ablative materials are likely to dominate future space missions because of the lightness of these materials and the cooling and protective effect of the ablation process. The ablation products absorb part of the radiation from the shock layer and carry this energy (189 s). Comparison between LBL, HSNB and LBL results from Bansal and Modest [START_REF] Bansal | Modeling of radiative heat transfer in carbonaceous atmospheres using k-distribution models[END_REF]. Contribution of the CN red and violet systems to the radiative source term along the stagnation line. downstream of the stagnation region. The relevant carbonaceous species for ablation studies have been included in the HTGR database and in the HSNB model. The HSNB model has been used to analyse coupled ablation-radiation phenomena for the peak heating trajectory point of the Apollo 4 command module entry [P6].

Another application of the HSNB model is the study of meteor entries. Monitoring radiation emission during meteor entry provides useful information on the mass and the composition of these objects. Meteors are composed of metallic species such as Iron, Magnesium, Potassium, Silicon or Sodium, which oxidize during the entry. Band parameters for the relevant metallic oxides have been computed and the HSNB model has been successfully used in luminosity calculations of meteor entry [START_REF] Dias | Luminosity calculation of meteor entry based on detailed flow simulations in the continuum regime[END_REF].

Backscattering modelling of soot fractal aggregate ensembles

References [A10], [C3], [T2] Soot particles are fine matter aerosols formed during the incomplete combustion of organic materials from natural or anthropogenic processes. They are primarily composed of carbon, therefore having a strong ability to absorb light in a wide spectrum. The molecular composition and submicronic size of these particles are problematic for human health. In addition, soot particles are found to induce a positive radiative forcing (warming effect) in the energy budget of the Earth, through both direct and indirect effects (increase in cloud lifetime, role as cloud condensation nuclei). The characterisation and monitoring of soot particles come up against the complexity of their morphology [START_REF] Okyay | Soot aggregate complex morphology: 3D geometry reconstruction by SEM tomography applied on soot issued from propane combustion[END_REF]. Freshly emitted soot particles are composed of clustered spherical primary particles, called monomers. The structures resulting from the aggregation of primary particles are complex, and are often referred to as fractal-like shape [START_REF] Köylü | Fractal and projected structure properties of soot aggregates[END_REF]. The composition of soot particles is also subject to uncertainties, as well as their dynamic evolution during their atmospheric lifetime.

Radiative property modelling of soot particles is of primary importance for studying the climate, for modelling fire propagation [START_REF] Parent | Spectral radiation emitted by kerosene pool fires[END_REF] or for optical diagnostics aiming at monitoring and characterising them [START_REF] Bouvier | A novel approach for in-situ soot size distribution measurement based on spectrally resolved light scattering[END_REF]. To study soot particles in atmospheric conditions, remote sensing techniques are especially of interest. Light Detection And Ranging (LIDAR) instruments rely on collecting and detecting the light backscattered by the particles and molecules inside the volume of an emitted laser pulse to study this interacting medium. Detecting this backscattered light and analysing the associated signal can provide information on the particles, such as their concentration, size distribution or morphological properties [START_REF] Müller | Microphysical particle parameters from extinction and backscatter lidar data by inversion with regularization: experiment[END_REF]. But using LIDAR signals requires an apriori knowledge of soot radiative properties, especially their backward scattering cross section. The Rayleigh-Debye-Gans theory for Fractal Aggregates (RDG-FA) [START_REF] Sorensen | Light scattering by fractal aggregates: A review[END_REF] provides good estimation of integrated cross sections but can lead to strong inaccuracies in estimating angle-dependent properties. The analysis of LIDAR signals is further complicated by multiple scattering effects, which can hardly be taken into account in inversion methods.

This work aims to provide a priori information on soot radiative properties and multiple scattering effects in LIDAR inversion methods. First, a sensitivity study of LIDAR-relevant radiative properties to morphological parameters has been carried out. Then, a Monte Carlo algorithm has been developed to accurately simulate LIDAR signals associated with clouds of soot particles in order to evaluate the potential impact of multiple scattering.

LIDAR framework

We consider spectrally-resolved and polarisation resolved LIDAR instruments which can record the co-polarised (∥) and cross-polarised (⊥) components of the backscattered light relative to the emitted beam polarisation plane in the wavelength range 300 to 1100 nm. The backscattered power P λ (z) ∥/⊥ recorded by the instrument, at a given range z from the receiver and wavelength λ, can be modelled by the LIDAR equation

P ∥/⊥ λ (z) = 1 z 2 K ∥/⊥ λ (z)O λ (z)U ∥/⊥ λ (z), (2.16)
where the superscripts refer to the state of polarisation of the received light relative to the emitted pulse, K is the LIDAR instrument function and O is the range-dependent overlap function. U is the attenuated backscatter function defined as

U ∥/⊥ λ (z) = σ ∥/⊥ λ,back (z) exp -2ξ ∥/⊥ λ (z) z 0 β λ (z ′ )dz ′ , (2.17)
where σ λ,back is the backscattering coefficient, β λ = κ λ + σ λ is the extinction coefficient and ξ λ is a correction factor introduced by Platt [START_REF] Platt | Remote sounding of high clouds. iii: Monte Carlo calculations of multiple scattered lidar returns[END_REF] to account for multiple-scattering. The single scattering approximation (ξ λ = 1) supposes that the electromagnetic wave scattered by a scatterer depends solely on the incoming laser electromagnetic wave. Hence, the contribution of the scattered wave by the other scatterers on the incident radiation is neglected. This neglected contribution is usually called multiple scattering, and can induce an increase in the received radiation by the LIDAR instrument (0 ≤ ξ λ ≤ 1). LIDAR inversion techniques require an a priori knowledge of the backward scattering coefficient (for both co-polarised and cross-polarised components) and of the extinction coefficient. These quantities are also referred to as the LIDAR products. In the following, we will only consider absorption and scattering by the particles and neglect the molecular contribution from the hosting medium (air). We will also consider the simple case where the medium contains a monodisperse distribution of aggregates meaning that all particles are described by the same morphological parameters. The LIDAR products are then expressed as

σ ∥/⊥ λ,back = n p C ∥/⊥ λ,back
(2.18)

β λ = n p C λ,ext (2.19)
where n p is the number density of the particles, C ∥/⊥ λ,back is the backscattering cross section and C λ,ext is the extinction cross section of the soot particles.

In addition, two quantities are of specific interest: the LIDAR ration (LR) and the LIDAR depolarisation ratio (LDR). The LR is a LIDAR parameter which is often used in LIDAR signal inversion methods. It quantifies the ratio of the extinguished light relatively to the backscattered light

LR λ = C λ,ext C λ,back . (2.20)
For the LDR, we will distinguish the particle LDR denoted δ p and the volume LDR denoted δ v . The particle LDR expresses the change of polarisation of the scattered wave by a particle and is defined as

δ p λ = C ⊥ λ,back C ∥ λ,back . (2.21)
It is a discriminating parameter to classify the scattering particles as δ p strongly depends on the morphology of the scatterer (δ p = 0 for spherical scatterers). The volume LDR expresses the change of polarisation of the backscattered power recorded by the LIDAR instrument and is defined as

δ v λ = P ⊥ λ (z) P ∥ λ (z) . ( 2

.22)

A non zero volume LDR results from a non zero particle LDR and/or from multiple scattering effects. For instance, an assembly of spherical particles may change the state of polarisation of the backscattered light through multiple scattering effects as shown experimentally in Ref. [START_REF] Carswell | Polarization anisotropy in LIDAR multiple scattering from atmospheric clouds[END_REF].

LIDAR-relevant radiative properties

Methodology The objective here is to assess the sensitivity of LIDAR-relevant radiative properties to the morphological parameters of the Soot Fractal Aggregates (SFA). To this end, we rely on a morphological fractal model with monodisperse monomers, where each spherical primary particle has the same size. This model relates the number of monomers N m in a particle with its overall size following the power-law

N m = k f R g r m D f , (2.23)
where R g is the radius of gyration, r m the radius of the monomers, D f the fractal dimension and k f the fractal prefactor. The radius of gyration is defined from the spatial position x i of each monomer with respect to the center of mass of the aggregate: R g = ( 1 Nm Nm 1 x 2 i ) 0.5 . The fractal aggregates have been generated using the tunable cluster-cluster aggregation algorithm developed by Mackowski [START_REF] Mackowski | A simplified model to predict the effects of aggregation on the absorption properties of soot particles[END_REF]. The algorithm generates monomer positions from given parameters N m , r m , k f and D f such that the radius of gyration of the generated aggregates satisfies Eq. (2.23) in a statistical sense. Twenty-seven ensembles of one hundred SFA have been generated, with the fractal prefactor always set to k f = 1.3 and the number of monomers to N m = 45, 125 or 450. The fractal dimension has been set to either D f = 1.6, 1.8 or 2.0 and the monomer radius to r m = 10 nm, 20 nm or 40 nm. These values are representative of those reported in the literature [START_REF] Köylü | Structure of overfire soot in buoyant turbulent diffusion flames at long residence times[END_REF].

The radiative properties of SFA are computed using the MSTM code [START_REF] Mackowski | A multiple sphere t-matrix fortran code for use on parallel computer clusters[END_REF] based on the Superposition T-Matrix method for aggregates made of non-overlapping touching spheres, which is the case here. The STM method relies on the expansion of the electric field on vector spherical wave functions and uses the superposition property of the scattered fields to account for the contribution of one sphere scattered field on another sphere. The radiative properties are computed over a large spectrum (300 nm ≤ λ ≤ 1100 nm, with a wavelength step λ = 20 nm). The soot optical index is modelled with the wavelength-dependent dispersion law from Chang and Charalampopoulos [START_REF] Chang | Determination of the wavelength dependence of refractive indices of flame soot[END_REF]. Results For each of the twenty-seven SFA ensembles, the mean and standard deviation of the lidar products (C ext , C back ), the LR (Eq. (2.20)) and the particle LDR (Eq. (2.21)) have been computed. Figure 2.6 shows the extinction and backscattering cross section of the different SFA ensembles of fractal dimension D f = 1.8. Cross-sections are normalised by the corresponding volume of the aggregates 4 3 πr 3 m N m in order to highlight the specific effect of the morphology. Both cross sections decrease with the wavelength. For a given wavelength, the normalised extinction cross section is little affected by the fractal parameters which means extinction is mostly driven by the volume of the scatterer. Still, non-negligible variations with respect to the monomer radius can be noticed (from 5 % to up to 20 %). The standard deviations of the extinction cross sections across all wavelengths and ensembles are below 0.1 %. However, the normalised backscattering cross section shows significant variations with monomer's numbers and radii: increasing twofold the monomer radius, leads to an increase of the normalised backscattering cross section by a factor from 1.5 to 6. The associated standard deviation is also more important, about 10 % over the whole spectrum for each ensemble of aggregates. This indicates that the link between backscattering and morphology can not be fully described by the volume of matter of an aggregate.

As the LR is the ratio between the extinction cross section and the backscattering cross section, both morphological and spectral dependence of the LR are closely related to those of the cross sections. Indeed, Fig. 2.7 shows an increasing LR with larger wavelengths, which is consistent with the wavelength dependence of the cross sections. Similarly, a larger monomer radius reduces the LR: LR values with different monomer radii show close values in the UV part of the spectrum, and diverge in the near infrared. Higher number of monomers also induces higher LR, with the noticeable exception of the ensemble of smallest aggregates (r m = 10 nm, N m = 45). The variation of the LR with SFA's number of monomers is still small, the associated standard deviations overlapping on many parts of the spectrum. Impacts of the fractal dimension are harder to evaluate, as several distinctive trends occur as a function of both wavelength and monomer radius.

Figure 2.8 shows a decrease of the particle LDR and its associated standard deviation with increasing wavelength, which means LIDAR signals will be more sensitive to the non-sphericity of SFA at small wavelengths. The particle LDR significantly increases with the monomer radius: for example, the LDR increases about tenfold when the monomer radius varies from r m = 10 nm 2.4: Radiative properties of the aggregates considered in the Monte Carlo simulation at λ = 355 nm using optical index of Ref. [START_REF] Chang | Determination of the wavelength dependence of refractive indices of flame soot[END_REF]. Optical index m, extinction and backscattering cross sections, albedo ω, LR and particle LDR.

to 40 nm at λ = 300 nm. This increase is proporionately smaller at high wavelengths. The particle LDR also increases with the number of monomers, although with a much less noticeable effect than that of the monomer radius. An impact of the fractal dimension is noticeable at small wavelength and large monomer radius.

To conclude, a strong influence of the monomer radius on all LIDAR-relevant properties has been observed, which means the processing of LIDAR signals might provide means to evaluate this parameter. However, a weak influence of the number of monomers and the fractal dimension has been detected. Still, these results do not exclude that very high (D f ≃ 3) or very low (D f ≃ 1) fractal dimensions, might have larger impacts. Even averaged over one hundred aggregates, the standard deviations of the backscattering cross section (and so the standard deviations of the LR and particle LDR) of ensembles of same morphological parameters are important. This indicates that the backscattering cross section is sensitive to the fine structure of the aggregates and that the chosen morphological description is not complete.

Monte Carlo modelling of multiple scattering

Problem description An accurate modelling of the backscatter function U ν (z) ∥/⊥ accounting for multiple scattering effects requires the solve of the RTE within the full spatial domain of interest instead of using the approximate Eq. (2.17). We rely here on the Monte Carlo method which is a stochastic algorithm to solve the RTE.

The propagating medium is modelled as a succession of homogeneous plane slabs, with definite thickness along the vertical axis z and infinite width along the horizontal axes. The emitter and receiver are modelled as circles of radius R em = 0.1 m and R rec = 0.045 m. They are positionned on the groung (z = 0) at the same central point. The Field Of View (FOV) of the emitter is set to θ em = 0.1 mrad while the FOV of the receiver is taken at either θ rec = 0.5 mrad or 5 mrad. Three wavelengths have been considered, namely λ = 355 nm, 532 nm and 1064 nm which correspond to Nd:YAG laser frequency tripled, doubled and main emission wavelengths respectively. For the sake of conciseness, we will only discuss the results obtained with the shortest wavelength as it is associated with the strongest multiple scattering effects. In what follows, the wavelength subscript will be omitted for all spectral quantities.

The medium is assumed to be composed of a cloud of aggregates with polydispersed monomers. Because of the cost of Monte Carlo simulation, the analysis is restricted to a single set of morphological parameters, corresponding to the mean values addressed in Sec. 2.3.2, namely D f = 1.8, k f = 1.3 and N m = 100. The mean monomer radius is set to r m = 20 nm and its standard deviation is set to σ rm = 1.1. A set of 100 aggregates is generated using the FracVAL code [START_REF] Morán | Fracval: An improved tunable algorithm of cluster-cluster aggregation for generation of fractal structures formed by polydisperse primary particles[END_REF] and the radiative properties are computed using the MSTM code [START_REF] Mackowski | A multiple sphere t-matrix fortran code for use on parallel computer clusters[END_REF]. The mean properties of this ensemble are reported in Tab. 2.4 using the optical index given in Ref. [START_REF] Chang | Determination of the wavelength dependence of refractive indices of flame soot[END_REF]. It should be noted that other optical index values from other literature sources have been considered to assess the sensitivity of simulated LIDAR signals to the optical index (see Ref. [T2]). Finally, a Gaussian particle concentration profile of these aggregates is considered and expressed as

n p (z) = n 0 exp - 1 2 z -z c σ z 2 (2.24)
where n p (z) is the particle number concentration as a function of the height from the receiver z, n 0 is the particle number concentration at the profile maximum, z c = 1.15 km is the profile center and σ z = 0.0375 km is the standard deviation. The particle concentration n 0 is varied to address different total cloud optical thicknesses ζ tot defined as

ζ tot = zc+4σz zc-4σz n p (z ′ )C ext dz ′ = n 0 C ext σ z √ 2π. (2.25)
Monte Carlo method The Monte Carlo method consists in simulating the emission of a large number of energy bundles (called photons hereafter for simplicity) to reproduce the laser pulse and track their propagation within the medium from a stochastic description of absorption and scattering phenomena. Counting the number of photons which travels back to the receiver allows us to estimate the backscatter function. The Stokes vector is transported here instead of the radiative intensity in order to model the polarisation state.

In this application, the standard Monte Carlo method would suffer from a very low convergence rate because the probability that an emitted photon is backscattered is typically of the order of 0.01 ‰. This implies that only one photon in one hundred thousand will actually contribute to the output signal. In order to optimise the method, we make use of three variance reduction techniques: (1) the absorption weighting which prevents a stochastic treatment of absorption; (2) the peel-off scattering [START_REF] Peest | Polarisation in Monte Carlo radiative transfer and dust scattering polarisation signatures of spiral galaxies[END_REF] which forces scattering events towards the receiver; (3) the scattering splitting which creates two photons at each scattering events, one in the forward direction and one in the backward direction. The final Monte Carlo algorithm is made of the following steps:

• initiation step A primary photon is initiated. The initial position is randomly drawn within the emitter disc and the initial propagation direction is randomly drawn within the emitter FOV. The Stokes vector is initiated with linear horizontal polarisation.

• transport step A propagation length is randomly drawn according to the extinction properties of the medium. The photon Stokes vector is weighted by the albedo to only account for scattering events.

• peel-off step It consists in forcing a scattering event towards the receiver, provided that the photon position is within the receiver FOV. The photon Stokes vector is appropriately modified according to the scattering event probability and to the extinction between the initial photon position and the peeled-off photon position on the receiver.

• scattering step It consists in sampling the scattering direction according to the current state of the photon and the medium properties. Scattering splitting can occur if the current photon is within the scattering splitting volume and if its current direction vector is directed outward this volume. Two photons are then created, one whose scattering angle is sampled within the forward hemisphere, and the second in the backward hemisphere. The Stokes vector of the scattered photons is calculated consistently to account for polarisation change induced by scattering and is appropriately weighted according to the scattering event probability given by the phase function.

Transport, peel-off and scattering steps are repeated for each photon issued from the primary photon until a stopping condition is triggered, either if a maximum scattering order is reached or if the Stokes vector becomes too small. This is referred to as the cascade loop as a large number of photons are being tracked from a single primary photon because of scattering splitting. Then, another primary photon is generated whithin the main loop until a large number of primary photons (here 10 9 ) is launched. The implementation of the present method has been validated against the literature results of Ref. [START_REF] Hu | Simple relation between lidar multiple scattering and depolarization for water clouds[END_REF]. In addition, the correction factor ξ(z) is calculated from

ζ tot = 10 ζ tot = 1
ξ(z) = 1 - 1 2 z 0 n p (z ′ )C ext dz ′ log P MS (z) P SS (z) , (2.26)
and the associated corrected signal is calculated from Eq. (2.17

MSF(z) = P MS (z) -P SS (z) P MS (z) . ( 2

.27)

Both are shown in Fig. 2.10. For a reliable results analysis, the only reported data points are those for which the computed relative errors on the MSF and on the LDR are below 10 % and those for which the received signal P MS (z) is greater than 1 % of the corresponding signal maximum. Overall, both MSF and volume LDR increases with the optical thickness and with the receiver FOV. Indeed, a larger optical thickness increases the probability of scattering events and a larger FOV increases the probability of multiple-scattering contributions to the return signal. The MSF and volume LDR reaches a maximum value of 50 % and 4 %, respectively. In all simulation cases with low optical thickness, both MSF and LDR decrease after reaching a maximum. This behaviour can be interpreted by an insufficient number of scattering events in the section of the cloud located after the maximum of particle number concentration. The MSF A HSNB model has been derived to model the radiative properties of atmospheric entry flows, in disequilibrium conditions. Compared to the LBL approach, the model is accurate and computationally efficient in order to address 3D geometries. The model is advantageous for taking into account thick molecular contributions and can easily include new radiating species and electronic systems, using the uncorrelation assumption inside narrow bands. However, solving the RTE in 3D spatial domains with the HSNB model would require deterministic or stochastic ray-tracing algorithms because of the model's formulation in terms of transmissivity. Also, a LBL description of atomic radiation is needed because of the weak spectral density of atomic lines.

A methodology has been proposed to simulate synthetic LIDAR signals from clouds of soot particles. It is based on a fractal description of the soot morphology, an accurate computation of the radiative properties of the aggregates with the MSTM method and a Monte Carlo algorithm to account for multiple scattering effects within the cloud. A sensitivity analysis of LIDAR-relevant radiative properties has shown that the monomer radius could be evaluated from LIDAR signal processing. Phenomenological relationships have been derived for the MSF and volume LDR to provide a priori information that can be used in LIDAR inversion algorithms. Nevertheless, the results were limited to a monodispersed distribution of SFA where each soot particle is described with the same morphological parameters.

Radiative transfer simulation

Challenges of the numerical resolution of the Radiative Transfer Equation (RTE) have been highlighted in introduction. Among them, a specific problem concerns the variation of the radiative intensity with respect to the propagation direction. This feature is found similarly in neutron transport and, in general, in any Boltzmann transport equation where particle densities vary in a position/velocity phase space.

The angular dependence of the radiation field is known to vary a lot with the optical thickness, which can be seen as the ratio between a macroscopic characteristic length and a mean free path of the photons resulting from interactions with the matter at a given wavelength. The optically thick regime (very large optical thicknesses) is relatively simple to model as the mean free path of the photons is so small that radiation transport becomes a local phenomenon. The Rosseland model reduces radiation transport to a conductive process, which combines with actual conduction in the energy balance of the matter. The truncated decomposition of the intensity on spherical harmonics (P N [START_REF] Modest | Radiative heat transfer[END_REF] or Simplified P N [START_REF] Larsen | Simplified pn approximations to the equations of radiative heat transfer and applications[END_REF] expansions) is also a well-suited method to describe optically thick media except near the boundaries. The optically thin regime (very low optical thicknesses) is harder to capture as the mean free path of the photons is large. Based on a discrete representation of angular variations, the Discrete Ordinates Method (DOM or S N expansion) is a popular technique to capture the optically thin regime but suffer from ray effects when the angular resolution is not sufficient [START_REF] Coelho | Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media, review[END_REF].

Actual absorbing gases are characterised by line spectrum with variations of the absorption coefficient (and thus of the optical thickness) of several orders of magnitude in the wavelength range of interest. Therefore, numerical methods need to capture both optically thin and thick regimes. Reference methods such as ray-tracing algorithms, either in their deterministic form or stochastic form (Monte Carlo), achieve very high accuracy whatever the optical thickness but require intensive computational resources. An original approach for decreasing computational requirements relies on multi-scale splitting. Roger et al. [START_REF] Roger | A hybrid transport-diffusion model for radiative transfer in absorbing and scattering media[END_REF] and Coelho et al. [START_REF] Coelho | Multi-scale methods for the solution of the radiative transfer equation[END_REF] developed a hybrid transport-diffusion model where the radiative intensity is split into a macroscopic part calculated by a diffusion equation and a mesoscopic part calculated by a Monte Carlo or a discrete ordinates method. Additionally, Luo et al. [START_REF] Luo | Multiscale solutions of radiative heat transfer by the discrete unified gas kinetic scheme[END_REF] developed a space-angle discretisation scheme, referred to as the discrete unified gas kinetic scheme (DUGKS), providing accurate and reliable numerical solutions for large, small, and in-between values of optical thickness.

In this section, we investigate two strategies to tackle the complexity associated with the angular dependence of the radiation field within a wide range of optical thicknesses and to reduce the number of angular degrees of freedom of the discretised RTE: (1) a goal-based angular adaptivity algorithm in Sec. 3.2 and (2) an angular reduced-order model in Sec. 3.3. But beforehand, the finite-element framework from which are built these two methods is described in Sec. 3.1.

Finite element framework

In this section is described the Boltzmann transport solver FETCH, based on the finite element method for both space and angular variables, and used as a framework for exploring novel angular discretisations. Emitting, absorbing but non-scattering medium are considered. The ADF model (see Sec. 2.1) is employed to account for the spectral dependence of the absorption coefficient.

Angular discretisation Following the finite element method, the angular dependence of the radiative intensity associated with the k-class of index i, is restricted to the function space spanned by a set of M i angular basis functions {G iq (Ω)} for q ∈ {1, 2,

• • • , M i } I i (Ω, r) ≃ M i q=1 G iq (Ω)I iq (r), (3.1)
I iq (r) being the coefficients of the projection onto the basis. It should be noted that the set of angular basis functions as well as the size of the basis can be different for each k-class of index i. The size of the basis may also vary in space M i (r) when angular adaptivity is performed (Sec. 3.2). Introducing the expansion (3.1) into Eq. (2.2) and applying Galerkin projection leads to the following angular linear system in Cartesian coordinates

A i,x ∂I i (r) ∂x + A i,y ∂I i (r) ∂y + A i,z ∂I i (r) ∂z + H i I i (r) = S i (r), (3.2)
where I i (r) is a vector of size M i containing the projection coefficients I iq (r). A i, x/y/z and H i are the angular matrices defined by

(A i, x/y/z ) p,q = Ω x/y/z G ip (Ω)G iq (Ω)dΩ, (3.3)
(H i ) p,q = k i G ip (Ω)G iq (Ω)dΩ, (3.4)
and S i (r) is the angular source term vector whose components are

S iq (r) = k i w i σ SB T 4 (r) π G iq (Ω)dΩ. (3.5)
This general framework holds for various angular discretisations such as discrete ordinates (S N ), spherical harmonics or wavelets. In the following, the considered discretisations will be derived from S N expansions. A finite element formulation of the S N method is implemented here which is not strictly equivalent to the standard quadrature-based S N method. It consists in dividing the angular domain into nonoverlapping solid angles ∆Ω p of equal size and associating to each angular interval a piece-wise constant angular basis function such that G p (Ω) = 1 if Ω ∈ ∆Ω p and 0 otherwise. Without scattering, it should be noted that the angular matrices associated to the discrete ordinates discretisation are diagonal.

Spatial discretisation A variational multiscale finite element method [START_REF] Buchan | The inner-element subgrid scale finite element method for the Boltzmann transport equation[END_REF] is used to perform the spatial discretisation. This approach combines a continuous finite element representation with a discontinuous subgrid scale (SGS) model. The spatial variations of the angular coefficients I i (r) are decomposed into a coarse component I i (r) and a subgrid component

∼ I i (r).
The coarse component is assumed to lie in a continuous finite element space, spanned by η N basis functions, while the subgrid component is assumed to lie in a discontinuous space, spanned by η Q basis functions, such that

I i (r) ≃ η N j=1 N j (r)I ij + η Q j=1 Q j (r) ∼ I ij . (3.6)
The continuous basis functions N j (r) and the discontinuous basis functions Q j (r) are associated with the same spatial mesh and are both piecewise linear across the spatial elements in this work.

The same spatial mesh is used for each k-class.

Using decomposition of Eq. (3.6) into Eq. (3.2) and applying Galerkin projection in space leads to a space-angle block linear system of size M i × (η N + η Q ) that has to be solved for each k-class i. The role of the subgrid component is to suppress any instability and nonphysical oscillations that may arise from a continuous finite element discretisation of the RTE.

However, the method enables the subgrid scale component of the solution to be eliminated from the resulting system of discretised equations so that a single linear system of size M i × η N is solved. The main advantage of the method is to provide accurate and numerically stable solutions, similar to those obtained from a discontinuous Galerkin formulation, but with much lower number of spatial degrees of freedom. More details on the implementation, the accuracy and the efficiency of the subgrid scale model can be found in Ref. [START_REF] Buchan | The inner-element subgrid scale finite element method for the Boltzmann transport equation[END_REF].

Boundary conditions are implemented using a Riemann decomposition method in order to separate the incoming and outgoing part of the radiative flux at the boundaries. This is a general technique that applies to any type of angular discretisation [START_REF] Buchan | Riemann boundary conditions for the Boltzmann transport equation using arbitrary angular approximations[END_REF]. The same set of angular basis functions is used for both coarse and subgrid scales. The space-angle linear system is solved using the generalised minimal residual method (GMRES), implemented with the PETSc library [START_REF] Balay | PETSc users manual[END_REF]. Either a spatial multigrid preconditioner [START_REF] Buchan | A sub-grid scale finite element agglomeration multigrid method with application to the Boltzmann transport equation[END_REF] or a successive over-relaxation (SOR) preconditioner are applied. Both GMRES and multigrid methods are implemented in a matrixfree fashion as the whole space-angle system is often too large to store in memory.

Goal-based angular adaptivity

References [A7], [C2]

A promising way to reduce the dimension of discretised space-angle systems in radiation transport is the use of adaptive discretisation methods that only refine the regions of interest. Several spatial adaptivity schemes for radiation transport have been developed during the past decade. An Adaptive Mesh Refinement (AMR) technique, based on a hierarchy of structured spatial meshes, was used by Ogando and Velarde [START_REF] Ogando | Development of a radiation transport fluid dynamic code under AMR scheme[END_REF], combined with a S N expansion in angle. Ragusa [START_REF] Ragusa | A simple Hessian-based 3D mesh adaptation technique with applications to the multigroup diffusion equations[END_REF] proposed a regular error measure, using the Hessian of the discretised spatial flux, while Lathouwers [START_REF] Lathouwers | Goal-oriented spatial adaptivity for the S N equations on unstructured triangular meshes[END_REF] and Goffin et al. [START_REF] Goffin | Minimising the error in eigenvalue calculations involving the Boltzmann transport equation using goal-based adaptivity on unstructured meshes[END_REF] experimented with goal-based error measures with various target functionals. More recently, other strategies based on h-refinement techniques [START_REF] Yang | h-Refinement for simple corner balance scheme of S N transport equation on distorted meshes[END_REF] or multilevel spatial meshes [START_REF] Liu | Multilevel mesh adaptivity for discrete ordinates transport calculation with spatial-moment-ratio indicators[END_REF] have been derived. Although all these space adaptivity techniques are successful in reducing the number of spatial degrees of freedom, their efficient implementation in parallel remains challenging.

A first attempt at adapting the angular resolution was proposed by Ackroyd and Wilson [START_REF] Ackroyd | Composite finite element solutions for neutron transport[END_REF] who used a variable order spherical harmonic expansion (P N ) across space. Using a P N expansion is advantageous for angular adaptivity because it forms a hierarchical basis and it does not need any angular interpolation in space. Goffin et al. took advantage of this property to develop goal-based angular adaptivity algorithms using a P N expansion [START_REF] Goffin | Goal-based angular adaptivity applied to the spherical harmonics discretisation of the neutral particle transport equation[END_REF]. Alternative hierarchical basis for angular discretisation can be formed using wavelet theory and multi-resolution analysis. Compared to spherical harmonics, wavelets have the advantage of adapting anisotropically due to their compact support. Watson [START_REF] Watson | The W N adaptive method for numerical solution of particle transport problems[END_REF] and Goffin et al. [START_REF] Goffin | Goal-based angular adaptivity applied to a wavelet-based discretisation of the neutral particle transport equation[END_REF] both made use of wavelet-based angular adaptivity methods with Haar wavelets and octahedral linear wavelets respectively. Other adaptive schemes have also been developed for non-hierarchical basis, such as the S N expansion [START_REF] Stone | Adaptive discrete-ordinates algorithms and strategies[END_REF]. Kópházy and Lathouwers [START_REF] Kópházi | A space-angle DGFEM approach for the Boltzmann radiation transport equation with local angular refinement[END_REF] proposed a generalised framework for local angular refinement with arbitrary angular basis functions and discontinuous Galerkin discretisation in space and angle, but the method involves complex algebra to compute the numerical fluxes. Favennec et al. [START_REF] Favennec | Ad hoc angular discretization of the radiative transfer equation[END_REF] allow unstructured adaptivity on their angular domain, however they use the same adapted angular discretisation across their entire spatial domain.

The objective of the present work is to develop an angular adaptivity method for thermal radiation, in the framework of coupled heat transfer problems, where the angular resolution is optimised to accurately compute the radiative power. A Haar wavelet angular expansion is used to perform anisotropic angular adaptivity. The angular resolution is allowed to vary in space and for the different absorption coefficient classes associated with the radiative properties of the medium.

Adaptivity algorithm

The goal-based adaptivity algorithm is build upon two main features: the use of a hierarchical angular basis and the definition of an error measure on the targeted goal.

Haar wavelets

The angular adaptivity algorithm relies on the use of a hierarchical functional basis which means that the set of basis functions at an order n is a sub-set of basis functions at an order n + 1. This allows us to use different sets of angular basis functions of different orders across space without the need to perform angular interpolation. The Haar wavelet expansion is a hierarchical version of the discrete ordinates (S N ) expansion. The Haar wavelets are piecewise constant functions built from a multiresolution analysis and associated with a partition scheme of the angular space [START_REF] Adam | Adaptive Haar wavelets for the angular discretisation of spectral wave models[END_REF]. An exact linear mapping relates the Haar wavelet and the discrete ordinates expansions built from the same angular partition. Thus, a transformation matrix can be applied to angular matrices in order to map from one discretisation to the other. This simplifies the implementation of Haar wavelets, allowing them to be built on existing S N -type frameworks. A 1D representation of the Haar wavelet basis functions over an arbitrary interval is given in Fig. 3.1 with the corresponding S N basis functions and the transformation matrix. The angular matrices are less sparse with a Haar discretisation than with a S N discretisation, making them more expensive to apply, but this is balanced by the hierarchical nature of the Haar wavelets with which angular adaptivity is simple. The Haar wavelet space will be noted W m,n and will be associated to the following 2D partition of the angular space: in each octant, the values of the cosine of the polar angle µ are divided into 2 m-1 equally-spaced intervals and the values of the azimuthal angle φ are divided into 2 n-1 equally-spaced intervals, leading to a total number of basis function of 8 × 2 m-1 × 2 n-1 in three dimensions of space. The 2D Haar wavelets are obtained from a tensor product of the 1D Haar wavelets in the polar and azimuthal directions.

Goal-based error measure We define the functional goal as the radiative power integrated in space F = r P rad (r)dr, as this is the physical quantity we want to control the accuracy in the framework of coupled heat transfer problems. Building error measures for this functional goal, requires solving the RTE defined by Eq. (1.1) (further referred to as the forward problem) as well as solving the adjoint RTE defined by

-Ω • ∇I * ν (r, Ω) = S * -κ ν I * ν (r, Ω), (3.7)
where I * ν (r, Ω) is the solution of the adjoint problem and S * = ∂f (I ν )/∂I ν = κ ν is the adjoint source term involving the targeted goal, with f being the integrand of the functional F .

The functional error measure could be expressed as the inner product in the phase space either (1) between the residual of the forward computation and the error in the adjoint computation -this will be referred to as the adjoint error measure -or (2) between the residual of the adjoint computation and the error in the forward computation -this will be referred to as the forward error measure. A local and discrete form of these adjoint and forward goal-based error measures are formally written as

e * ijq = | ϵ * ijq R ijq | ∆ F , e ijq = | ϵ ijq R * ijq | ∆ F , (3.8)
where ϵ and ϵ * are the forward and adjoint solution errors, R and R * are the forward and adjoint residuals and ∆ F is a normalising factor related to a user-defined tolerance for the functional goal. The indices i, j and q refer to the wavenumber, space and angular degrees of freedom and

vary for i = {1, • • • , N k }, j = {1, • • • , η N } and q ∈ M ij . Adjoint e *
ijq and forward e ijq error estimators are then combined to get a unique goal-based error measure e GB i,j,q such that e GB i,j,q = max{ e ijq , e * ijq }.

(3.9)

Implementation The angular adaptivity algorithm is an iterative process described in Algorithm 1. Starting from a uniform first order discretisation (W 1,1 ), the forward and adjoint problems are solved and the goal-based error measure e GB i,j,q is computed for each degree of freedom in the phase space. This error measure has been scaled such that it should be less than 1.0 if the tolerance on the functional has been met. If it is greater than 1.0, the angular resolution is locally increased by adding the basis functions of the next level of resolution that neighbour the basis function q under consideration. However, if the error measure is much smaller than 1.0, here smaller than an arbitrary threshold 0.01, the angular resolution is coarsened by removing the basis function q. The process is repeated until the angular discretisation has converged. In practice, we restrict the algorithm up to a maximum order M for the Haar wavelet expansion so that the adapted resolution does not significantly evolve after M iteration steps.

Algorithm 1 Goal-based angular adaptivity algorithm using M adaptive steps.

Set initial angular discretisation for m = 1, • • • , M do Solve forward (Eq. (1.1)) and adjoint (Eq. (3.7)) problems Calculate goal-based error measure (Eq. (3.9))

for i = 1, • • • , N k do for j = 1, • • • , η N do for q ∈ M ij do if e GB i,
j,q < 0.01 then Remove angular basis function q from node j else if e GB i,j,q > 1.0 then Add next level angular basis function to node j end if end for end for end for end for

In order to evaluate the relevance of the assigned goal, results of the goal-based adaptivity method will be compared with results of the regular adaptivity method in the next section. A regular adaptivity means that the error measure does not rely on a specific target but on the value of the radiative intensity itself. In that case, algorithm 1 is used but without solving the adjoint problem and replacing the goal-based error measure of Eq. (3.9) by the following

e REG ijq = | ϵ ijq |/(∆I max i ), (3.10)
where ϵ ijq is the foward angular error and ∆ is a user-defined tolerance relative to the maximum intensity I max i of the k-class i. This tolerance value has been selected such that the regular adaptive algorithm retains around the same number of angular basis functions as the goal-based, allowing for simple comparisons.

Accuracy and efficiency of the method

Street canyon test case In order to assess the accuracy and efficiency of the goal-based algorithm, we consider a 2D coupled flow/radiation problem representative of mixed convection in street canyons. The geometrical configuration and the boundary conditions of the problem are displayed in the left part of Fig. 3.2. The canyon has a reference length L and contains a humid air with a homogeneous water vapour molar fraction of 2 %. Its walls are rigid, black and isothermal at T 0 for the upper and lateral walls and T 0 + ∆T for the ground (T 0 = 294.2 K, ∆T = 1 K). An isothermal inlet flow at T 0 and u ref is imposed on the left open boundary while a zero conductive flux and and an outlet flow is imposed at the right open boundary. These two lateral open boundaries are assumed to be specular reflectors. On the top boundary, a zero normal flow and a zero conductive flux are imposed, as well as an incoming radiation flux corresponding to the flux emitted and transmitted by the atmosphere in the infrared (computed for a mid-lattitude summer atmosphere using the LBLRTM-LW radiation model [START_REF] Clough | Atmospheric radiative transfer modeling: a summary of the AER codes, short communication[END_REF]).

This coupled flow-radiation problem is solved for Rayleigh and Reynolds numbers of Ra = 10 8 and Re = 5 × 10 3 . The finite element fluid dynamics code Fluidity [START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF] is used to solve the flow equations. The radiation field is solved with a W 4,4 uniform expansion and with 12 kclasses to model the non-grey absorption of the humid air. For both flow and radiation, the spatial domain is discretised with an unstructured mesh made of 22,000 nodes, and a spatial domain decomposition was performed (48 partitions) in order to solve in parallel. A temperature snapshot is extracted from this coupled simulation once the transient features are removed and is shown in Fig. 3.2. The goal-based and regular adaptive algorithm are then applied to a single radiative transfer calculation using this coupled temperature snapshot Comparing uniform and goal-based adapt calculations, we can observe that the goal-based adaptive algorithm keeps almost all the angular basis functions up to the second order and then discriminates the higher order angular basis functions without compromising the accuracy: the W GB 4,4 resolution uses five times less angular basis functions than the W 4,4 resolution to reach the same accuracy in the functional (around 0.2 %). The goal-based adaptive algorithm is also quicker than the unifrom discretisation from the order 4: the W GB 4,4 resolution requires 5 times less CPU time than the W 4,4 resolution to reach the same accuracy in the functional, and even 6.5 times less CPU time if intermediate adaptive steps are not taken into account. Not including the time of the intermediate adaptive steps allows us to assess the actual computational cost of one radiation calculation with the adapted resolution which will be used in coupled calculations. If we now consider the regular adaptive calculation, Fig. 3.3 shows that the regular adaptive algorithm retains around the same number of angular basis functions than the goal-based one but with higher functional errors, so that there is no gain in accuracy for around the same computational cost compared with a uniform discretisation. It confirms that the goal-based error measure is more appropriate than the regular one to optimise the angular resolution according to the radiative power.

Distribution of adapted angles

The purpose here is to analyse which angular basis functions are retained by the goal-based adaptivity algorithm W GB 4,4 depending on the spatial position and the k-class. Figure 3.4 shows the number of angular basis functions, averaged in space, for each optical thickness class k i L. The distribution has a bell shape: the number of angular basis functions increases with the absorption coefficient up to a maximum reached for kL ≃ 1 and then decreases. The shape of this distribution is caused by two factors. First, in the radiative source term, the radiative intensity is weighted by the absorption coefficient which means that the higher the absorption coefficient is, the higher the contribution of the class to the source term is. This feature is captured by the source of the adjoint problem which is equal to the absorption coefficient. Secondly, the angular dependence of the radiative intensity vary with optical thickness, and we know that in thick media, the radiative transfer becomes local and the radiative intensity is close to isotropic. We can thus infer that a low number of angular basis functions is sufficient to accurately represent the radiative intensity for the thick classes. This explanation is confirmed by the spatial distribution of the number of angular basis functions, shown in Fig. 3.5 for three different classes. For the thickest class (k 10 L = 21), a small number of angular basis functions is used throughout the domain except near the nonreflecting boundaries, where the radiative intensity is more likely to deviate from isotropic. Although the number of angular basis functions is also quite low in average for the thinnest class (k 3 L = 2.2 × 10 -3 ), its spatial distribution strongly differs. Marks of ray effect can be observed that suggests that the radiative intensity is highly directional within the domain and probably under-resolved. Finally, the number of angular basis functions is the highest for the intermediate class (k 6 L = 0.11), and its spatial distribution is rather homogeneous except near the boundaries.

Coupled calculations So far, we demonstrated that the goal-based algorithm is useful in efficiently computing the radiative source term from a temperature snapshot. Numerical tests showed that the adapted resolution can be used in coupled calculations: the accuracy of the calculations is maintained over time compared to the uniform discretisation. Also, it was noticed that is was not necessary to change this adapted resolution over time by running the goal-based algorithm after a couple of time steps. Compared to keeping the initial adapted resolution, no significant gain in accuracy was obtained. Indeed, the spatial distribution of the adapted angles (Fig. 3.5) does not seem to be correlated with the temperature field and was thus not substantially varying over time. It seems the angular dependence of the radiation field is dominated by the boundary conditions, rather than determined by the local temperature variations in the gas.

Angular reduced-order modelling

References [A8] Reduced Order Models (ROM) are widely used in computational physics for decreasing the number of degrees of freedom of complex dynamical systems. They aim at finding a low order optimal basis able to represent a physical phenomenon or to condensate the response of a physical model to a set of parameters [START_REF] Benner | A survey of projection-based model reduction methods for parametric dynamical systems[END_REF]. Most of ROMs rely on the Proper Orthogonal Decomposition (POD) method which is based on a statistical analysis of high-fidelity numerical data or experimental data. This method extracts from the data a few dominant uncorrelated modes, that are able to capture the underlying physics. POD has been widely used in fluid mechanics, where spatial basis functions are extracted from high resolution spatial fields (snapshots) recorded through time [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF][START_REF] Sirovich | Turbulence and the dynamic of coherent structures. Part I: Coherent structures[END_REF], to model turbulent pipe flow [START_REF] Bakewell | Viscous sublayer and adjacent wall region in turbulent pipe flow[END_REF], natural convection [START_REF] Podvin | Low-order models for the flow in a differentially heated cavity[END_REF] or pollutant transport in street canyons [START_REF] Fang | Reduced order modelling of an unstructured mesh air pollution model and application in 2d/3d urban street canyons[END_REF].

In the field of radiative transfer, such model reduction techniques have been applied by Pinnau and Schultze [START_REF] Pinnau | Model reduction techniques for frequency averaging in radiative heat transfer[END_REF] for treating the spectral dependence of the radiative intensity. Buchan et al. [START_REF] Buchan | A POD reduced order model for resolving angular direction in neutron/photon transport problems[END_REF] proposed a POD reduced order model for resolving the angular dimension in neutron transport problems, where angular POD basis functions were derived from high order P N solutions. Tencer et al. [START_REF] Tencer | Reduced order modelling applied to the discrete ordinates method for radiation heat transfer in participating media[END_REF][START_REF] Tencer | Accelerated solution of discrete ordinates approximation to the Boltzmann transport equation for a gray absorbing-emitting medium via model reduction[END_REF] followed a similar approach to build optimal angular basis functions from high order discrete ordinates expansions, but only considered grey absorbing media. Another strategy has been followed by Behne et al. [START_REF] Behne | Minimally-invasive parametric model-order reduction for sweep-based radiation transport[END_REF], who introduced global POD modes in space and angle that are parametrised depending of absorption properties or source terms. In addition, other modal decompositions such as the Dynamic Mode Decomposition (DMD) have been recently employed to derive a ROM for radiation transport [START_REF] Sharak | A fast method based on Dynamic Mode Decomposition for radiative heat transfer in participating media[END_REF].

The goal of this study is to investigate the suitability of an angular ROM for thermal radiation involving actual participating gases and thus covering a large optical thickness range, associated with various angular behaviours of the radiation field. An angular ROM is derived from high order discrete ordinates expansions, based on the POD method, for each absorption coefficient k-class independently. The angular POD basis sets are truncated at different levels depending on the k-class in order to optimally compute the total radiative power.

POD-based angular discretisation

The Proper Orthogonal Decomposition (POD) method applied to the angular discretisation of the RTE aims at constructing an optimal angular basis based on a set of data representative of the angular variations of the radiative intensity. The method extracts from the data a hierarchical basis set that can be truncated to any size depending on the desired accuracy.

The method of snapshots is used here to build the POD angular basis functions. The method is applied independently to each k-class in order to form adapted POD basis functions to each optical thickness range. An accurate solution of the problem of interest is generated using a reference model for the angular discretisation. We use here a high order discrete ordinates method S N , with the same number of angular degrees of freedom M for each k-class. The reference angular coefficients at each spatial node j = 1, • • • , η are stored in column to form the following snapshot matrix S i of size M × η N for each k-class i

S i = [I i1 , • • • , I ij , • • • , I iη N ].
(3.11)

The POD method consists in searching for the eigenvalue decomposition of the symmetric matrix

S i S T i S i S T i = U i Λ i U T i , (3.12)
where Λ i is a diagonal matrix containing the eigenvalues of the decomposition, that are all positive and ranked in decreasing order. The matrix U i contains the eigenvectors of the decomposition: it represents a change-of-basis matrix between the S N angular basis functions and the POD angular basis functions. Therefore, the POD angular basis functions {G P iq (Ω)}, q ∈ {1, 2, • • • , M} are obtained by combining the POD eigenvectors (the columns of the matrix U i ) with the S N angular basis functions {G q (Ω)}, q ∈ {1, 2, • • • , M}, such that

G P iq (Ω) = M p=1 (U i ) p,q G p (Ω). (3.13)
The POD basis functions form a hierarchical basis set that can be truncated to any size. The aim of the reduced order model is to restrict this basis set to the first M P i basis functions (M P i << M) with the highest eigenvalues, which means they are the best for reconstructing the snapshot matrix in terms of the Frobenius norm. In the reduced order model, we will assume that the angular variations of the radiative intensity associated to the k-class i can be captured using a few numbers of POD angular basis functions

M P i I i (Ω, r) ≃ M P i q=1 G P iq (Ω)I P iq (r), (3.14)
where I P iq (r) are the angular coefficients in the POD space. The efficiency of the truncation for representing accurately the reference calculation will depend on the decay of the eigenvalue spectrum: the faster is the decay, the lower is the number of POD basis functions necessary.

Introducing the expansion (3.14) into the RTE and applying Galerkin projection leads to the angular linear system (3.2) using POD angular matrices A P i,x/y/z and H P i of size M P i × M P i (and POD angular source vector S P i (r) of size M P i ) which can be computed from S N angular matrices A x/y/z and H of size M × M (and S N angular source vector S(r) of size M) with the mapping

A P i,x/y/z = U T i A x/y/z U i , H P i = U T i HU i , S P i (r) = U T i S(r) (3.15)
Note that only the first M P i columns of the matrix U i are kept here, such that U i is of size M×M P i . This means the angular ROM can be easily implemented from the forward discretised angular system. It should be emphasized here that a unique POD angular basis set is formed per k-class by condensing the reference angular information throughout the spatial domain. Thus, the POD angular expansion does not change with the spatial position. In addition, the size of the POD truncation M P i is kept constant across space, although it can be varied with the k-class i.

Accuracy and efficiency of the method

Differentially heated cavity test case In order to assess the suitability of an angular POD model for radiative heat transfer calculations in non grey media, we consider the natural convection of a humid air mixture in a square differentially heated cavity. The cavity is made of two vertical walls which are black and isothermal and two horizontal walls which are adiabatic and perfectly diffuse reflecting. The Rayleigh number is set to Ra = 5 × 10 6 , the molar fraction of water vapour to 2 %, the mean temperature to T 0 = 294.2 K and the temperature difference between hot and cold walls to ∆T = 5 × 10 -2 K. The flow is solved, without radiation source term, with the finite element fluid dynamics code Fluidity [START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF] using an unstructured spatial grid of about 4400 elements. The resulting steady temperature field is shown in Fig. 3.6. From this temperature field, radiation calculations are carried out with nine k-classes (and one transparent class) to model the non-grey absorption of the humid air. The angular POD reduction is performed from a reference radiation calculation using a S 26 discrete ordinates method (364 angular basis functions) for the angular discretisation of each of the nine k-classes. The process detailed in Sec. 3.3.1 is applied to build the angular POD basis functions and derive the angular reduced order model.

POD spectrum and angular basis functions

The eigenvalue spectrum of the decomposition is presented in Fig. 3.7 for five of the nine k-classes. Larger eigenvalues are obtained for the thin classes because the weight w i (see Eq. (2.2)) associated with the radiation emission of the k-class i decreases with the optical thickness. Each spectrum is characterized by a dominant first eigenvalue, representing the mean angular distribution but the decay of the next eigenvalues changes with the optical thickness. For the thin classes, the decay is slower meaning that the angular information is more difficult to condensate and that a large number of POD basis function is necessary to reconstruct the data. At the optically thin limit, the propagation the same for each k-class i and that this error roughly scales with the square root of the eigenvalue spectrum of the POD. However, in the framework of coupled heat transfer, the angular POD-ROM is expected to be accurate in computing the radiative power

k 3 L = 3.88 × 10 -2 k 7 L = 14.2
P rad (r) = i (k i G i (r) -4σ SB w i k i T (r) 4
) rather than in computing the incident radiation.

In order to find a priori the best POD basis size M P i for each k-class i for optimally computing the radiative power, a simple error scaling is dervied based on the eigenvalue spectrum and the optical thickness. It consists in adding POD functions onto the angular basis until the following criteria is satisfied

λ i,m k i L < τ 2 , (3.16)
where τ is a given tolerance and λ i,m is the m th eigenvalue for k-class i, and fix the basis size M P i = m. The resulting distribution of angular basis size across the optical thickness range is presented in Fig. 3.9 for different tolerance values. This distribution has a bell shape and is similar to that obtained with the goal-based adaptivity algorithm (see Fig. 3.4): the maximum angular resolution is put for the fourth and fifth classes. A low angular resolution is put in the thinnest classes because they do not contribute much to the radiative power, although the associated eigenspectrum decay is slow. On the contrary the thickest classes contribute a lot to the radiative power but the angular information is efficiently reduced by the POD, and again only a few POD basis functions are retained. This optimised POD discretisation will be noted POD* in the following.

We will now compare the accuracy and efficiency of uniform POD (same number of basis functions for each k-class), optimised POD* (using the basis function distribution of Fig. 3.9) and uniform discrete ordinates using the following error measure e = (P rad (r) -P rad,ref (r))

2 dr P rad,ref (r) 2 dr , (3.17)
where the superscript "ref" refers to the reference S 26 solution. This error measure is given against the averaged angular basis size and the computational time in Fig. 3.10 for POD, POD* and S N models. For the optimised POD* model, the points in the graphs correspond to tolerance values in the range τ 2 ∈ [5 × 10 -5 ; 2 × 10 -6 ]. Compared to the discrete ordinates method, the uniform POD model has a faster convergence rate and uses less basis functions to reach a given error level but this angular reduction is not sufficient to obtain an actual gain in terms of computational time. However, the optimised POD* model allows for a much higher reduction of the averaged angular basis size leading to a significant time saving. Compared to the discrete ordinates method, the optimised POD* model uses about 10 times less angular basis The advantage of the angular adaptivity method is the ability to vary the angular resolution in space, although load balancing issues may arise for its efficient implementation in parallel.

The POD method provides a simple framework to build optimised angular basis functions but its predictive capabilities for unseen problems, that are outside the range of parameters used to generate the reference data, need to be clarified. However, when considering coupled interactions between thermal radiation and turbulent flows, the RTE has to be solved a large number of times for a given configuration, following the time evolution of temperature and concentration in absorbing species. It is then worthwhile to derive an optimised angular discretisation targeted for a specific problem.

Coupled convection-radiation interactions

Sections 2 and 3 introduced models tackling specific aspects of radiation transport such as its angular dependence or its spectral dependence. But radiation transport couples with the dynamic evolution of the state of the matter (energy, momentum, chermical composition), satisfying its own conservation laws.

Thermally-driven convection of radiating fluids is a typical example of a multiplysics problem involving radiation transport. Thermally-driven convection is a ubiquitous phenomenon, encountered in atmospheric physics, in astrophysics or in various engineering applications. These are flows generated by density variations caused by temperature differences. In radiating fluids, the emission and absorption of radiation affect the temperature, which in turn control the buoyant motion. Water vapour and carbon dioxide are the most common radiating gases in the infrared and are present in significant quantity in the atmosphere or in confined environments such as buildings or shallow caves.

Convective flows are usually turbulent, which means velocity and temperature fields vary chaotically over a wide range of time and length scales. Radiative transfer effects on turbulent convection have been first investigated by Spiegel [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]. Radiative dissipation of temperature fluctuations has been further confirmed by direct numerical simulations in turbulent channel flows [START_REF] Vicquelin | Effects of radiation in turbulent channel flow: analysis of coupled direct numerical simulations[END_REF]. In addition, in thermal convection, radiation also affects the mean temperature gradient and 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: for instance 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], 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]. Given that radiative transfer rather affects the large and intermediate scales of the flow, the question arises whether a low-order model could capture radiation effects on turbulent convection.

This section is dedicated to the understanding and modelling of coupled interactions between radiation and turbulent convection in confined environments. First, Direct Numerical Simulations (DNS) of coupled convection and radiation are performed and thoroughly analysed in Sec. 4.2 for a cubic Rayleigh-Bénard cell, whose size and content in absorbing species are representative of heat transfer in indoor environments. These results are then used in Sec. 4.3 to derive a low-order model able to capture radiation effects on the intermittent dynamics of convection. The last section (Sec. 4.4) investigates coupled convection-radiation interactions in shallow caves, that are humid and weakly ventilated cavities below ground surface, in order to better understand the climate within these environments. But beforehand, Sec. 4.1 presents the common modelling approaches and numerical methods used in these studies.

Numerical setup

Governing equations We consider the natural convection flow of air/CO 2 /H 2 O mixtures in enclosed cavities. The flow motion is solely caused by temperature gradients and potentially moisture gradients, which are small enough to make use of the Boussinesq approximation. Density variations with respect to temperature and moisture content are linearised around a reference state (T 0 , c w0 ) and only accounted for in the buoyant term of the momentum balance. Other thermophysics properties are assumed to be constant.

The parameters controlling the flow regime are the thermal Rayleigh number Ra and the solutal Rayleigh number Ra Φ defined by

Ra = gβ∆T L 3 ν f a , Ra Φ = gβ Φ ∆c w L 3 ν f D w , (4.1)
where g is the gravitational acceleration, β is the thermal expansion coefficient, β Φ is the solutal expansion coefficient, ∆T is a reference temperature difference, ∆c w is a reference water vapor concentration difference and L is a reference length. The flow fields are also affected by two other nondimensional parameters: the Prandtl number Pr = ν f /a and the Schmitt number Sc = ν f /D w which compare momentum diffusivity ν f with thermal diffisivity a and mass diffusivity D w of water vapour in air. Mass, momentum, energy and moisture balance are made dimensionless using the reference length L, the reference time L 2 /(a √ Ra), the reduced temperature θ = (T -T 0 )/∆T and the reduced concentration Φ = (c w -c w0 )/∆c w and write as

∇ • u = 0, (4.2)
∂u ∂t + u • ∇u = -∇p + Prθe z + Ra Φ Ra Pr 2 Sc Φe z + Pr √ Ra ∇ 2 u, (4.3)
∂θ ∂t + u • ∇θ = 1 √ Ra ∇ 2 θ + P rad , (4.4)
∂Φ ∂t + u • ∇Φ = Pr Sc √ Ra ∇ 2 Φ, (4.5)
where u and p are the dimensionless velocity and motion pressure and e z is a unit vector in the vertical direction pointing upward. Cavities of parallelepipedic shape with no slip walls are considered. Temperature or temperature gradients (as well as moisture or moisture gradients) are prescrined at cavity walls. Condensation and evaporation may occur at the boundaries but not within the volume so that there is no source term associated with phase change in transport equations (4.4)-(4.5).

The radiative transfer equation (RTE) (1.1) for a non scattering medium (σ ν = 0) has to be solved together with Eqs. (4.2)-(4.5) to compute the radiative power P rad (Eq. (1.3)) affecting the energy balance. In accordance with the Boussinesq approximation, the absorption coefficient κ ν is considered uniform. Cavity walls are assumed to be either black or perfectly diffuse reflecting. It is worth noting that it is not possible do define a limited set of nondimensional parameters to make the RTE dimensionless. This is due to the complex structure of absorption spectra of H 2 O and CO 2 in the infrared, and the resulting infinite number of optical thicknesses for a given length. The actual size of the cavity, the mean temperature level, and molar fractions of the absorbing species must be specified. For a given mixture and temperature level, the highresolution (line-by-line) absorption coefficient is computed using the HITRAN database [START_REF] Gordon | The HITRAN2016 molecular spectroscopic database[END_REF] and associated ADF model parameters (see Sec. 2.1) are determined for a practical implementation of spectral integrations.

Numerical methods Flow 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 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. Either a Direct Numerical Simulation (DNS) approach (in Secs. 4.2 and 4.3) or a Large-Eddy Simulation (LES) approach (in Sec. 4.4) is used. For the latter, the Spectral Vanishing Viscosity (SVV) method [START_REF] Tadmor | Convergence of spectral methods for nonlinear conservation laws[END_REF][START_REF] Karamanos | A spectral vanishing viscosity method for largeeddy simulations[END_REF][START_REF] Pasquetti | Spectral vanishing viscosity method for large-eddy simulation of turbulent flows[END_REF] is employed to model the effects of the unresolved subgrid scales. It consists in introducing an artificial dissipation term to ensure spectral convergence and dissipate the high modes of the Chebyshev polynomial development. Two parameters control the SVV modelling: a cut-off Chebyshev mode from which the artificial dissipation acts and a weight determining the intensity of the dissipation.

The RTE (1.1) is solved using a ray-tracing algorithm for each ADF k-class of absorption coefficient value for the associated mixture. The computation is made parallel by distributing the rays among the different processors. The 4π angular domain is uniformly discretized using Case radiation model solutal buoyancy and moisture transport C transparent gas (κ ν = 0, P rad = 0) disregarded (β c = 0, Φ = 0) CR radiating gas disregarded (β c = 0, Φ = 0) CRM radiating gas considered Table 4.1: Three considered cases with different levels of simplfication.

a large number of rays (typically 1000) from volume and boundary cell centres. At low Rayleigh numbers, the radiation mesh is built from the convection mesh, coarsened by a factor of two in each direction of space. At high Rayleigh numbers, the radiation mesh is coarsened by a factor of four/five in each direction of space compared to the convection mesh and we use a radiation subgrid model [A3] to account for the radiation of the spatial scales resolved by the flow mesh but filtered by the coarse radiation mesh. This subgrid model has been validated in various configurations and its accuracy is about a few percent on radiative power and wall fluxes [A5].

In the case of a flow LES (Sec. 4.4), the radiation of spatial scales unresolved by the flow mesh is ignored, which is justified by the uniformity of the absorption coefficient and the quasi-linearity of the Planck function in the small range of temperature considered. Flow calculations and radiation calculations are coupled explicitly and the radiative power is typically updated every 10 convection time steps. Indeed, the flow time step is imposed by numerical stability constraints and does not correspond to significant variations of the temperature field. In usteady regime, time integration is carried out over a period ∆t once the asymptotic regime (statistically steady) is reached. Validation and accuracy of numerical methods and coupled algorithm are not discussed here, but evidences can be found in my PhD thesis [Th1] and in references highlighted at the begining of each of the following subsections. 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). For the highest Rayleigh numbers considered, a coupled simulation can take up to 500,000 hours of CPU time on the scalar partition of the French supercomputer Jean-Zay.

In the following, different modelling assumptions will be used in the simulations in order to highlight specific effects associated with gas radiation and solutal buoyancy. Three cases as defined in Tab. 4.1 will be considered: the C case refers to the thermal convection of a transparent gas (Eqs. (4.2)-(4.4) with P rad = 0 and Φ = 0), the CR case refers to the thermal convection of a radiating gas (Eqs. (4.2)-(4.4) with Φ = 0 coupled with Eq. (1.1)), and the CRM case refers to the thermosolutal convection of a radiating gas (Eqs. (4.2)-(4.5) coupled with Eq. (1.1)).

Coupled heat transfer in Rayleigh-Bénard cells

References [A14] ,[P9], [P11], [T3]

Rayleigh-Bénard Convection (RBC), in which a fluid is confined between two horizontal plates, heated from the bottom and cooled from above, is a simple configuration for studying coupled radiation-convection interactions. First studies on the topic were focused on the convection onset of radiating fluids. Using linear stability analysis, Goody [START_REF] Goody | The influence of radiative transfer on cellular convection[END_REF], followed by Spiegel [START_REF] Spiegel | The convective instability of a radiating fluid layer[END_REF], have shown that radiative transfer delays the onset of convection and stabilizes the fluid layer. Two stabilizing physical mechanisms were highlighted: the decrease of the static temperature gradient in the core of the layer and the damping of thermal perturbations with radiative transfer. Though these studies were restricted to a gray fluid (radiative properties independent of the wavelength), the stabilizing effect of radiation has been further confirmed for real molecular gases [START_REF] Bdéoui | The onset of Rayleigh-Bénard instability in molecular radiating gases[END_REF] and supported by experiments conducted with ammonia [START_REF] Gille | Convection in a radiating gas[END_REF] or carbon dioxide [START_REF] Hutchinson | Effect of nongray gas radiation on thermal stability in carbon dioxide[END_REF].

While a substantial research effort has been devoted to understand turbulent RBC and predict the scaling laws for the Nusselt number [START_REF] Grossmann | Scaling in thermal convection: a unifying theory[END_REF], the large-scale organisation of the flow [START_REF] Brown | Reorientation of the large-scale circulation in turbulent 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] the turbulence properties [START_REF] Lohse | Small-scale properties of turbulent Rayleigh-Bénard convection[END_REF] or the plume dynamics [START_REF] Van Der Poel | Plume emission statistics in turbulent Rayleigh-Bénard convection[END_REF], the study of radiative transfer effects on RBC at high Rayleigh numbers has received little attention due to the computational cost associated with radiation calculations. Some authors have nevertheless attempted to account for radiative transfer in steady or weakly turbulent RBC. Lan et al. [START_REF] Lan | Stability analysis for threedimensional Rayleigh-Bénard convection with radiatively participating medium using spectral methods[END_REF] explored flow regimes above the onset of convection (Ra ∼ 10 3 ) using three-dimensional (3D) coupled calculations for a gray radiating fluid. Mishra et al. [START_REF] Mishra | Numerical analysis of Rayleigh-Bénard convection with and without volumetric radiation[END_REF] also performed two-dimensional (2D) coupled calculations for a gray participating medium for Rayleigh numbers around 10 4 and reported an increase of the number of convective cells with radiation. Radiation effects on the shape and on the number of large-scale convection rolls have also been observed in Ref. [START_REF] Sakurai | Radiation effects on mixed turbulent natural and forced convection in a horizontal channel using direct numerical simulation[END_REF] for mixed convection (Ra ∼ 10 6 ) using an optically thin model for radiation. In weakly turbulent regime (10 6 ≤ Ra ≤ 10 7 ), coupled numerical results for a real radiating gas in a a cubic cell [P3] show that radiative transfer significantly increases the kinetic energy of the mean and fluctuating flows, though the results were obtained within a limited integration time. This convection intensification with radiation is also noticed in experiments where a lightspot serves as heat source for an absorbing fluid [START_REF] Lepot | Radiative heating achieves the ultimate regime of thermal convection[END_REF][START_REF] Bouillaut | Transition to the ultimate regime in a radiatively driven convection experiment[END_REF]. Radiation is found to promote the mixing-length scaling regime as energy transfer is no longer restricted by the boundary-layers and can influence the flow field in the core of the cavity. It should be finally mentioned that surface-to-surface radiation can also affect RBC as shown for instance by [START_REF] Czarnota | Turbulent convection and thermal radiation in a cuboidal Rayleigh-Bénard cell with conductive plates[END_REF] who consider radiative and conductive horizontal plates. This brief literature review reveals a lack of reference results to understand the radiation effects on RBC and to derive and validate simple models able to capture these effects. To this aim, coupled DNS have been performed in a Rayleigh-Bénard cubic cell in the Rayleigh number range 10 6 ≤ Ra ≤ 10 9 . The cavity considered is shown in Fig. 4.1: it has horizontal isothermal walls (black) at T hot and T cold and adiabatic side walls (perfectly diffuse reflecting). Physical parameters are chosen to be representative of heat transfer in buildings: the reference temperature is taken at T 0 = 300 K, the Prandtl number is Pr = 0.707, the molar fraction of absorbing species are set to X H 2 O = 0.02 and X CO 2 = 0.02 and the size of the cavity L is set to either 1 m ot 3 m. Numerical parameters (mesh size, time step and time intervals) are given in Tab. 4.2. Water vapour gradients are disregarded and only C and CR cases will be discussed (the latter will be referred to as CR(1m) and CR(3m) depending on the size of the cavity). The analysis of the simulation results in this section will be focused on statistics. Statistical averaging combines time averaging and symmetry averaging over the 16 symmetries satisfied bv temperature and velocity fields in the cubic cell [START_REF] Puigjaner | Stability analysis of the flow in a cubical cavity heated from below[END_REF]. The temporal dynamics will be discussed in Sec. [START_REF] Babou | Spectroscopic data for the prediction of radiative transfer in CO 2 -N 2 plasmas[END_REF]

Flow fields and heat transfer

Figure 4.2 shows the temperature, the total kinetic energy and the radiative power, averaged statistically and over the horizontal plane, along the vertical z for cases C, CR(1m) and CR(3m).

In the C case, the temperature field displays thermal boundary layers near the walls, which shrink when the Rayleigh number increases, and a nearly isothermal core at θ ≃ 0. When gas radiation is taken into account (CR cases), radiation exchanges tend to break the temperatue homogeneity in the core and temperature gradients are no longer restricted to boundary layers. The kinetic energy is also significantly increased. However, gas radiation effects on temperature and velocity fields diminish when the Rayleigh number increases and when the cavity size is smaller. It can be explained by the following estimate of the radiative source term in the energy balance which roughly scales as

1 √ Ra P rad = O 1 √ Ra κ P σ SB T 3 0 L 2 λ , (4.6)
where κ P = κ ν I 0 ν (T 0 )dν × π/(σ SB T4 0 ) is the Planck mean absorption coefficient and λ is the thermal conductivity. Thus, the radiative source term varies in L 2 / √ Ra, whereas the order of magnitude of the convective source term u • ∇θ remains the same regardless of the Rayleigh number. The radiative power field shows an absorbing zone in the lower half of the cavity (z ≤ 0.5, P rad > 0) and an emitting zone in the upper half of the cavity (z ≥ 0.5, P rad < 0), except in a thin layer near the walls where this behaviour reverses. It can be checked in comparing CR(1m) and CR(3m) cases, that the magnitude of P rad varies in L 2 .

Figure 4.3 shows the three components of the total heat flux q tot q tot = -∂ θ s ∂z

q cond + √ Ra θw s q conv + L λ∆T ν 4π I ν Ω • e z dΩdν s q rad , ( 4.7)
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. The total flux q tot is constant with height as the side walls are adiabatic. In the C case, there is a balance between the conductive and convective fluxes. Variations are restricted to the boundary layer regions: in the cavity center, convection prevails and conduction vanishes.

In CR cases, an enhancement of the convective flux in the center of the cavity is observed. The variations of the convective flux are no longer restricted to the boundary layers due to a convection-radiation coupling in the quasi-isothermal part of the cavity. Though, the influence of gas radiation decreases as the Rayleigh number increases and the cavity size decreases.

Values of wall conductive and radiative fluxes are reported in Tab. 4.3. In the C case, the radiative flux represents the energy exchange between the isothermal black upper and lower

C CR(1m) CR(3m)

Ra q cond q rad (1m) q rad (3m) q cond q rad q cond q rad 8.30 walls. It is constant across the range of Rayleigh numbers but depends on the cavity size. Gas radiation tends to increase the wall conductive flux but not significantly. Comparing C and CR cases, it can also be noticed that the wall radiative flux decreases for a given cavity size because of a screening effect due to gas absorption.

Ra C CR(1m) CR(3m) terms ×10 3 ⟨τ u ⟩ v = ⟨ω θ ⟩ v ⟨τ u ⟩ v ⟨ω θ ⟩ v ⟨ω rad ⟩ v ⟨τ u ⟩ v ⟨ω θ ⟩ v ⟨ω rad ⟩ v 10

Energy budgets

We investigate in this section the effect of gas radiation on three averaged quantities: the kinetic energy e k = 1 2 u • u, the potential energy e p = -Prθ(z -0.5) and the thermal energy e θ = 1 2 θ 2 , where • denotes the statistical average. 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. 4.4. In the C case, the kinetic energy remains quite constant across the Rayleigh number range as the reference velocity

u ref = a √ Ra L
has been chosen to balance intertial and buoyant forces. Gas radiation clearly increases the kinetic energy but the relative effect of radiation is less important as the Rayleigh number increases and the cavity size decreases. Potential ernergy 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 and both are higher in CR cases because of a higher temperature gradient in the core.

The potential energy balance writes

∂e p ∂t + ∇ • (ue p ) = -Prθw - Pr √ Ra ∇ • ((z -0.
⟨ω θ ⟩ v - Pr √ Ra (z -0.5)P rad v ⟨ω rad ⟩ v = Pr θw v ⟨τu⟩ v , (4.9)
where Nu is the Nusselt number (conductive flux at the wall). Terms ⟨ω θ ⟩ v and ⟨ω rad ⟩ v correspond to sources of potential energy associated with conduction and radiation transport of temperature variations. These sources are balanced with the kinetic energy production term ⟨τ u ⟩ v associated with the work of buoyant forces. Table 4.4 reports values of these volumeaveraged source terms. In the C case, the kinetic energy production is solely balanced by the source of potential energy associated with conduction transport. Interestingly, this latter quantity ⟨ω θ ⟩ v is nearly constant between C, CR(1m) and CR(3m) cases at a given Rayleigh number. But in the CR cases, radiation acts as an additional source of potential energy (term ⟨ω rad ⟩ v ) resulting in larger kinetic energy production ⟨τ u ⟩ v . In the CR(3m) case at Ra = 10 6 and 10 7 , radiation and conduction sources of potential energy are of the same magnitude but the relative weight of the radiation source decreases with the Rayleigh number and with a cavity of smaller size (CR(1m) case). Balance of the total kinetic energy and the thermal energy will not be discussed here as we will rather focus on the balance of their fluctuating part in the folllowing section.

Turbulence budgets

This section analyses turbulence budgets of velocity and temperature fluctuating fields in the range 10 7 ⩽ Ra ⩽ 10 9 in the C and CR(1m) cases. The turbulent kinetic energy balance

e k (u ′ ) = 1 2 u ′ • u ′ writes ∂e k (u ′ ) ∂t + u • ∇e k (u ′ ) = -∇ • u ′ p ′ + 1 2 (u ′ • u ′ )u ′ - Pr √ Ra ∇e k (u ′ ) δ u ′ +Pr w ′ θ ′ τ u ′ -u ′ u ′ : ∇u ζ u ′ - Pr √ Ra ∇u ′ : ∇u ′ -ε u ′ . (4.10)
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 ′ is a diffusion term. Table 4.5 provides the volume average of these terms (except the diffusion one which vanishes for no slip boundaries). The integrated production term ⟨τ u ′ ⟩ v highly dominates the transfer term ⟨ζ u ′ ⟩ v (by a factor ∼ 10). 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 4.5 shows the vertical distribution of each term in the right hand side of Eq. (4.10), 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 θ (θ

′ ) = 1 2 θ ′ θ ′ writes ∂e θ (θ ′ ) ∂t + u • ∇e θ (θ ′ ) = -∇ • 1 2 θ ′ θ ′ u ′ - 1 √ Ra ∇e θ (θ ′ ) δ θ ′ -θ ′ u ′ • ∇θ ζ θ ′ - 1 √ Ra ∇θ ′ • ∇θ ′ -ε θ ′ - 1 √ Ra θ ′ (P rad ) ′ -ε rad ′ (4.11)
This equation shows 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 fluid, a radiative contribution to dissipation ε rad ′ is added. Table 4.5 provides the volume average of these terms (except the diffusion one which vanishes after integration for the considered boundary conditions). 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 can be explained by the ability of radiative transfer to dissipate large-scale thermal structures while conduction becomes more efficient for small-scale structures [START_REF] Soufiani | Temperature turbulence spectrum for high-temperature radiating gases[END_REF]. Figure 4.6 shows the vertical distribution of each term in the right hand side of Eq. (4.11), 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 average, 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 is very weak.

Reorientation dynamics in cubic Rayleigh-Bénard cells

References [A9], [A11], [A13], [C4]

One of the main features of interest of Rayleigh-Bénard convection is the large-scale motion from bottom hot wall to top cold wall, because it carries most of the kinetic energy and most of the heat. This large-scale motion, referred to as large-scale circulation (LSC) in the literature, potentially changes of orientation and shape over time, intermittently. 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] as represented in Fig. 4.7. 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 Rayleigh-Bénard convection[END_REF].

Several models have been proposed to explain the LSC dynamics in RBC for non-radiating gases. Brown and Ahlers [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] 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 in Ref. [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], Proper Orthogonal Decomposition (POD) modes [START_REF] Bailon-Cuba | Aspect ratio dependence 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] or spectral POD modes [START_REF] Thebault | Numerical study of the coherent structures in a transitional vertical channel natural convection flow[END_REF].

Up to the best of our knowledge, radiative transfer effects on the LSC dynamics have not been investigated in past research work. The goal of this section is to analyse these effects in cubic Rayleigh-Bénard cells, based on coupled DNS results discussed in Sec. 4.2. We restrict this study to coupled results obtained with the cavity size L = 1 m (CR(1m) case) and in the Rayleigh number 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]. Another objective of this section is to derive a reduced-order model able to capture the dynamics of large-scale coherent structures. To this aim we rely on POD to: (1) extract an optimal set of spatial basis functions called POD modes; (2) build a reduced-order model using Galerkin projection of the coupled system of equations onto the POD basis.

Radiation effects on reorientation dynamics

Rorientation dynamics can be monitored using the time evolution of the x and y components of the angular momentum L = (rr 0 ) × u dr with respect to the cavity centre r 0 . Figure 4.8 shows the time evolution of components L x and L y at the different Rayleigh numbers for the C and CR(1m) cases. Note that time integration is shorter for the highest Rayleigh numbers.

The C 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 characterised 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. 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 CR(1m) case, a significant change in dynamics compared to the C 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 characterised by a zero equilibrium value for one of the two components L x or L y . 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.

In order to quantify radiative transfer effects on the temporal dynamics, two characteristic frequency scales have been computed: 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. Frequencies are reported in Tab. 4.6. In the C case, the circulation frequency is nearly constant with the Rayleigh number as convective time units are used. In the CR(1m) case, the increase of the kinetic energy (discussed in Sec. 4.2) leads to an increase of the circulation frequency. This increase compared to the C case diminishes with the Rayleigh number. Data in Tab. 4.6 confirm the observations of Fig. 4.8: the reorientation frequency decreases with the Rayleigh number and, at a given Rayleigh number, 4.6: Circulation frequency f c and reorientation frequency f r 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 higher Rayleigh numbers, many rapid switches at Ra = 3 × 10 6 .

f c × 10 2 f r × 10 3 Ra C CR(1m) C CR(1m)
is higher when radiation is taken into account. It should be noted however that the uncertainty associated with reorientation frequencies is significant, owing to the few reorientations. If we assume a Poisson distribution of reorientation events [START_REF] Sreenivasan | Mean wind and its reversal in thermal convection[END_REF], the uncertainty on the reorientation frequency f r would be ± (f r /∆t), ∆t being the integration time. 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 seems to be statistically significant and has been reported in other works [START_REF] Sugiyama | Flow reversals in thermally driven turbulence[END_REF].

POD analysis

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) on a low-dimensional subspace. POD eigenfunctions are computed from the eigenvalue decomposition of the spatial correlation tensor of the data [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]:

3 j=1 U i (r, t)U j (r ′ , t)ϕ j n (r ′ )dr ′ = λ n ϕ i n (r). (4.12)
They form an orthonormal basis allowing the decomposition

U (r, t) = ∞ n=1 a n (t)ϕ n (r), (4.13)
where the projection coefficients a n (t) are statistically uncorrelated and their energy is equal to the eigenvalue λ n . The eigenvalue associated with a POD mode is thus a measure of its energy content and the objective is to restrict the decomposition in Eq. (4.13) 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). In thermal convection, we define the flow variable vector as U = {u, γθ}, where γ is a rescaling factor such that velocity and temperature fields have the same energy in the statistical sense [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 }. The eigenvalue problem (4.12) 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, covering the time sequences shown in Fig. 4.8, is extracted from each simulation. The four possible quasi-stable flow states (planar or diagonal) are not necessarily visited or not equally represented during the time sequence, although there is no physical evidence suggesting that these states are not equiprobable. In order to enforce an equal statistical weight for each flow state, we have built enlarged snapshot sets, obtained by the action of the symmetry group of the problem on the original snapshot sets [START_REF] Moin | Characteristic-eddy decomposition of turbulence in a channel[END_REF][START_REF] Holmes | Turbulence, Coherent Structures, Dynamical Systems and Symmetry[END_REF]. The Rayleigh-Bénard problem in the cubic cell 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]. When considering radiation, it 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 and work on an enlarged snapshot set made of 16,000 samples.

Results

The POD spectra obtained in the C case at different Rayleigh numbers are shown in the left part of Fig. 4.9. Each spectrum is normalized such that n λ n = 1 but since the total POD energy is nearly constant with the Rayleigh number in the C 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.

At high Rayleigh number, there is a different behaviour of large scales and intermediate scales, with a very fast decay of the low-order modes (n ≲ 10) 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, the energy contained in the intermediate scales (10 ≲ n ≲ 1000) is proportionately more important and a fast decay of the energy of the small scales is observed. In the right part of Fig. 4.9 is shown the POD spectrum ratio between CR(1m) and C cases data at each Rayleigh number. This ratio is always greater than one (except for one mode at Ra = 10 6 ) which confirms that the POD spectrum captures the energy increase associated with gas 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. The overall energy increase associated with radiation effects weaken when the Rayleigh number increases, as discussed in Sec. 4.2.

Figure 4.10 shows the first seven POD eigenfunctions at Ra = 10 7 in the coupled case and associated amplitude. 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. We briefly mention their properties and physical meaning:

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

L x L y Symmetries M - - S x S y S z S d L x X - S x AS y AS z - L y - X AS x S y AS z - D - - AS x AS y S z S d BL x X - S x AS y AS z - BL y - X AS x S y AS z - C - - S x S y S z AS d
Table 4.7: Name, angular momentum and symmetry properties of the POD modes in Fig. 4.10. 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 L x and L y modes form a pair of degenerate modes (the L y mode is the image of the L x mode 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 form a single large-scale diagonal roll.

• The D mode is an 8-roll mode that transports fluid from one corner to the other and strengthens the circulation along the diagonal. Its amplitude a D (t) changes sign at each reorientation events.

• 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.

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

To complete this description, the contribution to angular momentum and the symmetry properties of these modes are given in Tab. 4.7. 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 the connection between the LSC and the boundary layers. The D and C modes have opposite symmetry properties in the horizontal plane, highlighting their opposite effect on the LSC (stabilising for D, destabilising for C).

A remarkable feature is that these seven eigenfunctions are found at each Rayleigh number, for both C and CR(1m) cases, within the first twelve POD eigenfunctions (though with a different ordering sometimes). We also found that the other modes, which may appear within the first twelve eigenfunctions, interact weakly with this key 7-modes set and are not necessary to model the flow dynamics and reproduces correctly reorientation frequencies. When we project this key basis set from one condition to another (different Rayleigh number, different radiation condition), respecting the order, differences are smaller than 5 %. To conclude, POD eigenfunctions are globally preserved in the Rayleigh number range 10 6 ≤ Ra ≤ 10 8 , whether gas radiation is considered or not. However, POD eigenvalues significantly increase when gas radiation is taken into account and the distribution of the eigenvalue spectrum varies with the Rayleigh number. In the following, the set of 7-modes described above (modes M , L x/y , D, BL x/y and C) will serve as a basis to derive the reduced-order model.

POD modelling

Methodology POD low-order models are derived from Galerkin projection of Navier-Stokes equations onto the POD basis set. Using the decomposition (4.13), 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 + L R nm )a m (t) + Q nmp a m (t)a p (t) + T n (t), (4.14)
where L B nm , L D nm and L R nm are linear coefficients associated with buoyancy, diffusion and radiation and Q nmp are quadratic coefficients associated with advection. They can be directly determined from spatial modes ϕ n . T n (t) is a closure term, which models the effect of the truncation and corresponds to an equivalent dissipation term. Following Ref. [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.15)
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 σ.

Two reduced-order models, depending on whether gas radiation is considered or not, are derived: (1) an uncoupled model which aims at reproducing the dynamics of the C case in an a posteriori fashion; (2) a predicted coupled model which aims at reproducing the dynamics of the CR(1m) case in an a priori fashion, from DNS data of the C case ignoring radiation effects. For these two models, coefficients are computed as follow (see Tab. 4.8):

• uncoupled model Coefficients L B nm , L D nm and Q nmp are computed from POD eigenfunctions ϕ n extracted from the C case (L R nm = 0). Adjustable coefficients L A nm are determined for equilibria a n (t) = a eq n corresponding to diagonal states from the DNS (da eq n /dt = 0), equilibrium values a eq n being computed from the POD eigenvalues. Finally, the noise level σ is estimated from the ratio of the energy transferred from the low-order truncation to the next higher-order scales (computed using the POD eigenspectrum).

• predicted coupled model Coefficients L B nm , L D nm , Q nmp and L A nm are taken from the uncoupled model and coefficients L R nm are computed from POD eigenfunctions ϕ n extracted from the C case. This is justified by the overall preservation of POD spatial eigenfunctions between C and CR(1m) cases. However, the POD eigenvalues (energies) are not known a priori: we estimate them from equilibrium relations (da eq n /dt = 0), the adjustable coefficients L A nm being known. Finally, the noise level σ is determined from the energy increase of the large-scale modes compared with the uncoupled model (computed using estimated energies from equilibrium values a eq n ).

It is worth noting that radiative transfer effects can be considered as a linear contribution in the predicted coupled model owing to the weak temperature differences. The radiative emission (proportional to T 4 ) in the radiative transfer equation can be linearised around the mean reference temperature T 0 . Because the radiative transfer equation is linear, we can then use the superposition principle and write the radiative intensity as a sum of partial intensities associated with each term of the POD decomposition of the temperature field. It is therefore possible to define a modal-radiative power P rad,n , corresponding to the radiative response to a temperature eigenfunction ϕ θ n and such that P rad = n a n (t)P rad,n .

Results Uncoupled model and predicted coupled model are integrated at each Rayleigh number for a time period of 40,000 dimensionless time units (4 to 8 times longer than in the DNS). Results discussed below have been obtained using a 12-mode truncation (hereafter referred to as 11-D truncation as the mean mode amplitude a M is constant) but restriction to the key 7-mode set discussed in Sec. with the eigenvalues λ n and λ rad n found in the C and CR(1m) DNS. Except for the lowest Rayleigh number, corresponding to a large radiative increase which is overpredicted by the model, the total energy increase associated with radiation effects is correctly predicted by the model (within approximately 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 CR(1m) case is not always captured. Table 4.9 compares frequencies of zero crossings for modes L x , L y and D, obtained from the model and the DNS for both uncoupled and coupled cases. Frequency f D corresponds to the reorientation frequency f r , as the D mode changes sign at each reorientation. Uncertainties on model frequencies are about 4 to 8 times smaller than uncertainties on DNS frequencies, according to the respective integration times. Generally speaking, a good agreement is observed between model and DNS results (discrepancies of 10 % to 30 %, of the order of the uncertainty on DNS frequencies), except at Ra = 10 8 in the coupled case 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. At the highest Rayleigh numbers, it was 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.

Figure 4.12 shows histograms of modes L x and L y at different Rayleigh numbers obtained with the predicted coupled model and the coupled DNS (CR(1m) case). 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. 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 Fig. 4.12 that the system spends considerable time near the roll states (L x = 0 or L y = 0). 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 coupled 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.

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 . A linear stability analysis of the diagonal equilibria has been performed for both the uncoupled model and the 4.9: Average frequencies f n = 1/τ n in the DNS and reduced order model where τ n is the average time between zeros of a n (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 ). predicted coupled model at Ra = 10 6 and the associated eigenvalue spectrum is shown in Fig. 4.13 for 6-D and 11-D truncations. In the uncoupled model, the spectrum for both 6-D and 11-D truncations shows: (i) a pair of very stable eigenvalues associated with an eigenvector with a strong component along the D mode; (ii) two real eigenvalues associated with eigenvectors coupling L x/y and BL x/y modes; (iii) a pair of eigenvalues close to the stability limit associated with an eigenvector with a strong C component. This last point confirms the key role of the corner flows associated with the C mode in the reorientation process. In the predicted coupled model, the least stable pair of eigenvalues associated with mode C becomes markedly unstable. 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 corresponding to roll modes. 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. In summary, the predicted coupled model is able to foresee radiative transfer effects on the flow dynamics in an a priori fashion, that is from simulation data ignoring gas radiation. The model fairly captures the increase of energy of the mode amplitude, the reorientation frequencies and the change in dynamics observed at Ra = 10 6 .

Coupled transport phenomena in shallow caves

References [A15, A17], [P8] Shallow caves are partially or fully confined cavities near the surface (about 10 meters below ground), resulting from the chemical erosion of the surrounding rock, and sensitive to external climate conditions. The prediction of the climate in such caves is of primary interest for the analysis of speleothems and paleoenvironments [START_REF] Fairchild | Speleothem Science: From process to past environments[END_REF][START_REF] Jeannin | Assessing karst-hydraulic hazards in tunneling -the Brunnmühle spring system-Bernese Jura, Switzerland[END_REF][START_REF] Borsato | Temperature and altitudinal influence on karst dripwater chemistry: Implications for regional-scale palaeoclimate reconstructions from speleothems[END_REF], the study of karst flora and fauna [START_REF] Culver | The Biology of caves and other subterranean habitats[END_REF][START_REF] Perry | A review of factors affecting cave climate for hibernating bats in temperate North America[END_REF] or the conservation of parietal prehistoric paintings [START_REF] Quindos | Study of the environmental variables affecting the natural preservation of the Altamira Cave paintings located at Santillana del Mar, Spain[END_REF]107]. Painting damages may occur when this microclimate is disturbed. Human visits may result in significant climate perturbations in many different ways, directly (increase in temperature, humidity and CO 2 concentration, resulting in enhanced condensation and corrosion on cave walls [START_REF] Sánchez-Moral | Inorganic deterioration affecting the altamira cave, n spain: quantitative approach to wall-corrosion (solutional etching) processes induced by visitors[END_REF]) or indirectly (artificial ventilation, modification of cave entrance as in Lascaux [START_REF] Brunet | Re-establishing an underground climate appropriate for the conservation of the prehistoric paintings and engravings at Lascaux[END_REF] or Marsoulas [START_REF] Bourges | Conservation of prehistoric caves and stability of their inner climate: Lessons from Chauvet and other French caves[END_REF]). A deep understanding of the physical mechanisms driving heat transfer inside a cave is therefore necessary to improve the conservation of painting cave heritage.

Daily variations of external thermal conditions do not usually affect the climate of confined caves, as the heat conduction in the surrounding rock fully damps them. But seasonal variations generate temperature differences of the order of 0.1 K in shallow caves, between different zones in the cavity at a given time, which is sufficient to induce a convective airflow. For a cave height of about 5 m, the Rayleigh number is of the order of 10 9 , which means this flow is potentially turbulent, at least in some parts of the cave. In addition, convection is coupled with other heat transfer mechanisms such as the conduction in the rock, the radiative exchanges in the infrared, and the energy transfer associated with water phase change. Relative humidity is often close to 100 % in confined caves, due to thin liquid films of percolating water on the walls, and temperature variations may produce significant condensation and evaporation fluxes at the walls. Also, the presence of saturated water vapour and carbon dioxide (whose molar fraction can reach several percents in caves [START_REF] Houillon | Assessing cave internal aerology in understanding carbon dioxide (co2) dynamics: implications on calcite mass variation on the wall of Lascaux Cave (France)[END_REF]) makes the cave atmosphere radiating.

The numerical investigation of cave's climate is restricted by the complex topology of the cave and the convection time scale, which is of the order of one minute (typical velocities are of the order of 1 cm.s -1 ) while the time scale for climate studies is several decades. A simple approach to account for convection is to model the wall heat flux using a heat transfer coefficient, as it is done for instance by Li et al. [START_REF] Li | Behavior of an underground stone chamber with 3-d modeling to restrain water-related damage to mural paintings[END_REF] for modelling the microclimate of the Takamatsuzuka Tumulus chamber in Japan. Empirical estimates of heat transfer coefficients are also used by Guerrier et al. [START_REF] Guerrier | Climatology in shallow caves with negligible ventilation: heat and mass transfer[END_REF], who derivve a coupled model of the climate of a shallow cave over several decades including radiative exchanges between the opaque walls (though assuming a transparent atmosphere). In contrast, other authors rely on computational fluid dynamics to get better insights on heat transfer in caves [107,[START_REF] Lacanette | An Eulerian/Lagrangian method for the numerical simulation of incompressible convection flows interacting with complex obstacles: Application to the natural convection in the Lascaux Cave[END_REF]. Lacanette et al. [START_REF] Lacanette | An Eulerian/Lagrangian method for the numerical simulation of incompressible convection flows interacting with complex obstacles: Application to the natural convection in the Lascaux Cave[END_REF] develop specific numerical methods to solve for air velocity, temperature and moisture fields in an actual cave geometry corresponding to the Lascaux cave in France. The flow was assumed to be laminar, although the highest estimated Rayleigh number in the simulations reached about 10 8 .

The objective of this study is to improve the understanding of coupled heat transfer mechanisms within a cave and to assess the relative significance of turbulent convection, radiative transfer, and evaporation/condensation processes on cave's climate. Numerical simulations are performed in the simplified case of a parallelepiped cavity with prescribed thermal boundary conditions representative of external climate conditions at different times of the year. This approach provides the orders of magnitude of the various heat fluxes through the cave walls, including the time and location of condensation events, likely to damage painted walls.

Problem statement

We consider a confined parallelepipedal cavity embedded in a rock massif as shown in Fig. 4.14, whose size and position correspond to the entrance passage of the Lascaux cave in France (Hall of Bulls). Compared to previous sections, please note that x and z axes have been swaped and that the x axis is now along the vertical. In order to focus on cave's climate, the transport processes in the rock associated with large time scales (conduction time scale is of the order of the year) are decoupled from transport processes in the cave associated with small time scales (convection time scale is of the order of few minutes). Realistic wall temperatures are estimated using a simplified large-scale model of the conduction in the rock over a year, and then used as boundary conditions of a detailed LES model of coupled heat transfer processes in the cave for given times of the year.

Large-scale model

The large-scale model is a 3D extension of the 2D model defined in Ref. [START_REF] Guerrier | Climatology in shallow caves with negligible ventilation: heat and mass transfer[END_REF]. It consists in solving the heat equation in the rock (pure conduction) and neglecting convection in the cavity. A yearly periodic temperature is imposed at the upper surface of the massif T ex = T m + (δT /2) cos(2πt/τ ) with T m = 285.15 K, δT = 16 K. The lateral and bottom sides of the massif are adiabatic. At cave walls we impose that the conductive flux in the rock equals the radiative flux between the walls (walls are assumed to be black and air is assumed to be transparent). The contributions of conducto-convective fluxes and latent heat fluxes at cave walls are ignored. The model is solved using the finite element software Comsol Multiphysics until it converges towards a yearly periodic regime. From the periodic time evolution of the temperature field in the rock, six instantaneous temperature distributions at cave walls are extracted at t = 0.11τ (August), t = 0.20τ (September), t = 0.37τ (November), t = 0.61τ (February), t = 0.70τ (March) and t = 0.87τ (May). The periodic evolution of the temperature of the edge of the caves is shown in Fig. 4.14. The six cases considered are highlighted by filled circles. The temperature of the upper left edge evolves in significant and different proportions compared to that of the other edges, because it is the closest to the ground surface (the ground surface is inclined at 10 • from the horizontal direction). On the other hand, the temperature of the other edges is almost equivalent and differs from each other slightly.

Coupled LES Large-Eddy Simulation of flow equations (4.2)-(4.5) are solved together with the radiative transfer equation for each considered wall temperature fields, assumed to be timeindependent. Water vapour concentration at cave's walls is prescribed assuming saturation. The Rayleigh number, defined from the cavity height L x and the maximum temperature difference within cavity walls, varies in the range Ra ∈ [2.13 × 10 9 -8.43 × 10 9 ] (∆T ∈ [0.124 K -0.492 K]). Other dimensionless parameters are set to Pr = 0.712, Sc = 0.601 and Ra ϕ = 0.098 Ra. The molar fraction of absorbing species is set to X H 2 O = 0.014, X CO 2 = 0.02 to compute the radiative properties of the fluid. Flow simulations are carried out using a spatial mesh made of 241 × 241 × (32 × 21) collocation points. The Spectral Vanishing Viscosity (SVV) method is employed to model the effects of the unresolved flow subgrid scales. Radiation simulations are carried out using a spatial mesh of (48 × 48 × 128) points and the radiation subgrid model is used to account for the radiation of the spatial scales resolved by the flow mesh but filtered by the radiation mesh. Time integration is performed with a dimensionless convection time step of 0.001 and results are time-averaged over a dimensionless time period of 100, once the statistically steady state is reached (the reference time L 2

x /(a √ Ra) varies in the range ∈ [15 s -30 s]). Cases C, CR and CRM (see Tab. 4.1) are considered to highlight specific effects associated with gas radiation and solutal convection.

Climate regimes

We first focus on results obtained in the C case (ignoring radiation and solutal effects on convection) in order to describe the different convection regimes observed over the year.

Figure 4.15 shows the streamlines of the mean airflow coloured by the kinetic energy for each month studied, in the vertical y mid-plane (the mean flow and mean temperature fields are mostly bidimensional thanks to the symmetry of the equations and boundary conditions with respect to the y mid-plane). For all cases ascending or descending flows develop along the left and right vertical walls. They are connected through horizontal flows along the floor and the ceiling, giving rise to a primary rotating circulation in the cavity. However, we can define two distinct flow regimes depending on the period of the year: a one-cell regime in March, May and August (left panels in Fig. 4.15) characterised by a single large-scale circulation which extends over the entire cavity, and a multiple-cell regime in September, November and February (right panels in Fig. 4.15) for which the primary rotating flow near the walls is associated with more complex flow patterns within the core.

The different patterns are schematically classified in Fig. 4.16 according to the relative temperature of the four edges (given in Fig. 4.14) where rather cold edges are marked in blue and rather hot in red. The direction of rotation of the flow along the walls in the x-z plane is determined by the sign of the horizontal temperature gradient along the z axis. In March, May and February the right vertical wall is on average warmer than the left wall (the horizontal gradient is positive) giving rise to a counter-clockwise rotation. Indeed, due to the buoyancy forces, the hot wall and the cold wall drive the air flow upward and downward, respectively. Conversely, the horizontal temperature gradient is negative in August, September and November resulting in a clockwise rotation. On the other hand, when the floor is on average warmer than the ceiling (March, May, August), the vertical temperature gradient is negative resulting in an unstable thermal stratification in the core and this corresponds to the one-cell regime. When the floor is on average colder than the ceiling (September, November, February), the vertical temperature gradient is positive resulting in a stable thermal stratification and this corresponds to the multiple-cell regime. For example in November, two convection cells of weak intensity are observed. The air layer adjacent to the left wall is heated then rises but the outer part of this layer is cooled by the core then slows down and generates a horizontal current at mid-height. This current travels through the core to the right side then splits with a part incoming to the downward flow adjacent to the cold wall and another part, slightly warmer, moving upward thus creating recirculation cells.

Based on this simplified analysis of thermal boundary conditions, we can conclude that the global circulation along the vertical and horizontal walls is governed by the temperature variations between the vertical walls, i.e., by the horizontal temperature gradient, whereas the nature of the regime, one-cell or multiple-cell, is determined by the temperature variations between the horizontal walls, i.e., by the vertical temperature gradient.

Full coupled results

We compare in this section C, CR and CRM results for the months of May and November, associated with the largest temperature difference, and thus the largest potential heat fluxes.

Flow fields

The profiles of horizontal velocity component w, temperature T and water vapour concentration c w (CRM case only) along the vertical centerline of the cavity (z = L z /2 and y = L y /2) are displayed in Fig. 4.17. Velocity profiles in the C case confirm the flow patterns described in Fig. 4.15: (1) the large-scale circulation is weaker in November than in May, (2) the direction of rotation of this circulation changes between the two months (clockwise in November, counter-clockwise in May) and (3) a secondary-weak cell is noticeable in November as the sign of w varies in the cavity core. The comparison of velocity profiles between C and CR cases shows an intensification of convection when radiation is taken into account for both May and November. This intensification is similar to that observed in Rayleigh-Bénard configuration and can be explained by energy exchanges outside the boundary-layer through convection-radiation coupling. In November, this intensification seems to break the multiple-cell pattern. The vertical temperature profile is little affected by gas radiation in May while a moderate temperature increase in the core is observed in November.

The comparison between CRM and CR models does not show any noticeable effect of solutal convection on temperature and velocity fields. This weak effect of the moisture gradient can be readily explained by comparing the ratio between solutal and thermal buoyant forces given by the ration between solutal and thermal Rayleigh numbers Ra Φ /Ra = 0.098. In addition, the water vapour concentration profile shows a strong similarity with the temperature one, due to the analogy between heat and mass transfer. However, one can note non-monotonic variation of the concentration profile in November which are not observed in the temperature profile and may be explained by a homogenisation of the temperature field through gas-gas radiative exchanges. Table 4.10 gathers macroscopic flow quantities: the kinetic energy of the mean flow e k (u) = Wall heat fluxes The heat flux at cavity walls q wall include three components: a conductive one q cond , a radiative one q rad and a phase change contribution q lat associated with condensation and evaporation

q wall = -λ∇T • n q cond -LD w ∇c w • n q lat +σ SB T 4 - ν Ω•n<0 I ν (r b , Ω) |Ω • n| dΩ dν q rad , (4.16)
where L is the latent heat. Figure 4.18 shows the spatial distribution of the conductive and radiative fluxes in the center plane y = L y /2. We checked that CR and CRM models gave indistinguishable results because of the negligible effect of solutal buoyant forces. Differences between the results from C and CRM models is thus due to the presence of gas radiation in the latter case. Despite the strong convection enhancement associated with gas radiation, the conductive fluxes q cond (blue curves) are little affected by radiation coupling, as it has been observed in Rayleigh-Bénard cells. On the other hand, gas radiation significantly modifies the distribution of the radiative flux q rad along the walls (red curves). The absolute value of the wall radiative flux is rather decreased by gas radiation along the top and bottom walls due to a screening effect. For the left and right walls, the flux is rather increased but the analysis is made difficult by the complex structure of the temperature field in front of these walls.

The latent heat flux q lat is not shown in Fig. 4.18 but is found to be strongly correlated with the conductive flux. Figure 4.19 shows that the dimensionless conductive flux (Nusselt number Nu) almost equals the dimensionless diffusive mass flux (Sherwood number Sh) at cavity walls for both May and November, due to the strong similarities between temperature and water vapour concentration fields. Therefore, the latent heat flux can be directly estimated from the conductive one following

q lat ≃ LD w λ ∂c sat w (T ) ∂T T 0 q cond ≃ 1.7q cond , (4.17)
where c sat w (T ) is the water vapour concentration in air at saturation. The spatial distribution of the three components of the total wall heat flux for the whole cavity is displayed in Fig. 4.20 for May. It can be seen that radiation contribution to the total energy balance at the wall prevails over conductive and latent heat contributions. However, the added contributions of conductive and latent heat fluxes, both driven by convection and such that (q cond + q lat ) ≃ 2.7 q cond , are comparable with (and sometimes greater than) the radiation contribution in many places. The same remark applies for November (not shown). It means that more accurate predictions of wall heat and mass fluxes require to couple natural convection inside the cavity with heat conduction in the surrounded rock, which was not considered to generate wall temperatures. As stated in the introduction, condensation and evaporation are of primary interest for the conservation of painted caves. Condensation results in wall corrosion by limestone dissolution while evaporation induces limestone deposition. In May, condensation mostly takes place on the upper left wall (q lat < 0) with a maximum rate of 18 µm/day while in November (not shown) condensation is observed on the lower right walls with a maximum rate of 3 µm/day. Conversely, in May, evaporation mostly takes place on the lower right wall (q lat > 0) with a maximum rate of 6 µm/day while in November (not shown) evaporation is observed on the upper left walls with a maximum rate of 9 µm/day.

To conclude, the present simulations show that radiation of the cave's atmosphere leads to enhanced convection and mixing, significant change in the wall radiative flux but no effect on the wall conductive flux. Solutal convection can be neglected in front of thermal convection and mass fluxes and latent heat fluxes associated with moisture transport can be directly estimated from temperature gradients. This model gives estimate of period and location of condensation events, though more accurate predictions of condensation events would be obtained using realistic cave geometries and considering the coupling between convection in the cavity and conduction in the rock. * * * High-fidelity simulations of coupled heat transfer in Rayleigh-Bénard cells and shallow caves containing air/CO 2 /H 2 O mixtures have shown similar effects of radiation on turbulent thermal convection. Convection-radiation coupling allows for energy exchanges outside the thermal boundary layers and leads to an increase in the kinetic energy of the flow, an increase of the turbulence intensity and sometimes a change in the large-scale flow patterns. Gas radiation has little effect on the wall conductive flux but significantly alters the radiative flux exchanged between the opaque walls. These effects roughly scale with L 2 / √ Ra: they increase with the macroscopic reference scale L and decrease with the Rayleigh number Ra.

In the Rayleigh-Bénard cell, a predictive low-order model of radiative transfer effects have been derived from uncoupled simulation data. Radiative transfer effects are linear in the model owing to the weak temperature differences. This predicted coupled model is able to capture the energy increase associated with radiation, the characteristic time scale of the flow (such as the reorientation frequency) and potential changes in the large-scale flow dynamics. Given the high computational cost of radiation calculations, it is a valuable exploratory tool for investigating radiative transfer effects from uncoupled simulations. It could be potentially used to estimate the radiative power and the radiative flux in coupled shallow cave simulations, providing that a representative POD basis can be derived.

Conclusion

This section reviews the contributions presented in this professorial thesis to radiative property modelling, radiative transfer simulation and coupled convection-radiation interactions. Perspectives are drawn to guide future works and to deepen the knowledge gained and the results obtained. In addition, it is highlighted how this computational work relates to further experimental investigations and it is discussed how the numerical simulation of radiative transfer and coupled interactions could benefit from the emergence of machine learning and deep learning algorithms.

Radiative property modelling An HSNB model for radiative property modelling of hypersonic plasma flows was derived in Sec. 2.2. The formulation of the model in terms of transmissivity provides a flexible framework for adding many radiative contributions, whereas implementing alternative k-distribution models with a large number of radiative processes can be tedious. The approach is particularly efficient for treating thick molecular radiation and is compatible with disequilibrium conditions, using either multi-temperature or electronic-specific descriptions of the thermal state of the matter. Future work should focus on ablative thermal protection systems and new entry environments like the Martian atmosphere.

A contribution to the backscattering modelling of soot aggregates was presented for LIDAR applications in Sec. 2.3. Results show that the monomer radius of the aggregates could be a discriminating parameter in LIDAR signal processing. In addition, multiple-scattering effects may become significant at short-wavelengths and large optical thicknesses. Overall, the results highlight the need to better characterise the soot particle morphology and optical index. This requires the development of numerical models to simulate morphologies that account for complex features such as necking, overlapping and coating. Future work should also address the LIDAR processing of clouds made of polydispersed populations of aggregates, where morphological parameters vary in space.

Additional perspectives concern the development of specific models for mixtures of radiating gases and particles. A narrow-band model combining contributions of gas and soot particles was derived by Prasanna et al. [START_REF] Prasanna | Modelling radiative properties of participating species in a microwave plasma reactor for diamond deposition[END_REF] to model radiation in a microwave plasma reactor used for synthetic diamond deposition. In the field of combustion, an FSCK model was developed in Ref. [START_REF] Wang | Full-spectrum k-distribution look-up table for nonhomogeneous gas-soot mixtures[END_REF] to account for H 2 O and CO 2 radiation combined with soot radiation, whose properties are assumed to follow the Rayleigh scattering theory. An alternative approach for combustion systems relies on the CK narrow-band model and the RDG-FA theory for a more accurate description of soot morphology [START_REF] Torres-Monclard | A Quasi-Monte Carlo method to compute scattering effects in radiative heat transfer: Application to a sooted jet flame[END_REF]. The spectral dynamics of soot radiative properties is weak compared to line spectra associated with molecular or atomic radiation, so that soot contributions can be easily accounted for. However, an accurate description of soot properties with respect to the microphysical parameters characterising their morphology still requires further investigation.

Radiative transfer simulation Section 3 presented two strategies to address the angular dependence of the radiation field: one based on a Haar-wavelet adaptive angular discretisation and one based on a POD reduced order model. Both methods showed the advantage of using different types and numbers of angular basis functions across the optical thickness range. To efficiently calculate the radiative power, the computational effort has to be made on the intermediate optical thickness classes, as the thin classes do not contribute much to the radiative power and the thick classes behaviour is close to isotropic. This key result could be easily extended to any angular discretisation method, providing that an error estimate on the radiative power is derived.

The computational efficiency of the angular adaptivity could be further enhanced by coupling with spatial adaptivity, where the spatial mesh is only refined in the regions of interest. Space-angle adaptive algorithms have been developed in the literature [START_REF] Park | Coupled space-angle adaptivity for radiation transport calculations[END_REF][START_REF] Safarzadeh | Hybrid space-angle adaptivity for wholecore particle transport calculations[END_REF], but their efficient implementation in parallel would require further developments to address load balancing issues, especially in the case of coupled flow-radiation calculations. Regarding the POD reduced-order model, it could be noted that angular adaptivity (varying the angular basis size across space) has yet to be applied but would be possible given that POD angular functions form a hierarchical orthogonal basis. Finally, it would be worth investigating the extension of both methods (angular adaptivity and POD) to the case of inhomogeneous media. The radiative properties of such media could be still modelled with a global model, and it would still be possible to assign a different angular resolution to each optical thickness class.

Another strategy for future developments in the field would be to rely on the method of moments. In transport and kinetic theory, moment models are used to reduce the size of the phase space required to describe kinetic distributions. It has been widely applied to modelling polydispersed populations of droplets or particles, as it is done for instance in Ref. [START_REF] Vié | Size-velocity correlations in hybrid high order moment/multi-fluid methods for polydisperse evaporating sprays: Modeling and numerical issues[END_REF] for modelling evaporating sprays. In radiative transfer, the method of moments consists in deriving transport equations for a finite set of moments M n,m that are angular integrals of the radiative intensity defined as M n,m ν (r) = where µ and φ are the two variables tracking the propagation direction (respectively, the cosine of the polar angle and the azimuthal angle). However, the transport equations for moments are underdetermined and require closure. The closure must also fulfil realisability conditions, namely, ensuring physically-relevant solutions. For instance, the P N method (which is connected to the method of moments) does not fulfil such conditions as it does not necessarily lead to positive values of the radiative intensity when reconstructing. Entropy-based closures [START_REF] Hauck | High-order entropy-based closures for linear transport in slab geometry[END_REF] have been derived to this aim. Such method naturally retains fundamental properties from the kinetic formalism (such as hyperbolicity, entropy dissipation, and positivity) but their disadvantage is their high computational cost. Another promising closure for thermal radiation is the extended quadrature method of moments (EQMOM) closure [START_REF] Yuan | An extended quadrature method of moments for population balance equations[END_REF][START_REF] Vikas | Radiation transport modeling using extended quadrature method of moments[END_REF] which can capture both kinetic (directional) and diffuse (close to isotropic) regimes. Whatever the closure, a great attention must be paid to derive a suitable numerical scheme and appropriate boundary conditions for moment transport equations.

Coupled convection-radiation interactions Numerical results of Sec. 4 provided novel insights on radiative transfer effects on turbulent thermal convection. Convection enhancement through convection-radiation coupling in the bulk region was demonstrated in the idealised Rayleigh-Bénard configuration. The same trend was confirmed for cavities associated with complex wall temperature distributions, representing different climate conditions in shallow caves. Radiation leads to an increase in turbulent kinetic energy and an increase in temperature fluctuations but to a lesser extent because of radiative dissipation. Radiation may also alter the shape and orientation of the large-scale coherent structures and affect the flow regime. Nevertheless, the high computational cost of convection-radiation coupled simulations restrict these numerical studies and raise the question of efficient low-order descriptions of the essential features associated with radiation effects. It has been demonstrated that the POD method provides a robust description of the dynamics of large-scale structures in Rayleigh-Bénard cells. However, a drawback of the POD description is that POD modes, like Fourier modes, have global support and hardly capture localised coherent structures such as thermal plumes. Thermal plumes are an essential feature of thermal convection. They are coherent structures that carry most of the convective heat flux. The plumes are described as sheet-like structures near the boundary layers, progressively becoming stem-like and mushroom-like as they move into the bulk region. Many different criteria have been proposed in the literature to extract plumes from numerical or experimental temperature fields, but identifying such coherent structures consistently and objectively remains a long-running question.

Another limitation of the POD comes from its optimisation criterion based on the total energy of the fields (or variance of fluctuating fields if the mean is subtracted), which is suitable for extracting large-scale features but not for extracting small-scale features and their role in energy dissipation. Kinetic and thermal dissipation of velocity and temperature fields are key parameters of the scaling theory of thermal convection derived by Grossmann and Lohse [START_REF] Grossmann | Scaling in thermal convection: a unifying theory[END_REF]. Different regimes were identified depending on whether the bulk or boundary layer regions contribute the most to the dissipation. The theory was further refined by splitting the thermal dissipation rate into a plume region and a background region [START_REF] Grossmann | Fluctuations in turbulent Rayleigh-Bénard convection: The role of plumes[END_REF][START_REF] Vishnu | Statistics of thermal plumes and dissipation rates in turbulent Rayleigh-Bénard convection in a cubic cell[END_REF]. Deriving analysis tools capable of extracting coherent structures associated with thermal plumes and discriminating their contribution to dissipation should guide future works on thermal convection of radiating fluids. Especially, radiative transfer effects on plume dynamics and on the dissipation distribution need to be clarified. In the following, I highlight two promising approaches that could potentially address these goals: Latent Dirichlet Allocation and dissipation-optimised POD.

• Latent Dirichlet Allocation (LDA) Recently, Frihat et al. [START_REF] Frihat | Coherent structure identification in turbulent channel flow using latent Dirichlet allocation[END_REF] introduced a probabilistic method, LDA, that can extract localised latent factors from data. The method was originally developed in the context of natural language processing [START_REF] Blei | Latent Dirichlet allocation[END_REF], where it aims to extract topics from a collection of documents, each of which is represented by a non-ordered set of words taken from a vocabulary. The topics are defined by two distributions: the distribution of each document on the topics (each document is associated with a mixture of topics) and the distribution of each topic on the vocabulary words (a combination of words represents each topic). The method has been adapted to describe turbulent flows with the snapshots of the velocity field as documents and the cells (or pixels) as the vocabulary. The topics produced by the decomposition, renamed motifs, correspond to fixed (in the Eulerian sense), spatially coherent regions of the flow. The method was applied to the analysis of the turbulent Reynolds stress in wall turbulence [START_REF] Frihat | Coherent structure identification in turbulent channel flow using latent Dirichlet allocation[END_REF] and to pressure anomalies [START_REF] Fery | Learning a weather dictionary of atmospheric patterns using latent Dirichlet allocation[END_REF], where it was found to be well suited for the representation of sparse (intermittent) data. A first attempt to apply LDA to Rayleigh-Bénard convection has been carried out from 2D snapshots of the convective heat flux extracted from the database presented in Sec. 4.2 in the non-radiating case. For all Rayleigh numbers, the dominant motifs displayed in Fig. 5.1 consisted of elongated vertical structures located mostly within the vertical boundary layer, at a height of 0.25L. The width of these motifs scales with the boundary layer thickness and their prevalence depends on the large-scale organisation of the flow. The motifs can be classified according to their size, their prevalence and their position in either the boundary layer, the entrainment region or the bulk, which provides a valuable description of the organisation of the coherent structures of the flow [A16].

• Dissipation-optimised POD A novel POD approach based on an optimal representation of the gradient of the fields has been recently proposed by Olesen et al. [START_REF] Olesen | Dissipationoptimized proper orthogonal decomposition[END_REF]. This formalism provides a spectral decomposition of the strain rate tensor S ij and is called dissipation-optimised POD (d-POD), while the standard approach is referred to as kinetic energy-optimised POD (e-POD). The d-POD modes ψ ij n (r) are tensor modes, solution of the eigenvalue problem S ij (r, t)S kl (r ′ , t)ψ kl n (r ′ )dr ′ = λ n ψ ij n (r), (5.2) while e-POD modes are solution of Eq. (4.12). Applying this method to the study of a turbulent channel flow shows that d-POD modes display a richer small-scale structure compared to standard e-POD modes. In addition, the d-POD modes are better for field reconstruction in the boundary layer region, while e-POD modes are better for field reconstruction in the bulk region. Applying this method to thermal convection would provide a new classification of flow structures according to their contribution to the kinetic and thermal dissipation rate in radiating and non-radiating fluids.

Simulation and experiments

The research work presented in this professorial thesis is fully based on numerical simulation. Numerical simulations have been ubiquitous in science since the rapid growth of computers and aim at replacing experiments when they are too costly or not feasible. However, simulations cannot question or contradict the model they are built upon [START_REF] Jebeile | Computer simulation, experiment, and novelty[END_REF], contrary to experiments which keep a special place in scientific methodology. Experiments in the field of radiative heat transfer face the difficulty of defining a limited set of nondimensional parameters when working with real absorbing gases. It is due to the complex structure of absorption spectra of radiating molecules and the resulting infinite number of optical thicknesses for a given length. Extrapolating results obtained with small-scale experimental devices is then restricted. Experimental investigations are also challenging as non-intrusive measurement techniques are required when radiation comes into play.

In the field of atmospheric entries, ground experimental facilities attempt to reproduce inflight conditions undergone by space vehicles. Shock tubes are designed to replicate the shock wave ahead of the vehicle and the surrounding aerodynamic environment. The ESTHER facility [START_REF] Da Silva | Qualification of the european shock-tube for high enthalpy research[END_REF] is an example of a device that can reproduce any entry condition for all solar system planets (except Jupiter). Alternatively, plasma torches are used to recreate the chemical and thermal processes and study the response of ablative materials in recent studies: typical diagnostics are based on thermograms to measure the temperature of the material [START_REF] Fagnani | 3D infrared temperature maps measurements of ablative materials during plasma wind tunnel experiments[END_REF] or emission spectroscopy to validate radiation models [START_REF] Mcguire | Measurements and modeling of CO 4th positive (A-X) radiation[END_REF]. Results obtained with the HSNB model in Sec. 2.2 would benefit from experimental comparisons and validations.

Radiation effects on thermal convection in cavities have been investigated experimentally by several researchers. The insulation of the sides of the cavity, while mainting optical access, is challenging because of the low conductivity of radiating gases. Also, radiation of water vapour can not be studied experimentally at room temperature because of potential condensation and evaporation and researchers rather use carbon dioxide [START_REF] Hutchinson | Effect of nongray gas radiation on thermal stability in carbon dioxide[END_REF] or ammonia [START_REF] Prasanna | Convection induced by radiative cooling of a layer of participating medium[END_REF] as radiating fluids. Recently, the convection enhancement through convection-radiation coupling has been demonstrated by Kogawa et al. [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] in differentially heated cavities in the Rayleigh number range 10 8 -10 9 , using trifluoromethane as a radiating fluid. Still, experimental investigations of thermal radiation effects in Rayleigh-Bénard cells at high Rayleigh numbers would be required to confirm the numerical results discussed in Sec. 4.2. Deep learning applications to radiative transfer modelling have been attempted with different network architectures and for various purposes. A dense network (where neurons between two layers are fully inter-connected) was derived by Zhou et al. [START_REF] Zhou | A machine learning based full-spectrum correlated k-distribution model for nonhomogeneous gas-soot mixtures[END_REF] for predicting k-distributions associated with mixtures of soot and radiating gases such as CO, CO 2 and H 2 O, in a more efficient manner then pre-existing look-up tables. Dense networks were also used to predict the transmissivity of atmospheric layers given temperature, pressure and composition profiles [START_REF] Stegmann | A deep learning approach to fast radiative transfer[END_REF]. A widespread architecture in the literature is the Convolutional Neural Network (CNN), where filtering operations are applied between each network layer. Among the CNN architectures, the Unet [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] is a specialised one that makes predictions at every spatial grid point. For instance, Lagerquist et al. designed an Unet to emulate an atmospheric radiative transfer model that can be used as parametrisation in numerical weather predictions. CNN were also employed to simulate non-LTE populations of atomic hydrogen in astrophysical radiation models [START_REF] Chappell | SunnyNet: A neural network approach to 3D non-LTE radiative transfer[END_REF]. Another family of deep learning algorithms is the Physics-Informed Neural Networks (PINN) that take advantage of an a priori knowledge of the underlying physical equations behind the data. Mishra and Molinaro [START_REF] Mishra | Physics informed neural networks for simulating radiative transfer[END_REF] derived such an algorithm to solve the RTE. The network is trained using a gradient descent method to minimise a loss function that consists of the residual of the RTE computed at selected training points. Finally, a last application of neural networks is super-resolution, which is extrapolating results from a coarse grid to a fine one, as it is done for instance in turbulence [START_REF] Fukami | Machine-learning-based spatio-temporal super resolution reconstruction of turbulent flows[END_REF]. Royer et al. [START_REF] Royer | A new method for modeling radiative heat transfer based on Bayesian artificial neural networks and Monte Carlo method in participating media[END_REF] used this approach to extrapolate results obtained by a Monte Carlo method on a few probe points to the whole spatial domain.

For solving the RTE, neural networks could be employed for closing angular moment transport equations, as suggested in Ref. [START_REF] Scoggins | Machine learning moment closures for accurate and efficient simulation of polydisperse evaporating sprays[END_REF]. In order to speed up radiation calculations in coupled convection-radiation problems, autoencoder-like architecture could be investigated. Autoen-coders are generally used as generative models, as they provide a low-order representation of a dataset in an unsupervised way. They could be adapted to predict radiative power fields from temperature fields as illustrated in Fig. 5.2 for a given radiating mixture and boundary conditions. Finally, future work dedicated to the convection of radiating fluids could focus on applying deep learning to subgrid-scale modelling in Large Eddy Simulations or to the extraction of coherent structures [START_REF] Fonda | Deep learning in turbulent convection networks[END_REF].

So far, deep learning is driven by empirical developments rather than based on a rigorous mathematical theory. The difficulty in interpreting how neural networks work restricts their generalization in physics or engineering. Still, it is possible to relate their ability to contract and separate information with mathematical or physical concepts such as invariance, symmetries, or multi-scale analysis [START_REF] Mallat | Understanding deep convolutional networks[END_REF]. Another issue raised by these algorithms is their computational cost associated with data generation (when data come from numerical simulation) and training. Training costs are proportional to the size of the dataset and to the number of network parameters [START_REF] Schwartz | Green AI[END_REF] so that limiting the size of the networks remains a great concern. Addressing these computational cost and interpretability issues would benefit from a multi-disciplinary approach and cross-interactions between mathematics, physics and computational science. * * * This manuscript sums up ten years of research dedicated to radiative transfer modelling, exploring a diversity of applications related to environmental and energy issues. Numerical simulations are powerful means to improve knowledge of this complex phenomenon, with increased prediction capabilities but remaining challenges towards large-scale and multi-physics models.

5 Figure 2 . 1 :

 521 Figure 2.1: Left: Absorption spectrum of an air/H 2 O/CO 2 mixture (T 0 = 300 K, atmsopheric pressure, molar fractions X H 2 O = 0.02 and X CO 2 = 0.001) and Planck function at T 0 = 300 K. Right: Associated absorption distribution function F discretised in 16 intervals.

Figure 2 . 2 :

 22 Figure 2.2: List of molecular sytems and continuum processes included in the database.

Figure 2 . 3 :

 23 Figure 2.3: FIRE 2 (1634 s). Comparison between LBL, HSNB-Weak and HSNB models. Left: spectral and cumulated incident fluxes at the wall. Right: total radiative power along the stagnation line and differences with LBL calculations normalized by the maximum absolute value.

Figure 2 . 4 :

 24 Figure 2.4: Huygens(191 s). Comparison between LBL, HSNB-Weak and HSNB models. Left: spectral and cumulated incident fluxes at the wall. Right: total radiative power along the stagnation line and differences with LBL calculations normalized by the maximum absolute value.

Figure 2 . 5 :

 25 Figure 2.5: Huygens (189 s). Comparison between LBL, HSNB and LBL results from Bansal and Modest [9]. Contribution of the CN red and violet systems to the radiative source term along the stagnation line.

Figure 2 . 6 :

 26 Figure 2.6: Extinction (left) and backscattering (right) cross sections of SFA ensembles with respect to the volume of the aggregate 4/3πr 3 m N m . The three sets of coloured curves corresponds to SFA ensembles of monomer radius 10 nm (blue), 20 nm (red), and 40 nm (green). The full, dashed and dotted lines are associated with aggregates composed of 45, 125, and 450 monomers respectively

Figure 2 . 7 :

 27 Figure 2.7: LR of SFA ensembles. The three sets of coloured curves corresponds to SFA ensembles of monomer radius 10 nm (blue), 20 nm (red), and 40 nm (green). The full, dashed and dotted lines are associated with aggregates composed of 45, 125, and 450 monomers with D f = 1.8 (left) or associated with aggregates of fractal dimensions set to 1.6, 1.8 and 2.0 with N m = 125 (right).

Figure 2 . 8 :

 28 Figure 2.8: Particle LDR of SFA ensembles. The three sets of coloured curves corresponds to SFA ensembles of monomer radius 10 nm (blue), 20 nm (red), and 40 nm (green). The full, dashed and dotted lines are associated with aggregates composed of 45, 125, and 450 monomers with D f = 1.8 (left) or associated with aggregates of fractal dimensions set to 1.6, 1.8 and 2.0 with N m = 125 (right).

Figure 2 . 9 :

 29 Figure 2.9: Simulated LIDAR signals of the considered sooty plume with a total optical thickness ζ tot = 10 (left) and ζ tot = 1 (right). The receiver FOV is set to θ rec = 5 mrad.

  ) and also shown in Fig.2.9. The signals have a Gaussian shape whose peak is located below the concentration maximum (at z c = 1.15 km) and gets closer to the receiver when the total overall particle concentration (or total cloud optical thickness) increases. The effect of multiple-scattering is to increase the LIDAR signal especially at high optical thickness. The signal is stretched and the peak slightly shifts towards the concentration maximum. The smallest values of the correction factor ξ are associated with the highest multiple-scattering effects. The minimum ξ value is located ahead of the signal's peak and is equal to 0.87 at ζ max = 10 and 0.92 at ζ max = 1. Unphysical oscillations can be noticed in the correction factor profile in the region where the signals are close to zero because of the stochastic nature of the Monte Carlo calculations.The LIDAR signals obtained for different optical thicknesses and different receiver FOV values have been used to derive phenomenological relationships of two quantities that are useful in LIDAR signal inversion: the volume LDR and the multiple-scattering fraction MSF defined as

Figure 2 . 10 :

 210 Figure 2.10: MSF and volume LDR as a function of the soot plume optical depth.

Figure 3 . 1 :

 31 Figure 3.1: Left: representation of the Haar wavelet and S N basis functions over a 1D arbitrary interval at the order 3 (partition of the interval in 4). Right: the transformation matrix that relates the two spaces.

Figure 3 . 2 :

 32 Figure 3.2: Schematic diagram of the street canyon configuration (left) and snapshot of the temperature field (right).

Figure 3 . 3 :

 33 Figure 3.3: Functional error plot against averaged number of angular basis functions (left) and computational time (right).

Figure 3 . 4 : 4 .Figure 3 . 5 :

 34435 Figure 3.4: Averaged number of angular basis functions for each k-class of the global model for goal-based adapt W GB 4,4 .

Figure 3 . 6 :

 36 Figure 3.6: Snapshot of the reduced temperature field (T -T 0 )/∆T considered in the radiation calculations.

Figure 3 . 7 :

 37 Figure 3.7: Eigenvalue spectrum of the POD for five k-classes.

Figure 3 . 8 : 6 Figure 3 . 9 :

 38639 Figure 3.8: The first six POD basis functions, from left to right and from top to bottom, for the third and seventh k-class.

Figure 3 . 10 :

 310 Figure 3.10: Radiative power error norm (as defined in Eq. (3.17)) as a function of the averaged angular basis size (left) and computational time (right). Comparison between discrete ordinates (black plain lines), uniform POD (red dashed lines) and optimised POD* (blue squares) models.

  (a) C case | CR case.(b) Radiation subgrid model is used.

Figure 4 . 2 :

 42 Figure 4.2: Flow fields at different Rayleigh numbers. Temperature, kinetic energy and radiative power, averaged over the horizontal plane and statistically averaged. C case (red lines), CR(1m) case (black lines) and CR(3m) case (blue lines).

Figure 4 . 3 :Figure 4 . 4 :

 4344 Figure 4.3: Heat fluxes at different Rayleigh numbers. Conductive, convective and radiative fluxes, averaged over the horizontal plane and statistically averaged. C case (red lines), CR(1m) case (black lines) and CR(3m) case (blue lines). For radiative fluxes, the wall contribution (reported in Tab. 4.3) has been subtracted to focus on their variations within the cavity.

Figure 4 . 5 :

 45 Figure 4.5: Turbulent kinetic energy budget at different Rayleigh numbers. Production τ u ′ (green lines), transfer ζ u ′ (blue lines), dissipation ϵ u ′ (red lines) and diffusion δ u ′ (black lines) of turbulent kinetic energy for C and CR(1m) cases.

Figure 4 . 6 :

 46 Figure 4.6: Budget of temperature variance at different Rayleigh numbers. Transfer ζ θ ′ (blue lines), dissipation ϵ θ ′ (red lines), diffusion δ θ ′ (black lines) and radiative dissipation ϵ rad ′ (green lines) of temperature variance for for C and CR(1m) cases.

Figure 4 . 7 :

 47 Figure 4.7: Four quasi-stable diagonal states for the large-scale circulation in a cubic Rayleigh-Bénard cell.

Figure 4 . 8 :

 48 Figure 4.8: Time evolution of x and y components of angular momentum L x (blue lines) and L y (red line).

8 Figure 4 . 9 :

 849 Figure 4.9: POD eigenspectrum obtained from C case simulation data, normalized such that n λ n = 1 (left) and POD eigenspectrum ratio between CR(1m) case and C case simulation data (right).

Figure 4 . 10 :

 410 Figure 4.10: First seven POD modes ϕ n (r) at Ra=10 7 (CR(1m) case) with associated amplitude a n (t). Streamlines and isosurfaces of convective heat flux ϕ θ n ϕ w n = 0.25 coloured by mode temperature. Colour map for mode temperature ranges from +0.5 (blue) to -0.5 (red).

Figure 4 . 11 :

 411 Figure 4.11: Predicted energies compared with POD eigenvalues in the uncoupled (C case) and coupled (CR(1m) case) DNS. Energies are ranked according to the mode ordering in the uncoupled simulation.

Figure 4 . 12 :Figure 4 . 13 :

 412413 Figure 4.12: Histograms of modes L x and L y : Coupled DNS (a,c,e,g,i) and coupled model (b,d, f,h,j). Rayleigh numbers are: (a,b) Ra = 10 6 ; (c,d) Ra = 3 × 10 6 ; (e, f ) Ra = 10 7 ; (g,h) Ra = 3 × 10 7 ; (i,j) Ra = 10 8 .

Figure 4 .

 4 Figure 4.14: Top: Geometry setup of the rock massif and the cavity (not at scale). Bottom: Time evolution over a one-year period of the external temperature fluctuation T ex (t)-T m applied at the upper surface of the massif (left) and time evolution of the temperature fluctuations T (t)-T m averaged on four edges of the cave (right). Red, yellow, green and blue lines correspond to the left upper (z = 0 and x = L x ), left bottom (z = 0 and x = 0), right upper (z = L z and x = L x ) and right bottom (z = L z and x = 0) edges, respectively. The six months investigated are marked with filled circles.

Figure 4 .

 4 Figure 4.15: C case. Flow streamlines coloured by the kinetic energy of the mean flow. Streamlines are drawn from the y mid-plane and then projected onto the y mid-plane when they deviate from it.

Figure 4 . 16 :

 416 Figure 4.16: Simplified representation of thermal boundary conditions and associated flow regimes. Blue and red lines correspond to cold and hot edges, respectively.

Figure 4 . 17 :

 417 Figure 4.17: Vertical profiles of the horizontal velocity component w (left) and the temperature (right) along the vertical centerline of the cavity (z = L z /2 and y = L y /2) for the C (dasheddotted lines), CR (dotted lines) and CRM (solid lines) models in May (red curves) and November (blue curves).

Figure 4 .

 4 Figure 4.18: Conductive (blue lines) and radiative (red lines) heat fluxes at the intersection between the center plane y = L y /2 and the left (z = 0), right (z = L z ), bottom (x = 0) and top (x = L x ) walls, computed from the CRM (solid lines) and C (dashed-dotted lines) cases, in May and November.

Figure 4 . 19 :

 419 Figure 4.19: Cloud of points representing the local Nusselt number versus the local Sherwood number for May and November in the CRM case. Each point corresponds to a given spatial location at the wall.

Figure 4 . 20 :

 420 Figure 4.20: Spatial distribution of conductive q cond , latent heat (q lat ), radiative (q rad ) and total (q wall ) wall heat fluxes for the CRM model in May. The left panel corresponds to the top (X = L X ), left (Z = 0) and front (Y = L Y ) cave walls. The right panel corresponds to the bottom (X = 0), right (Z = L Z ) and rear (Y = 0) cave walls.

1 I

 1 ν (r, µ, φ) exp(inφ)µ m dφdµ,(5.1)

Figure 5 . 1 :

 51 Figure 5.1: Characteristic dominant LDA motif in a cubic Rayleigh-Bénard cell at different Rayleigh numbers. The black lines correspond to the boundary layer thickness.

Figure 5 . 2 :

 52 Figure 5.2: Sketch of an autoencoder that could be used for the prediction of radiative power fields from temperature fields.

Table 2 .

 2 1: Conditions for the trajectory points t = 1634 s of the FIRE 2 experiment and for the trajectory point t = 191 s of the Huygens probe entry. Radius of the vehicle R 0 , wall temperature T w and free stream conditions (temperature T ∞ , velocity u ∞ , total mass density ρ ∞ and mass fractions y ∞ ).

	%)

2 , O 2 , NO), photodissociation (O 2 ), photodetachment (N -, O -) Free-free processes N, O, N + , O + , N 2 , O 2 Titan entries Chemical species N, C, H, N 2 , C 2 , H 2 , CN, NH, CH, CH 2 , CH 3 , CH 4 , HCN Thick molecular systems N 2 (Birge-Hopfield 1 and 2, Worley-Jenkins, Worley, Caroll-Yoshino), CN violet, C 2 Swan Thin molecular systems N 2 (First and second positive), CN (red, LeBlanc), C 2 (Philips, Mulliken, Deslandres-D'Azambuja, Ballik and Ramsay, Fox-Herzberg)

Table 2 .

 2

2: Chemical species and radiative systems considered for Earth and Titan entries

Table 2 .

 2

	FIRE 2 (1634 s)	LBL	HSNB-Weak HSNB
	q rad,i w	(W/cm 2 )	146.78	151.94	150.85
	t CPU (s)	20480	41	242
	Huygens (191 s)	LBL	HSNB-Weak HSNB
	q rad,i w	(W/cm 2 )	82.68	104.26	86.24
	t CPU (s)	13158	5	105
	3: Comparison between LBL, HSNB-Weak and HSNB models for FIRE 2 (1634 s), FIRE
	2 (1642.66 s) and Huygens (191 s) cases. Incident radiative fluxes and computational times for
	one radiation calculation.

  mC ext (nm 2 ) C back (nm 2 .sr -1 )

			ω	LR (sr)	δ p
	1.663 + i0.715 74978 ± 763	349 ± 34	0.232 ± 0.006 217 ± 21 0.022 ± 0.003
	Table

  ,

			T cold	e=1	§ filled with air/H 2 O/CO 2
						mixture
	g	x	z	y	L	§ X H2O =0.02 / X CO2 =0.001 § L=1 or 3 m § side vertical walls
						adiabatic and reflecting
		T hot	e=1		e=0
	Figure 4.1: 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.

Table 4 .

 4 .3. 2: Simulation parameters: Convection mesh, radiation mesh, convection time step δt and integration time interval ∆t. 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.

	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 (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

Table 4 .

 4

		125.05	375.15	8.77 119.88	9.12 346.63
	16.24	125.05	375.15	16.66 120.82	16.90 352.55
	31.06	125.05	375.15	31.29 121.97	32.45 357.37
	61.38	125.05	375.15	61.82 122.36	62.51 360.92

3: Conductive and radiative wall heat fluxes at different Rayleigh numbers, averaged over the horizontal plane and statistically averaged.

Table 4 .

 4

	6	5.16	6.99	5.49	1.50	11.05 5.74	5.31
	10 7	3.41	4.23	3.50	0.73	7.44	3.55	3.88
	10 8	2.13	2.49	2.14	0.35	3.51	2.22	1.28
	10 9	1.35	1.47	1.36	0.09	2.02	1.37	0.62

4: Volume-averaged source terms of the potential energy balance.

Table 4 .

 4 5: Volume-averaged source terms of the turbulent kinetic energy balance and the balance of variance temperature fluctuations, in the C and CR(1m) cases at Rayleigh numbers Ra= 10 7 , 10 8 , and 10 9 . Variations due to radiation are provided in parentheses

	5)∇θ) +	Pr √ Ra	∂θ ∂z	-	Pr √ Ra	(z -0.5)P rad . (4.8)

  4.3.2 (hereafter referred to as 6-D truncation) leads to very similar results. A first output of the predicted coupled model are the estimated energies (or POD eigenvalues).Figure 4.11 compares for each mode and each Rayleigh number the predicted energy

		uncoupled model	predicted coupled model
	L B nm , L D nm , L R nm , Q nmp	computed from DNS data of the C case at each Rayleigh number (L R nm = 0)	computed from DNS data of the C case at each Rayleigh number, including L R nm
	L A nm	n /dt = 0) relations (da eq determined from equilibrium	taken from the uncoupled model
	a eq	determined from POD spectrum	determined from equilibrium
			relations (da eq n /dt = 0)
	σ	computed using POD spectrum	computed using equilibria a eq
		Table 4.8: Model parameter determination.

∼ λ rad n

1 2 u

 2 • u, the turbulent kinetic energy e k (u ′ ) = 1 2 u ′ • u ′ and the maximum in absolute value of the vertical velocity component |u| max (• denotes the time average). The kinetic energies are 2 to 3 times higher in the CR case compared to the C case for both May and November, and |u| max is increased by 11 % in May and by 14 % in November. Values do not differ much between CR and CRM case which confirms previous discussion.

			May			November
		C	CR	CRM	C	CR	CRM
	e k (u) × 10 4 (m 2 /s 2 )	2.80	6.80	6.90	0.31	1.10	1.10
	e k (u ′ ) × 10 4 (m 2 /s 2 )	0.93	2.70	3.20	0.07	0.14	0.14
	|u| max (m/s)	0.095 0.105 0.108	0.050 0.069 0.069
	Table 4.10: Macroscopic flow quantities

RADIATIVE PROPERTY MODELLING

COUPLED CONVECTION-RADIATION INTERACTIONS

CONCLUSION

Remerciements

L x < 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L y < 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L x > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L y > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > x < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > y < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > z < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L x > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L y < 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L x < 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > L y > 0 < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " ( n u l l ) " > ( n u l l ) < / l a t e x i t > A Curriculum Vitae