Steady-state operations of the next-generation fusion device ITER (International Thermonuclear Experimental Reactor) will require the development of reliable numerical tools to estimate key engineering parameters suitable for technological constraints at reasonable computational cost. For example, divertor monoblocks will have to be preserved from impinging heat loads exceeding 10 M W/m 2 to prevent possible mechanical damages.

So-called transport codes fulfill this requirement since they rely on two-dimensional fluid equations averaged over time fluctuations, similarly to Reynolds Averaged Navier-Stokes models commonly used for engineering applications in the neutral fluid community. Furthermore, transport codes can gather most of the physical ingredients ruling the edge plasma behavior, as well as realistic magnetic topology and wall geometry.

However, their predictability is limited by a crude description of turbulent (particle, momentum and energy) fluxes perpendicular to the magnetic field lines. In the plasma community, a special concern is devoted to acquire a detailed understanding and the predictive capability of cross-field transport, dominated by turbulence, since it strongly impacts the power extraction and the confinement of plasma over extended periods of time. In transport codes though, turbulent fluxes, which are assumed diffusive, are crudely determined by either homogeneous, or ad-hoc diffusive coefficients, or feedback-loop procedures applied a-posteriori on experimental data.

Motivated by these issues, in this work we introduce step-by-step a new approach with the aim to self-consistently estimate the distribution of turbulent fluxes in transport codes, when steady-state plasmas are concerned. The underlying strategy is inspired by the work done from the 60's in neutral turbulence and adapted here to plasma for fusion applications.

The first key concept is the Boussinesq assumption. It consists in assuming a colinearity between the Reynolds stress tensor -which represents the contribution of turbulence to the mean flow -and the mean rate of strain tensor -expressed by the gradient of the mean velocity through a coefficient: the so-called eddy-viscosity.

The second concept is to express this new eddy viscosity coefficient as a function of characteristic turbulence quantities. As a starting point, in this work, we have focused on the most popular in Computational Fluid Dynamics (CFD), the κ-ε model, where transport equations for the averaged kinetic turbulent energy and the turbulence dissipation rate are designed semi-empirically. Steady-state κ and ε allow for a self-consistent estimation of the eddy-viscosity coefficient, thus including the impact of turbulence in steady-state mean flows.

We propose a κ-ε-like model where two transport equations for turbulent kinetic energy and its dissipation rate are derived algebraically, including the physics of the linear interchange instability. Indeed, the interchange instability is known to drive most of plasma turbulence in the edge and in the scrape-off layer, particularly at the low-field side, where it leads to a ballooned distribution of the cross-field transport. Moreover, on the basis of dimensional analysis, transport coefficients are related to reference turbulent scales of velocity and time according to the so-called quasi-linear theory, widely used in the framework of reduced models for plasmas.

For the numerical implementation, we exploit the flexibility of the transport package SolEdge2D-EIRENE, developed for many years through the collaboration of the IRFM at the CEA and the laboratory M2P2 at Aix-Marseille University.

Since the new model is semi-empirical, it presents some free parameters to be closed. In this work, we have proposed different approaches. First, an automatic fitting procedure to optimize the numerical results with respect to a single experimental test-case. Then, to be more general and increase the predictive capabilities of the model, a reference scaling law for the width of the heat-flux profile in the scrape-off layer has been assumed, empirically determined from the experimental measurements of the outer target heat load in various machines.

The new model is integrated in SolEdge2D-EIRENE for simulations with diverted plasma in TCV and WEST-like geometries, for L-mode discharges. Steady-state results are discussed and shown to favorably compare with experimental data at both the outer mid-plane and the outer divertor. Moreover, self-consistent distributions of diffusivities are shown to exhibit poloidal asymmetries consistently with the ballooned distribution of cross-field transport observed at the same conditions in both first-principle codes and experiments.

Résumé en français

Le fonctionnement à l'équilibre du réacteur à fusion de prochaine génération, ITER, nécessitera le développement d'outils numériques fiables permettant d'estimer les paramètres d'ingénierie clés à un coût de calcul raisonnable. Les codes de transport répondent à cette exigence car ils reposent sur des équations fluides bidimensionnelles qui sont moyennées sur les fluctuations temporelles, de la même manière que les modèles Reynolds Averaged Navier-Stokes couramment utilisés dans la communauté des fluides neutres. De plus, les codes de transport intègrent la plupart des ingrédients physiques régissant le comportement du plasma de bord, ainsi que une topologie magnétique réaliste et la géométrie du mur.

Cependant, leur prévisibilité est limitée par une description inadéquate des écoulements turbulents perpendiculaires aux lignes de champ magnétique, qui affectent fortement le confinement du plasma sur de longues périodes. En effet les écoulements turbulents, supposés diffusifs, sont grossièrement déterminés par des coefficients de diffusion homogènes ou ad-hoc, ou par des procédures basées sur des boucles de rétroaction appliquées a-posteriori à des données expérimentales.

Motivés par ces questions, nous présentons dans ce travail un nouveau modèle pour estimer de manière cohérente la distribution des écoulements turbulents dans les codes de transport pour des plasmas en régime stationnaire. La stratégie consiste à introduire des modèles numériques efficaces largement utilisés dans la communauté de la turbulence neutre dans le contexte de la physique des plasmas. Deux concepts clés issu de la communauté des fluides neutres nous ont inspirés dans ce travail. Le premier est l'hypothèse de Boussinesq. Elle consiste à linéariser le tenseur de contraintes de Reynolds dans l'équation de Navier-Stokes moyennée dans le temps à partir d'une hypothèse de diffusion dans laquelle le terme de proportionnalité est appelé eddy viscosity. Le second concept est le modèle κ-ε, dans lequel les équations de transport pour l'énergie turbulente cinétique et le taux de dissipation de la turbulence sont obtenues de manière semi-empirique. A l'équilibre, κ et ε permettent une estimation auto-cohérente de l'eddy viscosity, intégrant ainsi l'impact de la turbulence sur les écoulements moyennés à l'état d'équilibre.

Ces concepts ne peuvent pas être appliqués directement pour enrichir la modélisation des flux perpendiculaires dans les plasmas en raison de différentes propriétés de turbulence. Par conséquent, nous proposons une adaptation du modèle κ-ε aux plasmas de fusion, où deux équations de transport pour l'énergie cinétique turbulente et son taux de dissipation sont dérivées algébriquement. Ces équations intégrent la physique de l'instabilité d'interchange linéaire, qui conduit à une distribution ballonnée du transport perpendiculaire dans le bord du plasma.

Différentes approches sont décrites afin de déterminer les paramètres libres: premièrement, une procédure de boucle de rétroaction pour optimiser les résultats numériques comparés avec un test expérimental. Ensuite, on utilise une loi d'échelle de référence pour la largeur du profil de flux de chaleur dans la SOL, déterminée empiriquement à partir des mesures expérimentales du flux de chaleur sur le divertor externe dans diverses machines.

Le nouveau modèle est intégré au package SolEdge2D-EIRENE, développé en collaboration entre le CEA et le laboratoire M2P2 de l'Université d'Aix-Marseille. Les résultats numériques à l'état d'équilibre sont discutés et on démontre qu'ils se comparent favorablement aux données expérimentales aussi bien à l'outer midplane qu'au divertor externe. De plus, on montre que les distributions de diffusivité présentent des asymétries poloïdales cohérentes avec la distribution ballonnée du transport perpendiculaire observée dans les mêmes conditions dans les codes de premier principe et les expériences.

1.1. Nuclear fusion

This chapter introduces the main concepts about thermonuclear fusion, magnetically confined plasmas and transport mechanisms. A brief introduction on the 2D transport code SolEdge2D-EIRENE is also provided.

Nuclear fusion

Thermonuclear fusion may have a crucial role in the energy economy of the future, due to almost limitless fuel reserves, negligible impact on the environment and inherent nuclear safety. The physical mechanism consists in a nuclear reaction between two light atoms which are compressed together until they overcome the electrostatic repulsion and combine into a more stable element. This process releases kinetic energy and subatomic particles (neutrons or protons). The amount of released energy depends on the mass defect ∆m = m R -m P between reacting and produced nuclei, according to the Einstein law:

E = ∆mc 2 (1.1)
where c is the speed of light in vacuum. A fusion process that produces a nucleus lighter than 56 Fe or 62 Ni is generally exothermic, namely it yields a net energy release. Indeed, when species with mass number lower than 56 F e or 62 N i fuse, the binding energy of the reaction increases leading to a favorable process to produce energy. In figure 1.1, the reaction rate σv of the three most efficient fusion reactions is displayed in m 3 /s against the temperature in 10 9 K. Temperatures required for a typical fusion reaction to be accessible are in the order of 10 8 K (the temperature of the solar surface is nearly 6000K). The deuterium-tritium fusion reaction, which produces one atom of 4 2 He at 3.5M eV and one neutron at 14.1M eV , is the best compromise between energy demand and gain, with a cross-section of σv = 10 -22 -10 -21 m 3 s -1 :

2 1 D + 3
1 T → 4 2 He(3.5M eV ) + 1 0 n(14.1M eV ) (1.2)

Turbulent transport in edge plasma for magnetic fusion 3

Deuterium is highly available in nature since it can be extracted from water, tritium has to be produced artificially from D but also 6 Li and 10 B.

The amount of energy produced per event by fusion is about 17.6M eV against nearly 200M eV in a typical fission reaction. However, while the fission of 1kg of 235 U supplies about 23 • 10 6 kW h, the equivalent amount of D-T fuel produces about 10 8 kW h, because the atomic mass is much lower in the latter case. Moreover, compared to fission, fusion does not imply chain reactions, thus it is intrinsically safe but needs continuous fueling.

We introduce the so-called amplification factor Q ≡ P f us P in as an indicator of the fusion energy gain, and eventually the viability of a fusion reactor. P f us is the total power produced by nuclear reactions, while P in is the power injected to enable the nuclear reactions.

At the extreme temperatures required by fusion reactions, the fuel is completely ionized, appearing as a globally neutral gas made of ions and electrons. This is the so-called forth-state of matter or thermonuclear plasma. Containing such a system in a solid vessel is unconceivable due to the prohibitive temperatures.

Extreme temperatures are not sufficient to ensure positive Q. In order to have a net production of energy, the Lawson criterion should be observed:

nT i τ E > 3 × 10 21 ][keV sm -3 ] (1.3)
where n is the density in m -3 , T i is the ion temperature in keV and tau E is the energy confinement time in s, namely the decay time of plasma energy when all the power sources are shut off. Inertial confinement and magnetic confinement are the main procedures to exploit fusion reactions in a thermonuclear facility. The former relies on high levels of heating and compression (n ∼ 10 31 m -3 ) of fuel pellets; the latter, which is the mechanism concerned in this work, on the magnetic confinement of charged particles composing the plasma. The next-generation fusion facility which is being built near Cadarache, in France, called International Thermonuclear Experimental Reactor (ITER), is designed to reach Q = 10 ([1]).

Magnetic confinement and tokamaks

By magnetic confinement we mean the constraining of plasma fuel, which is electrical conductive, into a defined reaction environment modeled by magnetic fields, where extreme temperatures are reached. In particular, a system of coils forces charged particles to follow helical orbits around the field according to the Lorenz force F L = qE + qv × B, where q, E and B are the particle charge, the electrostatic and the magnetic fields, respectively. When assuming a uniform magnetic field, the motion of a single charged particle is described by:

• a fast cyclotron rotation on the plane transverse (⊥) to the magnetic field lines, characterized by gyrofrequency Ω s = q s B/m s and radius ρ L,s =

Magnetic confinement and tokamaks

m s v ⊥ /q s B, also called Larmor radius. m s and q s are the mass and charge of species s, while v ⊥ is the average particle velocity in the ⊥ direction.

• a free motion parallel ( ) to the magnetic field lines

Assuming T = 10 8 K and B = 5T and v ⊥ ∼ v th = k B T /m (v th is the thermal velocity and k B is the Boltzmann constant), one finds ρ L,e = 4.5 × 10 -5 m and ρ L,D + = 2.8 × 10 -3 m, much lower than typical machine scales. The parallel length scale, at the contrary, relates to the free motion of charged particles and they are comparable with the machine ones. The solution adopted in tokamaks to confine plasma is thus to twist the magnetic field lines into a helical shape in order to constrain charged particles to move into closed trajectories.

Tokamaks represent the most common type of fusion machine in the framework of magnetic confinement. They consist in toroidal devices where ϕ is the axisymmetric angle around the main axis of the machine, while R and a are the major and minor radius respectively, as shown in figure 1.2. Finally, θ is the poloidal coordinate. In order to twist the magnetic field, a current (I plasma ) is driven in the plasma using the transformer principle. The primary inducing coil in placed at the center of the machine, while the secondary coil is the plasma itself. A poloidal magnetic field B P , induced by the plasma current, combines to the external, toroidal one B T , generated by the toroidal current I T circulating within the toroidal coils. The resulting twisted magnetic field B tot is shown in figure 1.2. B tot is helical, with each magnetic field line lying on one of the nested set of toroidal flux surfaces.

The distribution of the toroidal magnetic field B T varies proportionally to I T according to the Biot-Savart law integrated around the circle of radius R:

B T = µ 0 I T 2πR (1.4)
The shape of the poloidal magnetic field B P , on the other hand, is governed by the plasma current according to the same rule, this time integrated around the circuit of radius a:

B P = µ 0 I plasma 2πa (1.5)
The amplitude of B P is small compared to B T , thus the total magnetic field B tot scales approximately as 1/R, as represented in figure 1 The HFS is the region where the magnetic field is more intense, the LFS the one where it is the weaker.

be defined with respect to the magnetic axis (i.e. the direction of B T ) depending on the intensity of B tot : the so-called High Field Side (HFS) and Low Field Side (LFS). In particular, in the LFS both the magnetic field and plasma pressure decrease, therefore ∇p • ∇B > 0. Such a condition drives the plasma towards an unstable state, as discussed in 1.4.3.1.

Let us focus on the helical paths of magnetic field lines. At some toroidal angle ϕ, the field line has a given poloidal position, which it reaches again after a certain variation of the toroidal angle, ∆ϕ. We define the safety factor q:

q ≡ ∆ϕ 2π (1.6)
which stands for the number of toroidal turns that field lines undergo to complete a poloidal turn. For straight field lines and B P B T , this length can be estimated as:

q = 1 2π 1 R B T B P dθ (1.7)
where dθ is the infinitesimal variation in the poloidal direction while moving through a toroidal angle dϕ. When large values of the ratio R/a and circular cross-section tokamaks are considered, one obtains the approximation:

q ∼ aB T RB P (1.8) 6
1.2. Magnetic confinement and tokamaks

Limiter vs divertor configuration

Necessarily, in finite geometries, magnetic field lines can intercept the vessel and thus hot plasma interacts with wall materials: heavy sputtered atoms and molecules reach the plasma and cool it down. The simplest technical solution to control these heat and particle fluxes arising from imperfect confinement is the insertion of a solid component in the vessel design, with the specific aim of being the target where magnetic field lines can impinge on. These latter are called open field lines because the magnetic flux surface where they lie on close within the target. The last-closed field surface (LCFS), also called separatrix, separates the confined main plasma from the scrape-off layer (SOL), which is the boundary plasma. Usually, the distance between the toroidal magnetic axis and the separatrix corresponds to the minor radius a.

The most exploited target designs are the toroidal limiter and the divertor.

The toroidal limiter is a solid component tangent to the LCFS and extending all along the toroidal direction, as shown in 1.4. In this case, we say that the plasma is in limiter configuration.

𝜃 𝜑

Core plasma

Separatrix

Scrape-Off

Layer Limiter Nowadays, the divertor configuration, shown in figure 1.5, is preferred to the limiter one in most cases, since in the latter, a large amount of sputtered material from the limiter can reach the main plasma and cool it down, compromising the performance of the machine. Indeed, the divertor configuration allows for more distance between the plasma core and the solid targets, as one can see in figure 1.5. An additional current is driven in the toroidal direction to create a magnetic topology which in the poloidal plane is characterized by an X-point. At the X-point, the poloidal magnetic field is null and the safety factor q is theoretically infinite. In this configuration, the SOL region above the X-point is usually referred to as the main scrape-off layer while the one below is the divertor leg. Also, the private flux region is defined as a thin layer of plasma sustained by the transport of particles and heat from the main plasma. When just one X-point is induced, the magnetic equilibrium is indicated as single null and there are two divertor legs, as in figure 1.5. However, more exotic topologies can be induced, such as the double-null (with two X-points) and the snowflake (multiple X-points). In this thesis, the single-null diverted configuration is especially concerned.

Close flux surface

The sheath

The SOL plasma interacts with the solid target through the so-called sheath, a non-neutral region where the plasma potential drops down. The sheath forms because after the ignition of the plasma (within a few µs), electrons accelerate along the magnetic field lines and, due to their high mobility, they charge the wall negatively. Therefore, a thin sheath builds up to shield the bulk electrostatic potential, as shown in figure 1.6. This sheath is characterized by a positive electric field where ions are accelerated towards the wall. On the other hand, incoming electrons with kinetic energy lower than the potential barrier are repelled, implying that n e < n i . The equilibrium is reached when the electron and ion fluxes equal, namely n e v ,e = n i v ,i . The filtering effect of the sheath on electron energy is taken into account in the boundary condition on the energy fluxes are equal at the entrance of the sheath, which are convective:

Γ E,se, ,i = γ i n se,i T se,i v se, ,i + 1 2 m i n se,i v 3 se, ,i (1.9)
Γ E,se, ,e = γ e n se,e T se,e v se, ,e (1.10) γ i and γ e are called sheath transmission coefficients and weight the energy at the sheath according to the species. The estimation of both can be read in [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]. Usually, γ i = 5/2 and γ e = 5 are assumed.

The thickness of the sheath is in the order of magnitude of the Debye length:

λ D = 0 k B T e n e q 2 e
(1.11) where 0 is the vacuum permittivity, n e the electron density and q e the electron charge. For T e = 20eV and n e = 10 19 m -3 , one obtains λ D ∼ 10 -5 m.

Single particle motion

Sheath

In the sheath the electrostatic potential is provided by the Poisson equation:

∇ V = - q e 0 (n i -n e ) (1.12)
while the electrons can be described by the Boltzmann equation:

n ,e = n se exp(q e (V -V se )/k B T e ) (1.13)
where n se and V se are density and potential at the sheath limit.

The physics of the sheath roughly outlined here occurs at too small scale to be described by fluid codes, but is the base on which boundary conditions are defined.

Single particle motion

In this section, the motion of a single charge particle is investigated when it interacts with the electro-magnetic field.

In general, the force balance of such a particle (in non-relativistic cases) reads:

m dv ∂t = q E + v × B tot (1.14)
where B tot = B = B(x, t), E = E(x, t) and q and m are the charge and the mass respectively. The projection of 1.14 on the direction parallel ( ) to the total magnetic field B provides:

m dv dt = qE (1.15)
which states that the acceleration of the particle along B is governed by the parallel component of the electric field and its orientation is provided by the sign. On the other hand, the projection of equation 1.14 on the plane perpendicular (⊥) to the magnetic field lines, leads to different kinds of motion depending on both E and B.

• When E = 0 and B is uniform in space and constant in time, one recovers the gyro-motion in the ⊥ plane as the zero-order expansion of 1.14 around the quantity ω/ω c , where ω is the frequency of the motion and ω c is the gyro-frequency:

ω c = v ⊥ ρ L = |q|B m (1.16)
Indeed, one can assume that ω/ω c → 0 which implies a scale separation between the gyro-motion and the drift motion. This is called adiabatic assumption. In this case, the guiding center of such a particle moves in circles of radius ρ L .

• Assuming homogeneous magnetic and electric fields, the guiding center undergoes a drift of the velocity given by:

v (0) d = v E×B = E × B B 2 (1.17)
The E × B velocity is perpendicular to both the magnetic and electric fields and does not depend on the particle charge, thus it does not imply the formation of a net current.

• Finally, when B is not homogeneous as in equation 1.4, further drift mechanisms come into play. Indeed, according to the adiabatic assumption, when expanding 1.14 to the first order, one obtains for the drift velocity:

v (1) d = v E×B + v ∇B + v R = E × B B 2 + + 1 2 q |q| ρ L v ⊥ B × ∇B B 2 + m qB 2 v 2 R c × B R 2
where the second term is called ∇B drift velocity, while the third one is the curvature drift (R c is the curvature vector, defined from the axis of symmetry of the torus up to the field line). Both depend on the particle charge. Moreover, as shown in figure 1.7, v ∇B drifts the particle vertically, implying the generation of a vertical electric field which in turn causes an horizontal shift provided by the E × B drift.

Eventually, in the most general form, the velocity of a plasma particle has the structure below:

v = v B B + v g + v d (1.18)
where v g is the high-frequency gyro-motion of the charged particle.

Plasma transport

As mentioned above, charged particles undergo two kinds of motion as they are confined: a gyro-motion perpendicular to the magnetic field lines, and a free motion along these latter. However, further mechanisms come into play in toroidal plasmas that must be taken into account to assure the confinement over long periods of time. We refer to the mutual particle interactions, taken into account by the so-called classical description of transport and particle drifts, due to the interaction of the plasma with the electro-magnetic field. These latter are concerned in the neoclassical and turbulent description of transport.

Classical transport

Within the plasma, charged particles interact via Coulomb collisions. The displacement of their guiding centers can lead them to the neighboring magnetic flux surfaces, allowing for a local (i.e. short-range) transport of matter and energy. These interactions are random and isotropic, thus the transport can be modeled by diffusion, where the length scale is provided by the Larmor radius while the time scale is determined by the collision frequency ν C . ν C varies as nT -1.5 , ρ L as T 0.5 , thus the overall dependence of the diffusion coefficient

D C is T -0.5 . D C = ν C ρ 2 L (1.19)
For pure deuterium plasma, D C is of the order of 10 -2 -10 -4 m 2 /s, much lower than the values observed experimentally.

Neoclassical transport

When taking into account the profile of the magnetic field B ∝ 1/R, the neoclassical description of transport is invoked. This latter includes the effects of the plasma's toroidal shape to the diffusive transport of particles and energy due to collisions. Two main regimes are identified depending on the collisionality

ν * = ν C /(v ⊥ /qR) = λ C /L , where v ⊥ ∼ v th . L = 2πqR
is called connection length and it esteems of the total length of a magnetic field line. When ν * 1, the Pfirsch-Schlüter regime is concerned, characterized by small mean free path λ = v ⊥ /ν C , since the motion of a typical thermal particle is disturbed by collisions before an orbit has been completed. This is the usual picture of boundary plasma. The diffusion coefficient can be expressed as:

D P S = νλ 2 (1.20)
On the contrary, when ν * 1, we have the so-called banana regime, where particles are mostly trapped into banana orbits. The plasma core is usually in such a regime. In this case, the diffusion coefficient (here for electrons) can be estimated as in 1.21 ([4]):

D b ∼ q 2 3/2 ν e ρ 2
L,e (1.21) where = R/a is the aspect ratio, ν e is the electron collisional rate, q the safety factor and ρ L,e the Larmor radius for electrons. The magnitude of diffusive transport in the neoclassical framework is usually in the order of 10 -1 m 2 /s, at least one order of magnitude lower than experimental observations (see for example [START_REF] Wagner | Transport in toroidal devices-the experimentalist's view[END_REF] and [START_REF] Chankin | SOLPS modelling of asdex upgrade h-mode plasma[END_REF]).

Turbulent transport

For this reason, for quite some time the transport of particles and energy across the magnetic flux surfaces has been referred to as anomalous. Nowadays, it is well established within the plasma community that the most responsible for the discrepancy between neoclassical diffusivity and experimental observation is turbulence. Turbulence is a state of fluids where sudden changes in pressure and flow velocity occur, involving several temporal and spatial scales. In plasmas, it is caused by the interplay between plasma transport and self-consistent electro-magnetic fields due to linear and non-linear saturation of different instabilities ( [START_REF] Reiser | Introduction to drift wave turbulence models[END_REF]). The drive can also be triggered by geometrical effects (see for example [START_REF] Kinsey | The effect of plasma shaping on turbulent transport and e×b shear quenching in nonlinear gyrokinetic simulations[END_REF]). Due to its relevance in this work, a particular instability is outlined below, called the interchange instability.

Interchange instability

The curvature of the magnetic field allows plasma ions and electrons to experience the drifts discussed in 1.3. In analogy with the fluid Rayleigh-Bénard instability -where aligned buoyancy force and temperature stratification trigger turbulent flows above a threshold -in magnetized plasma an instability arises when ∇B is aligned with ∇p . This is called interchange instability and it is mostly present at the LFS of tokamaks, where ∇B • ∇p > 0. As a consequence, the transport of plasma results enhanced in such a region rather than on the HFS. This poloidal asymmetry is referred to as the ballooning effect, observed experimentally both in limiter ( [START_REF] Fedorczak | Experimental investigation of turbulent transport at the edge of a tokamak plasma[END_REF]) and diverted plasmas ( [START_REF] Smick | Transport and drift-driven plasma flow components in the alcator c-mod boundary plasma[END_REF]).

The word interchange refers to the exchange of energy among different flux tubes whenever the system departs from thermodynamical equilibrium in regions characterized by inhomogeneous magnetic field. This mechanism is shown in figure 1.8, left panel. A more detailed description is provided in the right panel, where the outer midplane region (i.e. the equatorial plane) is concerned.

Here indeed, a poloidal perturbation of the potential generates small convection cells with iso-potential contour line (in yellow in the right panel of figure 1.8), leading to local electric fields and thus radial E × B drift perpendicularly to the magnetic field. The ∇B (and the curvature) drift velocity are directed vertically and have opposite sign for ions and electrons. The combination of such drifts enable the mixing of the inner, warmer plasma, having a higher ∇B velocity, with the outer plasma, colder and relatively slower. Moreover, it fosters a charge separation leading to a perturbation of the plasma density distribution poloidally.

Assuming a sinusoidal shape of the potential and density distributions along the poloidal direction, the interchange instability occurs when the phases are in quadrature, namely when the phase shift is ∼ π/2. The expression of the linear interchange growth rate γ I is provided analytically in the framework of the linear stability analysis of a reduced model for turbulent plasma (see for example [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF], [START_REF] Cartier-Michaud | Cinétique Fluide Plasma Plasmas Plasmas, Théorie cinétique des Réduction de modèles Turbulence Vérification de codes Événements[END_REF]), resulting in the following expression:

γ I = p/ρ L R c p/∇p (1.22)
where R c is the curvature of the magnetic field and ρ L is the Larmor radius. Equation 1.22 can be arranged observing that R c ∼ B/∇B and p/ρ L ∼ c 2 s , thus becoming:

γ I = c s R • γ 0 R 2 ∇p i • ∇B ϕ p i B ϕ (1.23)

Turbulent transport in edge plasma for magnetic fusion 13

The amplitude coefficient γ 0 is usually set to 1.

Enhanced confinement

Instabilities-driven turbulent fluxes can be lowered significantly when the confinement of the plasma is improved. In 1982, the German tokamak ASDEX reached for the first time an High Confinement regime (called H-mode) in divertor configuration ( [START_REF] Wagner | Regime of improved confinement and high beta in neutral-beam-heated divertor discharges of the asdex tokamak[END_REF]). When transiting from the Low-Confinement (L-mode) to the H-mode, the plasma pressure gradients is observed to enhance and steep in the edge region, forming a pedestal and causing the shrinking of the SOL width at the outer midplane, extensively observed experimentally (see for instance [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF] and [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF]).

Modeling via the transport code SolEdge2D-EIRENE

As discussed in the previous sections, magnetically confined plasma is a complex system where several physical phenomena come into play and interact with each other: the magnetic confinement, the wall geometry, the transport of plasma and heat along both the magnetic field and the perpendicular direction, as much as the plasma-wall interaction and the transport of impurities within the vessel. Numerical simulations on such a complex system within reasonable computational time can only be afforded by so-called transport codes, relying on twodimensional fluid equations for main plasma quantities (for instance SOLPS [START_REF] Schneider | Plasma edge physics with b2-eirene[END_REF] or EDGE2D [START_REF] Simonini | Models and numerics in the multi-fluid 2-d edge plasma code edge2d/u[END_REF]). In this work, we have used the code called SolEdge2D-EIRENE.

SolEdge2D-EIRENE is a transport code package developed at CEA Cadarache, in France, in collaboration with the M2P2 laboratory of the Aix-Marseille University ( [START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code soledge2d-eirene[END_REF], [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF], [START_REF] Reiter | The eirene and b2-eirene codes[END_REF]). It was designed with the purpose of simulating the edge and the scrape-off layer of plasma in a toroidally axysimmetric spatial domain including realistic wall geometry and plasma-wall interactions. This latter is treated by coupling SolEdge2D with the Monte Carlo code EIRENE ( [START_REF] Reiter | The eirene and b2-eirene codes[END_REF]). Accounting for neutrals (both atoms and molecules) is crucial since they drive the source of particles, momentum and energy in the SOL.

Assumptions and equations

The model is derived from the fluid approach of Braginskii ([21]), which is assumed to be valid for high values of the collisionality, ν * = L /λ C , where the electron mean free path λ C is proportional to T 2 e /n e according to [START_REF] Richardson | NRL PLASMA FORMULARY[END_REF]. For a tokamak with major radius R = 2m and safety factor q = 4, when assuming T e = 100eV and n e = 10 19 m -3 , ν * ∼ 5. This result implies that the plasma is actually marginally collisional, and indeed kinetic corrections can be designed for the Braginskii equations. In the following, these latter are not taken into account. In particular, the description of the parallel transport follows a collisional approach, while the cross-field transport is usually modeled as a diffusive term.

Modeling via the transport code SolEdge2D-EIRENE

SolEdge2D solves the equations of n (density), T (temperature), v (velocity component aligned to the magnetic field) for both electrons and ions, assuming quasi-neutrality (n e = n i ) and ambipolarity (v ,e = v ,i ), which can be relaxed including an equation for the electric potential. The model is presented in detail in [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF] (equations 1 -4). n, T and v are the moments of the averaged kinetic distribution function.

Below, the density, momentum and energy equations are reported. The definition of the different terms can be found in [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF].

∂n ∂t + ∇ (nv ) = ∇ • (D n ∇ ⊥ n) + S n (1.24a) ∂ ∂t (m i nv ) + ∇ • (m i nv v) = -∇ p i + q e nE + R ei +∇ • ν ∇ v B B + m i nν∇ ⊥ v + S m,i (1.24b) ∂E i ∂t + ∇ • (E i v i + p i v i ) = q i n i v E + v R ei +∇ • κ i ∇ T i B B + nχ i ∇ ⊥ T i +∇ • ν v B B + m i n∇ ⊥ v +Q ei + S E,i (1.
In this system of equations, perpendicular transport terms rely on transport coefficients D n , χ e , χ i , ν. The main concern of this work, as discussed in section 1.6, is to develop a new modeling of these quantities to improve the predictive capabilities of the code.

In 1.9(a) and 1.9(b) examples of the code output are reported below for a WEST (W Environment in Steady-state Tokamak) test-case.

Boundary conditions

The boundary conditions of the code are derived from the physics of the sheath. In particular, at the entrance of the pre-sheath, the absolute value of plasma parallel velocity along the magnetic field lines reaches minimum sound speed. Indeed, due to the strong sink located at the wall and consequent loss of particles, the plasma behaves like a fluid which expands into vacuum after crossing the sheath entrance. This is the Bohm condition for the plasma parallel velocity, expressed by equation 1.25. The Bohm condition is used in SolEdge2D as immersed boundary condition in the penalization technique.

|v se, ,i | ≥ c s = q e (T e + T i ) m i (1.25)
Moreover, the parallel heat flux impinging on the solid target q ,BC reads:

q ,BC = (γnT v ) BC (1.26)
where γ is the total sheath transmission coefficient.

Numerical treatment

Finite volume for the spatial discretization and implicit/explicit Eulerian for the time integration are the numerical schemes exploited by SolEdge2D. The code can handle elaborate magnetic configurations thanks to its explicit domain decomposition technique, consisting in numerically solving each sub-domain in parallel. In particular, for numerical efficiency SolEdge2D generates a grid aligned to magnetic flux surfaces. Magnetic measurements at different locations surrounding the vacuum vessel are used as real-time inputs to achieve the numerical reconstruction of the plasma current density and the magnetic equilibrium (see for example [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF]).

To model complex and realistic vessel walls, a penalization technique has been implemented ( [START_REF] Paredes | A penalization technique to model plasma facing components in a tokamak with temperature variations[END_REF], [START_REF] Isoardi | Penalization modeling of a limiter in the tokamak edge plasma[END_REF]), inspired by the Computation Fluid Dynamics community. Its aim is to extend the mesh to the wall of the machine and apply a mask function which defines the cells to be solved (corresponding to the plasma domain) and the ones related to the solid material.

Perpendicular turbulent transport in SolEdge2D

The approach adopted by SolEdge2D to treat cross-field transport, as most transport codes for tokamak plasma, is assuming that the turbulent contribution to particle, momentum and energy flux is modeled by diffusion. As an example, we will consider the Braginskii equation for the density, however the same reasoning can also be applied to the momentum equation 1.24b and the energy equations 1.24c and 1.24d. We call Γ n,⊥ the perpendicular particle flux, which in general is defined by 1.27, where v ⊥ is the perpendicular velocity.

Γ n,⊥ = nv ⊥ (1.27)
The Braginskii equation for the density therefore reads:

∂n ∂t + ∇ (nv ) + ∇ • Γ n,⊥ = S n (1.28)
Γ n,⊥ models plasma turbulent phenomena occurring in the cross-field direction. Since transport codes cannot simulate high-frequency fluctuations of the flow, we average the continuity equation 1.24a, and so Γ n,⊥ , with respect to such frequencies. This is similar to the Reynolds average used for the simulation of turbulence in neutral fluids, as we will see in Chapter 2. The result is provided in equation 1.29:

Γ n,⊥ = n v ⊥ + ñṽ ⊥ (1.29)
where is the average operator and ñ and ṽ are the fluctuating component of the respective quantities. Here the vectorial nomenclature has been omitted for simplicity.

While v ⊥ can be related to the pinch velocity of plasma, ñṽ ⊥ is the transport of n due to high-frequency turbulence, which has to be modeled. A typical closure is to relate it to averaged quantities through a viscosity coefficient whose value has to be closed, as done in 1.30. Here, that coefficient is D n , while in the energy equation for electrons and ions the diffusion coefficients are χ e and χ i respectively (see equation 1.24). Note that when replacing 1.30 in equation 1.28, one recovers 1.24a.

Γ n,⊥ ∼ -D n ∇ ⊥ n (1.30)
In many cases, in transport simulations, density diffusion (or transport) coefficient is assumed constant and set to 1m 2 /s as a standard, while the work presented in this manuscript aims at providing a finer description of the cross-field transport.

The approach expressed by 1.30 is close to the Boussinesq assumption, adopted in the framework of the one-point closure models for neutral fluids, as discussed in Chapter 2.

Also, the approach expressed by 1.30 allows to estimate the characteristic lengths of the system, in particular the SOL width which, to a first approximation, can be determined from the equilibrium between parallel and perpendicular transport, as showed in 1.31. This expression can be obtained from equation 1.24a at steady state and when assuming S n = 0 (consistently with the so-called simple-SOL description of the SOL presented in [START_REF] Stangeby | The Plasma Boundary of Magnetic Fusion Devices[END_REF]).

1 L L 0 ∇ (nv )dl = ∇ • (D n ∇ ⊥ n) (1.31)
where D n is taken constant. Frequently, in cylindric coordinates, the SOL radial extension is smaller than the minor radius a, thus, to estimate its width, toroidal effects can be neglected and a plane geometry can be assumed. The first term on the left-hand side of equation 1.31 can be integrated recalling the Bohm condition for plasma at the solid target v = ±c s where the sign is opposite at l = 0 and at l = L :

1 L L 0 ∇ (nv )dl ∼ 2nc s L (1.32)
where as working assumption we set the average parallel velocity v ∼ c s . Finally, we perform a dimensional analysis on equation 1.32, namely we rewrite each term of the latter as a function of the reference dimensional quantities:

2nc s L = D n n λ 2 n (1.33)
Since here we assume that λ n ∼ λ SOL , an expression for the latter is readily obtained:

λ n = D n L 2c s (1.34)

Automatic fitting procedure

Adopting D n = 1m 2 /s, T = 100eV , L = 10m and deuterium plasma, one obtains λ n ∼ 10mm, much smaller than the size of the machine.

By the same reasoning on pressure (equation 1.24d), one finds:

λ p = χ e L 2γ e c s (1.35)
The e-folding length of the cross-field heat flux λ q in the SOL, on the other hand, is another useful estimation of the SOL width, since it allows to evaluate the surface region where the heat load impinges on. Therefore, in SolEdge2D simulations discussed in this thesis we refer to λ q as estimation fo λ SOL , or, alternatively, to the pressure e-folding length in the SOL (especially addressed in the framework of dimensional analysis).

Constant or ad-hoc diffusion coefficient, which can be easily exploited in transport codes, in most cases cannot provide but a rough description of the cross-field transport of plasma within the tokamak. Their distribution can be however traced basing on the data set available for a given experiment. This is the rationale behind the automatic fitting procedure outlined in the next section.

Automatic fitting procedure

In the current version of transport codes like SolEdge2D, SOLPS etc., ad-hoc effective transport coefficients are set as inputs, compared to experimental ones and adjusted with iterative schemes ( [START_REF] Chankin | SOLPS modelling of asdex upgrade h-mode plasma[END_REF]).

To address the issue of cross-field fluxes in transport codes, a more sophisticated fitting procedure can be invoked. Instead of providing transport coefficients as an input to the code, one directly provides midplane profiles of density and temperature. The code will automatically adjust the transport coefficients so as to match in the end the simulated midplane profiles with the experimental ones. More precisely, one uses the Proportional-Integral feedback loop described on Figure 1.11 to correct on the fly the transport coefficients while the temporal loop of the code is running.

The gain of the feedback loop τ i has been adjusted to ensure the code stability, while the difference of signals is given by the following formula:

≡ D ⊥ -D ⊥,0 = 1 τ i ∂ r n XP -∂ r n simu )dt + K(∂ r n XP -∂ r n simu ) (1.36)
where ∂ r n is the density gradient and the subscripts XP and simu stand for experimental and simulation respectively. Typically, the time constant associated with the integral correction must be larger than the parallel time τ = L /c s , with τ ∼ 10 -5 [s].

Using this feedback loop automatic procedure, the way of running SolEdge2D is made simple for experiments interpretation. As an input of the code, one directly The scheme of the feedback loop procedure implemented in SolEdge2D-EIRENE is represented as a control system.∂ r n Xp , provided as input, is the experimental profile while ∂ r n simu is the one provided by the simulation. is the error between these two. K and 1 τ i are the transfer functions of this scheme and finally D ⊥ is the output.

provides experimental measurements of density and temperature in the midplane.

As an output, one still obtains a 2D poloidal map of density and temperature but one additionally gets a 1D profile of the cross-field diffusivities that were adhoc before. The poloidal variation of these coefficients is, though, not taken into account. Since one sets 1D radial profiles of density and temperature at the outer midplane, one gets 1D radial profiles of diffusivities. These latter are assumed to be homogeneous on each flux surface. Once the transport coefficients are set by the feedback procedure, one knows that experimental and simulated midplane quantities are in rather good agreement.

Figure 1.12 shows the radial distributions of diffusivity and thermal conductivity for ions and electrons calculated with SolEdge2D-EIRENE by fitting Thomson scattering data for density and temperature at outer midplane, in the H-mode low-density ASDEX Upgrade shot discussed in [START_REF] Chankin | SOLPS modelling of asdex upgrade h-mode plasma[END_REF].

From this figure, it is manifest that diffusive transport coefficients have specific spatial distribution which must be taken into account for a proper estimation of plasma density and heat load in the poloidal plane.

Therefore, the next chapters are devoted to the discussion on possible strategies to model self-consistently transport coefficients and improve the predictive capabilities of transport codes concerning the estimation of plasma perpendicular fluxes. 20 1.7. Automatic fitting procedure Figure 1.12: Radial profiles of diffusivity D perp , in red color, and thermal conductivity for ions and electrons χ i and χ e , in blue and black color respectively. These profiles have been calculated by SolEdge2D-EIRENE by fitting experimental density and temperature data published in [START_REF] Chankin | SOLPS modelling of asdex upgrade h-mode plasma[END_REF].

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids

Introduction to turbulence concepts in neutral fluids

The aim of this chapter is to introduce the key features of neutral fluid turbulence theory which inspired the work presented in this manuscript. In particular, we introduce the concept of eddy viscosity and characterize one of the most popular models in the framework of Computational Fluid Dynamics (CFD): the κ-ε model.

Introduction to turbulence concepts in neutral fluids

Turbulence is an ubiquitous phenomenon in nature. It can be experienced every day, from the mixing of milk and coffee in a cappuccino to the annoying bumps occurring during a flight. It usually dominates all other flow phenomena and results in increased energy dissipation, mixing, heat transfer and drag. As a general concept, turbulence is a spatially and temporally complex regime of fluids characterized by random, unsteady nonlinear fluctuations due to the advection of a three-dimensional quantity called vorticity, defined as the curl of the velocity field u: ω ≡ ∇ × u(see appendix A). Turbulence arises as instability of the laminar flow, resulting from the competition between the non-linear inertial terms (destabilizing) and the linear (in Newtonian fluids) viscous term (stabilizing) contained in the Navier-Stokes equation, recalled below for an incompressible flow:

ρ ∂u i ∂t + ρ u j ∂ ∂x j u i = - ∂p ∂x i + ∂τ ij ∂x j -αu i + F i (2.1a
)

∂u i ∂x i = 0 (2.1b)
Here a tensorial notation has been used. The left-hand side is the material derivative of u i , where ρ ≡ mn is the volumetric mass density [kg/m 3 ]. p is the hydrostatic pressure while the viscous stress tensor τ ij accounts for (normal and tangential) stresses associated with dynamic viscosity µ[kg/ms] :

τ ij ≡ µ ∂u i ∂x j + ∂u j ∂x i (2.2)
Finally, -αu i is a large-scale drag force governed by an α coefficient and F i is a generic external force.

It is conventional to introduce the tensorial quantity called rate-of-strain tensor S ij , defined as:

S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i (2.3)
Thus a most compact form of the stress tensor reads:

τ ij = 2ρνS ij (2.4)
In figure 2.1, a cartoon of a sheared fluid element is shown to explain the origin of the stress tensor.

The ratio between inertial and viscous terms is called Reynolds number, defined as:

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 23

Re ≡ ρuL µ = uL ν (2.5)
where L is the system linear dimension in [m] and ν[m 2 /s] is the kinematic viscosity characteristic of the fluid. High values of such a number, depending on the given scale and on the fluid, imply unstable and eventually turbulent flow.

Energy cascade and Kolmogorov theory for isotropic homogeneous turbulence

Fluids with high values of Re are characterized by whirling structures (or eddies) of widely varying size and velocity, whose organization can be investigated by observing the turbulent energy spectrum of the flow. This spectrum displays in the Fourier's space the energy exchange occurring among eddies with different scales and it is usually referred to as energy cascade.

Since the concept of cascade is relevant further in this thesis, especially in the design of a new modeling for the cross-field fluxes in SolEdge2D-EIRENE, it is worth to present it in detail.

The first pioneering work on the derivation of the power law for the energy spectrum was made by Kolmogorov in 1941 in the framework of steady-state isotropic homogeneous three-dimensional turbulence. His theory relies on dimensional analysis to estimate the effective rate of energy transfer between different turbulent scales, depending on the characteristic scales of the system. Most important, this theory predicts the existence of a finite-size region in the space scales, called inertial range, where the mechanism of energy transfer does not change (this feature is called self-similarity) and does not depend on the viscosity.

The concept of scale hierarchy underlying the definition of energy cascade, on the other hand, was proposed by Richardson in 1922 for fully-developed, steady-state turbulence. According to this concept, the largest turbulent eddies, created by instabilities in the mean flow, are themselves subject to inertial instabilities and rapidly break-up to evolve into smaller swirls. Smallest eddies are in turn unstable and pass their energy and momentum onto even smaller structures, and so on. An energy cascade is thus established between neighing turbulent scales, evoking the idea of a cascading waterfall.

The direction of the energy cascade can actually be either direct (dominant in three-dimensional treatment of turbulence) and inverse. This distinction is related to the inviscid (ν = 0) invariants of the system. In hydrodynamical systems, an invariant is a quantity determined from the fluid variables (and possibly their spatial derivatives), advected by the flow and which holds "frozen" into this latter.

In three-dimensional turbulence, two global invariants are defined: the (kinetic) energy

E(t) = 1 2 u • u and the helicity H(t) = 1 2 u
• ω , which is a measure of linkage and/or knotting of vortex lines in the flow.

In two-dimensional turbulence, energy and enstrophy, defined as Z(t) = 1 2 w •w , are conserved, where Z is a measure of the strength of the vorticity field.

Direct cascade

The direct cascade is typical of three-dimensional turbulence, where helicity H(t) and energy E(t) are conserved. In particular, the conservation of energy leads the simple relation:

E in = E ν (t) + E α (t) (2.6)
where E in is the energy injected at a given scale k in ∝ -1 in , being the characteristic size of the eddy. E ν (t) is the energy dissipated by small (or "Kolmogorov") scales ( ν ), while E α is dissipated by drag mechanisms involving largest eddies. Equation 2.6 leads to the split of the spectrum displayed in figure 2.2. Here, an instability-driven source is injected in k in , the energy is transferred without dissipation in the inertial range and, finally, dissipated at the smallest length scales (k ≥ k ν ). On the other hand, the large length scales characterized by k < k in reach a statistical (thermal) equilibrium due to local and non-local energy diffusion across the wavenumbers, where all Fourier modes undergo a thermal bath described by a Gibbs-ensemble equipartition distribution E(k) ∝ k 2 (see [START_REF] Michel | Observation of thermal equilibrium in capillary wave turbulence[END_REF] for reference). Let us characterize a turbulent scale by its typical length , velocity u and time t = /u . Let κ = 1 2 u 2 be the turbulent energy per unit mass of scale , while ε is the power per unit mass coupled to the system, assumed to be scale-independent due to the conservation law, so that:

κ = t ε (2.7)
Given the definition of κ and (2.7), one obtains:

ε ∝ κ t ∝ u 3 (2.8)
Thus as the scale is reduced, the characteristic velocity decreases like 1/3 , and κ reads:

κ ∝ ε u ∝ ( ε) 1/3 = (ε ) 2/3
(2.9)

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 25

Thermal Equilibrium

E k ∝ 1 k ε k 2/3 = ε 2/3 k -5/3 (2.10)
This is the Kolmogorov scaling law of the spectral energy density as a power law in terms of the wavenumber k, with exponent equal to -5/3, where κ /t ∝ does not depend on . Using dimensional analysis, one can prove that ν /L ∝ Re -3/4 , where L = κ L ε is the integral scale and

ν = ν 3 ε 1/4 is the Kolmogorov scale,
showing that for high Reynolds numbers (and therefore for turbulent flows) the range of turbulent scales concerned widens.

For simplicity, we assume zero injection of helicity in the analysis above, however all experimental investigations, numerical simulations and phenomenological theories indicate that both energy and helicity have a simultaneous mean transfer to the small scales [START_REF] Alexakis | Cascades and transitions in turbulent flows[END_REF].

Inverse cascade

In the inverse cascade, typical of two-dimensional turbulence, both kinetic energy and enstrophy are conserved and a double cascade scenario can be outlined ( [START_REF] Robert | Inertial ranges in two-dimensional turbulence[END_REF]). The former undergoes an inverse cascade towards large scales, while the latter experiences a direct cascade to small scales, as sketched in figure 2.3.

This scenario implies the existence of two inertial ranges (for both energy and enstrophy) where the effects of the viscosity and the external forces are negligible. A forcing injection is localized in k in . The injected energy is then transferred to ever larger scales k α < k in in a quasi-stationary fashion. In the energy inertial range, E(k) depends only on the energy dissipation rate ε and the wavenumber k according to the power law E(k) ∝ ε 2/3 k -5/3 , retrieved both experimentally ( [START_REF] Shats | Spectral condensation of turbulence in plasmas and fluids and its role in low-to-high phase transitions in toroidal plasma[END_REF]) and numerically ( [START_REF] Smith | Bose condensation and small-scale structure generation in a random force driven 2d turbulence[END_REF]).

On the other hand, Z(t) is transferred to smaller scales through the enstrophy inertial range governed by the enstrophy dissipation rate η related to molecular viscosity. The power law proposed by Kraichnan in [START_REF] Robert | Inertial ranges in two-dimensional turbulence[END_REF] on the basis of dimensional analysis is Z(t) ∝ η 2/3 k -3 , numerically retrieved in [START_REF] Boffetta | Evidence for the double cascade scenario in twodimensional turbulence[END_REF] for Re = 32768 2 .

Both for direct and inverse cascades defined for neutral fluids, the interaction between scales is local, meaning that the transfer of energy occurs among neighbouring wave numbers.

Statistical approach and averaging techniques

We have learned that in turbulent flows, where the Reynolds number is sufficiently high, the range of scales describing the size of the eddies is wide, since ν /L ∝ Re -3/4 . A statistical approach can therefore be adopted aiming at modeling part of the scale motion instead of solving them directly. Such a model reduction can be achieved through different averaging techniques. In general, an averaging operator filters the variable f which therefore can be separated into mean and fluctuating components, f and f respectively:f = f + f .

• The ensemble average is one of the most commonly used in fluid mechanics.

Provided a number N of experiments, the averaging operator is defined as:

f ≡ 1 N N i=1 f i (2.11)
f i being the value of the variable for each experience.

• The time average is exploited for statistically stationary flows. In this case the averaged variable reads:

f ≡ lim T →∞ 1 T t+T t f (t )dt (2.12)

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 27

This averaging technique is used in the RANS numerical approach for turbulent flows, as we will see in the next section.

Starting from the second half of the 20 th century, the Favre density-weighted averaging has been introduced in the fluid mechanics for flows with varying density [START_REF] Favre | Statistical equations of turbulent gases in problems of hydrodynamic and continuum mechanics[END_REF]. The Favre-averaged component of the instantaneous flow variable f (x, t) is defined as f = ρf ρ , where ρ is the volumetric mass density and f and ρ are time-averaged.

In the present thesis, this approach has not been used despite its relevance with plasma flows.

Numerical approach for turbulent flows

The numerical resolution of the whole turbulent energy spectrum, achieved in the Direct Numerical Simulations (DNS) [START_REF] Orszag | Analytical theories of turbulence[END_REF], is usually too costly, therefore the Reynolds-averaged Navier-Stokes (RANS) statistical approach, based on the timeaveraging of the whole spectrum, is often exploited. An intermediate approach, adopted in the so-called Large Eddy Simulations (LES), consists in filtering only the dissipative scales [START_REF] Smagorinsky | Generale circulation experiments with the primitive equations[END_REF] [35] [START_REF] Fröhlich | Hybrid LES/RANS methods for the simulation of turbulent flows[END_REF]. For the sake of comparison, in figure 2.4 a comparison between solving capabilities in DNS, LES and RANS numerical approaches are showed, while in 2.5 the mixing between two streams with unequal scalar concentrations is numerically reproduced with DNS, LES and RANS simulations respectively [START_REF] Rasheed | Turbulence modeling: Large eddy simulation[END_REF]. Figure 2.4 shows that DNS simulations compute the entire spectrum of turbulent scale motion; LES simulations filer Kolmogorov scales out (which therefore have to be modeled) while solving the rest of the spectrum; finally RANS simulations apply an averaging procedure to the whole spectrum.

In figure 2.5, the comparison of results from DNS, LES and RANS simulations is provided, showing two streams with unequal scalar concentrations. Here, one can appreciate a definite structure of eddies in DNS case, where all the scales are solved, while turbulent structures tend to blur as far as LES and mostly RANS simulations are concerned.

The RANS averaging procedure

Experimental observations and numerical simulations of turbulent flows (extensively reported in many textbooks, such as in [START_REF] Davidson | turbulence[END_REF]) denote a strongly fluctuating behavior in time, as illustrated in figure 2.6, with a broad range in the frequency scale. Therefore, the RANS statistical averaging approach is often used. In the RANS approach, first introduced by Reynolds ( [START_REF] Reynolds | On the dynamical theory of incompressible viscous fluids and the determination of the criterion[END_REF]), the time average for the instantaneous flow variable f (x, t) is defined as:

f (x, t) ≡ lim T →∞ 1 T t+T t f (x, t )dt (2.13)

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 29

with f (x, t) = f (x, t) + f (x, t), with f (x, t) = 0 by construction. T → ∞ cannot be realized in any physical flow. However, when the time scale of the averaged (or mean) flow, T , is just much larger than the time scale of turbulent fluctuations, the scale separation between oscillating and mean contributions makes sense.

Time-averaging the non linear term leads to higher-order correlations in the equations governing the motion of the fluid, which are unknown a priori. The mean component of the Navier-Stokes equation indeed reads:

ρ ∂u i ∂t + ρ u i ∂ ∂x j u i = - ∂p ∂x i + ∂ ∂x j τ ij -ρu i u j -αu i + F i (2.14)
In equation 2.14, quantities are averaged with respect to turbulent fluctuations. Therefore, ∂u i ∂t does not vanish since the time evolution expressed by the partial derivative is much slower than the typical time scale of fluctuations, due to the scale separation. RS ij = -ρu i u j is the Reynolds stress tensor, which provides the momentum exchange between mean and fluctuating fields. Its explicit formulation

RS ij = -ρ    u x u x u x u y u x u z u y u x u y u y u y u z u z u x u z u y u z u z    (2.15)
shows that each element of RS ij is a quadratic velocity correlation function measuring the mutual, statistical dependence of the velocity components, unknown a priori. We refer to the trace of this tensor as the turbulent kinetic energy κ ≡ 1 2 u 2 i with i = x, y, z.

In the so-called Reynolds stress models (RSM ), transport equations are derived for each of the six distinct Reynolds stress terms by manipulating equation (2.1). However, this procedure lets another higher-order correlation term to arise for each equation. Moreover, it remains relatively expensive and not always straightforward to implement.

When the time evolution of moments of order n for some field f is coupled by the non-linearity to moments of order n + 1, a general approach is to apply a closurebased moment method. This latter is an approximation which provides an expression for the higher-order moments depending on lower-order ones. Schematically, an example of hierarchy composed by equations of increasing order is provided below up to order three:

∂ t f = f f + ν∆f ∂ t f f = f f f + ν∆f f ∂ t f f f = ν∆f f f (2.16)
Here f is the generic flow variable while the last term in the first and second equations is the non-linearity (governed by the dissipation coefficient ν), coupled

The eddy viscosity concept

to the moment of higher order. For the third equation, the closure is applied with a term of the same order, governed by coefficient ν. Most models describing turbulence in neutral fluids adopt this approach, assuming an algebraic relation for the turbulent stresses, which represent the non-linear contribution to be closed (see sections 2.3 and 2.4).

The statistical closure problem is a common feature in the moment-based description of both plasmas and neutral fluids. For instance, advection by the E ×B drift across the magnetic field surfaces is a remarkable example of plasma nonlinearity. In chapter 1 we have seen that transport codes designed for magnetically confined plasmas use the moment-based closure to estimate cross-field fluxes responsible for advection, which indeed are modeled as diffusive terms, namely as a function of lower-order moments.

The eddy viscosity concept

In equation (2.1), the statistical closure problem has been introduced for the Reynolds stress tensor RS ij when the average Navier-Stokes equation is concerned. In 1877, Boussinesq postulated the following shear-stress strain-rate relationship for one-dimensional time-averaged flows:

τ xy + RS xy = ρ(ν + ν t ) ∂u x ∂y (2.17)
where ν t is called turbulent or eddy viscosity. The rationale behind equation (2.17) is that the transfer of momentum caused by advecting fluctuations can be interpreted as an analogous phenomenon to the molecular transport momentum leading to laminar stress τ xy , where the effect of turbulence is to bump up the effective viscosity from ν to ν + ν t , with ν t ν.

The three-dimensional generalization of equation (2.17) is referred to as the Boussinesq's equation. It reads:

RS ij = -ρu i u j = ρν t ∂u i ∂x j + ∂u j ∂x i - 2 3
∂u k ∂x k δ ij - ρ 3 u k u k δ ij (2.18)
where δ ij is the Kronecker delta and ∂u k ∂x k = 0 for incompressible flow. The last term is the contribution of normal stresses, which acts like pressure forces.

Unlike the molecular viscosity ν, the eddy viscosity ν t is a property of the motion and not of the fluid.

The Boussinesq's equation on the Reynolds stress tensor represents an assumption and implies that turbulent and mean flow scales are locally proportional at any point of the flow. However, in principle there is no reason that the eddy viscosity should depend exclusively upon turbulence parameters, such as κ. In general, the ratio RS ij /S ij depends on both mean flow and turbulent scales.

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 31

Turbulence models based on one-point closure

Many turbulence models are based on the Boussinesq eddy viscosity approach (also called one-point closure) to determine the Reynolds stress tensor. The simplest one is a straight estimation of RS ij basing on dimensional arguments.

Prandtl mixing length theory

At the beginning of the 20 th century, Prandtl suggested a way to estimate the Reynolds stress tensor on the basis of dimensional analysis. He was inspired by the success of the kinetic theory of gases in predicting a macroscopic property such as the molecular viscosity from a simple relation with respect to the mean free-path length and the root-mean-square velocity of molecules. Indeed, in laminar flows, layers of fluid slide over each other and experience a mutual drag per unit area, exchanging momentum when they bounce around between the layers. In a similar fashion, in turbulent flows lumps of fluid are thrown around and jostled by turbulence.

Hence, Prandtl introduced a turbulent viscosity ν t = l mix V T , where l mix is the mixing length, namely the typical length where the fluid lump conserves its momentum, and V T is a measure of the amplitude of the fluctuating field. This expression is consistent with the evidence that the more energetic the turbulence, the greater the momentum exchange, and hence the greater the eddy viscosity.

With these assumptions, the Reynolds stress tensor can be modeled (here for a one-dimensional flow) by:

RS xy = ρV T l mix du x dy (2.19)
Furthermore, the following assumption holds:

V T l mix du x dy (2.20)
and thus:

ν t = l mix 2 du x dy (2.21)
The closure of zero-equation models consists in finding a proper prescription of l mix in order to close the equation (2.21). Depending on the geometry of the turbulent system, few values have been proposed (see [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]) in order to provide at least a qualitative prediction of the Reynolds stress tensor.

One-equation models

These models introduce one transport equation for the characteristic turbulent velocity scale V T , which is set proportional to the square root of the kinetic energy

Turbulence models based on one-point closure

transported by fluctuations, κ. The components of the Reynolds stress tensor are thus related to this scale. This choice seems reasonable since the more energetic the turbulence, the greater the momentum exchange, and so the eddy viscosity and the Reynold stress. As already mentioned, κ is the trace of RS ij :

κ ≡ 1 2 m u x u x + u y u y + u z u z (2.22)
An equation for κ can be derived analytically from the Navier-Stokes equations (see for instance [START_REF] George | Lectures in Turbulence for the 21st Century[END_REF] for further details), in the form:

ρ ∂κ ∂t + ρ ∂ ∂x j (u j κ) = RS ij ∂u i ∂x j -µ ∂u i ∂x k ∂u i ∂x k + ∂ ∂x k µ ∂κ ∂x j - 1 2 ρu i u i u j -p u j (2.23)
The left hand-side is the total derivative, composed by the time rate of κ and its advection.

The first term in the right-hand side is a source representing the specific kinetic energy per unit volume that an eddy gains per unit time due to the mean strain rate ∂u i ∂x j . It can be demonstrated that as far as the equation for the kinetic energy of the mean flow is considered, such a term appears as a sink, since the production of turbulence is a result of the mean flow loosing its kinetic energy. RS ij contains the Boussinesq eddy viscosity µ t .

The second term in the right-hand side is referred to as a dissipation, representing the mean rate at which the kinetic energy of the smallest turbulent eddies is transferred to thermal energy at the molecular scale. This term, which is always positive, is also called ε in the literature.

The third term is diffusion allowed by molecular motion, equally responsible for the diffusion of the mean flow momentum. Finally, the last terms refer to the rate at which turbulent energy is transported and redistributed though the flow via turbulent fluctuations.

The eddy viscosity is estimated as a function of some turbulent velocity scale V T , which is set as the square root of κ:

ν t = κ 1 2 (2.24)
This implies that we let a turbulent quantity such the eddy viscosity to scale as a mean quantity, here κ. Furthermore, equation (2.24) is an isotropic relation, since the same amount of momentum is assumed to be transported in all directions from one given point.

Two-equations models

These models, which are the most popular ones in a wide range of engineering analysis and research, provide independent transport equations for both the turbulent kinetic energy scale κ and another characteristic measure of the turbulent flow, depending on

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 33

the distinguishing type of model. They represent complete description of turbulence, in the sense that they take into account both time and space scales of turbulence. However, while the derivation of the κ equation is nearly analytical, a significant amount of uncertainty lies in the complementing transport equation.

In 1942, Kolmogorov indicated the specific dissipation rate ω as the most appropriate unknown to be coupled to κ. This quantity, measured in [s -1 ], can be interpreted as the inverse of the time scale of dissipation of turbulent eddies. Dimensional arguments lead to the following scales:

ν t ∼ ρκ/ω, ∼ κ 1/2 /ω, ε ∼ ωκ (2.25)
Three years later, Chou [START_REF] Chou | On the Velocity Correlations and the Solution of the Equations of Turbulent Fluctuation[END_REF] proposed to use ε, the dissipation rate measured in [m 2 /s 3 ], thus obtaining for the scales:

ν t ∼ ρκ 2 /ε, ∼ κ 3/2 /ε (2.26)
Further options came along in the following years, proposed by Rotta (1951) and Zeierman and Wolfshtein (1986), for instance. In all of them, dissipation, length scale and eddy viscosity are related on the basis of dimensional analysis. The κ -ω and the κ -models are so far the most used to deal with turbulent systems in the engineering framework.

The κ-ω model

The equations proposed for this model in [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] are:

ρ ∂κ ∂t + ρu j ∂κ ∂x j = ∂ ∂x j µ + µ t σ κ ∂κ ∂x j + τ ij ∂u i ∂x j -β * ρκω (2.27a) ρ ∂ω ∂t + ρu j ∂ω ∂x j = ∂ ∂x j µ + µ t σ ω ∂ω ∂x j + α ω κ τ ij ∂u i ∂x j -βρω 2 (2.27b)
where in 2.27b the last two terms are the production and the destruction (or dissipation) terms.

The closure coefficients are obtained here by considering two turbulent flow cases, one reproducing the decaying of homogeneous, isotropic turbulence, the other simulating a boundary layer.

α = 5/9, β = 3/40, β * = 9/100, σ ω = 1/2, σ κ = 1/2 (2.28)
The κ -ω model has been shown to give excellent predictions, within the 6% of DNS (Direct Numerical Simulation) ones, for pipe flows at Re ∼ 40000 [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF].

The κ -ε model

One of the most exploited models in CFD for the treatment of turbulence in neutral fluids is the κ-ε model, first introduced by Launder and Spalding in 1974 [START_REF] Launder | The numerical computation of turbulent flows[END_REF]. In

ε ≡ 2νS ij S ij (2.29)
when the fluid is incompressible, though, 2.29 can be simplified in:

ε ≡ ν ∂u i ∂x j ∂u i ∂x j (2.30)
A transport equation can be derived analytically for ε by performing the operation:

L(•) = ν ∂u i ∂x j ∂ ∂x j (•) (2.31)
The mathematical proceeding is described in detail in [START_REF] Celik | Introductory Turbulence Modeling[END_REF]. The complete model is described by the following system:

ρ ∂κ ∂t + ρu j ∂κ ∂x j = ∂ ∂x j µ + µ t σ κ ∂κ ∂x j + τ ij ∂u i ∂x j -β * ρκω (2.32a) ρ ∂ε ∂t + ρu j ∂ε ∂x j = -2µ u ik u jk + u ki u kj ∂u i ∂x j -2µu k u ij ∂ 2 u i ∂x k ∂x j -2µu ik u im u km -2µνu ikm u ikm + ∂ ∂x j µ ∂ε ∂x j -µu j u im u im -2νp m u jm (2.32b)
Let us focus on the 2.32b. As usual, the left-hand side represents the material derivative of ε. The right-hand side is composed essentially by three distinct contributions:

• The production term:

P ε = -2µ u ik u jk + u ki u kj ∂u i ∂x j -2µu k u ij ∂ 2 u i ∂x k ∂x j (2.

33)

P ε is modeled basing on the assumption of local equilibrium. Given: P ε ∝ P κ /t P , where P κ = τ ij ∂u i ∂x j is the production term of κ and t κε the characteristic time scale for the production of ε, then P ε can be modeled as:

P ε = τ ij ∂u i ∂x j (t κε ) -1 (2.34)
which is equivalent to state that the rate of change of ε must scale with κ, or:

t κε = κ ε (2.35)
thus eventually:

P ε = C ε1 κ ε P κ (2.36)

A brief introduction on the equations and the concepts of one-point closure models in neutral fluids 35

• The sink (or dissipation) term:

S ε = -2µu ik u im u km -2µνu ikm u ikm (2.37)
S ε can be modeled on the basis of dimensional arguments, recalling the definition of t κε :

S ε ∝ ε t κε = C ε2 ε 2 κ (2.38)
• The diffusion term:

D ε = ∂ ∂x j µ ∂ε ∂x j -µu j u im u im -2νp m u jm (2.39)
which can be modeled by means of a properly-scaled eddy viscosity µ t /σ ε :

D ε = ∂ ∂x j µ + µ t σ ε ∂ε ∂x j (2.40)
The combination of all the three modeled therms reads, for the complete κ -ε system:

ρ ∂κ ∂t + ρu j ∂κ ∂x j = ∂ ∂x j µ + µ t σ κ ∂κ ∂x j + τ ij ∂u i ∂x j -β * ρκω (2.41a) ρ ∂ε ∂t + ρu j ∂ε ∂x j = ∂ ∂x j µ + µ t σ ε ∂ε ∂x j + C ε1 ε κ τ ij ∂u i ∂x j -C ε2 ρ ε 2 κ (2.41b)
The eddy viscosity in this case reads:

µ t = ρC µ κ 2 ε (2.42)
Below the closure coefficients are reported. These have been obtained in a systematic fashion by empirical procedures specified in [START_REF] Rodi | Turbulent Buoyant Jets and Plumes[END_REF].

C ε1 = 1.44, C ε2 = 1.92, C µ = 0.09, σ κ = 1., σ ε = 1.3 (2.43)
It has been shown that the κ -ε model performs quite nicely for free shear flows given by jets, wakes and mixing layers [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF]. In general, despite its heuristic closure, it is exploited successfully in a broad range of applications ranging from industrial (notably airplanes design) to environmental flows. The concepts of eddy viscosity and κ-ε model are well-known modeling tools in the theory of neutral turbulence, however, following a new and ambitious approach, they may be thought to be cautiously adapted and exploited also in the framework of turbulent perpendicular transport in magnetically confined plasmas.

In the next chapter, a number of approaches is outlined for the modeling of turbulent transport in plasmas, which can be addressed to customize the concepts of eddy viscosity and κ-ε model.

3

A theoretical work on existing approaches related to the turbulence-driven transport

. The quasi-linear approach

This part of the thesis addresses some theoretical methods concerning the description of turbulence-driven transport in magnetically confined plasmas.

First, the quasi-linear approach is presented as a methodology to relate plasma diffusive cross-field fluxes to the velocity scale of turbulence. A qualitative analysis is also proposed to include the sink of turbulence in such an analysis.

In the second section, prey-predator models are shown to be a successful tool to describe the interplay between two competitive species, especially in transients, so they are found to fit the general picture of plasma turbulence as prey and its sink mechanisms as predators.

Finally, we recall the formalism of the Lorenz equations for chaotic systems, which benefit from several well-known properties, and sketch possible analogies with both prey-predator models in plasma and κ-ε models for neutral fluids.

The quasi-linear approach

In the limit of a electrostatic field characterized by small-amplitude oscillations, the Fourier transform of the fluid density equation (3.1) can be linearized to predict the transport associated to its fluctuations. The starting equation is the following one:

∂ t n + [ψ, n] = -ν ñ (3.1)
The transport of n is taken into account by the non-linear term in Poisson brackets [ψ, n] ≡ z • (∇ψ × ∇n), where z is the unit vector directed as the magnetic field B. Here ψ is the normalized electric potential. On the right-hand side, the linear oscillating restoring force ñ (called Krook's term) is characterized by frequency ν and relates to the effect of collisions on n.

The Poisson brackets approximate the transport of n due to the E × B drift:

∇ • (v E×B n). v E×B ≡ B×∇ψ B 2
is the E × B velocity. The transport of n thus relates:

∇(v E×B n) = ∇ • B × ∇ψ n B 2 = (∇ × B)(∇ψ) n B 2 -B(∇ × ∇ψ) n B 2 + (B × ∇ψ)∇ n B 2 = µ 0 J n B 2 ∇ψ + B • ∇ψ × ∇ n B 2 = µ 0 B 2 J • ∇ψ n + B 2 B B 2 ∇ψ × ∇ n B 2 ∼ ¨¨¨¨¨β -J E p n + B 2 ψ, n B 2 (3.2)
where we used the Maxwell-Ampère law ∇ × B = µ 0 J , µ 0 being the permeability in free space and J the current density. This term is further expanded by using

A theoretical work on existing approaches related to the turbulence-driven transport 39

the definition of β ≡ nT k B B 2 /2µ 0 , called plasma β and indicating the ratio between plasma and magnetic pressure. Therefore, µ 0 = βB 2 /p (p = nT k B ). As mentioned at the beginning of this section, we relate to the electrostatic plasma limit, thus β → 0. Finally, one expands the last term obtaining:

B 2 ψ, n B 2 = [ψ, n] + $ $ $ $ $ $ $ B 2 n ψ, 1 B 2 (3.3)
where the last term is neglected because we suppose that B is nearly uniform in space. This approximation is called straight field approximation.

We introduce the following averaging procedure:

n = n + ñ n = n y,t ñ y,t = 0 (3.4a)
where n = n y stands for the average of n against the magnetic flux surface (y). Let us average Eq.(3.1):

∂ t n + [ψ, n] + [ ψ, ñ] = 0 (3.5)
and subtract Eq.(3.5) from Eq.(3.1) to obtain:

∂ t ñ + [ ψ, n] + [ψ, ñ] + [ ψ, ñ] -[ ψ, ñ] = -ν ñ (3.6)
The quasi-linear approach consists in linearizing equation 3.6 under the assumption of small perturbations, thus [ ψ, ñ] → 0 and [ ψ, ñ] → 0.

First, equation 3.6 is Fourier-transformed to access the dispersion relation of n, namely n, against the spectrum of poloidal wavenumbers k y and frequencies ω, due to its transport:

ωn -k y n∂ x ψ + k y ψ∂ x n = -iν n (3.7)
The dispersion relation reads:

n = -k y ψ∂ x n Ω + iν (3.8)
where ψ is the Fourier-transformed of the potential ψ. Ω ≡ ω -k y ∂ x ψ is the total frequency accounting also for the Doppler shift, provided by k y ∂ x ψ.

Let us assume a simple notation for the drift velocity v E×B → U , where the bold character is omitted, and consider the hypothesis of incompressible flow, so that a stream function ψ can be defined with the same physical meaning as the electrostatic potential in plasma. Therefore, v E×B = -∇ × (ψ ẑ).

The non-linear term [ ψ, ñ] can be expressed as:

[ ψ, ñ] = ∂ y (ñ∂ x ψ) + ∂ x (ñ(-∂ y ψ)) = ∂ y (ñ Ũy ) + ∂ x (ñ Ũx ) (3.9) 40

The quasi-linear approach

which leads to:

[ ψ, ñ] = ∂ x ñ(-∂ y ψ) = ∂ x ñ Ũx (3.10)
where Ũx = -∂ y ψ and Ũy = ∂ x ψ are the component of Ũ in the x and y direction respectively.

In the Fourier space, the correlation between fluctuating fields reads:

Ũ ñ → k y ω k y ω
Ûk y ω nk y ω e -i(ω +ω )t+i(k y +k )x (3.11)

Applying the averaging operator Ũ ñ on the latter equation leads to k y + k y = 0 and ω + ω = 0, so that:

ñ Ũ → ky,ω Û-ky-ω nkyω = kyω nkyω Û * kyω (3.12)
where Û * kyω is the complex conjugate of Ûkyω . As a function of ψ, the non-linear term [ ψ, ñ] becomes:

[ ψ, ñ] = -∂ x ky,ω -ik y n ψ * = -∂ x ky,ω i|k y ψ| 2 Ω + iν ∂ x n (3.13)
Now the mean equation can be reformulated:

∂ t n + [ψ, n] -∂ x D QL ∂ x n = 0 (3.14)
where one defines:

D QL = ky,ω i|k y ψ| 2 Ω + iν (3.15)
Working out D QL :

D QL = ky,ω i|k y ψ| 2 Ω 2 + ν 2 (Ω -iν) ∼ i +∞ -∞ dΩ|k y ψ| 2 Ω Ω 2 + ν 2 + +∞ -∞ dΩ|k y ψ| 2 ν Ω 2 + ν 2 (3.16)
The first term in equation (3.16) is identically equal to zero for parity reasons (the argument of the integral is an odd function). It is also required that D QL is real. Concerning the second one, provided that ν → 0, it tends to π ky,ω |k y ψ| 2 δΩ. Note that this latter result is independent from ν.

A theoretical work on existing approaches related to the turbulence-driven transport 41

Finally, the quasi-linear transport coefficient reads:

D QL = ky,ω
|k y ψ| 2 πδΩ (3.17)

In equation (3.17), |k ψ| 2 refers to turbulent kinetic energy due to fluctuations of the advective velocity field. δΩ is the Dirac's δ function revealing that diffusion is localized in resonant modes. D QL can be expressed in a compact form which stresses the dimensional dependence: D QL = κτ c , where κ is the kinetic turbulent energy (this notation is consistent with the κ-ε model for neutral fluids), while τ c is the correlation time corresponding to the life time of turbulent structures.

Beyond the quasi-linear approach

The quasi-linear analysis does not provide a prediction for the correlation time. In this section, we investigate this contribution when approaching the quasi-linear procedure from a system of coupled equations, similar to the κ-ε one (see chapter 2). We use the same nomenclature to indicate the turbulent energy κ and its dissipation rate ε.

∂ t κ + [ψ, κ] = -ακ -ε (3.18a) ∂ t ε + ¨¨[ ψ, ε] = -βε - κ τ κ τ ε (3.18b) τ c = κ/ε (3.18c)
-ακ and -βε are the actual restoring forces for κ and ε respectively, as small as wanted, modeled to localize the resonant interaction. τ κ and τ ε are the restoring time scales needed for κ and ε to recover the state of equilibrium, while τ c provides the characteristic time of dissipation. Moreover, in 3.18b, the Poisson brackets are neglected since we suppose it is a sub-grid contribution.

We allow for a scale separation of κ into mean (κ) and fluctuating (κ) components, such that κ(x) = κ , thus κ = 0. The same notation is adopted for ε.

We average equations 3.18a and 3.18b and obtain:

∂ t κ + [ψ, κ] + [ ψ, κ] = 0 (3.19) ∂ε = - κ τ κ τ ε (3.20)
We substract 3.19 from 3.18a and 3.20 from 3.18b:

∂ t κ + [ψ, κ] + [ ψ, κ] + ¨¨[ ψ, κ] -$ $ $ $ [ ψ, κ] = -ακ -ε (3.21a) ∂ t ε = -β ε - κ τ κ τ ε (3.21b)
As usual, non-linear contributions with respect to fluctuations are neglected in the framework of a linear instability approach, thus [ ψ, f ] → 0 and [ ψ, f ] → 0.

Equations 3.21a and 3.21b are Fourier-transformed in space and time, κ and ε being the transformed functions:

(Ω + iα)κ + (k y ∂ x κ) ψ = -ε (3.22a) (Ω + iβ)ε = -i κ τ κ τ ε (3.22b)
The dispersion relation for ε thus reads:

ε = κ/(τ κ τ ε ) iΩ -β (3.23)
Comparing 3.18c with 3.23, one finds that:

τ c = τ ε τ κ (iΩ -β) (3.24)
where τ c ∝ τ κ τ ε . The dispersion relation for κ reads:

κ = -τ c (k y ψ) τ c (Ω + iα) + 1 ∂ x κ (3.25)
Therefore, the non-linear contribution to the transport of κ is provided by:

[ ψ, κ] → -∂ x ky,ω (-ik y ψ * )κ = -∂ x ky,ω i|k y ψ| 2 τ c (Ω + iα)τ c + 1 ∂ x κ (3.26)
The expression for D QL can thus be defined in the form:

D QL,κε = ky,ω i|k y ψ| 2 τ c (Ω + iα)τ c + 1 = i|k y ψ| 2 (Ωτ c +1) 2 (τ c ) 2 + α 2 Ω * (τ c ) -1 -iα ∼ i +∞ -∞ dX|k y ψ| 2 X X 2 + α 2 + +∞ -∞ dX|k y ψ| 2 α X 2 + α 2 (3.27)
where X = Ω + 1/τ c . We recall the same arguments as in the classical quasi-linear analysis to finally find the following expression for D QL,κε :

D QL,κε = kyω |k y ψ| 2 π (3.28)
Compared to the result obtained in 3.17, 3.28 shows that the total correlation time accounts for both frequency Ω and the inverse of turbulent energy dissipation time τ c :

τ c = πδ(Ω + 1/τ c ) (3.29)

A theoretical work on existing approaches related to the turbulence-driven transport 43

Eventually, the evolution of averaged κ is provided by equation 3.30.

∂ t κ + [ψ, κ] + ∂ x (D QL,κε ∂ x κ) = 0 (3.30)
In the asymptotic limit ε → 0, one has τ c → ∞ and τ → πδΩ, which is the result of the standard quasi-linear analysis.

These results show that, on the basis of a κ-ε-like system, it is possible to infer additional information in the correlation time of the quasi-linear diffusivity, concerning the dissipation of turbulent energy. Compared to the standard estimation of quasi-linear diffusivity, where the correlation time is a quantity generally fitted with gyrokinetic simulations [START_REF] Bourdelle | Turbulent transport in tokamak plasmas: bridging theory and experiment[END_REF], the method proposed above is a firststep qualitative approach towards a more predictable estimation of diffusivity in plasma turbulent fluxes.

Prey-predator models

Prey-predator models were first developed to describe the loss-win interactions between two actors in dynamical systems: a resource and a consumer.

In 1925, the US chemist Alfred J. Lotka proposed a basic set of equations to represent such as interplay in chemical reactions where the reagents concentrations oscillate [START_REF] Lotka | Elements of Physical Biology[END_REF]. In 1926, the Italian physicist Vito Volterra independently presented the same kind of mathematical system to outline the increase of predator fish (and corresponding decrease in prey fish) in the Adriatic Sea within a certain observation time [START_REF] Volterra | Variazioni e fluttuazioni del numero d'individui in specie animali conviventi[END_REF].

The so-called Lotka-Volterra model consists in two coupled differential equations:

∂ t x = x(α -βy) (3.31a) ∂ t y = y(γx -δ) (3.31b)
Here, as a general case, x = x(r, t) and y = y(r, t) are the prey and the predator species, respectively; ∂ t x and ∂ t y are their rate of change in time. α is the growth rate of the prey in absence of the predator, β measures the impact of predation ∂tx x , γ is the death (or emigration) rate of y when no interaction with species x is observed and δx denotes the rate of growth (or immigration) of the predator population as a function of the size of the prey. α, β, γ and δ are positive parameters. The characteristic prey-predator dynamics is shown in Figure 3.1(a).

This system exhibits a motion constant, namely a quantity which is constant along any solution (x(t), y(t)) and is determined by the initial conditions (x(r, 0), y(r, 0)):

C = α ln y(t) -βy(t) -γx(t) + δ ln x(t) (3.32)
In the xy-plane, these curves are closed, as illustrated in Figure 3.1(b), meaning that the solutions are periodic oscillations. The Lotka-Volterra system is just one flavor of prey-predator model. The Kermak and McKendrick model [START_REF] Kermack | Contributions to the mathematical theory of epidemics-i[END_REF], for example, studies the susceptile-infectious interactions in the framework of epidemics, the Galbraith one [START_REF] Galbraith | The Predator State[END_REF] investigates on the dynamics between the populace and predator institution occurring in economy. Although the biological origin of this kind of models, indeed, they appear to be very efficient in systems which show periodic loss-win behavior between two (or more) species.

Prey-predator models

An example of a prey-predator model in plasma

This is the reason why in the last years the mathematical formalism of the preypredator model has been found to well represent the picture of self-regulating turbulence via the E × B shear in the edge region of magnetically confined plasmas (see for instance [START_REF] Itoh | A mini-max principle for drift waves and mesoscale fluctuations[END_REF], [START_REF] Zhu | Robustness of predator-prey models for confinement regime transitions in fusion plasmas[END_REF], [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF]). The advantage of such models is the capability in capturing the main mechanisms undergoing the self-regulation of turbulence, while transport models, derived from the fundamental equations of plasma physics, often tend to be extremely detailed.

For instance, in the work proposed by Diamond and Miki in 2013 ( [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF]), the amplitude of turbulent fluctuations represents the prey, while two turbulence E × B shearing mechanisms are outlined as predators (called mean flows and zonal flows). In this model, which is 1D in the cross-field direction, the dynamics of Low-confinement → High-confinement transitions, the competition between shear flows in suppressing turbulence and finally the evolution of mean field density n and pressure p profiles are self-consistently reproduced using transport equations. For this reason, this work has represented a reference in the conception of the new

A theoretical work on existing approaches related to the turbulence-driven transport 45

modeling of cross-field fluxes in SOLEdge2D-EIRENE discussed further in this thesis (chapters 4 and 5).

Let us consider the evolution of turbulent energy I in the Miki and Diamond model:

∂ t I = I(γ I -∆ωI -α 0 E 0 -α V E V ) + χ N ∂ r (I∂ r I) (3.33)
Here the first term in the right-hand side is the source, the second one is the self-saturation of turbulence, regulated by the parameter ∆ω[s -1 ]. This is the nonlinear damping of I. α 0 and α V are coefficients regulating the coupling with the E × B shearing mechanisms and χ N is the turbulent thermal diffusivity.

The first term on the right-hand side represents the turbulence generation by the ion temperature gradient via linear instability, where

γ I = γ I0 (c s /R) (R/L T ) -(R/L T ) crit is its local growth rate.
Here γ I0 is the amplitude, c s the thermal velocity, R the major radius, L T the temperature gradient length and (R/L T ) crit a threshold parameter.

The diffusive components of density and pressure radial fluxes Γ n and Γ p , displayed in 3.34, are linked to turbulence dynamics by the transport coefficients D = χ.

Γ n ∝ D∂ r n Γ p ∝ χ∂ r p (3.34) D = τ c c 2 s I 1 + α t E V (3.35)
where τ c = a/c s , a is the minor radius and α t is a model parameter. The function describing 3.35 assures that the transport coefficient follows the dynamics of turbulence, due to the linear dependence to I (discussed in [START_REF] Diamond | Zonal flows in plasma-a review[END_REF]). D is also inversely proportional to the energy of the shearing flow (E V ). This flow locks the transition into an improved plasma regime, where cross-field turbulent transport is dropped, and so D and χ.

Figure 3.2 shows the time trace of turbulent energy (in blue) and E × B shearing flows (in green and red) on the left, and the phase-portrait of the predatorprey system on the right, highlighting the typical prey-predator dynamics, similar indeed to 3.1(a) and 3.1(b).

This behavior is expected because the coupling between predator and prey recalls the mathematical formalism of 3.31, which consists in the product xy and in 3.31 and in α 0 E 0 I and α V E V I in 3.33. From a mathematical point of view, the prey-predator models represent in general an interesting reference when modeling the competition between plasma turbulence and E × B shears. In this sense the prey-predator formalism can guide the improvement of transport models when the effects of turbulence on mean flows are included. Equation 3.33 keeps the original nomenclature of [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF]. It is worth to make clear that in the next chapters the turbulent kinetic energy (which in [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF] is denoted as I) will be indicated as κ to stress the parallelism with the κ-ε model exploited in Computational Fluid Dynamics. [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF]. On the left, time trace of turbulent energy (in blue), zonal flows (in green) and mean flows (in red). On the right, the phase-portrait of the predator-prey system. For both, a portion of the total time range is considered.

The Lorenz system

Another interesting mathematical reference for the study of plasma dynamics is the well-known Lorenz system, used to describe (the dynamics of) chaotic phenomena. The aim of this section is to display the main properties of such a model.

Equations and physical meaning

The Lorenz system is a mathematical model developed in 1963 by Edward Lorenz to study the atmospheric convection. It displays a set of non-linear deterministic equations relating the properties of a two-dimensional layer of depth H where a temperature difference ∆θ = θ H -θ C is imposed, the layer being uniformly warmed from the bottom side, where θ = θ H , and cooled from the upper one, where θ = θ C , as shown in figure 3

A theoretical work on existing approaches related to the turbulence-driven transport 47

viscosity ν, and undergoes a buoyancy force where g is the gravity acceleration.

In detail, the Lorenz equations describe the time evolution of the stream function ψ, the temperature θ and the heat source θ 2 :

∂ t θ = Raψ -θ -ψθ 2 (3.36a) ∂ t θ 2 = -α 2 θ 2 + θψ (3.36b) ∂ t ψ = Pr θ -ψ (3.36c) (3.36d)
These equations can be derived analytically from Navier-Stokes equations, see appendix C.

Fixed point analysis

The interest of system 3.36 consists in capturing many features about the physics of thermo-convection in a two-dimensional fluid layer. Standard procedure of fixed-point analysis allows to describe in detail such properties. The fixed points are the solutions of equations in 3.36:

0 = Raψ -θ -ψθ 2 (3.37a) 0 = -α 2 θ 2 + θψ (3.37b) 0 = Pr θ -ψ (3.37c)
where the over-line symbols denote the steady-state unknowns. Equation 3.37c readily yields θ = ψ and 3.37a leads to α 2 θ 2 = ψ 2 . Eventually, one has:

0 = ψ α 2 (Ra -1) -ψ 2 (3.38) 48

The Lorenz system

which admits three real solutions only for Ra > 1:

ψ = 0 ; ψ = ± α 2 (Ra -1) (3.39)
When Ra < 1, only the trivial solution ψ = θ = θ 2 = 0 is permitted, which corresponds to the one addressed in the linear analysis. Moreover, this solution describes the system when there is no convection .

When performing a stability analysis, we study the stability of the fixed point in response to a small perturbation, a scale separation being implied. Equation 3.40, shown below, is the fluctuating component of 3.36, where the familiar tilde (•) notation is assumed for the oscillating components of θ, θ 2 and ψ.

∂ t θ = -θ + Ra -θ 2 ψ -ψ θ2 - & & & ψ θ2
(3.40a)

∂ t θ2 = θψ + ψθ -α 2 θ2 + ψ θ (3.40b) ∂ t ψ = Pr θ -ψ (3.40c) (3.40d)
The eigenvalue γ is the solution of:

0 = -(γ + 1) Ra -θ 2 ψ Pr -(γ + Pr) 0 ψ θ -(γ + α 2 )
where:

-(γ + 1) Ra -θ 2 ψ Pr -(γ + Pr) 0 ψ θ -(γ + α 2 ) = -ψ Pr -(γ + Pr) ψ θ -(γ + α 2 ) -(γ + 1) Ra -θ 2 Pr -(γ + Pr) Therefore: -γ 3 -γ 2 (1 + α 2 + Pr) = -γ Pr(1 + α 2 ) + α 2 + ψ 2 -Pr(Ra -θ 2 ) = -α 2 Pr -Prψ 2 + α 2 Pr(Ra -θ 2 ) -Prψθ = 0 (3.41)
which is the dispersion relation. One finds that the fixed point θ = 0, θ 2 = 0, ψ = 0 is stable for Ra < 1 and unstable for Ra > 1, where Ra emerges like the natural control parameter of this analysis. Indeed, when this fixed point is concerned, 3.41 becomes:

(γ + α 2 ) γ 2 + γ(1 + Pr) + Pr(1 -Ra) = 0 (3.42)

A theoretical work on existing approaches related to the turbulence-driven transport 49

yielding to:

γ 1,± = -(1 + Pr) ± ((1 + Pr) 2 + 4Pr(Ra -1)) 1/2 ; γ 2 = -α 2 (3.43)
The condition for the stability of the fixed point is γ 1,± < 0 and thus Ra > 1.

The other fixed points, namely for ψ = θ and ψ 2 = α 2 θ 2 = α 2 (Ra-1) correspond to steady convection. In these cases, the dispersion relation relates:

γ 3 + γ 2 (1 + α 2 + Pr) + α 2 γ(Pr + Ra) + 2α 2 Pr(Ra -1) = 0 (3.44)
It is possible to demonstrate that there exist a critical value of the control parameter Ra = Ra * such that this pair of solutions bifurcates towards instability (see Appendix C). The critical value is given by Ra * ≡ Pr 3+α+Pr Pr-(1+α 2 ) . In particular, close to marginality, for Pr > (1 + α 2 ), ones has three cases:

• 0 < Ra < 1, there is only a single fixed point such that ψ = 0, θ = 0 and θ 2 = 0;

• 1 < Ra < Ra * , the fixed point ψ = 0 is unstable and there exist two stable fixed points with ψ = 0;

• Ra > Ra * : all fixed points are unstable.

The other regime is for Pr < (1 + α 2 ), such that:

• 0 < Ra < 1, there is only a single stable fixed point such that ψ = 0;

• Ra > 1, the fixed point ψ = 0 becomes unstable and there exists two fixed points with ψ = 0 that remain stable for all values of Ra.

In figure 3.5, a set of chaotic solutions for the Lorenz system is displayed, called Lorenz attractor, when α 2 = 8/3, Ra = 28 and Pr = 10, and for initial conditions (1, 1, 1). In this figure, Pr < (1 + α 2 ) Ra * ∼ 24.7 and all fixed points are unstable.

Motivations and strategy of the model

As mentioned in chapter 1, transport codes for magnetically confined plasmas usually rely on ad-hoc diffusion coefficients for the estimation of steady-state crossfield fluxes. This procedure is constraining because the spatial distribution of such quantities, net of time transients, is expected to vary significantly according to several factors, including plasma instabilities, absorbed power, magnetic equilibrium... As far as diffusive cross-field fluxes are modeled, a more generalized description of transport coefficients is thus required. This goal can be achieved by letting these coefficients scale with some property of the turbulent state of the system, as in the quasi-linear theory (see chapter 3). The latter describes a linear relation between diffusivity and turbulent kinetic energy, however more sophisticated dependencies should be taken into account both in transients and steady-state, due to the interplay of turbulent fluxes primarily with shear flows, observed both in experimental observations [START_REF] Cordey | Evolution of transport through the L-H transition in JET[END_REF], [START_REF] Colchin | Physics of slow L-H transitions in the DIII-D tokamak[END_REF] and first-principle simulations of plasma [START_REF] Diamond | Self-regulating shear flow turbulence: a paradigm for the L to H transition[END_REF].

In this chapter, a new modeling is presented for the calculation of (diffusive) crossfield fluxes in the edge and SOL plasmas of diverted tokamaks, when Braginskii's equations are concerned. This model is introduced into the two-dimensional transport code SolEdg2D-EIRENE and comes after the work begun in [START_REF] Bufferand | Interchange turbulence model for the edge plasma in soledge2deirene[END_REF] in the framework of limited plasmas. In its most general form, the model allows to calculate self-consistent transport coefficients when the system reaches equilibrium.

In the first section, we discuss the problematic related to this new modeling and we present the strategy of its design, conceived by steps of increasing complexity depending on both the closure of free parameters and the number of equations concerned. Each step of the modeling is presented exhaustively, from its characteristic equation(s) up to the implementation in SolEdge2D-EIRENE. Finally, for each modeling step, the numerical results are commented and compared to the experimental data of given L-mode discharges of TCV and WEST tokamaks.

Motivations and strategy of the model

Motivations

Chapters 1 and 3 raised two important issues, which represent the main motivations of this thesis:

• Being predictive when modeling cross-field fluxes?

Either homogeneous and ad-hoc distributions of transport coefficients are rough in most of the simulations and do not model consistently poloidal asymmetries. Moreover, the automatic fitting procedure discussed in chapter 1 is efficient a posteriori, meaning that it does not contain any method or information useful to improve the predictive capabilities of the simulation.

• Getting more information on the effects of steady-state turbulence on the main plasma in transport codes Transport of particles, momentum and energy in plasma can be related to the averaged turbulent energy at steady state, as argued in chapter 3, but in the two-dimensional Braginskii equations, where fluctuations due to turbulence are not present, the contributions to the turbulent energy have to be determined.

We address these issues drawing inspiration from the concepts of eddy viscosity and the κ-ε model described in chapter 2. Indeed they do not require a spectral treatment of turbulence, as in the quasi-linear approach. Instead, they rely on Reynolds-averaged Navier Stokes equations, similar to the Braginskii equations numerically solved by transport codes, turning to be appealing for the issue of modeling cross-field turbulent fluxes in plasma. However, the concepts of eddy viscosity and κ-ε model are based on the strong assumption of isotropic and homogeneous turbulence, therefore they cannot be applied to magnetized plasma, where turbulence is known to be anisotropic and non-homogeneous.

Since turbulence in plasmas shows some similitudes with 2D turbulence in neutral fluids ( [START_REF] Shats | Spectral condensation of turbulence in plasmas and fluids and its role in low-to-high phase transitions in toroidal plasma[END_REF]), the evolution of the quantity enstrophy (the integral of vorticity square), which in 2D turbulence is an inviscid invariant, may be addressed for plasmas too. In [START_REF] Coosemans | A new mean-field plasma edge transport model based on turbulent kinetic energy and enstrophy[END_REF], for instance, a set of transport equations is derived analytically for the evolution of turbulent kinetic energy and enstrophy in isothermal plasma.

As a working assumption, in this work we propose a heuristic approach to customize the semi-empirical transport equations for κ and ε for non-isothermal plasma, including the physics of interchange-driven turbulence. This system is coupled to transport equations of plasma particles, momentum and energy and it is closed employing two different approaches, based on experimental data on operating tokamaks. Other closures may be taken into account than an empirical one, such as integrating transport and first-principle fluid codes to allow for more accurate estimations of transport coefficients. In this work we focus on the former solution.

Heuristic approach means that while the structure of transport equations 2.32 is mostly preserved, some terms are adapted to model specific physical mechanism typical of magnetically confined plasmas and observed experimentally, like the enhanced cross-field transport at the low-field side. This process is mostly driven by the interchange instability, as mentioned in chapter 1. Such an adaptation has been inspired by the work of [START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF] in the framework of prey-predator models for plasma turbulence, already debated in chapter 3. In particular, we adapt the following terms:

• The drive. The source term, which in equation 2.32a is given by RS ij

∂u i ∂x i = µ t S ij 2
, is turned into γ κ κ, which is dimensionally consistent. γ κ is the growth rate of the instability driving the turbulence.

• The sink term. In the κ-ε model for neutral fluids, ε represents the dissipation rate of turbulent kinetic energy towards small, dissipative scales. Its meaning in the framework of magnetically confined plasmas cannot be hold, since the dissipation of turbulent energy is way more puzzling. Many kinds of instabilities may inject energy to different scales and several dissipation mechanisms may take place there at the same time (see for instance [START_REF] Shats | Spectral condensation of turbulence in plasmas and fluids and its role in low-to-high phase transitions in toroidal plasma[END_REF]). Therefore, ε does not keep its original meaning in the framework of this thesis: it should be rather interpreted like a general predator of turbulent kinetic energy κ.

Strategy

The model has been designed by steps concerning three different branches: (i) the equations; (ii) the closure; (iii) the experimental data set, as schematically shown in 4.1. In the first step, the model consists in one equation; it is closed by an improved automatic fitting procedure implemented in SolEdge2D-EIRENE, which fits the numerical results with diagnostics data from a single experiment; this diagnostics refer to the TCV tokamak.

In the second step, the one-equation model holds valid, but the closure involves the exploitation of an empirical scaling law. Three discharges from TCV are exploited as benchmark.

Finally, in the third step, two equations are solved for the model, using the same empirical scaling law as in the previous step, but the match with experimental data concerns the WEST machine.

First step: one-equation model and tuning with TCV data

The first step consists in implementing a one-equation model for turbulent kinetic energy κ in addition to the governing equations of SolEdge2D. Then, the diffusion of plasma density and temperature in the cross-field direction is set proportional to κ at steady-state, in the spirit of the quasi-linear approach. A feedback-loop procedure ensures that free parameters are optimized by tuning numerical results with one experimental test-case given by a TCV discharge. These results, displayed at the end of the section, have been published in [START_REF] Baschetti | Optimization of turbulence reduced model free parameters based on L-mode experiments and 2D transport simulations[END_REF].

Equation and dimensional analysis

The time evolution of κ is described by 4.1, which is a semi-empiric transport equation obtained with an heuristic procedure. Basically, κ is advected in space by both diffusion and convection like moments in the Braginskii model, while source and sink are algebraically designed. Moreover, it is semi-empirical because the dissipation of turbulent energy relies on free parameters, which have to be closed experimentally. Also, the κ equation undergoes the same averaging as density, momentum and energy in the Braginskii approach. Here the turbulent energy is defined as:

κ = 1 2 ũ2 m 2 /s 2 . ∂ t κ + ∇ (κu ) + ∇ ⊥ • (κu ⊥ ) -∇ ⊥ • (D κ ∇κ) = γ κ κ -∆ωκ 2 (4.1)
For the sake of simplicity, hereafter the bold notation for vectors will be omitted.

On the left-hand side we recognize the material derivative of κ, with separate components for contributions parallel and perpendicular to magnetic field lines. This latter is linked to mean field drift velocities u ⊥ , which so far are not addressed in this work, therefore hereafter ∇ ⊥ • (κu ⊥ ) is neglected. Indeed, we suppose that the transport of turbulent energy is purely diffusive, governed by diffusivity D κ . The first term on the right-hand side is the growth rate γ κ of leading instabilities. Here we assume γ κ = γ I + βT 2 e , as expressed in equation 4.2. γ I is the linear interchange instability growth rate, while βT 2 e , which is introduced heuristically, represents a turbulent background in hot plasma, where β is a free parameter whose dimension is eV -2 s -1 .

γ κ = γ 0 c s R R 2 ∇p i ∇B T p i B T + βT 2 e (4.2)
where we recall that R is the system size (tokamak major radius), c s the sound speed in plasma and B T the toroidal magnetic field. We also set γ 0 = 1. The last term on the right-hand side, ∆ωκ 2 , is defined as the self-saturation of turbulence: it keeps the system stable, by avoiding that κ grows indefinitely. The order of the exponent is the minimum to ensure such a feature. It is the second free parameter of the model and its dimensions are m -2 s.

To determine the leading terms in equation 4.1, dimensional analysis is mandatory. Hereafter, the subscript 0 denotes the reference dimensional quantities, while the index * describes non-dimensional variables, such that for a generic f variable, f * = f /f 0 .

Let us consider the normalized equation divided by κ 0 /t 0 , where t 0 = R/c 0 is the reference time and c 0 = T 0 /m i is the average thermal velocity of charged particles 56

First step: one-equation model and tuning with TCV data

in the parallel channel, T 0 being the reference temperature. Such a choice for the reference time addresses to relatively slow evolution scales (especially compared to the inverse of gyrokinetic frequency or even the drift wave time τ d ∼ c s /a, where c s is the thermal velocity and a is the tokamak minor radius). Therefore, it suits the modeling of turbulent cross-field fluxes at steady-state.

∂ t * κ * + c 0 t 0 L ∇ * (κ * u * ) = D 0 t 0 l 0,⊥ 2 ∇ * ⊥ • (D * κ ∇ * ⊥ κ) + γ * κ κ * -∆ω * κ 2 * (4.3)
On the right-hand side, D 0 and l 0,⊥ 2 are the reference diffusivity and perpendicular length. Moreover, basing on the relation between parallel and perpendicular channel, the following dimensional analysis can be outlined:

D 0 t 0 l 0,⊥ 2 ∼ c 0 t 0 L = 1 2πq (4.4)
where we used for the connection length: L = 2πqR, with q as safety factor. Therefore, the contributions related to the advection are at lower order of magnitude than the other terms in the equation.

As mentioned at the beginning of this section, 4.1 is a semi-empiric transport equation modeled exploiting heuristic arguments. In principle, one can approach the analytical derivation of a κ-like equation with the same rationale as in the derivation of the standard κ equation in Computational Fluid Dynamics, namely manipulating the fluid force balance equation. In the appendix B, however, which has to be intended as an exercise, the procedure is very general, with no in-depth expansion of either the pressure tensor, the external force (which is the electromagnetic field) and the collisional term.

Free parameters and model's closure

β and ∆ω are the free parameters of the model, the former referring to the modeling of a turbulent background in the core plasma, the latter relating to the selfsaturation of turbulent kinetic energy. In order to close the model, we tune them with diagnostic data available for a single experimental test case. First, an expression for the transport coefficient is provided proportional to κτ c , with τ c + R/c s , in the spirit of the quasi-linear approach. D n is the particle diffusivity, exploited hereafter as a reference transport coefficient. Indeed, as usually done in transport codes, we assume that electron and ion thermal conductivities are proportional to D n , thus in the following, conclusions are also valid for χ e and χ i .

κ is the non-trivial, steady-state solution of equation 4.1. It is easy to demonstrate that, when neglecting the higher-order terms (the advection contributions), κ = γ κ /∆ω, therefore: Experimental data on density and temperature (n exp , T e,exp and T i,exp ) are diagnosed at the outer midplane and each set of data is fitted by SolEdge2D (see figure 4.5 concerning density measurements). Here, we enforce the radial profile of the diffusivity (D ⊥,κ ) at the outer midplane to match the radial distribution inferred from the fitting procedure (D 1D ⊥,f b ), generating automatically a constraint for β and ∆ω. Hereafter, we refer to D 1D ⊥,f b as the feedback diffusivity.

D n = t c κ = t c γ κ ∆ω = R c s c s /R ∆ω R 2 ∇p i ∇B φ p i B φ + β ∆ω T 2 e ( 4
With such constraints on radial profiles of plasma quantities at the outer midplane, the code numerically solves their distributions in poloidal plane (n 2D , T 2D e and T 2D i ), which in turn are coupled to the κ model, resulting in a first-attempt transport coefficient D ⊥,κ . The expression first-attempt refers to the fact that the diffusivity is not optimized at this stage of the loop, and so are the free parameters.

Experimental setting

The experimental data set available for the fitting of the model is the L-mode, lower-single-null, ohmic, low-density, attached plasma discharge #51333 in TCV (Tokamak à configuration Variable). The sketch of the magnetic separatrix for this discharge is displayed in blue in figure 4.4. The main plasma parameters for this test-case are presented in Tab.4.1: B tor is the toroidal magnetic field, I pol is the poloidal electric current, P Ohm is the ohmic power injected in the plasma and n sep is the density at the separatrix. At this stage of the modeling, TCV has been chosen as reference machine since in general it benefits from a wide set of edge diagnostic (as shown in figure 4.4): high resolution Thomson scattering (HRTS) and reciprocating Langmuir probes (RCP) plunging at the outer midplane, well-embedded Langmuir probes (LP) to monitor the divertor and the infrared thermography system (IR) (see [START_REF] Coda | Overview of the tcv tokamak program: scientific progress and facility upgrades[END_REF], [START_REF] Maurizio | Divertor power load studies for attached l-mode single-null plasmas in tcv[END_REF] and [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF]).

Results

The model is implemented in the 2D transport code SolEdge2D-EIRENE assuming pure deuterium plasma, neglecting the mean field velocity drifts. Moreover, in the code, the contributions due to the advection of κ appearing in equation 4.1 are included and they contribute in smoothening the solution.

D n (𝜅) In figure 4.5, radial profiles of the model diffusivity (blue line) and the feedback diffusivity (black line) are displayed. Ion and electron thermal conductivities are expected to be proportional to D 1D ⊥,κ by construction. Despite the optimization procedure, a finite relative error persists between the two profiles. This relative error is provided by the fact that the interchange instability growth rate in equation poloidal plane. The model captures the poloidal asymmetry of the turbulent transport due to the interchange instability (ballooning). Similar qualitative results were obtained in [START_REF] Bufferand | Interchange turbulence model for the edge plasma in soledge2deirene[END_REF], where, also, the cross-field transport is driven by the interchange instability on the LFS for a limited plasma (MISTRAL test-case, see [START_REF] Paredes | A penalization technique to model plasma facing components in a tokamak with temperature variations[END_REF]).

Finally, the profiles of density, electron temperature, saturation current and perpendicular heat flux at the outer divertor are displayed in figure 4.8 as a function of the coordinate s along the wall. s SP denotes the position of the separatrix. The comparison with the Langmuir probes suggests a general overestimation of the peak value of a factor variable from 2 to 3, except for the electron temperature which shows a nice agreement with the peak value of the experimental temperature.

Z [cm] R [cm]

Second step: one-equation model, tuning with experimental scaling law

The model exploited above is interesting because it allows to calculate self-consistently the poloidal distribution of the cross-field transport. However, its free parameters still depend on a single, experimental test-case. When changing machine, or even discharge, they cannot be used anymore.

To enhance the predictive capabilities of the model and allow for the calculation of plasma mean quantities in the poloidal plane when changing geometry, magnetic configuration or input conditions, we introduce a new way to close the model [START_REF] Baschetti | A model for plasma anomalous transport in tokamaks: closure via the scaling of the global confinement[END_REF]. In this second step, the one-equation model is closed by the experimental scaling law proposed in [START_REF] Goldston | Heuristic drift-based model of the power scrape-off width in low-gaspuff h-mode tokamaks[END_REF] for the SOL width.

Equation

We recall the equation for κ seen in 4.1: term accounting for the background turbulence of the core plasma βT 2 e = 0, thus obtaining:

∂ t κ + ∇ (κu ) -∇ ⊥ • (D κ ∇κ) = γ κ κ -∆ωκ 2
γ κ = c s R R 2 ∇p i ∇B ϕ p i B ϕ (4.7)
The dimensional analysis holds as in equation 4.3.

Free parameters and model's closure

∆ω is the only free parameter of the model. To evaluate it, we first consider the non-trivial, steady-state solution of equation 4.6, provided that the transport contributions are higher-order terms:

κ = γ κ /∆ω (4.8)
Then, in the framework of dimensional analysis, we exploit the relation describing the equilibrium of the plasma pressure (1.35) discussed in chapter 1:

2γ e λ 2 SOL χ e L c s (4.9)
where γ e is the electron sheath transmission coefficient while χ e is the electron thermal conductivity. Hence, provided (4.8), the SOL width λ SOL can be estimated as:

λ 2 SOL L γ e c s χ e = L R 2γ e c 2 s γ κ ∆ω (4.10)
where the inverse relation returns the self-saturation coefficient ∆ω:

∆ω = L R 2γ e (λ scaling SOL ) 2 γ κ 1 c 2 s (4.11)
Here the closure comes into play by assuming for λ SOL the scaling law proposed in [START_REF] Goldston | Heuristic drift-based model of the power scrape-off width in low-gaspuff h-mode tokamaks[END_REF] and tested experimentally with nice good agreement in several machines ( [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] and [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF]).

This scaling law relies on a heuristic drift-based model for high-confinement and low-gas-puff plasma regime, stating that in the SOL vertical ∇B and curvature drifts are balanced against near-sonic parallel flows to the divertor plates, resulting in a SOL width whose expression is:

λ SOL ∼ 2qρ L,i (4.12)
q is the safety factor while ρ L,i is the ion Larmor radius. It is further assumed that the turbulent perpendicular electron thermal diffusivity is the main source of heat flux across the separatrix, flooding the SOL width with heat from the core plasma. Therefore, while the neoclassical interpretation is provided for the cross-field particle channel, turbulent arguments are used to describe the energy one in the SOL. Comparisons of the expected heat e-folding length λ q with scaling laws inferred from measurements in the JET and ASDEX tokamaks for both low and high-confinement regimes are shown to be in satisfactory agreement, as reported in [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF].

In this latter paper, λ q measured in L-mode discharges is found to be 2 -3 times larger than predicted by the scaling law in H-mode. We use such an empirical relation to infer the following expression to be exploited in (4.11):

λ SOL ∼ 4qρ L,i (4.13)
Indeed, in this stage of the modeling, we only refer to L-mode discharges. Recalling q = aB T /RB P and ρ L,i = m i v ⊥,i /q i B, equation 4.11 becomes:

∆ω = πA 2 16γ e q cyl ρ 2 * γ κ 1 c 2 s (4.14)
with A = R/a the aspect ratio and ρ * = ρ/a the normalized gyro-radius. Equation 4.14 shows that ∆ω is a non-trivial function of the magnetic field (through the term γκ q cyl ρ * ). For given tokamak and discharge conditions, the scan of the model on the outer divertor leg L div may represent a first step to test such closure and the model itself.

Experimental and numerical setting

For the experimental match, we chose a set of lower single null (LSN), Ohmic, L-mode discharges where the outer divertor leg length L div is varied while main plasma parameters are nearly unchanged. This experiment is widely described and investigated in [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF].

The three magnetic configurations consist in: the short-leg (L div = 21 cm) shot #51262, the medium-leg (L div = 36 cm) shot #51333, and the long-leg (L div = 64 cm) shot #51325. The following parameters are kept constant: major radius R = 89 cm, minor radius a = 22 cm, elongation k = 1.4, plasma current I p = 210 kA, poloidal magnetic field at the outer midplane B p,omp = 0.18 T. The latter are used as input parameters in SolEdge2D-EIRENE.

As boundary conditions in the code, the gas-puff is adjusted to obtain the same density at the separatrix for the three cases, while for the input power crossing the separatrix and entering the SOL P in , we set a reasonable value of 100kW .

Results

In this section, the numerical results from SolEdge2D-EIRENE (published in [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]) are displayed and discussed.

Outer mid-plane distributions of plasma are first investigated when varying L div . Since this procedure does not change the main plasma (input power, plasma current, gas puff, poloidal magnetic field B P are kept constant), density and temperature profiles are not expected to change meaningfully, in fact HRTS data related to the three configurations overlap. This is indeed retrieved in SolEdge2D outputs, see Figure 4.9. The overall trend is recovered thanks to the model and found to be in a very good agreement with experimental data from the HRTS diagnostics in the SOL, for both density and electron temperature, and also in the edge region for T e . Two-dimensional maps in figure 4.10 show the self-consistently calculated distribution of the cross-field particle diffusivity D n on the poloidal plane, at steady-state.

In all the three magnetic configurations, the ballooning feature is observed. The ballooning is responsible for the enhanced cross-field transport at the outer mid-plane downwards close to the X-point, due to the "bad" curvature of the LFS (∇p•∇B > 0). Simulations show that the model, which cannot count on the calculation of fluctuations, is also able to capture the fact that this enhanced transport broadens up to the common flux region (where the SOL plasma is constrained by the separatrix below the X-point), but not in the private flux region. This picture is consistent with TOKAM3X isothermal simulations in TCV-like geometry [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF] and Compass-like geometry [START_REF] Galassi | Drive of parallel flows by turbulence and large-scalee×btransverse transport in divertor geometry[END_REF], when neutrals are not modeled. This picture is consistent also with experimental observations in the MAST tokamak [70]. We verify whether the distribution of the cross-field particle flux Γ n,⊥ = -D n ∇n is also concerned by the interesting features seen for the particle diffusivity, in figure 4.11: Qualitatively, both the ballooning at the LFS and the enhanced cross-filed transport in the common flux region are conserved. In particular, when compared to Γ n,⊥ from TOKAM3X turbulent simulations (in figure 4.12), the distribution of particle flux estimated by our model denotes a sharper difference in transport between LFS and HFS, where turbulent fluxes are almost absent, probably due to the lack in the modeling of instability mechanisms other than interchange.

𝑅 -𝑅 𝑠𝑒𝑝 [𝑚]

Finally, the heat flux experimental profiles measured for each discharge by the IR system are compared to numerical profiles in figure 4. [START_REF] Wagner | Regime of improved confinement and high beta in neutral-beam-heated divertor discharges of the asdex tokamak[END_REF], where a nice overlap in both absolute value and shape of the profiles can be observed. In the SOL region, the slope of heat flux decay is overestimated for the short and medium configurations, while in the long one the estimation of the SOL width is very good.

We use the standard parametrization described by the Eich's formula ( [START_REF] Asdex Upgrade Team | Inter-elm power decay length for jet and asdex upgrade: Measurement and comparison with heuristic drift-based model[END_REF]) to estimate λ q as the best approximation for λ SOL . We find that the value of λ q for each discharge is comprised in the experimental range of [5 -10]mm (see [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF]): The trend of λ q against the outer divertor leg length is shown in figure 4.14, left panel. The linear trend observed in the experiment and reported in [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF] is not retrieved for the three configurations. This result may be expected since these discharges are

51333 51325

λ SOL = 7.8mm λ SOL = 10.3mm λ SOL = 9.1mm not following the scaling law for the λ SOL used to close the model.

Third step: two-equation model, tuning with experimental scaling law

The third step consists in coupling the transport equations of κ and ε, which is a quantity related to the dissipation of turbulent energy. The advantage is adding in the description further information on turbulence thanks to a second reference scale provided by steady-state ε (the first one is related to κ). We also hold the same closure on the scaling law of λ SOL as in the second step modeling. against the outer divertor leg length L div . In the right panel, the empty squares stand for target infrared data and the green diamond stands for the outer midplane reciprocating probe data (see [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF] for further details).

𝑅 -𝑅 𝑠𝑒𝑝

[𝑚] 𝑅 -𝑅 𝑠𝑒𝑝 [𝑚] 𝑅 -𝑅 𝑠𝑒𝑝 [𝑚] 𝑞 ⊥ [𝑊𝑚 -2 ] 𝑞 ⊥ [𝑊𝑚 -2 ] 𝑞 ⊥ [𝑊𝑚 -2 ]

Equations and dimensional analysis

The equations describing the evolution of κ and ε are displayed below:

∂ t κ + ∇ (κu ) -∇ ⊥ • (D κ ∇κ) = γ κ κ -∆ωκ 2 -ε (4.15a) ∂ t ε + ∇ • (εu ) -∇ ⊥ • (D ε ∇ε) = γ ε ε - P κ 3/2 ε 2 (4.15b)
where the equation for κ is the same as 4.6 except for the coupling term ε. In analogy with equation 4.15a, in 4.15b the left-hand side contains the material derivative and the diffusive term governed by diffusivity D ε . On the right hand side, the source is driven by the growth rate γ ε , while the last contribution relates to an intricate dissipation term which is a function of κ, ε, and a new parameter, P . As working assumption, the structure of this latter is designed such that at steady state (and net of transport contributions ∇ • (εu ) and ∇ ⊥ • (D ε ∇ε)) 4.15b reads:

∝ κ 3/2 (4.16)
In the theory of neutral fluids this scaling relation occurs in the energy inertial range of inverse energy cascade when two-dimensional turbulence is addressed (see chapter 2). In such a range, energy is transferred from the injection scale to scales with lower wavenumber, corresponding to larger characteristic lengths. Also in magnetized plasmas, large-scale flows called zonal flows, which are responsible for the shear of turbulent fluxes, can be excited or damped due to the non-linear coupling of smaller scales (see [START_REF] Diamond | Zonal flows in plasma-a review[END_REF]). For this similarity between inverse energy cascade in neutral and plasma turbulence, 4.16 is enchosen as constraint for the heuristic design of equation 4.15b.

Free parameters and model's closure

At this stage of the modeling, ∆ω is only assumed to assure that κ does not diverge to infinity when the dissipation rate ε goes to zero. To determine ∆ω, we consider equation 4.15a, net of transport terms and at steady-state, when ε = 0 and therefore κ reaches a prescribed maximum value κ max :

γ κ -∆ωκ max = 0 (4.17)
Thus one readily obtains: ∆ω = γ κ /κ max . P (m/s) is a dimensional free parameter governing the saturation of turbulence and that is closed by the scaling law for λ SOL from [START_REF] Goldston | Heuristic drift-based model of the power scrape-off width in low-gaspuff h-mode tokamaks[END_REF] used in the previous section.

To determine its expression, we start recalling the dimensional analysis for the equilibrium of the electron pressure (equation 4.18):

2γ e λ 2 SOL χ e L c s (4.18)
and therefore:

λ SOL = L χ e 2γ e c s ( 4.19)
A new functional dependence of the transport coefficient is now introduced. In analogy with the κ-ε in neutral fluid, we set:

χ e = C µ κ 2 ε (4.20)
where C µ is prescribed. To determine κ and ε, we solve the system described in 4.15 at the equilibrium, net of the contributions provided by the transport of both κ and ε: The non-trivial solution is plugged in 4.20 which reads:

γ κ κ -ε = 0 (4.21) γ ε ε - P κ 3 2 ε 2 = 0 (4.22)
χ e = C µ P γ ε 2 γ κ (4.23)
The inverse relation provides:

P = λ SOL γ ε 2γ e c s C µ γ κ L (4.24)
Finally, we exploit the scaling law given by 4.13 for L-mode plasma: λ SOL ∼ 4qρ L,i and make each term explicit against plasma and machine parameters, obtaining:

P = γ ε γ κ c s ρ * γ e q cyl C µ π γ κ R c s α S A (4.25)

Experimental and numerical setting

In the framework of a multi-machine application, the full κ-ε model is implemented in the SolEdge2D-EIRENE code assuming as reference test-case an early lowersingle-null discharge from WEST tokamak during its second experimental campaign. For this latter, the diagnostics available are: interfero-polarimetry, Reciprocating Langmuir Probes (RCP) at the outer midplane, and Langmuir probes (LP) at the outer divertor.

A WEST discharge has been chosen as a reference experiment to validate the new modeling of the cross-field transport against a different machine from TCV (multi-machine approach). However, for the sake of completeness, the full κ-ε should also be tested and validated against TCV data exploited in the previous sections.

In table 4.4 the operation parameters are reported.

I p [kA] ne [m -2 ] LH1/LH2[M W ] 700 1.8•10 19 1.6/0.5
The simulation is run assuming pure deuterium plasma, no drifts and the boundary conditions of table 4.5. n BC is the density of plasma at the entrance of n BC [m -3 ] P in [kW ] 1.9•10 19 400 the SOL, P in is the rate of energy which flows from the main plasma crossing the separatrix.

Results

In this section very preliminary results of the κ-ε model implemented in the SolEdge2D-EIRENE are displayed.

At the outer midplane, the numerical profile of the density in the SOL region falls within the experimental data range defined by the outputs of the interferopolarimetry and the reciprocating Langmuir probes (RLP) (see figure 4.15), with a certain degree of incertitude. Thus only a qualitative comparison can be provided. The agreement is good at the separatrix, since the density profile is constrained only at the core boundary.

𝑅 -𝑅 𝑠𝑒𝑝 [𝑚] 𝑛 [𝑚 -3 ]

Model

Interf.-pol. RCP The poloidal distribution of the particle diffusion coefficient is represented in 4. [START_REF] Schneider | Plasma edge physics with b2-eirene[END_REF]. In this two-equation version of the modeling, the ballooning at the LFS is retrieved and appears as a sharp region of enhanced cross-field transport, which is radially confined in the edge region.The X-point is affected by enhanced cross-field transport as well, which in future analysis should be further investigated.

Electron density, parallel density current and parallel heat flux at the outer divertor and remapped at the outer midplane are displayed in figure 4.17 as a function of the normalized magnetic flux. They are compared to Langmuir probes (LP) denoting a nice agreement. Finally, figure 4.18 shows the heat flux perpendicular to the outer divertor plate and remapped at the outer midplane. The estimation of the heat flux width, taken as best approximation for λ SOL , provides λ q 9mm. This is in quantitative agreement with measures of heat flux in the main plasma obtained with reciprocating Langmuir probes across multiple shots, shown in the right figure of 4.18 (courtesy of [START_REF] Ciraolo | First modeling of strongly radiating west plasmas with soledgeeirene[END_REF]). Further analysis is required for this latter modeling step of the κ-ε-like approach to turbulent fluxes in SolEdge2D-EIRENE.

𝑅 [𝑐𝑚] 𝑍 [𝑐𝑚] 𝐷 𝑛 [𝑚 2 /𝑠]

5

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 76

System of equations

A new semi-empirical modeling of the cross-field fluxes in diverted plasma has been outlined in chapter 4 and it has been showd that the model, implemented in SolEdge2D-EIRENE, increases the predictive capabilities of the code, since numerical results favorably compare to experimental data.

However, a parametric analysis should be implemented, consisting in testing the sensitivity of the model against plasma and machine parameters, to compare numerical results against empirical scaling law and disclose the range of validity of the model itself.

This kind of analysis is demanding when exploiting two-dimensional transport codes, due to the physics (a detailed description of the charged particle motion along the magnetic field lines or the plasma-wall interactions), and the numerics (two-dimensional grid). Therefore, fast simulations based on one-dimensional reduced models may be appealing.

In this chapter, we project the new semi-empirical model presented in chapter 4 on the radial direction of the poloidal plane. Such a simplification aims to build a sort of toy-model, where the main physical mechanisms driving plasma quantities are preserved, in the same spirit of prey-predator models introduced in chapter 3. The equations and the numerical setting are outlined. Numerical scaling laws for both the width of the SOL and the energy confinement time are compared to the experimental ones provided by [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF] and [START_REF] Transport | Plasma confinement and transport[END_REF] respectively and commented. Finally, some conclusive perspectives are outlined.

System of equations

A system of partial differential equations describing the evolution of plasma density, energy, vorticity, κ, ε and neutrals density in the core up to the scrape-off layer is integrated using the integral operator:

1 L L 0 dl
(5.1)

The resulting system reads:

∂ t n + 2χ m nc s L -∇ r (D n ∇ r n) = S n (5.2a) ∂ t E α + 2χ m nT α c s L -∇ r (D n T α ∇ r n + χ α n∇ r T α ) = S E (5.2b) ∂ t Ω - 2χ m nc s L φ T e -Λ -∇ r (ν∇ r Ω) = 0 (5.2c) ∂ t κ = γ κ κ -ζκ 2 -ε + ∇ r (D κ ∇ r κ) (5.2d) ∂ t ε = γ ε ε - P κ 3/2 ε 2 + ∇ r (D ε ∇ r ε) (5.2e) ∂ t n 0 = D n0 ∇ 2 r n 0 + I n0 -S n (5.2f)
where n is the density, E α = 3 2 nT α is the energy (per unit volume), where α = e, i for electrons and ions respectively), Ω is the vorticity, κ and ε the turbulent energy

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 77

and its rate of dissipation, n 0 the neutral density. All of them are functions of the radial coordinate r and time t. Equations 5.2a and 5.2c are integrated from 1.24a, 1.24c and 1.24d in chapter 1. Equations 5.2e and 5.2f for κ and ε are integrated from 4.1 and 4.15b in chapter 4. Finally, a vorticity equation and a diffusive-like equation for neutrals have been added. The former is derived by manipulating the force balance in the radial direction (see appendix A), Ω being defined as:

Ω ≡ m i ∇ r q e n∇ r φ B 2 + ∇ r p i B 2 (5.3)
where φ is the electrostatic potential, q e is the electron charge and p i is the ion pressure. In most results presented in this chapter, equation 5.2c is not coupled to the others, except in last section (5.5). Equation 5.2f is a simple model to take into account the particle source due to neutrals ionization, in particular the rearrangement of the ionization front depending to the plasma temperature: S n = σv iz nn 0 , where σv iz is the ionization rate coefficient measured in m 3 s -1 , and E i is the ionization potential for the hydrogen, measured in eV . An expression for σv iz can be found in [START_REF] Richardson | NRL PLASMA FORMULARY[END_REF]:

σv iz = 10 -5 (T e /E i ) 0.5 E 1.5 i (6 + T e /E i ) exp(-E i /T e ) (5.4)
As the input power, and so the plasma temperature, increases, the peak of σv iz moves outward namely towards the wall. In figure 5.1, the steady source rate S n is . Both the SOL region and the core are simulated by this model, which represents a distinctive feature compared to other predator-prey models designed for magnetized plasma. The distinction between these two is hold by the mask function χ m in the

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 79

parallel contributions of equations 5.2a, 5.2c and 5.2c, which describe the Bohm condition for the respective variables near the wall. χ m = 1 in the SOL while χ m = 0 in the core.

S E is the external heat source rate, localized in the core region and expressed by the following relation:

S E = P in exp - 1 1 -(r/a

The numerical setting

The system described in equation 5.2 is solved numerically after discretization in space and time according to the finite volume method and the backwards Eulerian method respectively. The time step is assumed constant, while the mesh is modeled to be denser in the separatrix region.

On the inner and outer boundaries, we impose zero-gradient for each unknown. A diffusive operator is defined and recalled for each equation to build the component of the characteristic matrix which depends on transport coefficients. To test such an operator, we compare the analytical and numerical profiles of density, in the simple case of constant particle diffusivity (D n = 1m 2 /s) isothermal plasma (T α = cst but also S n = cst) and zero-density at the outer boundary (therefore, the Bohm condition is not addressed here and there is no physical distinction between core and SOL).

Let us consider the density equation (5.2a) at steady state:

1 r ∂ r (D n r∂ r ) = -S n (5.6)
We integrate 5.6 twice obtaining:

n(r) = n(r = 0) - S n r 2 4D n (5.7)
We also exploit the zero-density condition at the outer boundary (where r = r max ):

n(r) = S n (r 2 max -r 2 ) 4D n (5.8)
which indeed provides n(r max = 0). The comparison between numerical and analytical results is shown in figure 5.6. The quality of this diffusive operator is proved by estimating the L 2 norm for increasing mesh refining. Figure 5.7 shows that the relative error calculated as L 2 norm increases with the square of the average space step h 2 ave , resulting in a reliable diffusion operator.

Results for a scan against the input power

Equations 5.2a-5.2f are solved for given value of the input power: P in is increased within the range [0.1M W -6M W ] for a WEST-like set of data (see 5.1). The power released from the core region, P in , is a key parameter to scan the solutions of equation 5.2. Indeed, the rate of energy entering the system affects plasma and neutrals balance by varying the ionization rate coefficient, regulates the equilibrium of ion and electron energy an therefore the drive of turbulent energy, as illustrated in figure 5.8. Moreover, P in has an engineering relevance since it can be related to the auxiliary power injected to the plasma (net of the radiative power which so far is not concerned in this model).

Results for a scan against the input power

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 83

Comparison with experimental scaling laws

Scan of the SOL width

Fast scans of plasma main quantities with respect to machine parameters allow for direct comparison with experimental scaling laws. We compare engineering quantities calculated from the model (notably the SOL width and the energy confinement time) to experimental scaling laws widely used in plasma community, in order to identify the parametric limits of validity of the model and, likely, regimes where these engineering quantities bifurcate due to different operational regimes of plasma. λ SOL is the first quantity addressed due to its crucial role in governing the power exhaust in tokamaks. The numerical scaling of the plasma pressure width, assumed as the best approximation for the λ SOL , is plotted against the input power.

Then it is compared to the empirical scaling law for the heat flux width from [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF].

Figure 5.9(a) and 5.9(b) display the scaling of λ SOL (continuous black line) against increasing input power P in ∈ [0.1M W, 6M W ] and poloidal magnetic field B P . This trend is compared to the scaling law for the heat flux width at the divertor (remapped at the outer midplane) from [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF], in red in the figure. This scaling law reads: λ q = 1.44 ± 0.67B -0.8±0.32 T q 1.4±0.67 P 0.22±0.1 in R -0.03±0. 28 (5.9)

As P in increases, the ion Larmor radius ρ = mv qB is also expected to increase due to the thermal velocity, and so the λ SOL , while as B P rises, the safety factor q cyl = aB T RB P drops because it is inversely proportional to this control parameter, and so the SOL width.

λ q is estimated by the slope of the exponential regression of pressure in the first region of the SOL, near the separatrix (see figures 5.10(a) and 5.10(b)). The width of this region covers about 20% of the SOL width. The grey-shaded area is the 95% confidential interval corresponding to such regression. Figure 5.9(a) and 5.9(b) show that the numerical values obtained for the λ SOL follow the same trend of the scaling law, with a maximum relative error of nearly 15% within the range of input power explored, and nearly 11% within the range investigated for B P .

It is worth to stress that the nice agreement between numerical results and empirical data shown in figure 5.9(a) cannot be retrieved when assuming a constant transport coefficient, as usually accepted in the framework of transport coefficients. In this case, indeed, λ .

Finally, we address the sensitivity analysis of numerical pressure width against the major radius R and the toroidal magnetic field B T , retrieving that λ SOL decreases when the major radius R drops down and there is almost no dependence with respect to B T , consistently with 5.9.

Scan of the energy confinement time τ E

The same exercise is proposed for the thermal energy confinement time scaling against P in and B P and compared to the empirical law for L-mode plasma discussed in [START_REF] Transport | Plasma confinement and transport[END_REF], usually called ITER96-th. In the previous section, a nice agreement between the λ SOL calculated by the model and the scaling law has been verified. Physically, there exists a link between the SOL width and the energy confinement time, which has been investigated experimentally ( [START_REF] Brunner | High-resolution heat flux width measurements at reactor-level magnetic fields and observation of a unified width scaling across confinement regimes in the alcator c-mod tokamak[END_REF]). We wonder whether this 1D model is able to capture such a link, by comparing the numerical output of the energy confinement time to the ITER96-th scaling law.

τ L E,th = 0.023I 0.96 B 0.03 T P -0.73 in n 0.4 M 0.2 R 1.83 -0.06 k 0.64 (5.10) where I is the plasma current in [M A], B T is the toroidal magnetic field in [T ], P is the loss power in [M W ], n is the line-averaged density, in [m -2 ], M is the average ion mass in [AM U ], R the major radius in [m] and finally and k are the inverse aspect ratio and the elongation, respectively. τ E,th is shown to decrease as P in increases, probably due to the rise of free energy available to instabilities, driving a confinement degradation. On the contrary, τ E,th is directly proportional to the plasma current and thus to B P .

In the model, we define τ E as:

τ E = W P in -dW/dt (5.11)
where W is the global plasma energy content, P in is the total heating power applied to the plasma and dW/dt is the time variation of W at steady state (dW/dt = 0).

Hence, equation 5.11 is reduced to:

τ E = 3/2 nT dV P in
, where at the numerator the local energy content is integrated over the spatial domain. 5.12(b) show the comparison between the trend of τ E obtained numerically (black line) and the scaling law (red line), when varying P in and B P respectively. The red-shaded area corresponds to the root mean square error related to the power regression of the scaling law, equal to 15.8%.

The qualitative trend of the empirical scaling law is retrieved, even though it does not constrain the model by construction.

Case with shear: preliminary modeling and results

So far, the modeling proposed in this work has referred to a low-confinement regime for plasma, meaning that the saturation of turbulent energy κ is regulated by the scaling laws which close the model itself.

In the perspective of capturing also the physical mechanisms behind the L-H transition, we introduce in the model a new element aiming at competing with turbulent energy as in prey-predator models. This element is the E × B shear, defined by:

ω E×B ≡ ∇ r v θ E×B θ = 1 B ∆ r φ (5.12)

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 87

Despite detailed origins and physical mechanisms on the L-H transition are still unclear within the fusion community, several tokamak experiments showed that it occurs when the radial electric field E r = -∇φ rapidly changes near the separatrix, causing a shearing of radial turbulent transport and finally an improvement of the confinement. During the transition, the E × B shear would overcome the rate of generation of turbulent energy. This is usually known as E × B shear suppression paradigm ([75]).

The E × B shear, which has the dimension of a frequency, is introduced in the model via the growth rates γ κ and γ ε in the source terms for κ and ε respectively, such that:

γ κ = γ I + rβ E×B ω 2 E×B ; γ ε = γ I + β E×B ω 2 E×B (5.13)
where r is a parameter of the system. The definitions provided in 5.13 represent a working assumption. Here the coupling of the E × B shear with the equations governing the plasma (equation 5.2) is provided by the vorticity equation 5.2c which estimates Ω and therefore the electric potential φ.

The idea is that in L-mode, the dominant mechanism of drive is the interchange instability γ κ/ε ∼ γ I , while approaching to the H-mode, γ κ ∼ rβ E×B ω 2 E×B and γ ε ∼ β E×B ω 2 E×B (β E×B regulates the amplitude of the shear). The transition is led by the input power P in .

We solve the one-dimensional system varying P in from 0.1M W to 8M W . In these simulations, as first step, we assumed quasi-linear-like diffusivities: D n ∝ χ e ∝ χ i ∼ (R/c s )κ. In this case one has:

χ e = R c s κ (5.14)
Since the steady-state solutions of the κ-ε system are: (0; 0) and P γκ γε 2 ; γ 3

χ e = R c s P γ κ γ ε 2 (5.15)
We recall the relation describing the equilibrium of plasma pressure in the framework of dimensional analysis (1.35) and introduce 5.15:

λ 2 SOL = L 2γ e c s χ e = L 2γ e c s R c s P γ κ γ ε 2 (5.16)
Finally, the inverse relation provides the expression of the parameter P :

P = λ scaling,L SOL γ ε γ κ c s 2γ e L R (5.17)
where λ scaling,L SOL is the empirical scaling law for the low-confinement regime discussed in [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF].

When the plasma tends to the low confinement (for relatively small input power P in ), the ratio γε γκ → γ I γ I → 1, while, increasing the input power and approaching to the H-mode, γε γκ → 1 r . In this latter case:

P = λ scaling,L SOL 1 r c s 2γ e L R (5.18)
According to [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF] and [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF], the empirical scaling law of the SOL width in highconfinement regime λ scaling,H SOL is smaller than the one in low-confinement by a factor 2. Therefore, by setting r = 2, we retrieve immediately that P scales as 1 2 λ scaling,L SOL = λ scaling,H SOL . Since the diffusivity is governed by the parameter P containing the scaling law, we conclude that the distribution of diffusivity in our model is sensitive to a possible change of plasma regime. This is indeed observed in numerical results, showed below.

𝑟 [𝑚]

(a) Radial profile of electron temperature for increasing input power. These simulations include the physics of the E × B shear.

𝑟 [𝑚]

(b) Radial profile of plasma diffusivity for increasing input power. These simulations include the physics of the E × B shear. Some key features can be highlighted:

• In figure 5.13(a) the radial profile of electron temperature at steady state denotes a pedestal (namely a sharp increase localized at the separatrix) from the threshold power P in = 2.5M W on.

• The radial profile of diffusivity, on the other hand, displays a graduallyincreasing drop for corresponding values of input (see figure 5.13(b)). In particular, at P in = 6M W , a deep drop develops at the separatrix, representing a region where the transport of plasma is highly hampered, since here the transport coefficient undergoes a reduction of more than two orders of magnitude. This feature recalls the mechanism played by transport barriers,

1D parametric study of the new model and comparison with empirical scaling laws for the SOL width and the energy confinement time 89

regions where the radial turbulent transport is hindered by the E × B shear ( [START_REF] Diamond | Zonal flows in plasma-a review[END_REF]).

Interestingly, the trend of τ E against input power compared to the ITER96-th scaling law, represented in figure 5.14, bifurcates at the threshold power P in = 2.5M W corresponding to the arising of the pedestal and the drop of D n , suggesting a transition to a high confinement regime. Future goals require to upgrade the model with E × B shear, by using the eddyviscosity-like diffusivities. Therefore, the SOL width and the energy confinement time should be investigated, as done in the case with no-shear, and compared against H-mode scaling laws.

Conclusions

This Ph.D. work contributes in developing a new modeling of turbulent cross-field fluxes in edge diverted plasma with the goal to improve the predictive capabilities of the transport code SolEdge2D-EIRENE. To achieve such a result, we bridged the concepts of eddy viscosity and κ-ε model -widely used in Computational Fluid Dynamics to model turbulence semi-empirically in neutral flows -to the quasi-linear theory and the physics of interchange-instability-driven turbulence exploited in the framework of reduced models for turbulent plasma.

We recall that two-dimensional transport codes based on Braginskii fluid equations are crucial to support steady-state operations in the next-generation fusion device ITER because of their reasonable computational cost. Among them, we exploited the flexibility of SolEdge2D-EIRENE, which allows to simulate versatile realistic magnetic configurations and plasma-wall interactions.

Transport codes rely on ad-hoc diffusion coefficients for the estimation of steadystate cross-field fluxes. These latter can be provided a posteriori from experimental data on plasma density and temperature acquired from the machine diagnostics after each discharge. We showed in chapter 1 that this proceeding can be optimized by means of an automatic fitting procedure implemented in SolEdge2D-EIRENE, which however is constrained by the specific experimental observations, thus degrading the predictive capabilities of the code.

Motivated by this issue, in chapter 4 we proposed a step-by-step modeling of plasma turbulent cross-field fluxes where two transport equations are derived algebraically, one for the turbulent kinetic energy and the other for its dissipation rate, inspired by the κ-ε model adopted in Computational Fluid Dynamics. We included a production term based on the linear interchange instability. These additional equations are coupled to the Braginskii equations solved by SolEdge2D-EIRENE through the modeling of the cross-field fluxes of density, momentum and energy. We assumed for these latter a diffusive expression where the transport coefficients vary with the local value of turbulent kinetic energy and its dissipation rate at the equilibrium, consistently with the eddy viscosity concept outlined in chapter 2.

We showed in chapter 3 that similarly to the eddy viscosity assumption, also in the framework of plasma physics transport coefficients can be related to reference turbulent scales of velocity and time according to the so-called quasi-linear theory, leading to D QL = κτ c . While the standard approach allows to find a linear dependence between transport coefficients and the turbulent energy κ, thus providing the reference velocity scale, it does not predict a closure for the time scale τ c . We proposed a qualitative discussion where we argued that it is possible to infer additional information in the expression of such a time scale, concerning the dissipation of turbulent energy, using the dispersion relations of a κ-ε-like system. In particular, we found that τ c ∝ ( κ ε ) -1 and for ε → 0 one recovers the standard quasi-linear analysis. Compared to this latter, where the time scale is generally fitted with gyrokinetic simulations [START_REF] Bourdelle | Turbulent transport in tokamak plasmas: bridging theory and experiment[END_REF], the method that we propose is a first-step towards a more predictable estimation of diffusivity in plasma turbulent fluxes.

In Chapter 4, the step-by-step modeling has been formally outlined. In the first step, the model consists in one equation for the turbulent energy κ, which is coupled with the Braginskii set of equations. A steady state solution for κ allows to estimate the transport coefficient according to the standard quasi-linear approach. The free parameters are closed by an improved automatic fitting procedure implemented in SolEdge2D-EIRENE, which fits the numerical plasma profiles at the outer midplane with diagnostics data of a single discharge. The comparison between numerical and experimental profiles of the transport coefficient at the outer midplane shows a good accordance. The model's free parameters still need to be closed by experimental data available from a single discharge.

To add predictive capabilities to the model, in the second step the one-equation model is closed by the empirical scaling law for λ SOL proposed by [START_REF] Goldston | Heuristic drift-based model of the power scrape-off width in low-gaspuff h-mode tokamaks[END_REF] and [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF]. Transport coefficients are evaluated self-consistently according to the standard quasi-linear approach. Three discharges from the TCV machine are referred to as benchmark, where the main plasma is kept in the same conditions while the outer divertor leg length is varied. We showed that the kind of closure that we adopted improves the predictive capabilities of the model. Indeed, the model has been able to reproduce two-dimensional distributions of plasma quantities in the edge and near SOL, just providing the boundary conditions of plasma in the core. At the outer midplane, radial profiles of density and temperature nearly overlap, successfully reproducing the fact that the main plasma is unchanged when the magnetic configuration is varied. Moreover, these profiles compare very well with the experimental data in the edge region.

The two-dimensional distribution of density transport coefficient and flux is also calculated self-consistently for the three magnetic configurations (assuming that the thermal conductivity and the viscosity are proportional to the particle diffusivity). For all of them, the ballooning feature is observed at the low-field side as a consequence of the interchange-instability-driven turbulence on plasma mean fields.

Simulations show that the model is also able to capture that this enhanced transport broadens up to the common flux region, but not in the private flux region, consistently with turbulent three-dimensional isothermal simulations made with TOKAM3X in TCV-like geometry and Compass-like geometry. This picture is consistent also with experimental observations in the MAST tokamak.

The estimation of the SOL width is also consistent with the observed range of values, however, λ SOL does not scale as the outer divertor leg length as found experimentally. The reason is that the three TCV discharges are not fitted by the scaling law used to close the model. Indeed, numerical results show λ SOL is nearly unchanged when varying the outer divertor leg length, while experimentally it linearly increases with the latter.

Finally, in the third step of the new modeling, two equations are outlined describing the evolution of the turbulent energy and its dissipation rate, exploiting the same empirical scaling law as closure. This time, the expression of the transport coefficient scales as both the turbulent energy and its dissipation rate, consistently with the eddy-viscosity concept and the qualitative observations on the extended quasi-linear analysis discussed in chapter 3. Numerical simulations refer to an early WEST discharge where experimental observations exhibit a certain level of incertitude. Preliminary results show a good accordance between the plasma density profile and diagnostics observations at the separatrix. Also, at the divertor the estimation of the SOL width is shown to be in quantitative agreement with measures of heat flux in the main plasma obtained with reciprocating Langmuir probes across multiple shots. The two-dimensional distribution of the particle transport coefficient denotes a sharp ballooning at the LFS, especially in the edge region. Such a result needs however a further investigation, due to the evidence of enhanced transport in the X-point.

We showed in chapter 5 that a one-dimensional reduction of the new model allows to perform fast simulations where we scan the SOL width λ SOL and the energy confinement time τ E at steady-state against crucial machine parameters such as the input power and the poloidal magnetic field. The strength of such approach is a reliable and fast verification of numerical results provided by the new model across multi-machine simulations against widely-used scaling laws for λ SOL and τ E . The comparison of such scans with the predictions of respective scaling laws is showed to be favorable.

Future improvements of the model can be implemented by exploiting its flexibility. Indeed, the expression of the turbulence growth rate can be easily improved to take into account different kinds of instability, provided the relation between its growth rate and mean plasma quantities. Moreover, in future simulations we should perform multi-machine simulations and account for the multi-fluid capability of SolEdge2D-EIRENE by simulating impurities in plasma discharges.

We think that this new modeling may also be adapted to reproduce improved confinement plasma regimes by introducing a predator-like mechanism for turbulent energy within the equations, provided by the E × B shear flow. First promising efforts have been done to this aim, as outlined in chapter 5. Here, simulations with increasing input power show that when the λ SOL bifurcates towards the scaling law proposed by [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] for H-mode plasma, temperature radial profile exhibits a pedestal-like behavior at the separatrix, where the radial transport coefficient is found to decrease significantly.

B B.1 Derivation of a general "κ" equation for magnetized plasma

The derivation of a "κ" equation for magnetized plasmas is similar to the one adopted for neutral fluids, namely the momentum equation is manipulated in order to obtain transport equation for the turbulent kinetic energy κ.

In the most general form, when assuming incompressible flow and the Boussinesq's assumption for the density, the momentum equation reads:

mn 0 ∂u i ∂t + mn 0 u l ∂u i ∂x = - ∂P ,i ∂x + F ,i + R ,i (B.1)
where i, j and are the space component of the velocity vector u, P ,i is the pressure tensor and finally F ,i and R ,i are generic terms representing an external force (averaged over the distribution function f ) and the collisional term respectively.

In neutral fluids, F ,i = F i = mn 0 g i is the buoyancy force, while R ,i = ∂τ ,i ∂x , where τ ,i is the stress tensor heuristically introduced in the Euler's equation in order to take into account for irreversible exchange of momentum due to mutual friction of fluid tubes.

In plasmas, F ,i is the electro-magnetic force and R ,i is the collisional operator. Hereafter, we will keep handling the both in a generic fashion.

First We introduce Π ij = P ij /(mn 0 ), mn 0 Z i = R i + F i f and U ij = u i Z j so that:

∂ t u i + u ∂u i ∂x + ∂Π i ∂x = Z i (B.3a) ∂ t u i u j + u ∂ u i u j ∂x + u j ∂Π i ∂x + u i ∂Π j ∂x = U ij + U ji (B.3b)
The mean part for this system is then given by:

∂ t u i + u ∂u i ∂x + ∂ u u i ∂x + ∂Π i ∂x = Z i (B.4a) ∂ t u i u j + u ∂ u i u j ∂x + u j ∂Π i ∂x + u i ∂Π j ∂x + ∂ t u i u j + ∂ u u i u j ∂x + u ∂ u i u j ∂x + ∂ u u j u i ∂x + ∂ u u j u i ∂x + u j ∂ Π i ∂x + u i ∂ Π j ∂x = U ij + U ji (B.4b)
The equations for the fluctuating part are then obtained.

∂ t u i + ∂ u u i -u u i ∂x + u ∂u i ∂x + u ∂ u i ∂x + ∂ Π i ∂x = Z i (B.5a) ∂ t u i u j -u i u j + ∂ u u i u j -u u i u j ∂x + u ∂ u i u j -u i u j ∂x + ∂ u u i u j -u u i u j ∂x + ∂ u u i u j -u u i u j ∂x + u j ∂ Π i ∂x -u j ∂ Π i ∂x + u i ∂ Π j ∂x -u i ∂ Π j ∂x = U ij + U ji (B.5b)
As done in the derivation of the fluid equations the term ∂ t (u i u j ) in B.4b can be determined using B.4a. In the Favre averaging approach let us define the average nu = n 0 V f b . Here V f b is the velocity governed by the force balance, therefore:

u
∂ t V f b + 1 m i n 0 ∇ p + 1 n 0 ∇• nu ⊗ u -ν∇ 2 V f b -Θ g x = 0 (C.2b)
For the sake of simplicity in the present work we shall consider the Boussinesq approximation, namely that one can approximate n by n 0 , so that u = V f b and u = V f b + u, and consequently:

∂ t V f b + ∇• V f b ⊗ V f b -ν∇ 2 V f b = - 1 m i n 0 ∇ p + Θ g x -∇• u ⊗ u (C.2c)
Given the relationship between the flow and the stream function one has V f b = z × ∇ ψ . Since ψ only depends on x, one finds therefore that: V f b = y∇ x ψ . There are thus constraints governed by the underlying symmetry as readily found with the projections in the x and y directions respectively.

∂ t V f b -ν∇ 2 V f b = -∇• u u y (C.3a) ∇• u u x = - 1 m i n 0 ∇ x p + Θ g (C.3b)
To obtain equation: y-projection average momentum balance v1, we have taken into account that V f b •∇ = V f b ∇ y and that the averaged fields V f b and p only depends on x. As a consequence, the only source term for V f b is the kinematic pressure tensor governed by the velocity fluctuations u. Furthermore, the velocity V f b does not contribute to the force balance in the x direction, equation: x-projection average momentum balance v1. To complete equation: average momentum balance v1, one must consider the equation governing the fluctuating velocity u = z × ∇ ψ:

∂ t u + 1 m i n 0 ∇ p + ∇• u ⊗ u -ν∇ 2 u -Θg x + V f b ∇ y u + u x ∇ x V f b = ∇• u ⊗ u (C.4)
The kinematic pressure tensor is the non-linear contribution to these equations that will govern the mode coupling and the development of turbulence.

C.1.2.2 Mean and fluctuating Rayleigh-Bénard equations

We now consider the separation in mean and fluctuating part in the system equation: momentum conservation v5 coupled to equation: temperature balance v5

∂ t Θ + ψ, Θ + ∇ x ψ ∂ y Θ -∇ x Θ ∂ y ψ -κ∇ 2 Θ = S (C.5a) ∂ t W + ψ, W -∇ x W ∂ y ψ + ∇ x ψ∂ y W -ν∇ 2 W + g ∂ y Θ = 0 (C.5b) C.
107

In this system we have allowed the temperature, stream function and vorticity fields to be split into the previous terms and a linear contribution in x, hence replacing Θ by Θ + x∇ x Θ where ∇ x Θ is a constant parameter and ψ by ψ + x∇ x ψ, ∇ x ψ being a constant parameter and W by W + x∇ x W , ∇ x W being a constant parameter. As such the system is not quite consistent since we have ignored here the relationship between ψ and W , W = ∇ 2 x ψ. These terms provide source terms to the system that allow one removing S from equation: temperature field v0. In this framework one can split the fields into their average and fluctuations with respect to the y-average. ∂ t Θ + ψ, Θ -κ∇ 2 Θ = 0 (C.6a)

∂ t Θ + ψ, Θ -ψ, Θ -κ∇ 2 Θ + ∇ x ψ + ∇ x ψ ∂ y Θ -∇ x Θ + ∇ x Θ ∂ y ψ = 0 (C.6b) ∂ t W + ψ, W -ν∇ 2 W = 0 (C.6c) ∂ t W + ψ, W -ψ, W -ν∇ 2 W + g ∂ y Θ + ∇ x ψ + ∇ x ψ ∂ y W -∇ x W + ∇ x W ∂ y ψ = 0 (C.6d)
On finds here that the terms ∇ x f and ∇ x f have similar contributions. At this stage their separation can appear to be be artificial. However, when stepping to the Fourier modes, the former will govern a mode coupling while the latter does not. Introducing a difference between these two terms is therefore justified in that particular framework.

C.1.2.3 Kinematic pressure tensor

We consider both the Eulerian equation: momentum conservation and Lagrangian forms of the evolution equation of u.

∂ t m i nu + ∇• pI + m i nu ⊗ u -ν∇ 2 m i nu -m i n 0 Θg x = 0 (C.7a) m i n ∂ t u + u•∇ u + ∇p -ν∇ 2 m i nu -m i n 0 Θg x = 0 (C.7b)
Given the identity: The latter can also be written as:

m i ∂ t nu ⊗ u = ∂ t nu ⊗ u + nu ⊗ ∂ t u (C.
m i n∂ t u ⊗ u + ∇p + m i n u•∇ u -ν∇ 2 m i nu -m i n 0 Θg x ⊗ u + u ⊗ ∇p + m i n u•∇ u -ν∇ 2 m i nu -m i n 0 Θg x = 0 (C.8c)
where we have used the particle balance equation ∂ t n + ∇• nu = 0.

C.1.3 Equations in Fourier space

Let us now introduce the projection on Fourier modes, both in the x and y directions, using the position vector r determined by (x, y) and the wave vector k given by (k x , k y ). The structure of these equations is therefore an evolution of a field, either Θ or W , which includes damping and drift, a linear coupling term to the other field, and non-linear coupling to the stream function ψ, keeping in mind that stream function and vorticity are coupled by the linear Poisson relation equation: Fourier vorticity potential v0. This system can be simplified in the case ∇ x ψ = 0 and ∇ x W = 0, a case that is mostly considered in the literature.

T Θ,Ψ (k, k ) = z • k × k (C.12a) T Ψ,Ψ (k, k ) = k -k 2 -k 2 k 2 z • k × k (C.12b)
where we have used the identity k x k y -k y k x = z • (k × k ).

C.1.4 Stepping to the Lorenz model

We reduce the model by considering a reference mode Θ (kx,ky) and a correction to the mean temperature gradient of the form Θ (2kx,0) , hence a given mode and its impact on transport determined by the second term.

∂ t Θ (kx,ky) = ik y ∇ x Θ Ψ (kx,ky) -κk 2 Θ (kx,ky) -2k x k y Ψ (-kx,ky) Θ (2kx,0) (C.13a)

This equation introduces a coupling to both Ψ (kx,ky) and Ψ (-kx,ky) . Neglecting higher harmonics for Ψ reduces the evolution of these two modes to the linear contributions.

∂ t Ψ (kx,ky) = i k y k 2 gΘ (kx,ky) -νk 2 Ψ (kx,ky) (C.13b)

In the framework of the Lorenz model a symmetric boundary condition with respect to x = 0 is imposed so that Θ (-kx,ky) = Θ (kx,ky) and Ψ (-kx,ky) = Ψ (kx,ky) . Furthermore, the structure of the mode coupling highlighted by the two previous equations, then ensures that one can consider the modes Θ (kx,ky) to be imaginary while the modes Ψ (kx,ky) are real, hence: Ψ (kx,ky) = Ψ (-kx,-ky) = Ψ (-kx,ky) = Ψ (kx,-ky) Θ (kx,ky) = -Θ (-kx,-ky) = Θ (-kx,ky) = -Θ (kx,-ky)

One can then modify equation: temperature main mode kx ky v1 to obtain:

∂ t θ 1 = -k y ∇ x Θ ψ 1 -κk 2 θ 1 -2k x k y ψ 1 θ 2 (C.15a) ∂ t ψ 1 = k y k 2 gθ 1 -νk 2 ψ 1 (C.15b)
with the definitions θ 1 = iΘ (kx,ky) , θ 2 = iΘ (2kx,0) and ψ 1 = Ψ (kx,ky) . This system is only complete provided one considers the evolution equation of Θ (2kx,0) .

∂ t Θ (2kx,0) = -κ 4k 2 x Θ (2kx,0) + 2k x k y Ψ (kx,-ky) Θ (kx,ky) (C.16a) and therefore the evolution equation for θ 2 :

∂ t θ 2 = -κ 4k 2 x θ 2 + 2k x k y ψ 1 θ 1 (C.16b)
We then consider as reference time the inverse of κk 2 , normalise the stream function by κk 2 /(2k x k y ) an the temperature field by Θ such that Θ = (νk 2 /(k y g), one can 110 C.1. The 2D Rayleigh-Bénard system then introduce the Rayleigh dimensionless parameter R a = -k 2 y g∇ x Θ/(k 2 νk 2 κk 2 ) and the Prandtl dimensionless parameter P r = ν/κ, one then obtains the normalised system:

∂ t θ 1 = R a ψ 1 -θ 1 -ψ 1 θ 2 (C.17a) ∂ t ψ 1 = P r θ 1 -ψ 1 (C.17b) ∂ t θ 2 = - 4k 2 x k 2 θ 2 + ψ 1 θ 1 (C.17c)
In this system the drive is due to the term proportional to R a , the Rayleigh parameter being proportional to the external drive by the imposed temperature gradient, The Prandtl number is a measure of the relative importance of the damping terms, viscosity versus heat diffusion, and the third parameter characterises the geometry according to the choice of the relevant k 2 x and k 2 .

C.1.5 Bifurcation and thermoconvection

In this subsection, the details about the stability analysis concerning the Lorenz model are discussed, in particular when the non-trivial fixed points are concerned. As mentioned in chapter 3, for ψ = θ and ψ 2 = α 2 θ 2 = α 2 (Ra -1), the dispersion relation relates: Since the first three coefficients of the eigen value equation C.18 are real and positive, this function intersects the real axes once and there is a real root γ R and two complex conjugate roots γ C and γ C * . Furthermore, one finds the product of the three roots γ R |γ C | 2 is eqaul to -2α 2 Pr(Ra -1) < 0 which implies that γ R < 0. The stability of the fixed point then depends on on the sign of the real part of the imaginary roots, given that these are complex conjugate. A change of sign of this real part implies that there exist a value of the control parameter Ra such that γ = iω + ε with ω and ε real. We consider the case close to marginality, hence such that ε ω, in order to expand equation C.18 for small ε. To facilitate the calculation, it is convenient to reformulate the dispersion relation as:

γ 3 + γ 2
γ 3 + Aγ 2 + Bγ + C = 0 (C.19)
The equation for ω and ε is then:

-iω 3 -3ω 2 ε -Aω 2 + 2Aiωε + iBω + Bε + C = 0 (C.20)

We can split this equation into two independent equations, one for the real part and one for the imaginary part, obtaining: Finally, for Pr > (1 + α 2 ), ε > 0, if Ra > Ra * , three cases arise:

• 0 < Ra < 1, there is only a single fixed point such that ψ = 0, θ = 0 and θ 2 = 0;

• 1 < Ra < R * , the fixed point ψ = 0 is unstable and there exist two stable fixed points with ψ = 0;

• Ra > R * : all fixed points are unstable.

The other regime is for Pr < (1 + α 2 ), such that:

• 0 < Ra < 1, there is only a single stable fixed point such that ψ = 0;

• Ra > 1, the fixed point ψ = 0 becomes unstable and there exists two fixed points with ψ = 0 that remain stable for all values of Ra.

The case ω = 0, which is the second solution of equation C.21, leads to:

ε = C B = 2Pr(Ra -1) Pr + Ra (C.29)
Thus it describes the marginality state of the solution ψ = 0, hence ε > 0 for Ra > 1.

Figure 1 . 1 :

 11 Figure 1.1: The reaction rate against temperature for the three most efficient fusion reactions (D-T,D-D,D-He3). Figure from Wikipedia web page.

Figure 1 . 2 :

 12 Figure 1.2: Cartoon of the axisymmetric geometry of a tokamak. The major and minor radius R and a are displayed in red, ϕ and θ are the toroidal and poloidal angles while B T (in blue), B P (in yellow) and B tot (in green) are the toroidal, poloidal and total magnetic fields respectively. Finally, I plasma is the plasma current induced by the central solenoid, flowing in the toroidal direction, while I T is the toroidal current flowing in the poloidal plane.

. 3 .Figure 1 . 3 :

 313 Figure 1.3: Shape of B T and B P varying with R. The HFS is the region where the magnetic field is more intense, the LFS the one where it is the weaker.

Figure 1 . 4 :

 14 Figure 1.4: Cartoon of a cross-section of the plasma in a toroidal limiter magnetic configuration. The separatrix is outlined by a dotted line.

Figure 1 . 5 :

 15 Figure 1.5: Cartoon of the divertor magnetic configuration in the single-null variant. Picture from [2].

Figure 1 . 6 :

 16 Figure 1.6: Cartoon representing the sheath region.

Figure 1 . 7 :

 17 Figure 1.7: E ×B drift and ∇B drift causing horizontal and vertical motion respectively.

Figure 1 . 8 :

 18 Figure 1.8: Cartoon of the interchange instability mechanism driving turbulence at the LFS of the tokamaks, left panel. In the right panel, a zoom of the LFS region, where convection cells with iso-potential contour lines (in yellow) are created by the perturbation of the potential. The black arrows indicate the direction of the E × B and ∇B drifts.

Figure 1 . 9 :

 19 Figure 1.9: Poloidal maps of plasma density and electron temperature calculated by SolEdge2D-EIRENE for a WEST test case.

1. 6 .

 6 Perpendicular turbulent transport in SolEdge2D

Figure 1 . 10 :

 110 Figure 1.10: Example of mesh decomposition for a WEST magnetic equilibrium. Each sub-domain is characterized by a different color.

Figure 1 .

 1 Figure 1.11:The scheme of the feedback loop procedure implemented in SolEdge2D-EIRENE is represented as a control system.∂ r n Xp , provided as input, is the experimental profile while ∂ r n simu is the one provided by the simulation. is the error between these two. K and 1 τ i are the transfer functions of this scheme and finally D ⊥ is the output.

22

 22

Figure 2 . 1 :

 21 Figure 2.1: Distortion of a fluid element in a parallel shear flow.

Figure 2 . 2 :

 22 Figure 2.2: Log-log sketch of the energy spectrum in direct cascade. Figure from [27].

Figure 2 . 3 :

 23 Figure 2.3: Log-log sketch of the energy and enstrophy spectrum in inverse cascade. Figure from [27].

Figure 2 . 4 :

 24 Figure 2.4: Comparison of solving capabilities in RANS, LES and DNS approaches [38].

2. 2 .Figure 2 . 5 :

 225 Figure 2.5: Comparison of DNS, LES and RANS simulations for two streams with unequal scalar concentrations [37].

Figure 2 . 6 :

 26 Figure 2.6: Time history of the fluctuating velocity at a point in a turbulent flow. Picture from [40].

34 2. 4 .

 344 Turbulence models based on one-point closure this model, the second variable is the energy dissipation rate, ε[m 2 /s 3 ]:

Contents 3 . 1 38 3. 1

 31381 The quasi-linear approach . . . . . . . . . . . . . . . . . 38 3.1.1 Beyond the quasi-linear approach . . . . . . . . . . . . . 41 3.2 Prey-predator models . . . . . . . . . . . . . . . . . . . . 43 3.2.1 An example of a prey-predator model in plasma . . . . 44 3.3 The Lorenz system . . . . . . . . . . . . . . . . . . . . . 46 3.3.1 Equations and physical meaning . . . . . . . . . . . . . 46 3.3.2 Fixed point analysis . . . . . . . . . . . . . . . . . . . . 47

  (a) The periodic loss-win dynamics between the prey(in blue) and the predator (in green). (b) Level curves of the conserved quantity. In this zero-dimensional simulation, α = 1, δ = 5 and β = γ = 5. Equilibrium points are marked by a star in (0, 0) and (1, 0.5).

Figure 3 . 1 :

 31 Figure 3.1: The prey-predator model.

Figure 3 . 2 :

 32 Figure 3.2:Plots from[START_REF] Miki | Spatio-temporal evolution of the L-I-H transition[END_REF]. On the left, time trace of turbulent energy (in blue), zonal flows (in green) and mean flows (in red). On the right, the phase-portrait of the predator-prey system. For both, a portion of the total time range is considered.

. 3 .Figure 3 . 3 :

 333 Figure 3.3: 2D fluid layer of width H, uniformly warmed from below, where temperature θ = θ H , and cooled from above, where θ = θ C .

2 Figure 3 . 4 :

 234 Figure 3.4: Time evolution of ψ, Θ and Θ 2 in a given time range.

Figure 3 .Figure 3 . 5 :

 335 Figure 3.5: Lorenz attractor when α 2 = 8/3, Ra = 28 and P r = 10. Moreover, x, y and z correspond to ψ, θ and θ 2 in the notation used in 3.36.

  Lorenz attractor when α 2 = 8/3, 1 < Ra = 10 < Ra * and P r = 10. The fixed point ψ = 0 is unstable and there exist two stable fixed points with ψ = 0.

  Lorenz attractor when α 2 = 8/3, 0 < Ra = 0.1 < 1 and P r = 10. There is only a single fixed point.

Figure 3 . 6 A 4 . 1

 3641 Figure 3.6

Figure 4 . 1 :

 41 Figure 4.1: Strategy of the modeling. Choice of the model (one or two equations, blue boxes), the experimental test-cases to use as benchmarks for the numerical results (pink boxes) and the closure (single experiment or scaling law, green boxes).

. 5 ) 4 .

 54 A new modeling for cross-field turbulent transport of diverted plasma 57 Equation 4.5 allows to determine the spatial steady-state distribution of the cross-field transport in the poloidal plane.β and ∆ω are closed iteratively, by tuning the radial profile provided by 4.5 to the one inferred by experimental data using the automatic fitting procedure discussed in chapter 1. Such a method is well outlined in the scheme of figure 4.2.Model D ⊥,𝜅

Figure 4 . 2 :

 42 Figure 4.2: Scheme of the one-equation model implemented in the feedback-loopprocedure. The green loop is the feedback loop analyzed in chapter 1, while the red one is the new loop integrating the one-equation model.

  The radial profile of D ⊥,κ is compared to D 1D ⊥,f b in figure4.[START_REF] Wagner | Transport in toroidal devices-the experimentalist's view[END_REF], where α 0 = ∆ω -1 and ζ = β/∆ω are more convenient definitions for the free parameters. The errorErr(α 0 , ζ) = D ⊥,κ -D 1D⊥,f b is a function of α 0 and ζ; therefore the minimization of Err(α 0 , ζ) leads to explicit values of the free parameters and finally to close the model. The contour plot in 4.3 shows the result of this optimization: the lighter the color of contour regions, the smaller the error.Optimized α 0 and ζ re-enter the model which in turn defines the full twodimensional distribution of the cross-field transport. The advantage of such an approach compared to the feedback-loop procedure seen in chapter 1 consists in the gain of the poloidal dependence of the diffusivity calculated by the model, provided by the modeling of the interchange instability.

Figure 4 . 3 :

 43 Figure 4.3: The optimization of the error Err(α 0 , ζ) allows to find uniquely the values of the free parameters α 0 and ζ and close the model.

Figure 4 . 4 :

 44 Figure 4.4: Cross section sketch of the TCV vacuum vessel with main diagnostics. Also, the magnetic separatrix for the three outer-divertor leg configurations presented in section 4.3: the long leg, in orange, the medium leg, in blue and the short leg, in green.

Figure 4 . 5 :

 45 Figure 4.5: Numerical profile of plasma diffusivity at the outer midplane (continuous blue line) compared to the feedback-loop diffusivity (continuous black line).

60 4. 2 .

 602 First step: one-equation model and tuning with TCV data 4.5 enforces a precise profile of the radial diffusivity, while ∆ω only adjusts its amplitude. Despite the error, the advantage of such modeling is to provide a first-attempt profile of the transport coefficient calculated self-consistently.The profiles of density and electron temperature at the outer midplane are compared to the fit of the HRTS data (see figure4.6), resulting in a decent agreement of the electron temperature at the separatrix. Although the free parameters of the model are closed to match as much as possible the experimental data at the outer midplane, the model profiles of plasma density and temperature do not overlap the experimental fits, due to the relative error between diffusivity profiles mentioned for figure 4.5.

Figure 4 .Figure 4 . 6 :

 446 Figure 4.6: Radial profiles of density and temperature at the outer midplane, compared to HRTS data. In the density plot (a), the black continuous line denotes the fit of experimental acquisitions. Note that in (b) the agreement between numerical and fitted profiles of the electron temperature is nice in the neighborhood of the separatrix (R = R sep ).

Figure 4 . 7 :

 47 Figure 4.7: Distribution of self-consistently calculated plasma diffusivity in the poloidal plane in TCV geometry.The high value of diffusivity in the common flux region is thought to be due to numerical diffusion. In the far SOL, where the gradient of pressure is expected to be low, the diffusivity is large because here the effect of the model closure is dominant.

(4. 6 )Figure 4 . 8 :

 648 Figure 4.8: Outer target profiles of density, electron temperature, saturation current and heat flux perpendicular to the divertor.

  Density [𝑚 -3 ] 𝑅 -𝑅 𝑠𝑒𝑝 [𝑚] 𝑇 𝑒 [𝑒𝑉] * HRTS short * HRTS medium * HRTS long

Figure 4 . 9 :

 49 Figure 4.9: Radial profiles of plasma density (left panel) and electron temperature (right panel) at the outer midplane, outlined by continuous lines, for the three magnetic configurations. These numerical results, obtained with SolEdge2D-EIRENE, are compared to high-resolution Thomson scattering data (sketched in the small cartoon of TCV vessel).

  𝑚 2 /𝑠] 𝐷 𝑛 [𝑚 2 /𝑠] 𝐷 𝑛 [𝑚 2 /𝑠]

Figure 4 . 10 :

 410 Figure 4.10: Particle diffusivity D n self-consistently calculated on the poloidal plane for the three magnetic configurations. Note the ballooning and the existence of enhanced transport in the common flux region.

Figure 4 . 11 :

 411 Figure 4.11: Particle flux Γ n,⊥ = -D n ∇n self-consistently calculated on the poloidal plane by SolEdge2D-EIRENE for the three magnetic configurations.

Figure 4 . 12 :

 412 Figure 4.12: Estimation of the 2D map of Γ n,⊥ from TOKAM3X turbulent simulations (i.e. with fluctuating components).

Figure 4 . 13 :Figure 4 . 14 :

 413414 Figure 4.13: Perpendicular heat flux remapped at the outer midplane for the three magnetic configurations and compared to infrared data. Note that both the absolute magnitude and the scaling are in really good agreement with experimental measurements.

70 4. 4 . 2

 7042 Third step: two-equation model, tuning with experimental scaling law System 4.22 exhibits two solutions: the trivial one (0; 0) and P γκ γε

Figure 4 . 15 :

 415 Figure 4.15: Radial profile of density at the outer midplane (blue curve), compared to two sets of experimental data: the interfero-polarimentry (balck crosses) and the reciprocating Langmuir probes (black circles).

Figure 4 . 16 :

 416 Figure 4.16: Poloidal distribution of the particle diffusivity self-consistently calculated by the κ-ε model coupled to SolEdge2D-EIRENE.

Figure 4 . 17 :

 417 Figure 4.17: Electron density n, parallel density current j sat, and parallel heat flux (q ) at the outer divertor and remapped at the outer midplane (blue lines) compared to experimental measurements from Langmuir Probes (red dots).

  𝑞 ⊥ [𝑀𝑊/𝑚 2 ] 𝑅 -𝑅 𝑠𝑒𝑝 [𝑚] 𝑞 ∥ [𝑀𝑊/𝑚 2 ] 𝝀 𝒒 ∼ 𝟗 𝒎𝒎 𝑍 [𝑚]

Figure 4 . 18 :

 418 Figure 4.18: Heat flux perpendicular to the outer divertor plate and remapped at the outer midplane (left panel). In the right, measures of heat flux in the main plasma obtained with reciprocating Langmuir probes across multiple shots (picture from [72]).

Figure 5 . 1 :

 51 Figure 5.1: Radial profiles of plasma source S n (in magenta) and neutral source I n0 (in grey) at steady state, for a simulation run according to the reference parameters of 5.1 and input power P in = 1M W .

3 ]Figure 5 . 2 :

 352 Figure 5.2: Radial profiles of neutrals density for given plasma parameters.

) 2 ( 5 . 5 )Figure 5 . 3 : 2 ε κ 2

 2555322 Figure 5.3: Plot of the energy source S E against the radial coordinate for arbitrary amplitude P in . Note that it is localized in the core.

Figure 5 . 4 :Figure 5 . 5 :

 5455 Figure 5.4: Cartoon of the cylindrical system of coordinates used for 5.2. The red line represents the spatial domain of the system .

5 .Figure 5 . 6 :

 556 Figure 5.6: Radial profiles of the numerical density (continuous red line) vs the analytical one (blue squares).

.

Figure 5 . 7 :

 57 Figure5.7: L 2 norm relative error between numerical and analytical density is in the order of 2% in the worst case (higher average space step h ave ).

.

Figure 5 . 8 :

 58 Figure 5.8: One-dimensional profiles for increasing input power.

2

 2 SOL ∝ χαL √ Tsep would rather decrease, diverging from the trend of the scaling law.

5. 4 .

 4 Comparison with experimental scaling laws𝑃 [𝑀𝑊] Scan of the λ SOL calculated numerically vs. the input power, compared to the experimental scaling law for L-mode discharges discussed in[START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF]. Scan of the λ SOL calculated numerically vs. B P ,compared to the experimental scaling law for L-mode discharges.

Figure 5 . 9 :

 59 Figure 5.9: Results from the model compared to the scaling laws when varying two different control parameters.

  Zoom of pressure profile in the neighborhood of the separatrix. The grey-shaded area is the 95% confidential interval of the regression.

Figure 5 . 10 :

 510 Figure 5.10: Radial pressure profile for the simulation with P in = 1M W .

5 .Figure 5 . 11 :

 5511 Figure 5.11: Study of the sensitivity of the plasma pressure width in the SOL compared to the major radiu R (left panel) and the toroidal magnetic field B T . Note that for this latter, the model denotes nearly no dependence, as expected in the scaling law of equation 5.9

  Figures 5.12

  Scan of the energy confinement time τ E vs the input power (continuous black line), compared to the experimental scaling law for the thermal energy confinement time (continuous red line) discussed in[START_REF] Transport | Plasma confinement and transport[END_REF]. Scan of the energy confinement time τ E vs B P (continuous black line), compared to the experimental scaling law for the thermal energy confinement time (continuous red line).

Figure 5 . 12 :

 512 Figure 5.12: Numerical results compared to scaling laws when varying two different control parameters.

2 (see 4 .

 24 21 and 4.22), the resulting diffusivity reads:

Figure 5 . 13 :

 513 Figure 5.13: Radial profiles of T e and D n for increasing P in .

Figure 5 . 14 :

 514 Figure 5.14: Energy confinement time calculated numerically (continuous black line) and compared to the ITER96-th scaling law (continuous red line). A bifurcation arises in P in = 2.5M W , which may suggest a transition towards an improved confinement regime of plasma.

  8a) one can obtain the evolution equation for the tensor u ⊗ u:m i ∂ t nu ⊗ u = -∇p + ν∇ 2 m i nu + m i n 0 Θg x ⊗ u -m i n u•∇ u + m i u∇• nu ⊗ u + u ⊗ -m i n u•∇ u -∇p + ν∇ 2 m i nu + m i n 0 Θg x (C.8b) 108 C.1. The 2D Rayleigh-Bénard system

Ψ

  k e ik•r (C.9b) Since the two fields are real one has the standard relations Θ-k = Θ * k , Ψ -k = Ψ * k .Depending on the needs the notation k can be replaced by the equivalent (k x , k y ). The Rayleigh-Bénard equation then take the form: ∂ t Θ (kx,ky) + ik y ∇ x ψΘ (kx,ky) -∇ x ΘΨ (kx,ky) + κk 2 Θ (kx,ky)+ k Ψ (k x ,k y ) Θ (kx-k x ,ky-k y ) k x k y -k y k x = 0 (C.10a) ∂ t W (kx,ky) + ik y ∇ x ψW (kx,ky) -∇ x W Ψ (kx,ky) + ik y gΘ (kx,ky) + νk 2 W (kx,ky)+ k Ψ (k x ,k y ) W (kx-k x ,ky-k y ) k x k y -k y k x = 0 (C.10b) W (kx,ky) = -k 2 Ψ (kx,ky) (C.10c)

∂

  t Θ (kx,ky) -ik y ∇ x Θ Ψ (kx,ky) + κk 2 Θ (kx,ky)+ k T Θ,Ψ (k, k )Ψ (k x ,k y ) Θ (kx-k x ,ky-k y ) = 0 (C.11a) ∂ t Ψ (kx,ky) -i k y k 2 gΘ (kx,ky) + νk 2 Ψ (kx,ky) + k T Ψ,Ψ (k, k )Ψ (k x ,k y ) Ψ (kx-k x ,ky-k y ) = 0 (C.11b)

( 1 +

 1 α 2 + Pr) + α 2 γ(Pr + Ra) + 2α 2 Pr(Ra -1) = 0 (C.[START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code soledge2d-eirene[END_REF]

ω 3 =

 3 ω(2Aε + B) (C.21) Aω 2 = -3ω 2 ε + Bε + C (C.22) C. 111 For ω = 0, equation C.21 reads: ω 2 = B + 2Aε (C.23) and (C.22) becomes:AB + 2A 2 ε = -3Bε + Bε + C (C.24)Therefore:2ε(A 2 + B) = C -AB (C.25)For ε = 0 (at marginality), one has ω 2 = B and C = AB, providing:ω 2 = α 2 (Pr + Ra) (C.26)and:Ra * = Pr 3 + α + Pr Pr -(1 + α 2 ) = 1 + 1 + α 2 + (2 + α 2 )Pr + Pr 2 Pr -(1 + α 2 ) (C.27)Ra * > 0 only for Pr > α 2 + 1. In this latter case, S * > 1.When ε = 0, the sign of ε is given by the sign of C -AB. Since α 2 = 0, C > AB requires: Ra(Pr -(1 + α 2 )) > Pr(3 + α 2 + Pr) (C.28)

Contents 2.1 Introduction to turbulence concepts in neutral fluids . 22

	2.1.1 Energy cascade and Kolmogorov theory for isotropic
	homogeneous turbulence . . . . . . . . . . . . . . . . . . 23
	2.1.2 Statistical approach and averaging techniques . . . . . . 26
	2

.2 Numerical approach for turbulent flows . . . . . . . . . 27 2

  .2.1 The RANS averaging procedure . . . . . . . . . . . . . .

28 2.3 The eddy viscosity concept . . . . . . . . . . . . . . . . . 30 2.4 Turbulence models based on one-point closure . . . . . 31

  One-equation models . . . . . . . . . . . . . . . . . . . . 31 2.4.3 Two-equations models . . . . . . . . . . . . . . . . . . . 32

	2.4.1 Prandtl mixing length theory . . . . . . . . . . . . . . .	31
	2.4.2

58 4.3 Second step: one-equation model, tuning with exper- imental scaling law

  . . 54 4.2.1 Equation and dimensional analysis . . . . . . . . . . . . 55 4.2.2 Free parameters and model's closure . . . . . . . . . . . 56 4.2.3 Experimental setting . . . . . . . . . . . . . . . . . . . . 58 4.2.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.3.2 Free parameters and model's closure . . . . . . . . . . . 62 4.3.3 Experimental and numerical setting . . . . . . . . . . . 64 4.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.

. . . . . . . . . . . . . . . . . . . . . . 61 4.3.1

4 Third step: two-equation model, tuning with experi- mental scaling law

  Equations and dimensional analysis . . . . . . . . . . . 68 4.4.2 Free parameters and model's closure . . . . . . . . . . . 69 4.4.3 Experimental and numerical setting . . . . . . . . . . . 70 4.4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . 67 4.4.1 71 51 52

Table 4 . 1 :

 41 Main plasma parameters for TCV shot 51333.

	-3 ]

B tor [T ] I pol [kA] P Ohm [M W ] n e,sep [m

Table 4 . 2 :

 42 Boundary conditions for the three simulations n seo,OM P [m -3 ] P in [kW ] 0.7•10 19 100

Table 4 . 3

 43

: Value of the λ SOL at the divertor, remapped at the outer midplane, estimated numerically.

Table 4 . 4 :

 44 Main plasma parameters for WEST shot 52698.

Table 4 . 5 :

 45 Boundary conditions for density and input power.

Table 5 .

 5

			1: WEST data
	R[m] a[m] B T [T ] B P [T ] k
	2.5	0.5	3.7	0.2	1.5 5

  , equation (B.1) is multiplied by u :

					mn 0	∂(u i u j ) ∂t	+ mn 0 u	∂(u i u j ) ∂x	=
	-u j	∂P ,i ∂x	-u i	∂P ,j ∂x	+ u

j (R ,i + F ,i ) + u i (R ,j + F ,j ) (B.2) 102 B.

1. Derivation of a general "κ" equation for magnetized plasma

  j ∂ t u i + u j uTaking this expression into account one can then simplify B.4b.Z j + u j Z i (B.7b)The trace defines the energy K of the fluctuating field in the κ model.We start here from the momentum balance equation equation: momentum conservation taking into account the simplifications introduced in the buoyancy force.∂ t m i nu + ∇• pI + m i nu ⊗ u -ν∇ 2 m i nu -m i n 0 Θg x = 0 (C.2a)

	106	B.							C.1. The 2D Rayleigh-Bénard system	103
	C.1.2.1 Momentum balance
		∂ t u i u j +	∂ u u i u j ∂x
		+ u	∂ u i u j ∂x	+ u u j	∂u i ∂x	+ u u i	∂u j ∂x
		+ u j = u i K = 1 ∂ Π i ∂x + u i ∂ Π j ∂x 2 mn 0 u i u i	(B.8a)
		∂ t K + u	∂K ∂x	+ 1 2 mn 0	∂ u u i u j ∂x
			+ mn 0 u u i	∂u i ∂x	+ u i	∂ P i ∂x	= mn 0 u i Z i	(B.8b)
				∂u i ∂x	+ u j	u u i ∂x	+ u j	∂Π i ∂x	= u j Z i	(B.6a)
		u i ∂ t u j + u i u	∂u j ∂x		+ u i	∂ u u j ∂x	+ u i	∂Π j ∂x	= u i Z j	(B.6b)
	So that one obtains:
		∂ t u i u j + u	∂ u i u j ∂x		+ u j	∂ u u i ∂x	+ u i	∂ u u j ∂x
			+ u j	∂Π i ∂x	+ u i	∂Π j ∂x	= u j Z i + u i Z j	(B.7a)

Acknowledgements

Although this analogy may seem overused, I must admit that in my experience

A

A.1 Vorticity equation for neutral fluids

Let us consider the momentum balance for neutral fluids:

where d Θ (m i n) is the dilatation coefficient depending on the temperature variation Θ, G is the buoyancy force aligned along the î direction. k is the axis of symmetry of the system. For the sake of simplicity let us the Boussinesqu's approximation, namely the density is assumed constant in time and space, the only effect of density variation being present in the buoyancy force driven by dilatation. The curl of equation (A.1) reads:

where ω ≡ ∇ × u and α Θ ≡ d Θ (log(n)). We use the following relation to treat the convection term:

where the second term vanishes when applying the curl operator and the third term cancels out due to the assumption of nearly incompressible flow. When projecting equation (A.2) along the k direction, one obtains:

where we define W ≡ ω • k. We recall another useful relation:

A.2. Vorticity equation for magnetically confined plasma

the first two terms vanishing by assumption and by construction respectively. Moreover, one has:

where u = u • k. Now one can rewrite equation (A.4):

where

) can be further manipulated if one explicit the velocity field with respect to its stream function ψ: u = k × ∇ψ. Then one has:

with • standing for Poisson's brackets. So finally equation (A.7) becomes:

where the third term in the left-hand side takes into account for the stratification of density due to temperature variation while the second term in the right-hand side may be further developed according to the system once the equation is integrated along the axis of symmetry.

A.2 Vorticity equation for magnetically confined plasma

We start from the conservation of total electric charge:

where b = B/B and j E×B , j p and j * are total E ×B, polarization and diamagnetic charge densities respectively. j E×B = en i u E×B -en e u E×B = 0.

Assuming that ∇ × B ∼ 0, reasonable in an electrostatic framework, it can be demonstrated for the divergence of the diamagnetic current [A/m 3 ]:

While for the divergence of the polarization current [A/m 3 ], one has:

For convenience, we define the vorticity the following quantity:

where the kinetic pressure is in [eV ] and the potential in [V ].

One can develop ∇ • j p and consider the first-order approximation of the ion drift velocity u ⊥,i ∼ u

⊥,i , namely u i ∼ u ,i + u E×B,i + u ∇B,i , and obtain:

The parallel current density ∇ • (j b) can be computed from the generalized Ohm's law:

where σ [S/m] is the electric conductivity (j = σE).

Let us stress some important assumptions about the magnetic field B that fit the 1D κ-ε model:

2. ∇ × B ∼ 0 acceptable for low β ≡ nT B 2 /2µ 0 namely in the electrostatic framework.

3. B is homogeneous and constant, so b × ∇B = 0, ∇ 1 B = 0 and ∇ • u E×B = 0 For these reasons, it results: u ∇B,i = u ∇B,e = 0 and equation (A.2) reads:

Rayleigh-Bénard equations

Let us consider equation A.9 in appendix A when (ω • ∇)u = 0. The Rayleigh-Bénard set of equations is completed by a transport equation of the temperature difference Θ of the same form as that for the vorticity W :

This expression characterises the evolution of the field Θ given a source S and the divergence of the flux Q defined by:

This heat flux is governed by a convective contribution with the incompressible velocity u depending on the stream function ψ, leading to the Poisson bracket in equation: temperature balance v5 and on a diffusive flux -κ∇Θ.

C.1.2 Averaging and fluctuations

In the chosen example of this Section, the z-direction is a direction of symmetry, i.e. there is no dependence on z and the x-direction is that of stratification governed by the buoyancy force gΘx. We shall further assume that the direction y is periodic which allows one introducing an average with respect to this direction and consequently to split the fields f into their average f and fluctuations f , namely the departure from the average, f = f + f . One can then apply this projection on the momentum balance equation equation: momentum conservation as well as to the vorticity equation: momentum conservation v5 and temperature equations equation: temperature balance v5. We first consider the momentum balance and in particular the non-linear tensor.