S Baschetti

H Bufferand

G Ciraolo

Ph Ghendrih

E Serre
email: eric.serre@univ-amu.fr

P Tamain

West Team

Self-consistent cross-field transport model for core and edge plasma transport

A two-equation model to self-consistently determine cross-field fluxes in the edge and scrape-off layer region of diverted plasma is used to complete 2D mean-field edge transport description of plasma wall interaction. Inspired by the Reynolds Average Navier-Stokes simulations for neutral fluids, this model is based on the local evolution of the turbulent kinetic energy κ and its dissipation rate ε. These two equations are algebraically derived for RANS modeling and are very slightly modified and adapted to describe self-consistent plasma turbulent transport. The general features of the model are discussed and bridged to the well-known predatorprey and quasilinear models commonly used to investigate plasma transport. Specific closures are proposed based on the interchange turbulence. Results of the 1D model are confronted to experimental evidence by analyzing the computed SOL width and comparing the results to the existing scaling law for L-mode plasmas. Introducing a dependence on the shear of large scale flows, typically the zonal flows, 1D simulations can exhibit an H-mode like transition when increasing the input power, generating an increased stored energy thanks to an interface barrier located at the separatrix. Further 2D plasma-wall interaction simulations for WEST are analyzed that show a good match with the experimental profiles, as well as a ballooned transport driving turbulent transport in the divertor SOL and nearly no transport in the private flux region. The SOL width of WEST is also recovered. These results show the remarkable capability of the κ-ε model to capture key aspects of the physics of turbulent transport throughout the plasma knowing that a unique scalar free parameter is available to tune cross field transport in the whole 2D cross section of the plasma.

Introduction

The control of plasma-wall interaction in next step devices aiming at burning plasma operation is presently understood to be a key research topic bridging the physics of advanced divertor scenarios and technological constraints governing the heat flux exhaust capability of the wall components [START_REF]ITER Physics Expert Group on Divertor Modelling and Database and ITER Physics Basis Editors. Chapter 4: Power and particle control[END_REF][START_REF] Loarte | the ITPA Scrape-off Layer, and Diver Group. Chapter 4: Power and particle control[END_REF]. The ITER divertor design is based on a long standing effort of transport simulations of the edge plasma in axisymmetric geometry [START_REF] Kukushkin | Effect of neutral transport on ITER divertor performance[END_REF][START_REF] Pitts | Status and physics basis of the ITER divertor[END_REF]. This simulation effort has played a key role in defining guidelines for the divertor design. Recent effort has been made towards design optimization [START_REF] Dekeyser | Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method[END_REF] and determination of reliable transport coefficients [START_REF] Baelmans | Efficient parameter estimation in 2D transport models based on an adjoint formalism[END_REF]. However, stepping towards high performance experiments near the operational limits, taking into account the aging of the components and experimental feedback, will require improved reliability of the numerical tools.

A suite of models is being developed, ranging from simplified models for optimization and uncertainty propagation to state-of-the-art first principle models of plasma turbulence transport in relevant plasma conditions. In that respect, fullf gyrokinetic simulations of the core and edge plasma are being used in the fusion community, but remain extremely costly from the computational point of view. In particular, the near-wall region requires addressing particle transport, hence electron and ion dynamics on the same footing, taking into account ionization particle sources and given magnetic as well as boundary condition geometries that are much more complex than those considered in the core [START_REF] Chang | Integrated particle simulation of neoclassical and turbulence physics in the tokamak pedestal/edge region using XGC[END_REF][START_REF] Churchill | Pedestal and edge electrostatic turbulence characteristics from an XGC1 gyrokinetic simulation[END_REF][START_REF] Caschera | Global confinement properties of Tokamak plasmas in global, flux-driven, gyrokinetic simulations[END_REF]. Performing such simulations with 3D fluid codes that handle self-consistently all scales of the flow, from the grid spacing to device size, are generally still restricted to rather simplified geometries and only take into account a fraction of the atomic physics at play in plasma wall interaction, see [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF][START_REF] Bufferand | Numerical modelling for divertor design of the WEST device with a focus on plasma-wall interactions[END_REF] and references therein.

Thus, in view of ITER operation or reactor design, the reduced models implemented in the transport codes are well adapted to provide relevant information on appropriate return times, similarly to the Reynolds Averaged Navier-Stokes (RANS) codes commonly used for engineering applications in the neutral fluid community [START_REF] Pope | Turbulent Flows[END_REF]. These transport codes are based on reduced fluid models, generally assuming axisymmetry of the plasma and accordingly using axisymmetric averaged equations. The reduction of the number of degrees of freedom allows one to take into account additional equations describing surface physics processes and atomic physics in realistic tokamak geometries. There exists in the fusion community a certain number of such 2D state-of-the-art transport codes. The reference effort started with B2 [START_REF] Schneider | B2-EIRENE simulation of ASDEX and ASDEX-Upgrade scrape-off layer plasmas[END_REF], then leading to SOLPS [START_REF] Schneider | B2-solps5.0: Sol transport code with drifts and currents[END_REF][START_REF] Schneider | Plasma edge physics with B2-Eirene[END_REF] and lately to SOLPS-ITER [START_REF] Wiesen | The new SOLPS-ITER code package[END_REF]. These three plasma transport modules are coupled to the kinetic code Eirene [START_REF] Reiter | The EIRENE and B2-EIRENE codes[END_REF] for the neutral particle transport. The other major codes are UEDGE [START_REF] Rognlien | A fully implicit, time dependent 2-D fluid code for modeling tokamak edge plasmas[END_REF] (coupled to the DEGAS code [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF] for neutral particles), SONIC [START_REF] Kawashima | Simulation of radiative divertor plasmas by Ar seeding with the full W-wall in JT-60SA[END_REF], EDGE2D [START_REF] Simonini | Models and numerics in the multi-fluid 2-D edge plasma code EDGE2D/U. Contrib[END_REF] (first coupled to the NIMBUS [START_REF] Cupini | Nimbus-Monte Carlo simulation of neutral particle transport in fusion devices[END_REF] for neutral particles and later to Eirene), and SolEdge2D-Eirene [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][START_REF] Bufferand | Near wall plasma simulation using penalization technique with the transport code SolEdge2D-Eirene[END_REF]. The latter, operated with configurations of most existing tokamaks in Europe and used in the present work, simulates the plasma edge and scrape-off layer (SOL) in a toroidally axisymmetric spatial domain including realistic wall geometry and detailed plasma-wall interaction. It is coupled to the Monte Carlo code Eirene [START_REF] Reiter | The EIRENE and B2-EIRENE codes[END_REF], which generates the particle source from the various neutral ionization processes, for both atoms and molecules.

A major mathematical change in the structure of the transport equations is the step from microscopic convective turbulent transport to a larger scale diffusive transport. Models that are presently used for plasma-wall interaction assume that the transverse fluxes are driven by local gradients, and consequently defined by diffusion or conductivity coefficients (equivalent to the turbulent eddy viscosity introduced in neutral fluids). The apparent universality of diffusive transport models is connected to the fact that it is in practice the simplest local transport model that can be implemented in equations. The effective particle diffusion D n stemming from this process, and similarly the effective viscosity ν n for the parallel momentum and effective heat conductivity χ e and χ i , respectively for the electrons and the ions, must be determined at each point of the mesh. As an example the perpendicular turbulent flux of particles ñṽ ⊥ , namely the convective transport of the density fluctuations n by the cross-field drifts ṽ⊥ is closed as:

ñṽ ⊥ cg ∼ -D n ∇ ⊥ n cg (1)
where cg stands for the coarse graining procedure. This step is most important since it determines in fact a projection operator and consequently the physics that one retains and that which is orthogonal [START_REF] Zwanzig | Memory effects in irreversible thermodynamics[END_REF][START_REF] Mori | Transport, Collective Motion, and Brownian Motion*)[END_REF]. Ensemble averaging is often presented as the backbone for coarse graining. However, this elegant argument is not practical and not really informative regarding the rules that govern the coarse graining process. In practice an averaging operator over the high-frequencies is used; D n is then an effective diffusion coefficient and ∇ ⊥ n cg the local transverse gradient of the coarse grained density field.

For ITER, the global behavior of the plasma must be analyzed with time scales of a fraction of the energy confinement time in the core and a fraction of the thermal time scale of the plasma facing components, both in the range of seconds. The chosen coarse graining time scale τ cg is therefore in the range of a second while the turbulence time scales τ turb are usually assumed to range between microseconds and milliseconds. It is to be noted however that in flux driven turbulence simulations, no spectral gap is observed between this turbulence time range and the macroscopic time τ cg [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF]. This can be an issue in the coarse graining procedure [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF]. Regarding space, two length scales are observed, the system size and the turbulence scale ṽ⊥ τ turb . For magnetic fusion plasmas, the length ṽ⊥ τ cg exceeds the system size. In the plasma core, with system size given by the minor radius a, the regime ṽ⊥ τ turb a can be assumed to hold. For plasma-wall interaction, characterized by the SOL width λ n , the regime ṽ⊥ τ turb ≈ λ is more likely. Defining and implementing a coarse graining procedure is then less straightforward [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF].

In present transport simulations, the particle, momentum and heat diffusion coefficients are not determined self-consistently but are tuned to match experimental radial profiles usually known at a single poloidal location, typically in the midplane. This is done either by providing transport coefficients as input to the code, having checked a good match between the empirical profiles and the reconstructed ones, either in a more sophisticated manner by adjusting automatically these values according to the midplane profiles as part of the initial step of the run. In SolEdge2D-Eirene, the automatic fitting procedure based on a proportional-Integral feedback loop is implemented [START_REF] Baschetti | Optimization of turbulence reduced model free parameters based on L-mode experiments and 2D transport simulations[END_REF]. In both cases, radial profiles of the effective transport coefficients are determined and it is assumed that these do not exhibit any other dependence that strongly limits the predictive capability of such codes. We remind that accurate values of the transport coefficients are required to determine the SOL width, which is the crucial parameter that governs plasma-wall interaction. Indeed, the competition between the SOL parallel conductive/convective transport and cross field transport combining large scale drifts, and turbulent and collisional transport [START_REF] Stangeby | The plasma boundary of magnetic fusion devices / Peter C. Stangeby[END_REF] determines the width of the heat and particle channels impinging onto the target plates, and consequently the operation constraints governed by the exhaust limitations. It is also important to underline that in most transport codes used to investigate plasma-wall interaction, the local gradients, for example the gradient of the thermal energy ∇T are estimated in eV /m, irrespective of the magnetic geometry, while the transport coefficients are assumed to be constant on a given magnetic surface. Since the minimum distance between flux surfaces is found in the midplane on the low-field side, where the transport coefficients are determined, away from this location, the gradients then appear to be smaller due the magnetic field flux expansion. Consequently they yield a reduced flux. A ballooned transport is thus induced but not as a consequence of the known properties of micro-turbulence but because the diffusion and conductivity coefficients are assumed to be homogeneous on each flux surface. In such a framework, the induced ballooning is not governed by the physics of turbulent transport as reported in all tokamak configurations [START_REF] Labombard | Transport-driven scrape-off-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF][START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF][START_REF] Dif-Pradalier | The Mistral base case to validate kinetic and fluid turbulence transport codes of the edge and SOL plasmas[END_REF], but stems directly from the geometry of the magnetic surfaces.

Thus, these codes require stepping towards predictive transport modeling, freezing on a physics basis all the free parameters that account for cross-field plasma transport induced by the drifts [START_REF] Baelmans | Efficient parameter estimation in 2D transport models based on an adjoint formalism[END_REF]. Following our recent work in Ref. [START_REF] Bufferand | Interchange turbulence model for the edge plasma in SOLEDGE2D-Eirene[END_REF], we propose in this paper a model for the self-consistent estimation of cross-field fluxes in the edge and scrape-off layer regions of diverted plasma that has been implemented in SolEdge2D-Eirene. The present effort is inspired by the work done from the 70's in hydrodynamics [START_REF] Launder | The numerical computation of turbulent flows[END_REF], neutral fluid turbulence in contrast to plasma turbulence in the fluid framework, and adapted here to magnetically confined plasmas for fusion applications. In neutral fluids, the Boussinesq assumption that the turbulent stresses and the deformation speeds of the mean flow are proportional, defines the eddy viscosity ν t . This key quantity defining the transport properties of these fluids is then related to the turbulence kinetic energy κ ≡ 1 2 ṽ2 and turbulence dissipation rate ε. We proposes here to follow a similar procedure, thus introducing evolution equations for κ and ε and from these deriving the dependence on space and time of the plasma transport. We will model this transport by a diffusive process, as done for neutral fluids: convective transport will also be discussed since since turbulent transport in the edge appears to be governed by bursts of ballistic events [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF][START_REF] Sarazin | Intermittent particle transport in two-dimensional edge turbulence[END_REF][START_REF] Fedorczak | Flow generation and intermittent transport in the scrape-off-layer of the Tore Supra tokamak[END_REF]. An approach similar to that addressed here was recently proposed by [START_REF] Coosemans | A new mean-field plasma edge transport model based on turbulent kinetic energy and enstrophy[END_REF][START_REF] Coosemans | Turbulent kinetic energy in 2D isothermal interchange-dominated scrape-off layer E × B drift turbulence: Governing equation and relation to particle transport[END_REF], but in a simpler configuration for isothermal plasma and 2D closed field lines domain. In this latter work, the perpendicular transport coefficients are determined from the turbulent kinetic energy κ and the enstrophy which is an invariant in 2D turbulence in neutral incompressible fluids. Their equations are analytically derived from the interchange turbulence model implemented in the TOKAM2D code [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF][START_REF] Ghendrih | Generation and dynamics of SOL corrugated profiles[END_REF]. This work suggests that the connection between the κ-ε model and fundamental properties of turbulence, such as the energy circulation between fluctuations and mean field [START_REF] Scott | The character of transport caused by E×B drift turbulence[END_REF], could be investigated based on these results from plasma physics. It is not the scope of this paper to step in that direction, but rather to present the capability of such a model to describe the edge and SOL transport with self-organization of turbulent transport. The κ-ε model can be seen as an intermediate model reduction between that inspired by the quasilinear results [START_REF] Kinsey | The first transport code simulations using the trapped gyro-Landau-fluid model[END_REF][START_REF] Bourdelle | Core turbulent transport in tokamak plasmas: bridging theory and experiment with QuaLiKiz[END_REF] and direct coupling of turbulence to mean-field transport [START_REF] Nishimura | Coupling of perpendicular transport in turbulence and divertor codes[END_REF][START_REF] Zhang | Self-consistent simulation of transport and turbulence in tokamak edge plasma by coupling SOLPS-ITER and BOUT++[END_REF].

The paper is organized as follows. We first discuss the general features of the κ-ε framework and bridge this approach to predator-prey and quasilinear models that are used to investigate plasma transport. We then introduce the SolEdge2D-Eirene suite of codes in Section 2. The general advection diffusion transport equation for κ and ε, implemented in SolEdge2D-Eirene, are detailed and the rules used to close the system are presented. In the fourth section, a one-dimensional reduction is proposed to carry out fast scans of plasma parameters and determine the main trends of the solutions. In the fifth section, the numerical results are confronted to experimental data in two Lmode plasma discharges of the WEST tokamak [START_REF] Bourdelle | WEST physics basis[END_REF]. Finally, a discussion and conclusion Section closes the paper.

Local κ-ε evolution, comparison to existing turbulent transport models

The κ-ε transport model

It this work, we consider the standard κ-ε model used in neutral fluids [START_REF] Launder | The numerical computation of turbulent flows[END_REF] where κ is the kinetic energy per unit mass of the fluctuating transverse velocity, typically ≈ v 2 E and ε is a damping process acting on κ. One thus retains the same evolution equations.

∂ t κ + ∇•K = S κ -ε (2a) ∂ t ε + ∇•E = S ε (2b)
Here ∇ • K is the divergence of the flux of κ and S κ -ε is the local evolution term of κ. The evolution equation for ε (2b) is built to be similar to that of κ. The term ∇•E accounts for the divergence of the flux of ε and S ε is the local evolution term of ε. The coupling between these equations that is made explicit by the term -ε on the right hand side of (2a) identifies the field ε as a rate of turbulent energy dissipation and determines the relative normalization of ε with respect to κ. The source terms S κ and S ε are expanded up to order 2:

S κ = γ κ κ-β κ κ 2 , S ε = γ ε ε-β ε ε 2 .
The linear part depends on the growth rates γ κ and γ ε . The non-linear parts are chosen to be restoring terms with β κ and β ε positive. The evolution equations (3) for two fields κ and ε in R + exhibit therefore a generic local form. Assuming diffusive transport in the transverse direction and the expansion of S κ , S ε described above, the system [START_REF] Loarte | the ITPA Scrape-off Layer, and Diver Group. Chapter 4: Power and particle control[END_REF] takes the following form:

∂ t κ -∇ ⊥ • D κ ∇ ⊥ κ + ∇ •K = γ κ κ -β κ κ 2 -ε (3a) ∂ t ε -∇ ⊥ • D ε ∇ ⊥ ε + ∇ •E = γ ε ε -β ε ε 2 (3b)
It can also be shown that equation (3b) can be rewritten as a Ginzburg-Landau amplitude equation [START_REF] Wim Van Saarloos | Fronts, pulses, sources and sinks in generalized complex Ginzburg-Landau equations[END_REF][START_REF] Cross | Pattern formation outside of equilibrium[END_REF] for the amplitude A such that ε = |A| 2 . This further underlines the generic structure of the κ-ε system (3). In the two following Sections, we discuss the connection of the κ-ε approach to models presently used to investigate plasma turbulence: the predator-prey and quasilinear transport models.

Local κ-ε dynamics and comparison to predator-prey models

Let us consider here the local evolution equations for the fields κ and ε, hence system (3) without cross-field and parallel transport.

∂ t κ = γ κ κ -β κ κ 2 -ε (4a) ∂ t ε = γ ε ε -β ε ε 2 (4b)
One can consider the system (4) as a predator-prey system and more generally any kind of reservoir system with internal dynamics and coupling [START_REF] Grisolia | Plasma wall particle balance in Tore Supra[END_REF][START_REF] Loarer | Particle balance modelling in ergodic divertor experiments on Tore Supra[END_REF]. The predatorprey dynamics [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF] have been introduced in plasma-physics to investigate bifurcationlike phenomena of transport properties [START_REF] Diamond | Self-regulating shear flow turbulence: A paradigm for the L to H transition[END_REF][START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF]. However, rather than adjusting the predator-prey model to match a particular behavior, we only consider in this paper the local evolution terms defined by the κ-ε model (3) as chosen for neutral fluids [START_REF] Launder | The numerical computation of turbulent flows[END_REF]. The field ε is an effective predator of the prey κ. Indeed the two species exhibit growth rates that are similar and do not depend on either species κ or ε. The effective predator ε does not feed directly on κ. Rather the existence of ε inhibits the prey development as if modifying the ecosystem in an unfavorable way for κ. The fixed points of the predator equation are ε = 0 and ε = γ ε /β ε . For β ε > 0, one thus finds a phase transition behavior such that for γ ε < 0, ε = 0 is the only possible fixed point while for γ ε > 0 two fixed point exist. One obtains similarly that for γ κ < 0 the only fixed point is κ = 0. This transition for the coupled equations is reminiscent of the bifurcation in the Ginzburg-Landau amplitude equation such that for negative growth rate the stable fixed point is zero that becomes unstable when the growth rate is positive. The comparison to the predator-prey models provides a first understanding of the field ε as the predator of the turbulent energy. One can also remark that ε is akin to the time derivative of κ. The set of equations can then be understood as the first two of a hierarchy of equations coupling increasing order time derivatives. The closure at ε is performed at the same order as the Newtonian closure in mechanics. The fact that κ ε belong to such a hierarchy of time derivatives justifies the choice to retain similar equations for both κ and ε. One can conjecture that this structure is the key point that drives the efficiency of the κ-ε model. Let us now address the properties of the non-linear saturation terms. In neutral fluids [START_REF] Launder | The numerical computation of turbulent flows[END_REF] the definitions of β κ and β ε are guided by dimensionality arguments. Hence β κ = 1/D ω where D ω is a diffusion coefficient. Regarding β ε we note that ε is dimensionaly equivalent to the rate of energy transfer between scales that governs the cascade paradigm of steady-state turbulent spectra. Compared to the choice made for neutral fluids, we then make one step towards interpreting the κ-ε model. The field κ stands for the turbulent energy of the modes that govern turbulent transport. In the framework of turbulent cascades, ε is a dissipation mechanism that can be understood as the rate of energy transfer by direct and inverse energy cascade [START_REF] Pope | Turbulent Flows[END_REF] to regions of the spectrum that weakly contribute to transport. In the generic case of a direct cascade ε couples this energy to the Kolmogorov region of the spectrum where energy is dissipated. This interpretation in terms of a direct energy cascade is not mandatory for the present work, we believe that it provides an understanding of the physics taken into account by ε. It is to be noted that it is often argued that 2D turbulence leads to an inverse energy cascades. This argument is too loose and does not hold when considering turbulent transport, namely turbulent convective transport of density, pressure, etc. coupling the turbulent velocity to fluctuations of the transported field at all scales. In this case the existence of a second quadratic invariant akin to the enstrophy is possible [START_REF] Suzana | Resistive drift-wave turbulence[END_REF] but not guaranteed. Lastly, the parallel transport in the third dimension cannot be ignored, especially the parallel losses that characterize the SOL physics, and can break the symmetry conditions required in neutral fluid 2D turbulence for a second quadratic invariant to exist. For steady state conditions at a given scale , a dimensional analysis readily yields the relationship between the rate of turbulent energy transfer ε and the turbulent energy κ.

κ ∝ 2/3 ε 2/3 (5)
This universal feature drives both direct and inverse energy cascades. Stepping to determining the power law dependence in the region of the spectrum governed by the energy cascade (direct as well as inverse), one recovers the well known power law 5/3 ε -1/4 of the Kolmogorov theory [START_REF] Kolmogorov | The Local Structure of Turbulence in Incompressible Viscous Fluid for Very Large Reynolds' Numbers[END_REF].

To enforce the steady-state relation of turbulent spectra [START_REF] Dekeyser | Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method[END_REF] in the κ-ε model, we set

β ε ∝ κ -3/2
. This is a slight change with respect to the standard neutral fluid κ-ε model where β ε ∝ κ -1 [START_REF] Launder | The numerical computation of turbulent flows[END_REF]. For the sake of generality, we shall also consider β ε ∝ κ -η , which depending on η encompasses the choice made for neutral fluids η = 1 and that made in [START_REF] Dekeyser | Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method[END_REF], η = 3/2. It will be shown that all values of η ≥ 1 exhibit comparable properties, the choice of η ≥ 1 does not introduce critical changes to the model. The fixed points of the local evolution for ε is then zero and ε = γ ε /β ε = γκ * (κ/κ * ) 3/2 , γ > 0 is a proportionality factor, dimensionaly the inverse of a time, and κ * ∈ R + is a convenient normalization of κ to be chosen hereafter. The control parameter γ governs the amplitude of the dissipation via the energy cascade for a given turbulent energy.

The effect of varying this control parameter is addressed in the following. The steady state non-trivial fixed point thus leads to κ ∝ ε 2/3 as in [START_REF] Dekeyser | Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method[END_REF]. Given ε = γκ * (κ/κ * ) 3/2 , the fixed point for κ is then determined by the equation:

γ κ κ - 1 D ω κ 2 - γ κ 1/2 * κ 3/2 = 0 (6a)
Solutions other than κ = 0 are thus determined by:

κ κ * + κ κ * 1/2 - γ κ γ γD ω κ * = 0 (6b)
More generally, we set ε = γκ * (κ/κ * ) η , with η = 1 for standard neutral fluids, and η = 3/2 for (5). We define κ * = |γ κ |D ω > 0 so that the fixed point is either κ/κ * = 0 or the positive root of:

κ |κ * | + γ |γ κ | κ κ * (η-1) - γ κ |γ κ | = 0 (6c)
For the two reference values of η, η = 1 and η = 3/2, analytical solutions other than κ = 0 are obtained:

η =1 ; κ = γ κ D ω (7a) η = 3 2 ; κ 1/2 = - γ 2|γ κ | + γ 2|γ κ | 2 + γ κ |γ κ | 1/2 |γ κ |D ω 1/2 (7b)
For fixed γ > 0, and varying the control parameter γ κ , one finds, as with the standard Ginzburg-Landau wave amplitude equation, that for γ κ < 0 the only positive fixed point is κ = 0. Conversely, for γ κ ≥ 0 the stable fixed point is either κ = 0 or the solution given by (6c). Let us now assume that γ κ > 0 is fixed and vary γ that governs the coupling between the two fields κ and ε via the damping process at play in the spectrum, Figure 1. For γ = 0, ε = 0 and the two fields are decoupled, κ/κ * = 1. For the relevant values η ≥ 1, as γ is increased, κ decreases monotonically from its maximum value κ * towards zero, Figure 1 left hand side, while ε first increases with γ before decaying towards zero when γ → +∞, the maximum being achieved for γ/γ κ = 0.5 2-η and ε/ε * = 0.25 with ε * = γ κ κ * , Figure 1 right hand side. Compared to the case without coupling, readily obtained by setting γ = 0 so that ε = 0, adding the damping by ε in the equation for the turbulent energy, introduces a regulation with a possible time delay because this damping is governed by an independent equation with different time scales. However, the latter effect will not occur if the time scales that govern the evolution of ε are small, hence γ → +∞ so that ε exhibits an adiabatic response. Regarding the dependence on η, there is a change in behavior between η < 1 open circles, and η ≥ 1 closed symbols, as shown on Figure 1. In particular κ → 0 when γ → 0 for η < 1, while κ/κ * → 1 for η ≥ 1. Furthermore, instead of obtaining an asymptotic convergence towards zero when the control parameter γ is increased as with η > 1, for η = 1 both κ and ε switch to zero above the critical value γ = γ κ with a discontinuity of their derivative. The value η = 1 used in the standard κ-ε model is a marginal value with respect to the behavior for η > 1. It also appears that the cases η ≥ 1 are qualitatively comparable so that the particular choice of η > 1 does not determine specific properties of the fixed points of the local dynamics. It is to be underlined that the present version of the predator-prey model implemented in this κ-ε model is not chosen to exhibit particular physics in terms of bifurcations, limit cycles etc. Although it is convenient to see in the loss term ε that governed by the losses at the Kolmogorov dissipation scale via the energy cascade this is not the unique loss channel. Indeed, in the fusion plasma literature, an important loss path controlling turbulent transport is governed by a coupling to large scale flows [START_REF] Busse | Generation of mean flows by thermal convection[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF][START_REF] Itoh | Physics of zonal flows[END_REF]. This result fits with the picture of an inverse energy cascade although disparate wave coupling is often considered to be in fact the main drive. The latter can be understood as an alternative loss channel described by ε. In the following such a feature is taken into account by a modification of the ε local evolution. Indeed, it is understood that either the nonlinear coupling of the source term to the free energy, hence the dependence of γ κ and γ ε on the various gradients [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF], or the shearing effect via the self generated zonal flow [START_REF] Busse | Generation of mean flows by thermal convection[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF][START_REF] Itoh | Physics of zonal flows[END_REF] are the main players in the turbulence energy evolution, something of the form of ε but governed by different dynamics [START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF]. In this context, it is also important to mention the approaches based on prey-predator models [START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF] and used to investigate the H-mode transport barrier, its onset as well as its dynamics in the vicinity of the threshold. The aspect that is considered there, as well as in [START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF], is that the control mechanism of the kinetic energy κ of the quasi-2D plasma turbulence is governed by large scale and meso-scale processes such as the E × B shearing of turbulence driven by zonal flows and mean flows [START_REF] Busse | Generation of mean flows by thermal convection[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF][START_REF] Itoh | Physics of zonal flows[END_REF]. In this framework, ε, the predator for turbulence [START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF], stands for the stabilizing shearing effect of these large scale flows. The more complex local dynamics for κ-ε govern the features that are reminiscent of the L-H transition, as reported in [START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF][START_REF] Norscini | Self-organized turbulent transport in fusion plasmas[END_REF]. The comparison of the κ-ε model to the predatorprey models used to describe plasma self-organization opens the way to further adapt the standard κ-ε model by modifying the local evolution. We leave this investigation to future work. However, we show in this paper that including a shearing effect in the standard framework, without complex and highly nonlinear dynamics, already provides a scheme to self-consistently generate an interface barrier, Section 3.4.

κ-ε transport model and quasilinear framework

In the framework of reduced models, in particular for real time control, a renewed interest has been given to the quasilinear theory [START_REF] Vedenov | Quasi-linear theory of a plasma[END_REF][START_REF] Drummond | Non-linear stability of plasma oscillations[END_REF][START_REF] Frieman | Kinetic theory of a weakly unstable plasma[END_REF] by using a large data base of first principle simulations and experiments to constrain the quasilinear model [START_REF] Bourdelle | A new gyrokinetic quasilinear transport model applied to particle transport in tokamak plasmas[END_REF][START_REF] Bourdelle | Core turbulent transport in tokamak plasmas: bridging theory and experiment with QuaLiKiz[END_REF][START_REF] Citrin | Real-time capable first principle based modelling of tokamak turbulent transport[END_REF].

We show here that the κ -ε model can be seen as an extension of the quasilinear theory, determining consistently, at each position and each time, the turbulent kinetic energy, typically κ, and the characteristic time governing the width of the resonance. Let us now consider a simplified quasilinear model for plasma transport governed by the following generic equation:

∂ t f + Lf = S (8)
We assume that one can define an averaging procedure such that f 0 is the projection of f with this averaging procedure, and such that the projection of f 1 = f -f 0 is zero, hence < f >= f 0 and < f 1 >= 0. The averaging projection can also be applied on the evolution operator L = L 0 + L 1 with < L >= L 0 and < L 1 >= 0. For simplicity we assume < S >= S. Equation (8) can then be split into a set of two equations, for the mean and the fluctuations:

∂ t f 0 + L 0 f 0 + L 1 f 1 = S (9a) ∂ t f 1 + L 0 f 1 + L 1 f 0 + L 1 f 1 -L 1 f 1 = 0 (9b)
The quasilinear theory simplifies (9) by dropping the non-linear fluctuation terms in (9b) hence removing L 1 f 1 -L 1 f 1 while retaining L 1 f 1 in (9a). One could justify this step by considering some ordering such that the mean is of order 0 while fluctuations are of order 1, however, this does not imply that L 1 applied to f 1 is of order 2. This is clear with the example L 1 = v 1 ∇ ⊥ , v 1 and f 1 being of order 1, but one cannot state that ∇ ⊥ f 1 is of order 1 and not order 0. Furthermore, should L 1 f 1 be of order 2, it can be questionable to neglect terms of order 2 in an equation where the other terms are of order 1 as in (9b) while retaining an order 2 term in (9a) where the other terms are of order 0. It is more interesting here to consider symmetries. Indeed the term L 1 f 1 -L 1 f 1 can be seen as a correction to L 0 f 1 , and considering that f 1 does not belong to the kernel of L 0 , this correction does not change the structure of (9b). The latter can then be rewritten as:

∂ t f 1 + L * 0 f 1 + L 1 f 0 = 0 (10a) L * 0 f 1 = L 0 f 1 + L 1 f 1 -L 1 f 1 (10b) f 1 = - 1 ∂ t + L * 0 L 1 f 0 (10c)
The formal solution (10c) is obtained by assuming that the operator ∂ t + L * 0 can be inverted implying in particular that f 1 does not belong to the kernel of ∂ t + L * 0 . One can then consider the transport equation for the mean f 0 :

∂ t f 0 + L 0 f 0 -L 1 1 ∂ t + L * 0 L 1 f 0 = S (11)
Let us now discuss this result in the framework of transport in a magnetized plasma, typically setting L 0 = v ∇ and L 1 = v 1 ∇ ⊥ where the subscript and ⊥ refer to the direction of the magnetic field. We thus focus in this discussion on the position dependence of f 0 with respect to the magnetic surfaces; ∇ depending on the variation within a magnetic surface, and, ∇ ⊥ depending on the variation between magnetic surfaces, hence ∇ ⊥ ≈ ∇ψ∂ ψ where ψ is a magnetic surface label. Setting L 1 = 0 and S = 0 in [START_REF] Bufferand | Numerical modelling for divertor design of the WEST device with a focus on plasma-wall interactions[END_REF], and given the symmetry ∇ ψ = 0, f 0 is a steady state solution when f 0 belongs to the kernel of L 0 , hence f 0 only depends on ψ. In this case perfect confinement is obtained, the source term S is not needed. Breaking the symmetry with fluctuations, L 1 = 0 then governs the need for a source term to achieve a steady state.

With the proposed operators, we have implicitly assumed the simplification ∇ ⊥ v 1 = 0, furthermore, approximating L * 0 by L 0 as done in the quasilinear framework, one obtains:

∂ t f 0 + L 0 f 0 -∇ ⊥ v 2 1 ∂ t + v ∇ ∇ ⊥ f 0 = S (12a)
D ql = 1 2 v 2 1 τ ql = κ τ ql (12b) τ ql = 2 ∂ t + v ∇ (12c)
The quasilinear analysis thus yields three important features, (i) one obtains a diffusive like transport in the cross-field direction, (ii) the diffusion coefficient is proportional to κ, and (iii) a characteristic time is required to completely determine the diffusion coefficient D QL (12b). It is interesting to note that the result can be extended to the general case hence without approximating L * 0 by L 0 . This can modify some aspects of the commutation between v 1 and the operator τ ql but does not change the main features of the result, in particular the diffusive structure. Two key assumptions constrain the validity of this result, first the possibility of defining the averaging procedure in line with the symmetries that govern the operator L and second that the inversion of either ∂ t +L 0 or ∂ t + L * 0 is possible, and, when possible, does not govern a change in the structure of the final result such as the explicit dependence on ∇ ⊥ f 0 . Closing the quasilinear approach then requires to determine both the proper κ and τ ql . In the most recent effort, this is achieved using a large data base of local gyrokinetic simulations completed by setting a proportionality factor using experimental evidence [START_REF] Bourdelle | Core turbulent transport in tokamak plasmas: bridging theory and experiment with QuaLiKiz[END_REF][START_REF] Citrin | Real-time capable first principle based modelling of tokamak turbulent transport[END_REF]. The dimensional argument used for the κ-ε model can be considered to determine the time τ ql since the ratio κ/ε has the dimension of a time, one can then set τ ql ∝ κ/ε and one obtains therefore:

D ql = C ql κ 2 ε ( 13
)
where C ql is a proportionality factor. It is also interesting to note that the time κ/ε also appears in the evolution of κ (3a).

1 κ ∂ t κ - 1 κ ∇ ⊥ • D κ ∇ ⊥ κ + 1 κ ∇ •K = γ κ -β κ κ - ε κ (14)
Indeed, the local evolution, right hand side of (14) depends on three characteristic time scales, the linear drive 1/γ κ , the non-linear restoring time D ω /κ and κ/ε. The relations we have obtained at that stage appear to be comparatively complicated because we have chosen the proportionality factor γ as free parameter. When considering as free parameter the characteristic time τ ql = κ/ε, one obtains:

γ κ τ ql = γ κ κ ε = γ κ γ κ κ * -(η-1) (15)
Most of the complexity of this system is thus governed by the relationship between the time τ ql and κ. Conversely, the steady state solutions exhibit simple dependencies in terms of γ κ τ ql , Figure 2:

κ κ * = 1 - 1 γ κ τ ql (16a) ε ε * = 1 γ κ τ ql 1 - 1 γ κ τ ql (16b)
One finds here that 1/(γ κ τ ql ) ≤ 1, hence τ ql ≥ 1/γ κ so that τ ql is the longest time scale in the system. In this framework, setting D = κ 2 /ε, one finds a very simple expression for the diffusion coefficient:

D D * = κ 2 ε * εκ 2 * = γ κ τ ql -1 (16c)
The form [START_REF] Schneider | B2-EIRENE simulation of ASDEX and ASDEX-Upgrade scrape-off layer plasmas[END_REF] is the reference one for the diffusion transport coefficients in the κε framework. Therefore, one can consider that the κ-ε model is an extension of the quasilinear framework such that the fields κ and ε are determined within the local transport model, and evolved self-consistently rather that being imported using other tools, either fitting of experimental evidence or determined using a data base of local gyrokinetic simulations. Regarding plasma-wall interaction, which exhibits at least a 2D dependence in space, 2D measurements of the gradients are out of reach, the κ-ε approach then appears as an extension of the quasilinear theory that is better suited to address the complexity and variability in geometry and parameter space of the problem. Furthermore, rather that fitting the characteristic time, the κ-ε model allows one determining self-consistently this time and in particular the dynamics with the possible occurrence of time delays in the response of ε.

2. Closure of the κ-ε model for plasma transport

The κ-ε model with diffusive cross-field transport

In Section 1.2 we have discussed the form taken by the non-linear restoring terms of the local evolution. For convenience we now write the restoring term as β ε = V /κ 3/2 where V is free parameter, dimensionaly equivalent to a velocity. Furthermore, as in (3) we neglect the large scale flows that arise from the non-homogeneity of the magnetic field as well as from large scale electric fields. This is in line with standard modeling of plasmawall interaction but will fall short whenever a finer description of cross-field transport is required. In particular when contemplating the κ-ε model as a tool to investigate the properties of turbulent transport described by first principle fluid or gyrokinetic simulations. To reduce the number of free parameters, we impose in the present work

D κ = D ε = D where D = κ 2 /ε.
Further work and comparison to transport data bases could lead to introduce further control parameters to ensure a better match with experimental and simulation evidence. Taking into account these different changes, one then obtains:

∂ t κ -∇ ⊥ • D∇ ⊥ κ + ∇ •K = γ κ κ - κ 2 D ω -ε (17a) ∂ t ε -∇ ⊥ • D∇ ⊥ ε + ∇ •E = γ ε ε -V ε 2 κ 3/2 (17b)
The transverse transport that governs the plasma fields is also assumed to be diffusive, and the coefficients proportional to D = κ 2 /ε. Various boundary conditions have been implemented for κ and ε. For the simulations presented in this paper the radial boundary conditions are set by vanishing second derivatives. A similar boundary condition has been used in the parallel direction for the 2D SolEdge simulations. This choice does not constrain the outflux and boundary value. Several 1D simulations have been run with vanishing boundary values for both κ and ε. The change only modifies a thin layer in the outer SOL, in agreement with the importance of the local dynamics for κ and ε compared to transport.

A first closure issue consists in determining these various diffusion coefficients. we first assume that the turbulent (subscript t) fluid viscosity ν t is proportional to D with proportionality factor C ν . The relevant turbulent particle diffusion D n is then defined by the turbulent Schmidt number Sc t = ν t /D n , while the electron and ion turbulent heat diffusivity are determined by a single turbulent Prandtl number P r t , hence:

ν t = C ν κ 2 ε ; D n = ν t Sc t ≈ ν t ; χ e = χ i = ν t P r t ≈ 2ν t (18)
The constant parameter C ν is to be determined with an appropriate constraint, typically a fitting procedure of either empirical or simulation data. In the present work, we further assume, as rather standard in neutral fluids Ref. [START_REF] Gualtieri | On the values for the turbulent schmidt number in environmental flows[END_REF] and references therein, that the turbulent Schmidt number Sc t = ν t /D n is of order unity. We shall also consider an empirical observation from midplane profiles which have been matched by steady state transport simulations with the constraint χ e,i = 2 × D n [START_REF] Chankin | SOLPS modelling of ASDEX upgrade H-mode plasma[END_REF]. Accordingly, we shall set the turbulent Prandtl number P r t = ν t /χ to be equal to 1/2. Again future work with fitting procedures using empirical or simulation data should be used to provide a firm basis to set these closure parameters, possibly different from those chosen in this paper.

Growth rates closure for the plasma κ-ε model

Since κ and ε are determined locally in space and at each time step by the model (3), the diffusion coefficient D determined from κ and ε is also defined locally and exhibits an evolution in time. Consistently, we consider therefore that the growth rates in (3) or ( 17) are local and time dependent. This is a major particularity of the κ-ε model in neutral fluids that we transpose to plasmas. Plasma turbulence in magnetically confined plasmas mostly arises from instabilities usually identified according to their linear drive. For electrostatic turbulence two main classes can be identified: (i) interchangelike instabilities governed by the magnetic field inhomogeneity that drives a current across the magnetic surfaces, (ii) the drift wave instability mostly driven by parallel currents. The interchange-like mechanism drives Ion Temperature Gradients [START_REF] Guzdar | Ion-temperature-gradient instability in toroidal plasmas[END_REF] and SOL-interchange [START_REF] Garbet | A model for the turbulence in the scrape-off layer of tokamaks[END_REF]. It is characterized by a poloidal asymmetry of the drive, leading to turbulence ballooning on the low field side of the torus. Conversely the drift wave instability [START_REF] Horton | Drift waves and transport[END_REF][START_REF] Scott | Drift wave versus interchange turbulence in tokamak geometry: Linear versus nonlinear mode structure[END_REF], including the sheath negative resistivity mechanism [START_REF] Berk | Temperature-gradient instability induced by conducting end walls[END_REF], does not exhibit a poloidal asymmetry in the drive. However, in the non-linear stage working back to a particular mechanism is questionable [START_REF] Scott | Drift wave versus interchange turbulence in tokamak geometry: Linear versus nonlinear mode structure[END_REF]. Furthermore, it has been found that linearly unstable mechanisms can "cooperate" [START_REF] Maeyama | Cross-scale interactions between electron and ion scale turbulence in a tokamak plasma[END_REF], or compete, one instability inhibiting the other [START_REF] Maeyama | Suppression of ion-scale microtearing modes by electron-scale turbulence via cross-scale nonlinear interactions in tokamak plasmas[END_REF]. Such predator-prey behavior could be transposed to κ-ε models. This is left for future work. For the present purpose of using the κ-ε model to investigate plasma-wall interaction, we have decided to single out the interchange-like instability driven by the plasma pressure gradient. As discussed, such a mechanism is understood to drive the poloidal asymmetry observed in the turbulent transport between the low field side and the high field side, the so-called ballooning of the turbulence, which is not consistently taken into account in 2D transport codes addressing plasma-wall interaction. The motivation is to recover this imbalance and does not mean that we regard the specific interchange instability as the only player. As in the standard κ-ε model, we shall furthermore consider the simplification γ κ = γ ε , so that γ κ = γ ε = γ I where γ I is the linear interchange instability growth rate [START_REF] Garbet | A model for the turbulence in the scrape-off layer of tokamaks[END_REF][START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF][START_REF] Ghendrih | Generation and dynamics of SOL corrugated profiles[END_REF]. The case with γ κ = γ ε is also of interest and addressed elsewhere [START_REF] Auroux | Limit cycles and spreading with k-ε plasma transport[END_REF].

γ I = c s R R 2 ∇p•∇B T pB T (19)
The variation in space of the axisymmetric toroidal magnetic field B T is governed by the tokamak major radius R. The time scale R/c s , where c s is the sound velocity, is chosen consistently with a drive depending on the pressure gradient ∇p. The pressure p and its gradient play a key role since they vary in space and time while RB T is constant when assuming axisymmetry. More complete models can be addressed, either by assuming competing growth rate of the different instabilities that add up to drive the same turbulent energy κ, or competing κ-ε systems, one for each instability, such that one would then add the various diffusion coefficients. However, the implicit rule in [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF], namely that γ I = 0 when (∇p/p) • (∇B T /B T ) ≤ 0 governs a ballooning of this local drive. It is then possible that a class of instabilities characterized by such a ballooning feature can be addressed with this particular growth rate. As will be shown in Section 2.6, a vanishing drive combined to the other transport properties will govern an effective damping of the turbulent energy. The overall effect is then reminiscent of the growth rates derived in Refs. [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF][START_REF] Ghendrih | Generation and dynamics of SOL corrugated profiles[END_REF][START_REF] Garbet | A model for the turbulence in the scrape-off layer of tokamaks[END_REF] that take into account the transport features.

Ballooned transport in the poloidal plane

In order to check that the growth rate ( 19) allows one to recover the salient features of the ballooned transport patterns in the poloidal plane, we analyze various simulations of a particular TCV experiments [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF][START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] performed with the 3D turbulence code TOKAM3X [START_REF] Tamain | The TOKAM3X code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF] and SolEdge2D-EIRENE. Depending on the simulations, ballooning effects are either included self-consistently, as with TOKAM3X and SolEdge2D with the κ-only transport model, or enforced in the transport model of SolEdge2D, or ignored as with standard SolEdge2D transport model. Let us first analyze the SolEdge2D-EIRENE simulations without self-consistent turbulent transport from Ref. [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF], Figure 3. On panel (a), we have reproduced the particle flux as obtained with the procedure of midplane profile tuning of the transport coefficients [START_REF] Baschetti | Optimization of turbulence reduced model free parameters based on L-mode experiments and 2D transport simulations[END_REF], and assuming that these coefficients can be extrapolated to the whole poloidal plane only taking into account the radial dependence. This is the standard procedure used in 2D transport modeling of plasma wall interaction. The patterns of strong outward particle flux are located on the one hand at the core to SOL interface, in both the high and low field regions, and on the other hand in the SOL region of the divertor for both the high and low field side legs. Finally, the region with moderate particle outflux exhibits poloidal symmetry around the core plasma and along the divertor legs. This very homogeneous transport pattern did not provide agreement with the experimentally estimated transport properties based on the measurement of the SOL widths. Furthermore, such symmetric transport pattern is known to disagree with experimental evidence that suggests ballooned transport [START_REF] Labombard | Transport-driven scrape-off-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF][START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF][START_REF] Dif-Pradalier | The Mistral base case to validate kinetic and fluid turbulence transport codes of the edge and SOL plasmas[END_REF]. Another simulation was performed with a different choice of transport coefficients, tuning the radial profile to fit the midplane profiles [START_REF] Baschetti | Optimization of turbulence reduced model free parameters based on L-mode experiments and 2D transport simulations[END_REF] and then extrapolating these to the poloidal plane by enforcing strongly ballooned transport, Figure 3, panel (b). A modest improvement when comparing to the TCV data is reported [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF].

Let us now analyze TOKAM3X and SolEdge2D-EIRENE with the κ-only turbulent transport model simulations, Figure 4. The first observation that can be made is the clear imbalance between the high field side region with no significant turbulent transport and the low field side region where cross-field transport is localized. The transport [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. Regions with largest particle diffusion coefficient are labeled "high", those with still significant particle diffusion are labeled with "medium".

coefficient is observed to increase from the inner core boundary towards the separatrix, following a universal trend of the fluctuation profile measurements [START_REF] Fonck | Long-wavelength density turbulence in the TFTR tokamak[END_REF][START_REF] Hennequin | Scaling laws of density fluctuations in tokamak plasmas[END_REF][START_REF] Hornung | Turbulence correlation properties measured with ultrafast sweeping reflectometry on Tore Supra[END_REF], recovered in some gyro-kinetic simulations [START_REF] Dif-Pradalier | Evidence for global edgecore interplay in fusion plasmas[END_REF]. Regarding the ballooning aspect, one finds that the poloidal extent is typically ±80 • , therefore broader than the ±45 • reported for the Tore Supra limiter experiment [START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF][START_REF] Dif-Pradalier | The Mistral base case to validate kinetic and fluid turbulence transport codes of the edge and SOL plasmas[END_REF][START_REF] Fedorczak | Experimental investigation on the poloidal extent of the turbulent radial flux in tokamak scrape-off layer[END_REF]. One can notice also that the ballooned aspect is superimposed to a near constant D n pattern extending poloidally in the low field side SOL and following approximately the separatrix. Comparing panels (a) and (b) indicates that qualitative agreement is achieved. The ballooned feature agrees with experimental observations Refs. [START_REF] Labombard | Transport-driven scrape-off-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF][START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF] and the localization along the low field side divertor leg as reported with visible imaging of turbulent filaments as in the MAST tokamak [START_REF] Harrison | Filamentary transport in the private flux region in MAST[END_REF]. The particular choice made for the growth rate ( 19) therefore generates a ballooned structure of turbulent transport that exhibits key similarities with micro-turbulence simulations [START_REF] Galassi | Drive of parallel flows by turbulence and large-scale E × B transverse transport in divertor geometry[END_REF][START_REF] Dif-Pradalier | Evidence for global edgecore interplay in fusion plasmas[END_REF][START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF]. It clearly indicates that using such a growth rate of the form ( 19) is quite effective in providing a means to capture self-consistently a major property of turbulent transport throughout the plasma volume [START_REF] Labombard | Transport-driven scrape-off-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF][START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF][START_REF] Fedorczak | Experimental investigation on the poloidal extent of the turbulent radial flux in tokamak scrape-off layer[END_REF].

Diffusive transport driven by the local dynamics: closure constraint for D ω

To address the closure of the κ-ε model with respect to D ω and V , we consider the local evolution equation for κ and ε (17):

∂ t κ = γ κ κ - 1 D ω κ 2 -ε (20a) ∂ t ε = γ ε ε - V κ 3/2 ε 2 (20b)
The control parameter of the quadratic term of the equation for κ, D ω in (20a), ensures that κ has no fixed point at infinity and remains therefore bounded. Such a property is most important in predator-preys models when all fixed points are lost so that the system must then enter a limit cycle. The parameter D ω is a free parameter, with dimension m 2 s -1 akin to a diffusion coefficient. It can also be written as

D ω = κ ω /∆ω.
The coefficient D ω is split into the frequency width ∆ω of the turbulent energy spectrum and a characteristic value of the latter κ ω , providing a connection to published predator prey models [START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF]. As shown in Section 1, the steady state solution is obtained as the solution of a second order equation. It can be recast into an equation for the diffusion coefficient D = κ 2 /ε:

D 2 + DD ω -D V D ω = 0 (21a) D V = γ κ V 2 γ 2 ε (21b)
Note that for D V to be positive one assumes here that γ κ > 0. this is the condition to have solutions different from D = 0, therefore one positive solution determined by:

D = -1 2 D ω + 1 2 D 2 ω + 4D V D ω (21c)
Depending on V , one finds therefore two regimes, the weak turbulence regime for D V D ω , such that the expansion of the solution for D in (21c) yields D ≈ D V , which scales like V 2 , and a Bohm strong turbulence regime such that D ≈ √ D V D ω and is therefore linear in V [START_REF] Pettini | Chaotic diffusion across a magnetic field in a model of electrostatic turbulent plasma[END_REF]. The ratio D V /D ω , which can be interpreted as a Kubo number [START_REF] Balescu | Anisotropic electrostatic turbulence and zonal flow generation[END_REF][START_REF] Qin | The role of the kubo number in two-component turbulence[END_REF], governs this transition. As a possible closure of the system, we enforce the diffusion coefficient D to be in the weak turbulence regime. In such a framework, we therefore consider that the transition to the Bohm regime for D will occur for much larger values of D than required to match the experiments. Consequently, we choose the parameter D ω such that D ω D n , where D n is the typical diffusion coefficient consistent with the SOL width. The parameter D ω then only plays a role during transients with large κ amplitude. In the simulations D ω is defined in terms of κ max = 10 10 , D ω = κ max /γ κ , where γ κ is a characteristic value of the growth rate. For the sake of simplicity, D ω is also assumed to be constant in the whole simulation volume. In this work we have therefore set the value of the parameter D ω so that is has a minor role on the steady-state behavior of the κ-ε system.

SOL width driven by the local dynamics: closure constraint for V

Given the choice made for the growth rate γ κ = γ ε = γ I (19), the local fixed point yields κ = V 2 and consequently D ≈ D V ≈ V 2 /γ I . The choice for the closure of the control parameter V is a key element of the model since it characterizes the magnitude of the turbulence energy and consequently that of turbulent transport. We have chosen here to define V with the aim of improving the predictive capability of simulations of plasma-wall interaction, in particular divertor operation. The control parameter V is then chosen such that the local evolution properties allow one to recover the SOL width λ SOL as determined by the L-mode empirical scaling [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF]. Since the scaling law is consistent with data from several devices, the choice of V should be appropriate to describe this class of devices with no further free parameter tuned to ensure the match. For the simulations, we have therefore enforced the relationship:

D n = C ν D ≈ λ 2 SOL τ (22a)
where τ = L /c s is the parallel connection characteristic time, c s the sound velocity and L = qR the connection length. The recent scaling laws for the SOL width are based on a proportionality to the so-called poloidal Larmor radius qρ 0 [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF][START_REF] Eich | Inter-ELM power decay length for JET and ASDEX Upgrade: Measurement and comparison with heuristic drift-based model[END_REF][START_REF] Eich | Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER[END_REF]. This scaling property is readily recovered with simple arguments. First the SOL width is a length and consequently must scale as a length. Second, it is understood that microturbulence determines the SOL width so that the characteristic length associated to micro-turbulence, the Larmor radius ρ 0 , is the reference length of interest, consequently λ SOL /a ∝ ρ * where ρ * = ρ 0 /a. However, dimensionless parameters can still play a role. Using the long standing empirical dependence on the plasma current [START_REF] Guilhem | Actively cooled pump limiters and power scrape-off length measurements in Tore Supra[END_REF], hence on the poloidal magnetic field rather than the toroidal magnetic field, one can then readily conjecture:

λ SOL a = α s qρ * (22b)
where α s is a matching scalar parameter. Alternatively, as shown in Ref. [START_REF] Norscini | Self-organized turbulent transport in fusion plasmas[END_REF] and sketched here, one can consider the quasi-linear framework for electrostatic turbulencealthough questionable given the large amplitude fluctuations in the SOL-with

λ SOL /a = (D ql τ /a 2 ) 1/2 . Then for D ql ∝ (ρ * v 0 ) 2 τ , one obtains λ SOL /a ∝ qρ * A, where A = R/a
is the aspect ratio. Such a result holds provided one assumes that the characteristic time τ ql (12c) is identified to be proportional to τ . The question of the SOL width scaling has been a long standing issue and the matter of considerable work as for the ITER design [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF]. Only recently [START_REF] Eich | Inter-ELM power decay length for JET and ASDEX Upgrade: Measurement and comparison with heuristic drift-based model[END_REF] has a particular empirical scaling law been considered to be generic. However, one could argue that the scaling law is assumed and that the significant error bars on the exponents [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF] accommodate for the experimental dispersion. Note that data base used for the empirical analysis [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF] does not account for aspect ratio A variation. The dependence on such a parameter would be useful to discriminate various models. The results discussed here are also tightly related to choosing τ as characteristic time. Given the closure (18) yielding D n = C ν D and that introduced here (22a), we then obtain D ≈ D V :

D = λ 2 SOL τ 2 τ C ν = λ 2 SOL /a 2 τ 2 /a 2 τ C ν = ρ 2 * c 2 s τ α 2 s /A 2 C ν (22c)
One can then determine the control parameter V where the combination of proportionality factors defines the factor

V 0 = α s / √ C ν . V = c s ρ * V 0 A τ γ 2 ε γ κ (22d)
The typical value we thus obtain for the parameter V is therefore ρ * c s , the order of magnitude of the drift velocities. For high-confinement and low-gas-puff plasma regime, the H-mode scaling law leads to λ SOL /a ∼ 2qρ * [START_REF] Eich | Inter-ELM power decay length for JET and ASDEX Upgrade: Measurement and comparison with heuristic drift-based model[END_REF], hence α s ∼ 2. In L-mode, the SOL width is typically 2 -3 times larger with the same scaling [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF], α s ≈ 4 -6. Since the growth rates γ κ and γ ε are defined as local parameters, one readily finds that V (22d) is also a local parameter. The latter appears as the local effective velocity of micro-turbulent transport. For magnetized plasmas, the magnitude of these velocity fluctuations is typically c s ρ * , which is the characteristic amplitude of the E × B electric drift velocity. The transport cross-field velocity will depend on this drift velocity as well as some form of repetition rate of the transport events. Indeed the latter can occur in bursts followed by short quiescent phases [START_REF] Ghendrih | Scaling intermittent cross-field particle flux to ITER[END_REF]. The effective velocity V must account for this behavior, typically multiplying the drift velocity by the ratio of the burst duration divided by the time delay between two bursts of turbulence. We therefore use global features such as the interchange growth rate and the empirical scaling of the SOL width to constrain the local free parameters of the κ-ε model. The latter introduces new possibilities for transport self-organization, but these do not enforce that the global properties of the simulations will exhibit the empirical SOL width scaling. To obtain the expression of V , we have used the characteristic time scale τ = qAa/c s . In fact two time scales appear in the model, τ and the linear time scale associated to the growth rate of the interchange instability, namely 1/γ I [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF].

γ I τ = q R 2 ∇p•∇B T pB T ∝ qA 1/2 (23a)
We assume here that ∇B ≈ B/R while that ∇p ≈ p/a, this scaling property is commonly considered in linear approximations. To take this scaling into account while retaining the dependence on the departure from equilibrium Θ we then define:

γ I τ = qA 1/2 G(Θ) ; Θ = a∇p p • R∇B T B T (23b)
The function G(Θ) is null for Θ < 0 and otherwise equal to √ Θ. One then obtains for V τ the value of V for τ /τ = 1:

V τ = c s ρ * V 0 q 1/2 A 3/4 G(Θ) (24a)
Given a characteristic time τ and the velocity amplitude V one can define the length scales λ b = V τ for ballistic transport and λ d = √ Dτ for diffusion. In the first part of this Section, the time scale τ is chosen to be τ , leading to different ballistic and diffusive expressions of the SOL width, when setting setting γ κ = γ ε = γ I one obtains λ b /λ d ≈ (τ γ I ). However, when retaining τ = 1/γ I the ballistic and diffusive length scales defined above are identical. With the latter assumption, one obtains V γ :

V γ = c s ρ * V 0 q A 1/2 G(Θ) (24b)
For both cases, one has

κ = γ 2 κ V 2 /γ 2 ε , ε = γ κ κ and D = γ κ V 2 /γ 2
ε , so that one has:

D τ = ac s V 2 0 q ρ 2 * A -1 ; τ /τ = 1 (25a) D γ = ac s V 2 0 q 2 ρ 2 * A -1/2 G(Θ) ; τ γ I = 1 (25b)
The choice of the characteristic time scale has an impact of the final result. To account for this fact, we introduce three undetermined exponents ρ , q and A in the definition of D.

D ∝ ac s ρ 2 * ρ q q A A (26)
For τ γ I = 1, q = 2, A = -0.5 while for τ = τ q = 1, A = -1. The exponent of ρ 2 * allows one taking into account a possible ρ * dependence of the characteristic time scale. A gyroBohm scaling, such that τ does not depend on ρ * , is recovered with ρ = 1, and a Bohm scaling with τ depending on ρ * is obtained with ρ = 1/2. Expression [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF] provides a more general framework than has been used in the simulation work presented in this paper. It will be used when comparing empirical scaling laws and that derived from the κ-ε model for the energy confinement time and SOL width. Regarding the simulations, with both the 1D model and SolEdge2D, we have used the closure V τ , slightly modified. Indeed, when G(Θ) → 0, for instance on the high field side, both fixed points values of κ and ε tend to zero. This occurrence is difficult to handle numerically and the expression for the diffusion coefficient is modified to be robust against numerical overflow with D = κ 2 /(ε + ε f ). The parameter ε f is small and is introduced to avoid any spurious divergence when ε → 0. As a consequence, at the limit G(Θ) → 0, one obtains that D → 0 rather than the limit given by (25a). Recovering the latter would be numerically demanding as well as being nonphysical. Further investigation of this point has been achieved and is to be published [START_REF] Auroux | Limit cycles and spreading with k-ε plasma transport[END_REF].

Linear analysis of the κ-ε model

When investigating dispersion equations for plasma instabilities one obtains most often a growth rate that exhibits a threshold. This is a key property of self-organization that is not taken into account by the chosen form of the growth rate [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF] for the present κ-ε model. By analyzing the linearized fluctuation equations we show that the transport terms of the evolution equations (3b) generate a threshold for the effective growth rate. Let us consider the linearized κ-ε model in Fourier space:

∂ t κ + ik u κ + k 2 ⊥ D κ * κ = γ κ κ - 2 D ω κ * κ -ε (27a) ∂ t ε + ik u ε + k 2 ⊥ D ε * ε = γ ε ε - 2V κ 3/2 * ε * ε + 3V κ 5/2 * ε 2 * κ (27b)
where κ * , ε * are the steady state solution homogeneous in space, hence solutions of the local system [START_REF] Kawashima | Simulation of radiative divertor plasmas by Ar seeding with the full W-wall in JT-60SA[END_REF] and where κ and ε are the Fourier modes with wave vectors k and k ⊥ . The parallel velocities are chosen to be constants and the diffusion coefficients D κ * and D ε * are computed with the steady state solution κ * , ε * . The system [START_REF] Baschetti | Optimization of turbulence reduced model free parameters based on L-mode experiments and 2D transport simulations[END_REF] describes the evolution of the system disturbed from the steady state solution by a perturbation varying in space. The dispersion relation is then given by setting the determinant of the linear system to zero and takes the generic form:

γ -A γ -B + C = 0 (28a)
The coefficients A, B, C of this dispersion relation are:

A = γ κ -k 2 ⊥ D κ * - 2 D ω κ * -ik u (28b) B = γ ε -k 2 ⊥ D ε * - 2V κ 3/2 * ε * -ik u (28c)
C = 3V κ 5/2 * ε 2 * (28d)
Stabilizing effects are governed by the diffusion coefficients of the form k 2 ⊥ D * and by parallel transport, the latter depending on the physics controlling the Doppler shift k u . The contribution of the transport terms thus generate damping processes that inhibit the drive terms and consequently govern the occurrence of threshold effects. One then recovers effectively a more complex form of the growth rate, including in particular a threshold in pressure gradient as well as damping processes governed by parallel transport. These are usually taken into account when computing γ I as in Refs. [START_REF] Sarazin | Etude de la turbulence de bord dans les plasmas de tokamaks[END_REF][START_REF] Ghendrih | Generation and dynamics of SOL corrugated profiles[END_REF][START_REF] Garbet | A model for the turbulence in the scrape-off layer of tokamaks[END_REF]. However, since the physics of parallel transport and diffusive crossfield transport appear independently in the evolution equations of κ-ε (17), the present approach consists of using the growth rate proportional to γ I [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF] and let the physics included in the κ-ε system (17) modify this drive. The closure based on the SOL width involves the choice of the characteristic time but also that of the dominant transport process, either ballistic with characteristic velocity V b or diffusive with diffusion coefficient D. The standard approximation for transport is to consider a diffusive process such that the flux depends on a gradient and consequently on the departure from thermodynamic equilibrium. Introducing ballistic transport, typically of the form V b n for density yields transport even without departure from thermodynamic equilibrium, which simply makes no sense unless V b itself depends on a gradient. In the present case, both expressions of V τ (24a) and V γ (24b) depend on the departure from equilibrium Θ (23b). It would therefore be possible to add to the diffusive transport an outward ballistic transport proportional to κ 1/2 since κ = 0 is governed by a departure from equilibrium, either directly via the growth rate of the form 19 or via spreading from a region with positive growth rate. In the present work, for the sake of simplicity, we make the choice of only retaining diffusive transport as done in neutral fluids. Still, the non-linearity of the diffusion coefficient can generate transport properties that are quite similar to ballistic transport. Indeed, the diffusive transport of any field F with diffusion coefficient of the form D = κ 2 /ε yields a divergence of the flux of the form:

-∇ ⊥ D∇ ⊥ F = -D∇ 2 ⊥ F -D∇ ⊥ F 2 ∇ ⊥ κ κ - ∇ ⊥ ε ε (29a)
The divergence is simplified by considering here a slab geometry. Such a structure is comparable to that induced by convective transport, hence driven by a flux of the form -D pinch ∇ ⊥ F + V pinch F , where V pinch must stand for the coupling to another field and be proportional to the gradient of that other field. The divergence of this flux is then:

∇ ⊥ -D pinch ∇ ⊥ F + V pinch F = -D pinch ∇ 2 ⊥ F -D pinch ∇ ⊥ F ∇ ⊥ D pinch D pinch - V pinch D pinch + ∇ ⊥ V pinch F (29b)
Identifying the latter form of the flux divergence with expression (29a) for constant D pinch and V pinch , one obtains

V pinch = D(∇ ⊥ Log(ε) -2∇ ⊥ Log(κ)
). The structure of these two transport models departs when the variation of V pinch , ∇ ⊥ V pinch contributes significantly to the transport. This seems to be a rather particular transport property. But for this specific case, the fact that κ and ε can vary in space generates transport features that are reminiscent of convection although the mathematical structure of the transverse transport term is diffusive. The non-linear dependence of D in κ and ε can therefore generate complex transport properties that depart significantly from the standard case with constant diffusion parameter, such as spreading [START_REF] Singh | When does turbulence spreading matter?[END_REF]. Given these properties, the reduced κ-ε model of cross-field transport can be relevant to model slow transients and not only steady states. This is particularly important for a non-linear system, such as that governing plasma-wall interaction, since one cannot assume that a steady state does exist. Indeed, clear oscillatory behavior has been reported [START_REF] Loarte | Self-sustained divertor plasma oscillations in the JET tokamak[END_REF].

Confrontation of κ-ε transport with scaling laws

We first test the model in a 1D transport model. The results can then be compared to the available empirical scaling laws characterizing global properties such as the energy confinement time and the SOL width. One can also test means to extend the capability of the model towards self-consistent generation of transport barriers Section 3.4.

The 1D κ-ε transport model

The 1D model is obtained by averaging the transport equations in the poloidal and toroidal directions. Transport is taken into account by diffusion that replace all convective transport processes at the microscopic level. The diffusion coefficients are chosen as discussed in Section 2.1 equation [START_REF] Rognlien | A fully implicit, time dependent 2-D fluid code for modeling tokamak edge plasmas[END_REF]. The plasma transport equations for the particle density, and the density of electron and ion thermal energy are presented in Appendix A. In these equations a mask function defines the SOL region for r ≥ a where parallel losses to the wall apply. The plasma equations are completed by a vorticity equation (A.1b). It will be used when addressing the impact of velocity shear on the transport properties in Section 3.4. A simplified transport model for neutrals is also introduced so that the particle source by ionization, and in particular its location, changes consistently with the plasma parameters. A feed-back controlled gas injection Φ N 0 is set to ensure that the separatrix density is fixed, equal to 1. 10 19 m -3 . The cross-field diffusion coefficients of plasma transport are determined by the κ-ε coupled equations ( 17), restricted to the 1-D radial direction (A.1d). A core energy source with given radial profile is used. The scan in additional heating power is performed by varying the source amplitude (A.2b).

Impact of a parameters scan on the SOL width

The SOL width estimated by the numerical model is compared to the empirical scaling law for the heat flux width λ q given by ref. [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF] in mm for L-mode discharges:

λ q = (1.44 ± 0.67)B -0.8±0.32 T q 1.4±0.67 P 0.22±0.1 in R -0.03±0.28 (30a)
As shown in Appendix B, this scaling law is not dimensionaly correct since the adimensional ratio λ q /a still depends on B T as well as on the dimensionless parameters q, A, β, ν * and ρ * .

λ q a ∝ q 2.15 A 0.03 β 0.46 ν -0.17 * ρ 0.85 * B 0.156

T (30b)
Within the error bars of (30a), the only modification of a single exponent that allows one recovering a dimensionaly correct expression is to change the exponent of the power law on B T from -0.8 to -0.8 -0.156 = -0.956 so that the same expression (30b) is obtained but with an exponent for B T equal to zero. In order to determine the dependence of λ q /a on the dimensionless parameters, we use the dimensionless scaling law for the energy confinement time ITER96-th [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF]. It is to be remarked that also for this empirical scaling law, a correction on the exponent of R is made to achieve a dimensionaly correct expression [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF], see also Appendix B. In both cases, the modification of the quality of the scaling law by the dimensionality constraint, the so-called Kadomtsev constraint [START_REF] Lacina | Similarity rules in plasma physics[END_REF][START_REF] Kadomtsev | Tokamaks and dimensional analysis[END_REF][START_REF] Connor | Scaling laws for plasma confinement[END_REF], is lacking. Finally, it is to be underlined that the assumption λ SOL /a ∝ qρ * used to define the free parameter V , Section 2.5 (22d) is not recovered here since one obtains typically λ q /a ∝ q 2.15 ρ 0.85 * . There seems to be a consistency issue between the assumed dependence λ q /a ∝ qρ * as claimed in the literature and the effective scaling law and that stemming from the proposed scaling law.

Since the decay rate of the energy flux in the SOL region depends on the parallel heat conductivity, which is not taken into account by the 1D model, we determine here the e-folding length of the pressure profile λ p and use it as a proxy for λ q . The fall-off region of the pressure profile used to determine λ p is restricted to the vicinity of the separatrix, typically the first 20% of the SOL width. Based on the oversimplified use of e-folding lengths of density and temperature, one can predict that λ q /λ p ≈ 0.75 ‡, this ratio being a constant depending on the sheath transmission factor and respective heat and particle cross-field transport. The SOL width determined by this method could overestimate the actual width but should exhibit the same parameter dependence.

In a first series of simulations analyzed here, we perform a scan of the input power

P in ∈ [0.1, 0.5, 1, 2, 3, 4, 6] M W .
The lowest values of P in of 0.1 and 0.5 M W are not realistic but useful when analyzing the results in terms of scaling laws. The particle source is adjusted to impose the same density at the separatrix n a = 10 19 m -3 . One can then expect that the neutral penetration into the plasma will be roughly constant so that the ionization source will have the same shape but varying amplitude to match the particle outflux. Experimentally this would correspond to a feedback process on a density measurement at the midplane separatrix, a control scenario particularly suitable to investigate divertor physics [START_REF] Baschetti | Radiation Driven Bifurcations in Fusion Plasmas[END_REF]. With such a scenario for the power scan one readily expects an increase of the plasma thermal energies T e and T i with P in as well as an increased SOL width given the scaling law (30a). Based on a straightforward analysis [START_REF] Engelhardt | Influence of an ergodic magnetic limiter on the impurity content in a tokamak[END_REF], one can show that the ratio between core and separatrix density n 0 /n a is expected to be n 0 /n a = 1 + λ I /λ SOL . The neutral ionization mean free path λ I determines the radial distance between the particle source and the separatrix. One thus finds that the source is all the more effective that it is closest to the plasma center. This effect is balanced by particle transport characterized by λ SOL . The main change when scanning the power is then the impact on the transport features characterized by λ SOL . This effect drives a reduction of the core density with increasing power. These trends are recovered in the simulations, Figure 5. On the top-left figure labeled (a), are plotted the outer density profiles, these being characterized by flat core profiles. As the power is increased the core density decreases while the density profile at the separatrix becomes flatter. The ratio n 0 /n a is observed to decrease -n a the separatrix density is feedback controlled to be constant-when increasing P in , Figure 6 left hand side. However, one can observe that the density decrease is less pronounced for P in ≥ 1M W compared to the very low power cases. Conversely, the thermal energies, Figure 5 top right labeled (b), increase with the injected power. Analyzing the core T 0 and separatrix T a thermal energy, one finds that the temperature profile scales with the injected power, T a ∝ P 0.537 in and T 0 ∝ P 0.513 in , 6 ‡ This value is estimated using an e-folding length for both density λ n and and temperature λ T , and given equal ion and electron temperatures. Assuming that the profile of the heat flux is governed by the sheath conditions, one then obtains λ q /λ p = (1 + λ T /λ n )/( 32 + λ T /λ n ) equal to 2 3 for λ T /λ n → 0 and increasing to 1 for λ T /λ n → +∞. The quoted value λ q /λ p = 0.75 is obtained for λ T /λ n → 0.5. right hand side. The overall effect on the plasma pressure p = nT is an increase. In the core, Figure 6 left hand side, the pressure increase p 0 does not exhibit a power law because of the particular response of the core density. At the separatrix, where the density is maintained constant by the chosen feedback, the pressure is governed by the behavior of the thermal energy, and exhibits the same power law with respect to the injected power. The pressure gradient at the separatrix drives the turbulent energy κ as well as the turbulent energy transfer rate ε, these determine the plasma transport and consequently the pressure gradient at prescribed fluxes: feedback controlled for the particles and externally controlled for the energy. The particle radial diffusion D n governed by the fluctuating convection process is a key element of this non-linear loop: see Figure 5 bottom left labeled (c) where the profiles of D n are plotted for different powers P in in the vicinity of the separatrix. The profiles are similar in the region r ≤ a, and exhibit differences in the SOL r > a and in particular in the far SOL region at lowest power in regions where the plasma is weakly ionizing. Considering the values at the separatrix, one finds that D n exhibits a power law dependence on the injected power, D n ∝ P 0.644 in , Figure 7 left hand side. Furthermore, from the density profiles Figure 5 top-left, one can determine the density gradient length L n = (∇ r n/n) -1 . The latter is also observed to scale with the injected power, L n ∝ P 0.263 in , Figure 7 left hand side. Assuming that this scaling is consistent with the choice made to determine the parameter V (22d), one expects L n /a ∝ qρ * and thus L n ∝ T 1/2 a , yielding the power law exponent 0.269 in good agreement with that obtained for L n . One can also compare the scaling law for L n to that observed experimentally as described by the scaling law of λ q (30a). Balancing the cross-field diffusive flux and the parallel losses, one expects the exponent of L n to be 0.19, half of the exponent of D n minus one fourth of the exponent of the thermal energy T a §. Both exponents stand within the error bars of the experimental scaling law (30a). Their difference stems from the particle source in the SOL that must be taken into account for the steady state particle flux balance together with the parallel and cross-field fluxes. One can expect this particle source to exhibit another power law dependence on P in than the two flux divergence terms. One thus finds that the SOL behavior is governed by the balance of three terms that are characterized by different underlying scaling laws. The overall behavior can only exhibit a power law, governed by a single exponent, on a restricted range of P in . The amplitude of the ionization source term can be addressed by analyzing the cross-field particle flux at the separatrix that balances in steady state neutral influx. The flux of neutrals that are ionized in the edge plasma is therefore determined by Φ n = D n n/L n at the separatrix, Figure 7 right hand side. This flux, and consequently the particle source, is governed by the feedback loop that maintains the separatrix density n a = 10 19 m -3 . Owing to the enhanced particle transport with the injected power, the neutral particle influx must increase with the latter and indeed one finds Φ n ∝ P 0.386 in . A last observation from the simulation is the plasma electric potential profiles shown on Figure 5 bottom-left labeled (d). One finds that a potential well develops in the core, typically governed by the ion temperature profile, while the SOL values are positive and slaved to the electron temperature in the SOL, the outer wall being grounded and defining the reference electric potential. As a consequence of this drive, a maximum E × B shear region is located in the vicinity of the separatrix where the electric field reverses.

Let us now address the scaling of λ SOL determined by the pressure e-folding length λ p and compare it to the experimental scaling law (30a). On Figures 8 and9 are plotted the variation of λ p from the present κ-ε 1D model -black curves open squaresand power law variations: the experimental scaling law -dash dot black curve, and a fit of the results of the 1D model, dashed blue curve. The gray-shaded area stands for the 95% confidence interval of the regression used to determine λ p .

parameter 1. Exponents of the power law dependence on the engineering control parameter for the fit of the κ-ε simulations and the experimental scaling law. The second value of the scaling law coefficient for B T is that proposed to recover a dimensionally correct experimental scaling law.

P in B pol R B T κ-ε
Table 1 summarizes the various results for the main parameter scans. One can notice that with the chosen closure, the κ-ε model exhibits the same trend as the experimental § One considers here L 2 n = D n L /c s and therefore

L n ∝ D 1/2 n T -1/4 .
scaling law since all the power law exponents have the same sign. However, the magnitude of these exponents observed in the κ-ε simulations are systematically smaller. We first discuss the results for the power scan, Figure 8 left hand side. The values obtained with the simulations are typically a factor 3 larger than the values given by the scaling law. Furthermore, there is a factor two in the exponents. For this parameter, the match between the experimental scaling law and the κ-ε simulations is rather poor.

The dependence of λ q on the magnitude of the poloidal magnetic field, hence on the plasma current, is governed by that on q in (30a) when the other engineering parameters are kept constant, hence λ q ∝ B -1.4 pol . In the experimental findings on a single machine this is the main parameter dependence as evidenced in [START_REF] Guilhem | Actively cooled pump limiters and power scrape-off length measurements in Tore Supra[END_REF][START_REF] Gunn | Scrape-off layer power flux measurements in the Tore Supra tokamak[END_REF]. The latter trend, typically λ q ∝ B -0.5 pol reported in [START_REF] Guilhem | Actively cooled pump limiters and power scrape-off length measurements in Tore Supra[END_REF][START_REF] Gunn | Scrape-off layer power flux measurements in the Tore Supra tokamak[END_REF], or the direct behavior stemming from λ q ∝ q, hence λ q ∝ B -1

pol is approximately recovered with the κ-ε simulations, Figure 8 right hand side. However, this exponents departs by nearly a factor 1.5 from that in (30a). The dependence on R, keeping all other engineering parameters constant, hence minor radius, poloidal and toroidal magnetic field, density and input power, Figure 9 left hand side, is mainly governed in (30a) by the dependence of q on 1/R, leading to λ q ∝ R -1. 43 . Again the trend of the scaling law and the simulation agree, see Figure 9 left hand side, but the power law behavior is clearly different. Finally, when varying the toroidal magnetic field alone, Figure 9 right hand side, the κ-ε model yields a rather weak dependence while that of the scaling law is a factor two larger. The latter combines the explicit dependence on B T as well as that embedded in the safety factor leading to the effective power law B 0.6

T . With the proposed correction such that the experimental scaling law fulfills the Kadomtsev constraint, the exponent for the B T scaling is decreased to 0.444 reducing the mismatch with the κ-ε result.

When addressing the comparison between the κ-ε results and the scaling law (30a) several issues are raised which can readily explain the observed differences. First, as we have just recalled, the reference expression is not dimensionally correct, second, when written in terms of the standard dimensionless parameters, one does not recover the assumed underlying form, and, third, the error bars on the coefficients of the power law are considerable. These issues can explain the mismatch, or when taking into account the error bars, enable finding some form of agreement despite significant differences in behavior and predictions. The κ-ε results are further discussed in the following Section considering the core energy confinement together with the SOL width.

Physics background of the scaling law for the SOL width and global energy confinement time

Considering the two definitions τ E ≈ a 2 /D core and λ 2 q = τ D SOL one can readily determine the relationship between the energy confinement time and the SOL width, τ E = τ (a/λ q ) 2 D SOL /D core , therefore:

τ E λ 2 q ∝ τ a 2 D SOL D core (31a)
When assuming the same scaling properties for D SOL and D core , one can then estimate λ q given the scaling of τ E or vice-versa that of τ E given that of λ q .

τ E ∝ τ a 2 λ 2 q ≈ q 2 R Ω 0 λ q R 2 A 2 λ 2 q = (qA -1 ) 2 (Rλ -1 q ) 3 Ω -1 0 (31b)
Such a connection between the global energy confinement time and the width of the boundary layer has been investigated experimentally [101]. With the approximation λ q ≈ qρ 0 with ρ 0 = c s /Ω 0 , Ω 0 being the reference ion Larmor gyration pulsation, one obtains the confinement time τ

(qρ 0 ) E,λ : Ω 0 τ (qρ 0 ) E,λ ∝ q -1 A ρ -3 * (31c)
Using the scaling of λ q (30a) obtained for L-mode operation, one thus obtains a scaling law for

τ L E τ L E,λ ∝ q -2.2 A -2 R 3.09 B 1.4 P -0.66 in (32a)
Recalling the approximation made for the safety factor 1/q ≈ A 2 I p /(RBk), where I p is the plasma current and k the elongation, see Appendix B, the scaling for τ E derived from that of the SOL width is:

τ L E,λ ∝ I 2.2 p B -0.8 P -0.66 in R 0.89 A 2.2 k -2.2 (32b)
The latter expression is to be compared to a scaling law for the energy confinement time, such as the ITER96-th scaling law for L-mode [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF]:

τ L E,th = 0.023 I 0.96 p B 0.03 P -0.73 in R 1.83 A 0.06 k 0.64 n 0.4 M 0.2 (33a)
where τ E is obtained in Apart from a similarity regarding the scaling on the injected power P in , the agreement between the two expressions is poor, a major difference being the difference in scaling with the plasma current governed by the strong dependence of the λ q -scaling on the safety factor. Using the scaling of τ E , one can also determine that of λ q : λ E,L q = q 0.99 P 0.26 n -0.13 M -0.07 R 0.07 A -0.05 k -0.63 (33b)

As readily expected, this scaling law exhibits a weaker dependence on the safety factor q, and a similar scaling on the injected power P in . The other engineering scaling factors appear to have a small impact but for the elongation k at given safety factor. Plasma shaping appears to enhance confinement (33a) which consequently drives a reduction of the SOL width (33b).

In the framework of the κ-ε approach, one can also estimate the energy confinement time given the diffusion coefficient D (25) stemming from κ and ε. The latter are governed by the the closure constraint on the SOL width that is used to determine V .

For the chosen scaling relationship of D V [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF], one can estimate the confinement time τ

(κ-) E considering that τ (κ-) E ∝ a 2 /D V : Ω 0 τ (κ-) E ∝ ρ -1-2 ρ * q -q A -A (34a)
The Larmor pulsation Ω 0 = eB/m is used to enforce the standard time normalization. The values obtained for the different exponents ρ , q and A for the different scaling laws and for the κ-ε model are summarized in Table 2. As shown in Appendix B, the adimensional parameter form of the τ L E,th scaling is:

Ω 0 τ L E,th ∝ ρ -1.85 * q -3.74 A 0.37 β -1.41 ν 0.19 * k 3.22 (34b)
The expression of ν * used here is the collision frequency normalized by the passing particle transit frequency 1/τ . Therefore, it does not contain the aspect ratio dependence governed by the trapped particles physics so that ν * is proportional to A. Because of this difference in definition, the scaling exponent on the aspect ratio is different from that of the ITER physics basis [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF]. Comparing (34a) and (34b) then yields the condition ρ = 0.425, q = 3.74 and A = -0.37 to match the power law dependencies. For the two characteristic times addressed above, one finds that the choice τ ∝ 1/γ I discussed at the beginning of this Section departs less from the present results than τ ∝ τ used in Section 2.2. While the transport governed by the local diffusion coefficient derived in the κ-ε framework is typically GyroBohm ρ = 1, the value obtained to match the L-mode scaling law is not even Bohm-like ρ < 1/2. The β dependence also leads to a significant difference. However, it is to be mentioned that the dependence on β is presently being questioned [START_REF] Sips | Assessment of the baseline scenario at q 95 ≈ 3 for ITER[END_REF]. Finally, the ν * dependence is relatively small and has therefore a rather weak effect. Given the differences, in particular regarding the β and ρ * scaling, it is interesting to compare the κ-ε expression of τ E (34a) to the latest analysis of the empirical data base of the energy confinement time [START_REF] Sips | Assessment of the baseline scenario at q 95 ≈ 3 for ITER[END_REF].

Ω 0 τ H E,th ∝ ρ -3.0 * β 0.0 ν -0.14 * q -1.7 k 3.22 A 0.9 (34c)
This expression is slightly different from that in Ref. [START_REF] Sips | Assessment of the baseline scenario at q 95 ≈ 3 for ITER[END_REF] because in this reference ν * includes the aspect ratio dependence of the trapped particle physics but does not appear to include the dependence on the safety factor. We have thus modified the scaling law to take into account the differences in the definitions of ν * . This scaling law being GyroBohm, one readily finds ρ = 1. One can also identify q = 1.7 and A = -0.9 that are somewhat closer to the values obtained with τ ∝ 1/γ I , q = 2, A = -0.5, than with τ ∝ τ , q = 1, A = -1, see Table 2. This is when having in mind that the dependence on the safety factor is based on a more complete experimental evidence than that on the aspect ratio. The fact that the new scaling law of the energy confinement time does not depend on β, makes the connection between energy confinement time and SOL width more relevant since the SOL width does not appear to depend on β.

The 1D model allows one comparing the predicted energy confinement time to that scaling law

ρ q A IT ER -L [91] 0.425 3.74 0.37 IT ER -H [102] 1 1.7 -0.9 κ -ε, τ ∝ 1/γ I 1 2 -0.5 κ -ε, τ ∝ τ 1 1 -1
Table 2. Given the definition D = ac s ρ 2 ρ * q -q A -A , see [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF], leading to (34a), values of the exponents ρ , q , A for the ITER confinement time scaling laws compared to values obtained with the κ-ε model assuming λ/a ∝ qρ * (22b).

achieved in the steady state 1D κ-ε simulations. The appropriate trend of the dependence on the input power P in is recovered, Figure 10 left hand side. The power law fitting the numerical results, black open squares, is P -0.64 in blue short-dash line, comparable to P -0.73 in given by the ITER L-mode scaling (33a), dashed black line. One finds that the simulation data also compares well with the power law (32b) stemming from the empirical SOL width scaling P -0.66 in , dash-dot black line. The gyroBohm scaling of transport used for the κ-ε model, hence for ρ = 1 in (34a), yields P -0.6 in blue longdash line. The scan in B pol , hence in I p appears to qualitatively right, Figure 10 right hand side but the quantitative power law fitting leads to more discrepancy The model data, open black squares, is fitted by the power law B 0.98 pol short dash blue line, which is quite close to B 0.96 pol given by the ITER L-mode scaling (33a). However, the power law computed using various closures do not give such a good match. Using the SOL width scaling law one obtains B 2.2 pol (32b) black dash-dot line. The theory based κ-ε models (34a) yield B 0.8 pol for q = 2 long-dash blue line B 0.4 pol for q = 1 dotted blue line. For such models the perfect fit would be obtained with q = 2.4. One finds therefore that the characteristic time that must be considered to recover the turbulent transport properties appears to be more difficult to define and understand than the properties governed by the characteristic amplitude of the fluctuating velocity.

Including turbulence shearing by large scale flows

The proposed κ-ε model for plasma turbulent transport is closely linked to that used in neutral fluids. In particular, the field ε governs the energy cascade providing conditions for steady state turbulent spectra. However, a key physics mechanism that controls plasma turbulent transport is eddy shearing by large scale flows, typically the zonal flows [103]. In particular, the regimes of improved confinement reported in magnetically confined plasmas are observed to be governed by turbulence shearing by large scale flows, see Ref. [START_REF] Burrell | Role of sheared E x B flow in self-organized, improved confinement states in magnetized plasmas[END_REF] and references therein. A means to include such a mechanism is to introduce a second predator of the energy of the turbulence, typically the zonal flows [START_REF] Diamond | Self-regulating shear flow turbulence: A paradigm for the L to H transition[END_REF]103,[START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF][START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF]. Alternatively, we propose to take into account the shear stabilization in a simplified way. The previous Sections indicate that the physics enforced by the κ-ε model of this paper can be addressed in terms of the effective microscopic velocity V that governs 32b) is plotted with dash-dot black lines while the theory based κ-ε models (34a) are plotted with blue long-dash lines for q = 2 and blue dotted line for q = 1. Note that for the injected power scaling law, these two cases yield the same behavior.

the convective transport at the smallest scales. As discussed, this velocity is found to be of the order of magnitude as the E × B drift velocity, thus proportional to c s ρ * . We now change the definition of V to take into account the large scale velocity shear effect. Let ω z be the radial shear of the zonal flows and let τ z be the characteristic time to enforce the shearing effect, the new dimensionless control parameter is then (ω z τ z ) 2 since the shearing stabilization does not depend on the direction of the large scale flow, and therefore does not depend on the sign of ω z . We then write the modified expression of V as:

V = c s ρ * V 0 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 q ηq A η A G(Θ) 1/2 (35a)
The parameter r governs the ratio between the two asymptotic limits, the small shear limit (ω z τ z ) 2 → 0, such that the term standing for the shearing effect is equal to 1, and the large shear limit, (ω z τ z ) 2 → +∞, such that this factor is equal to 1/r. The values r > 1 govern a reduction of the turbulent energy amplitude. The proposed dependence of V on the shearing efficiency (ω z τ z ) 2 is devised as a simple expression that governs a reduction of V by the factor r between the no-shear asymptotic limit and the strong shear asymptotic limit. The fixed point for κ-ε is such that κ = V 2 and ε = γ I κ so that D = κ/γ I . Given 1/γ I = a √ A/c s , the modified diffusion coefficient is then:

D = V 2 a √ A c s = c 2 s ρ 2 * a √ A c s G(Θ) 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 2 q 2ηq A 2η A = ac s ρ 2 * V 2 0 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 2 q 2ηq A 2η A -0.5 (35b) λ = a qρ * 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 q ηq-1/2 A η A -3/4 (35c)
In the 1D simulation case the parameter ω z is defined given the electric potential profile solution of the vorticity equation (A.1b). Two free parameters then remain to be fixed, τ z that governs the sensitivity to the shearing mechanism and r that governs the amplitude of the transition from the low to the high shear limits. For a factor two change in the diffusion coefficient, one must choose r ≈ √ 2. For the more important variation between the L and H mode characterized by a decrease of the SOL width by a factor 2 to 3, as discussed previously for the parameter α s , one must consider r ≥ 2. However, the modification of the ionization source in the core plasma, and consequently modified density gradients in the edge and SOL region also contribute to the changes of κ and ε making the overall effect more complicated than the argument solely based on the value of the correction of V due to the flow shear. The modified model is found to undergo a transition from the no-shear regime to the strong shear regime when the injected power is increased Figure 11. The effect of shear appears to be localized just inside the separatrix, close to the inversion point of the electric field Figure 11 left hand side. Indeed, the profiles appear to undergo a transition between the low power regime, closed symbols, and the high power regime, open symbols. Steepening of the ion thermal energy gradient is observed in the vicinity of the separatrix, a behavior that is reminiscent of an interface barrier [START_REF] Norscini | Self-organized turbulent transport in fusion plasmas[END_REF]. The mechanism driving the shear layer is the transition from the core where the electric potential is typically proportional to -T i /e to the SOL region where the sheath boundary condition lead to an electric potential proportional to T e /e. The increase of the electric field in the vicinity of the separatrix is then governed by enhanced temperature gradients in the edge when increasing P in . When plotting the variation of the plasma diffusion coefficient D n with injected power P in , Figure 11 right hand side, One can observe that for the inner most radial position ρ = 0.96, D n increases monotonically with P in and as the position gets closer to the separatrix, a decrease of D n sets in for 2 M W ≤ P in ≤ 2.5 M W . The enhanced gradient is driven both by the heating source and by the decrease of the turbulent transport D n driven by the velocity shear. The latter depending on the gradients as well as on turbulent transport provides a mechanism for a nonlinear feedback process, and therefore a transition to improved confinement driven by the heating power. The modification of turbulent transport governed by the flow shear we have presented here is an ad hoc way of generating a feedback on turbulence and governing a further process of self-organization of plasma transport. The transition to barrier formation could then be modeled self-consistently provided the free parameters can be determined by matching the simulations to experimental evidence. It is to be underlined that the chosen expression of V determines the impact of the flow shear on transport. It is devised to govern a reduction of the cross-field diffusive transport. The observed transition follows the standard feedback loop where the reduction of transport drives steeper gradients that in turn enhance the shear stabilization [START_REF] Burrell | Role of sheared E x B flow in self-organized, improved confinement states in magnetized plasmas[END_REF]. However, the reduced model used here does not account for an initial step enforced by Reynolds's stress generation of zonal flows [START_REF] Manz | Spatial redistribution of turbulent and mean kinetic energy[END_REF][START_REF] Held | Non-Oberbeck-Boussinesq zonal flow generation[END_REF]. Alternative reduced models governing the local behavior of κ and ε driven by this physics have been proposed to investigate the transition to transport barriers [START_REF] Diamond | Self-regulating shear flow turbulence: A paradigm for the L to H transition[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF][START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF]. It is to be noted that in the ad hoc model proposed here, a modest change in the magnitude of V with a monotonic dependence on the shear, is enough to drive a bifurcation behavior. This result must be confronted to experimental evidence and matching parameters tuned to obtain a universal behavior. If this can be achieved, this minimal model could be used to investigate the compatibility of high performance divertor regimes preserving the H-mode confinement in the core.

SolEdge2D-Eirene modeling of a WEST discharge

This Section introduces the first results provided by the 2D κ-ε model of a WEST experiment [START_REF] Bufferand | Density regimes and heat flux deposition in the WEST shallow divertor configuration[END_REF][START_REF] Bourdelle | WEST physics basis[END_REF] in the lower-single-null configuration, discharge #55049. The transport code SolEdge2D-Eirene [START_REF] Bufferand | Numerical modelling for divertor design of the WEST device with a focus on plasma-wall interactions[END_REF] is used for pure deuterium plasma (no impurity physics), and without the mean field velocity drifts. The equations are recalled in Appendix C. The coefficients of the diffusive and conductive perpendicular fluxes (D n , χ e , χ i and ν) are self-consistently estimated using the κ-ε model, Section 2. The equations are discretized using a second-order finite volume scheme associated to a volume penalization technique [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] on a structured mesh based on a grid aligned on the magnetic flux surfaces for numerical efficiency, Fig. 12. Typically 10 6 iterations Figure 12. Left hand side: Example of mesh decomposition for a WEST magnetic equilibrium with double X-point [START_REF] Bufferand | Density regimes and heat flux deposition in the WEST shallow divertor configuration[END_REF][START_REF] Bourdelle | WEST physics basis[END_REF]. Each sub-domain is characterized by a different color. The penalization technique allows one to add an axisymmetric object, such as a baffle or a toroidal secondary limiter within this computational domain and investigate its impact on plasma-wall interaction. Right-hand side: corresponding mesh for Eirene, overlapping the plasma domain in black, and in regions shielded from the plasma in red within the circular vacuum vessel.

are needed to reach steady state conditions with, as well as without, the κ-ε model. The plasma parameters for this WEST discharge as well as the control parameters of the simulation are presented in Table 3 andTable 4, respectively.

As expected from the chosen dependence of the growth rate [START_REF] Heifetz | A monte-carlo model of neutralparticle transport in diverted plasmas[END_REF], one finds that the SolEdge2D-Eirene simulation of the WEST shot 55049 exhibits a clearly ballooned pattern for the field κ localized on the low-field side of the torus, Figure 13 left hand side. Using a logarithmic scale for κ normalized by its maximum value, allows one to highlight the difference of magnitude and to observe spreading governed by parallel transport in the near SOL region of the turbulence energy κ. The far SOL region on the high field side, with near zero plasma pressure, exhibits spurious values of κ, also highlighted by the logarithmic scale. These have very weak effect on the plasma properties. The ballooned pattern is recovered for the diffusion coefficient D, Figure 13 right hand side. In the radial direction, D becomes small in the outer SOL, beyond the secondary separatrix and up to the wall. One can also notice a spreading of D in the Table 3. Main parameters of the shot W EST -55049. Top line particle and energy control: line averaged density obtained by feedback control ne , additional heating power P add by Lower Hybrid, split onto two antennas, core P core rad and edge P edge rad radiated power. Middle line: magnetic equilibrium parameters, plasma current I p , toroidal magnetic field B T , major R and minor a radius. Lower line: dimensionless control parameters, safety factor q 95 , elongation k, aspect ratio A and X-point height from the divertor floor h X normalized by the elongation times the minor radius k a.

ne [m -2 ] P add [M W ] P core rad [M W ] P edge rad [M W ] 4.0 10 19
4.0 2.0 0.6 T e are shown on Figure 14 at the outer midplane and at the divertor. A logarithmic scale is used to highlight the SOL part of the profiles. The profiles indicate that the thermal energy and density only weakly vary from the midplane to the divertor that corresponds to the so-called sheath limited or hot divertor regime. Let us notice that in the experiment however one observes a factor two drop in the thermal energy from midplane to the divertor, this is closer to the so-called high recycling regime. The difference appears to be related to the particle recycling pattern and via the boundary condition to a change in the parallel transport behavior. This issue will be further analyzed elsewhere. For the present comparison, the shape of the divertor profiles, which is governed by the cross-field transport and thus the κ-ε model, is more relevant than the actual values. One can notice the sharp drop of the divertor profiles towards the private flux region and slightly stronger gradients in the outer SOL, R -R sep ≥ 0.03 m.

I p [M A] B T [T ] R[m] a[m] 0.5 3.6 2.50 0.5 q 95 k A h X /(k a)[m]

Midplane profiles

The radial density and electron thermal energy profiles are compared to the experimental profiles at the outer midplane, Figure 15. One can notice the agreement between the measurement points and the simulation output in the near SOL region 0.01 ≤ R -R sep ≤ 0.03 m. The SolEdge2D-Eirene profile of the electron thermal energy exhibits a too strong gradient, in particular in the outer SOL 0.03 ≤ R -R sep ≤ 0.05 m, driving a departure between the simulation profile and the measurements. For the density, the mismatch is restricted to the outer SOL 0.03 ≤ R-R sep ≤ 0.05 m where the simulation profile exhibits a too small decay length. These differences are highlighted by the logarithmic scales, and are hardly noticeable with linear scales. Given the fact that there is a single free parameter to tune all crossfield transport in the whole simulation domain the overall agreement is remarkable. 4.0.2. Divertor profiles Figure 16 and Figure 17 show radial divertor profiles, remapped to the midplane in terms of the distance to the separatrix, the SOL profile connected to the midplane for R -R sep ≥ 0 and for R -R sep ≤ 0 the profile towards the private flux region. The figures use left and right axis presentation, with on the left hand side axis the scale for the WEST experimental data labeled "WEST", closed blue point dashed line, and on the right hand side axis the scale for the SolEdge2D-Eirene output, labeled "SE2D", black plain line. The difference between these two scales underlines the shortfall regarding the magnitude, while splitting the results between left (experimental) and right (simulation) scales highlights the profile similarities. On Figure 16, are presented the parallel fluxes to the divertor, the saturation current j sat proportional to the particle flux, left hand side panel, and the parallel heat flux Q onto the divertor target plate, right hand side panel. The heat flux pattern exhibits comparable magnitude with a peak value for the simulation typically 10 % smaller than the experimental measurement. Agreement regarding the shape is also found with close to exponential fall-off in the SOL region and sharp decay towards Regarding the ion saturation current, the experimental peak value is typically 1.6 times the value achieved in SolEdge2D-Eirene. The simulation is thus characterized by a too large fueling efficiency, and consequently, smaller particle flux to the target plate for the same core density. Since the electron thermal energy T e is proportional to the divertor heat flux divided by the particle flux, one expects a factor 1.9 drop of T e between simulation and experiment, Figure 17 right hand side. In practice one finds a ratio of the order of 2.2, therefore Table 5. SOL width for WEST: λ SE2D pe from the pressure e-folding length given by SolEdge2D-Eirene in the plasma midplane, λ W EST pe from the pressure e-folding length given by the reciprocating Langmuir probe data in the plasma midplane, from λ SE2D SOL the heat flux channel width in SolEdge2D-Eirene simulation, λ W EST SOL the heat flux channel width from Langmuir probe embedded in the divertor. SOL width prediction for WEST according to available scaling laws, λ H according to [START_REF] Eich | Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER[END_REF], λ L = 2 λ H , λ q according to [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF], λ max q , λ min q maximum and minimum value according to [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF] taking the error bars into account. All values in [mm]. 16 left hand side should behave as the inverse of the thermal energy, and one finds indeed that the density is about 2.3 times larger in the experiment compared to the simulation together with a difference in shape that is the opposite to that seen for the electron thermal energy. Although, there are differences, one finds consistency most likely governed by the conservation of the total plasma pressure along the field lines and by the agreement between simulation and experiment of the plasma energy flux onto the divertor.

λ SE2D pe λ W EST pe λ SE2D SOL λ W EST SOL | λ H λ L λ q λ max q λ min q 8.
Regarding the various ways of determining the SOL width, from e-folding length of the plasma pressure in the midplane, given by SolEdge2D-Eirene λ SE2D SOL , or by the reciprocating Langmuir probe λ W EST pe , or from the heat flux deposition pattern, λ SE2D SOL given by SolEdge2D-Eirene, and λ W EST SOL given by the probe data, Table 5 yield quite comparable values 8 ± 0.8 [mm] in the SolEdge2D-Eirene simulations and 8.2 ± 0.9 [mm] for the WEST experiment. This value agrees with the scaling law λ H = 0.64 B -1.15 pol [START_REF] Eich | Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER[END_REF] for the SOL width in H-mode, with a factor 2 increase from H-mode to L-mode, one obtains λ L ≈ 2λ H = 8.2 [mm]. The latter value agrees with both simulation and experimental results. Conversely, the scaling law for L-mode [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF] leads to λ q which is about 2.8 times smaller than the observed value. Furthermore, the considerable uncertainty governed the error bars on the exponents leads to a range of values from 0.3 [mm] to 18 [mm]. One finds that the SolEdge2D-Eirene simulation with the κ-ε transport model has better predictive capability than the scaling law built for L-mode data. The latter appears to exhibit a too strong dependence on the safety factor. The uncertainty on the exponents and that for the safety factor is also an issue. Finally, the WEST experiment analyzed here is characterized by a large aspect ratio. The dependence on that parameter is not well known because of the lack of appropriate data in the databases that have been used. As discussed in Section 3.3, a complete picture of cross-field dependence on aspect ratio and safety factor is still missing. This issue is partly made more complicated because of the dependence of the safety factor on the aspect ratio A, q cyl A = B T /B pol , where q cyl is the cylindrical approximation of the safety factor. The open issue appears to be therefore the balance between magnetic field properties B T /B pol and geometry A = R/a.

With the confrontation to the WEST experiment we thus find that the κ-ε cross-field model of transport is robust and compares well with experimental evidence despite the fact that the simulation exhibits differences in the divertor regime, most likely governed by a reduced particle screening capability in the simulation compared to the experiment.

Discussion and conclusion

In the present paper we have developed the theoretical background of the κ-ε model, and we have shown its relevance to quite significantly improve the predictive capability of 2D plasma transport models based on averaged fluid equations. Assuming turbulent diffusive transport, the model allows to self-consistently estimate the cross-field fluxes in the edge and SOL regions of diverted plasmas. Instead of tuning by hand the coefficients of perpendicular transport for particles, momentum and heat with some kind of trial process to match as well as possible a given experiment, they are computed at each time and each position from two fields κ and ε governed by two additional transport equations algebraically derived. κ is the turbulence kinetic energy, and ε is interpreted in this paper as the loss via the energy cascade to the energy dissipation scale. One can then define a diffusion coefficient D = κ 2 /ε and a characteristic time τ = κ/ε.

While inspired from neutral fluid dynamics, we show that both predator-prey models [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF][START_REF] Diamond | Self-regulating shear flow turbulence: A paradigm for the L to H transition[END_REF][START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF] and quasilinear [START_REF] Vedenov | Quasi-linear theory of a plasma[END_REF][START_REF] Drummond | Non-linear stability of plasma oscillations[END_REF][START_REF] Frieman | Kinetic theory of a weakly unstable plasma[END_REF] frameworks used to describe plasma turbulence provide a theoretical background to the κ-ε description of turbulent transport.

The local evolution of each field is governed by a linear term driving exponential transients and a quadratic saturation term. Non vanishing fields are obtained with an interchange like growth rate γ I . The dissipation contribution in the κ equation is of the form κ 2 /D ω + ε. Consistently with the neutral fluid approach we have used the asymptotic limit D ω 1 in the simulations so that only ε acts as dissipation of the kinetic energy of the velocity fluctuations. The quadratic dissipation term in the ε equation, of the form βε 2 , is defined with β = V /κ 3/2 , where V is a characteristic velocity. The choice made for the velocity V is governed by the dimensionless scaling of the width of the plasma boundary layer, the SOL width λ SOL that is supposedly the background of well known empirical scaling laws of the SOL width [START_REF] Eich | Inter-ELM power decay length for JET and ASDEX Upgrade: Measurement and comparison with heuristic drift-based model[END_REF][START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF]. We have then obtained V /c s ∝ ρ * √ γ I τ , c s being the sound velocity, ρ * the ratio of a characteristic Larmor radius and plasma minor radius a, and τ ∝ qAa/c s the parallel connection time, see (22d).

Results of the 1D model are confronted to experimental evidence by analyzing the computed SOL width and comparing the results to the existing scaling law for L-mode plasmas [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF]. Since the chosen closure for V accounts for the latter, one could expect close agreement between the empirical scaling law and the 1D model. It is not the case. This has led us to analyze the scaling law and several shortcomings have been identified, the most striking being that its dimensionless form departs significantly from the reference qρ * dependence and that in fact it is not even dimensionaly correct. One also finds that the dependence on the safety factor appears to have a too big exponent leading to a too strong dependence on the plasma current. Since the transport model also governs the energy confinement time, we have also compared the 1D simulations with scaling law of the energy confinement time, addressing both the dependence in terms of the engineering parameters and the dimensionless aspects. Modification of the scaling exponents of ρ * , A and q in the expression of V /c s have been calculated to match the L-mode and H-mode confinement time dependencies. Finally the 1D model has been used to illustrate one of the possible problems that can be addressed in the κ-ε framework. Introducing in the definition of the velocity V a dependence on the shear of large scale flows, typically the zonal flows [START_REF] Diamond | Zonal flows in plasma-a review[END_REF], so that V exhibits a monotonic decrease from the low shear limit to the large shear limit, the ratio between these values ranging from ≈ √ 2 to ≈ 2, we have reported in 1D simulations an H-mode like transition when increasing the input power, generating an increased stored energy thanks to an interface barrier located at the separatrix [START_REF] Norscini | Self-organized turbulent transport in fusion plasmas[END_REF]. This result is not quite a surprise since the predator-prey models have been introduced in fusion plasmas to investigate such barrier formation [START_REF] Diamond | Self-regulating shear flow turbulence: A paradigm for the L to H transition[END_REF]103,[START_REF] Miki | Spatio-temporal evolution of the L → I → H transition[END_REF][START_REF] Floriani | Self-regulation of turbulence bursts and transport barriers[END_REF]. Furthermore, the model has been devised to exhibit the appropriate trend, namely a decrease of V with flow shear. However, the bifurcation aspect is not introduced specifically in the model. The transition thus appears via a monotonic dependence of V on the shear completed by a feedback mechanism which governs a further increase of the flow shear as transport is reduced. This shows that the κ-ε model has the capability to generate changes of the transport regimes in a selfconsistent manner.

Stepping to the 2D plasma-wall simulations in L-mode, we find that the κ-ε model yields a good match with the experimental profiles, both at the divertor and at the midplane. The transport is shown to be ballooned, as expected for an interchange like instability, driving turbulent transport in the divertor SOL and nearly no transport in the private flux region. The SOL width of WEST is recovered far more precisely than when using the L-mode scaling for the SOL width [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF]. It appears in fact that twice the H-mode scaling [START_REF] Eich | Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER[END_REF] does a better job in predicting the SOL width. These first simulations with the κ-ε model demonstrate the ability of the model to predict equilibrium profiles at midplane and at the outer divertor target as well as providing a 2D distribution of turbulent transport with radial and poloidal dependencies determined self-consistently, in particular the ballooning aspect reported in experiments [START_REF] Labombard | Transport-driven scrape-off-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF][START_REF] Gunn | Evidence for a poloidally localized enhancement of radial transport in the scrape-off layer of the Tore Supra tokamak[END_REF][START_REF] Dif-Pradalier | The Mistral base case to validate kinetic and fluid turbulence transport codes of the edge and SOL plasmas[END_REF].

We have thus shown the remarkable capability of the κ-ε model to capture key aspects of the physics of turbulent transport throughout the plasma. The present simulations, which are first tests of the model, use a single scalar as tuning parameter to describe the whole 2D dependence of the turbulent diffusion coefficient. The confrontation to experimental data is quite convincing, both for midplane and divertor profiles. The width of the energy exhaust channel is also recovered in an L-mode simulation of WEST.

Further adjustment of the model can be achieved using experiments and/or microturbulence simulations. The model can be used to model the core plasma as well as self consistent transitions to improved confinement regimes. Indeed, the non-linear dependence of the diffusion coefficients could allow one extending the use of such transport modeling to transients. This modeling effort would no longer be restricted to steady state as the usual diffusive ansatz would imply. Similarly, because the diffusive coefficients governing the transport of the two fields κ and ε are defined to depend non-linearly on their values, typically like κ 2 /ε, the transport dynamics can depart significantly from that stemming from fixed diffusion processes. Finally, in the present implementation of the model a single instability drives the turbulent fluctuations on the ion characteristic length scale but the model could be completed using other instability growth rates, and one could also split the model according to fluctuations on either ion or electron scales. including the polarization drift current for the ions.

∂ t Ω - 1 r ∇ r rν∇ r Ω = H(r -a) 1 τ eφ T e -Λ (A.1b)
The constraint governed by current loss at the sheath is modified to take the linear form of a restoring force of the electric potential φ towards ΛT e /e. The steady state solution in the core plasma thus governs φ ≈ -T i /e while in the SOL one has φ ∝ T e /e, provided |∇T e | |∇p i |/n. At the separatrix the viscosity ν will then bridge these two asymptotic behaviors. In the standard use of this 1D model, as well as for the SolEdge2D simulation addressed in the paper, large scale flows are not taken into account in the transport. On the same footing the non-linear drive of the vorticity via the Reynolds's stress, mandatory for micro-turbulence simulations is not included in (A.1b). Solving the vorticity equation (A.1b) is only useful when considering an ad hoc modification of the control parameter V governed by the E × B shear. A simplified transport model for neutrals is also introduced so that the particle source by ionization S n , and in particular its location, changes consistently with the plasma parameters. Only the particle balance equation is used for the neutrals with density N 0 , diffusive transport with constant diffusion coefficient D N 0 and a feed-back controlled gas injection Φ N 0 . In this 1D model, we address an ionizing plasma regime such that the ionization and charge exchange cross sections are comparable. The diffusive nature of neutral transport is not particularly well suited since ionization occurs in a couple of steps, charge exchange being subdominant.

∂ t N 0 - 1 r ∇ r rD N 0 ∇ r N 0 = Φ N 0 H(r -a) -S n (A.1c)
The source term Φ N 0 should be localized at the outer wall as well as peaked towards r = a to describe recycling via the private flux region of a divertor or in the vicinity of a limiter tip. We simplify these aspects by taking Φ N 0 constant throughout the SOL region, hence of the form Φ N 0 H(r -a). Finally, the cross-field diffusion terms of the plasma transport equations are determined by the κ-ε coupled equations.

     ∂ t κ - 1 r ∇ r rD∇ r κ = γ κ κ -ζκ 2 -ε ∂ t ε - 1 r ∇ r rD∇ r ε = γ ε ε - V κ 3/2 ε 2 (A.1d)
It is to be noted that no parallel convective loss term is retained here to account for the SOL transport to the wall. These would reduce the drive governed by the local dynamics, for instance substituting in the κ equation γ κ by γ κ -σ κ /τ where σ κ is a constant of order unity. However, a similar parallel loss term occurs on the closed magnetic surfaces for parallel currents, which have a stabilizing effect on the interchange drive. Consequently, a similar correction should be applied to γ κ in the region r < a. The control parameters that appear in the 1D model (A.1d) should therefore be understood as effective parameters, as discussed for plasma transport (A.1a) and (A.1b). Although, the control parameters in (A.1d) should be modified to account for the parallel transport, in practice we have assumed that this effect is small for both fields κ and ε.

Appendix A.2. The sources

The particle source distribution in space depends mostly on the neutral penetration into the plasma. The crude model for neutral transport is tuned to obtain realistic profiles of S n = nN 0 σv i , where σv i is the ionization rate. The latter is computed using an expression of the form σv i ∝ x 1/2 /(x i + x) exp(-1/x) [START_REF] Richardson | NRL PLASMA FORMULARY[END_REF], where x = T e /E i and E i = 13.6 eV is the ionization energy. With x i = 6 the maximum of the ionization rate occurs at x = 10. Overall, this analytical expression yields a dependence of the ionization rate with the electron temperature that is roughly comparable to published data [START_REF] Janev | Survey of atomic processes in edge plasmas[END_REF][START_REF] Post | A review of recent developments in atomic processesfor divertors and edge plasmas[END_REF]. The neutral source Φ N 0 is set to ensure that the density at the separatrix is maintained equal to 1.10 19 m -3 . This can be understood as a particle injection rate with feedback on the plasma separatrix density. A typical radial profile for the neutral density in this model is displayed on Figure A1 left hand side for the following reference plasma parameters typical of a WEST experiment: R = 2.5 m, a = 0.5 m, B T = 3.7 T , B P = 0.2 T , P in = 1 M W . One finds that the neutral density decreases from the wall into the plasma. The typical e-folding length is approximately 0.05 m, hence a tenth of the chosen minor radius a. This decay rate is a combined effect of the diffusive transport and of the ionization sink, the latter prevailing in the confined plasma region. The This generates a constraint, the so-called Kadomtsev similarity constraint [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF].

ν * = n a T 2 qA A sν ν * 0 ; ν * 0 = e 2 /(4π 0 ) 2 n 0 a 0 T 2 0 (B.2a) β = nT B 2 β 0 ; β 0 = n 0 T 0 B 2 0 /(2µ 0 ) (B.2b) ρ * = T 1/2 aB ρ * 0 ; ρ * 0 = m p T 0 eB 0 a 0 (B.2c) q = F q (k, δ)q cyl ; q cyl = aB T RB pol ; B pol = µ 0 2π I p akF I (k, δ) (B.2d)
In these expressions the quantities with 0 subscript stand for the parameter values used for normalization, most often done implicitly by specifying the units. When ν * is governed by the physics of the trapped particles s ν = 3/2, otherwise, as more readily done for the edge region s ν = 0. The parameter ρ * is defined with the proton mass m p , however, regarding the physics of cross-field transport, the effective ion Larmor radius is of interest. One must therefore take this effect into account when addressing the dependence on the effective ion mass M . The safety factor is split into the cylindrical contribution q cyl and two dimensionless functions of the magnetic equilibrium F q and F I , that depend typically on dimensionless parameters such as the ellipticity k, the triangularity δ, etc. Given these 3 dimensionless parameters one can express the density n, the product aB = RB T /A, an the square root of the thermal energy T

1/2 . n = T 2 a B ν * /A sν qA B 1 ν * 0 (B.3a) a B = T 1/2 1 ρ * ρ * 0 (B.3b) T 1/2 = B 1/5 qAβ ρ * ν * /A sν 1/5 ρ * 0 ν * 0 β 0 1/5 (B.3c)
Given (B.3c), one thus obtains:

n = n 0 B 8/5 ρ * ν * /A sν qA 2/5 β 3/5 (B.4a) a B = (aB) 0 B 1/5 qAβ ρ 6 * ν * /A sν 1/5 (B.4b) T 1/2 = T 1/2 0 B 1/5 qAβ ρ * ν * /A sν 1/5 (B.4c)
The various constants that appear in (B.4) lead to the definitions of n 0 , (aB) 0 and T 0 . In a similar way one can express the plasma current from the definition (B.2d) of the safety factor:

I p = 2π µ 0 k A q cyl aB = I 0 k A q cyl qAβ ρ 6 * ν * /A sν 1/5 B 1/5 (B.5a)
The four dimensional parameters n, (aB), T , I p are found to depend on the dimensionless parameters q, ν * , β and ρ * as well as on the dimensional parameter B. Finally, one must express the heating power. This engineering parameter is specific insofar that it will depend on the internal energy and on the confinement time for steady state conditions as follows:

P = n T V τ E = n T B 2 /(2µ 0 ) V Ωτ E B 2 /(2µ 0 )Ω 0 = Akβ Ωτ E a 3 B 3 F V 2π 2 B 2 0 eB 0 2µ 0 m = P 0 Akβ Ωτ E a B 3 = P 0 Akβ Ωτ E qAβ ρ 6 * ν * /A sν 3/5 B 3/5 (B.5b)
The dimensionless function F V (k, δ) is introduced when writing explicitly the plasma volume V = πa(ka)2πRF V . The confinement time is normalized here by Ω, where Ω is a reference Larmor pulsation depending on B T .

Ωτ E ∝ I α I p B α B +1 n αn P -αp M α M a α R A α A +α R k α k (B.6a) Ωτ E ∝ k a k q aq M a M A a A β a β ρ aρ * ν aν * (B.6b)
One can then identify the exponents of the power law dependence of the confinement time on the dimensionless parameters. In particular one obtains a B : a B = (α I + 8α n -3α p + 5α B -4α R + 5) 5(1 -α p ) (B.7a)

With the coefficients of the ITER1996th L-mode scaling, Table B1 one obtains a B ≈ -0.15 instead of zero. The scaling law is not homogeneous. One also recovers the expressions for the other exponents published in Ref. [START_REF] Sarazin | Impact of scaling laws on tokamak reactor dimensioning[END_REF]. The three exponents a β , a ν and a ρ characterize the plasma: and the main ion species effective mass dependence is:

a M = α M (1 -α p ) (B.7h)
In the ITER physics basis [START_REF]ITER Physics Basis Editors ITER Physics Basis Expert Groups on Confinement and Transport and Confinement Modelling and Database. Chapter 2: Plasma confinement and transport[END_REF], a correction on the exponent α R is proposed to recover the constraint a B = 0. In practice any of the coefficients contributing to a B could be corrected as well as any combination of corrections yielding a B = 0. In fact, the scaling law and the error minimization procedure using the data base should be performed with the constraint a B = 0. In the ITER physics basis this point is not discussed and no information is given on the impact of a change of α R on the accuracy of the scaling law.

Another possibility is to modify directly the exponent α B that only contributes to a B . This would imply a change from α B = 0.03 to α B = -0.01. This can appear as a small change, however it might be clearly inconsistent with the data, in particular regarding the change of sign.

Appendix B.2. Adimensional scaling of the empirical SOL width

Let us now consider the scaling law proposed for λ q in Ref. [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF]:

λ q = λ * ,0 B u B T q uq P up R u R (B.8a)
As for the energy confinement time, this scaling law for λ q in millimeters, assumes that B T is given in Tesla, P in in M W , R in meters. In a first step, the expression is rearranged to step towards the dimensionless expression. We then introduce r B = u B + 1 -u R , r q = u q , r p = u p , r a = u R -1, r A = u R and therefore:

λ q a = λ * ,0 (B T ) r B q rq P rp (aB) ra A r A (B.8b)

The coefficients u of the scaling law (B.8a) are given in Table B2.

λ * ,0 u q u B u p u R 1.44 ± 0.67 -0.8 ± 0.32 1.4 ± 0.67 0.22 ± 0.1 -0.03 ± 0.28 Table B2. Exponents of the engineering scaling law of the SOL width proposed in Ref. [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in L-mode plasmas at JET and AUG[END_REF].

λ q a = (1.44 ± 0.67)B r B q rq A r A P rp 0 (Ωτ E ) -rp q 3rp/5 A 8rp/5 A 3sν rp/5 β 8rp/5 ν -3rp/5 * ρ -18rp/5 * B 3rp/5

(a B) ra 0 q ra/5 A ra/5 A sν ra/5 β ra/5 ν -ra/5 * ρ -6ra/5 * B ra/5

(B.9a)

The dependence on Ωτ E can be expressed directly in dimensionless parameters. It seems more appropriate when doing to use the expression that satisfies the Kadomtsev constraint, hence with a B = 0.

λ q a ∝ A r A +8rp/5+ra/5 A sν (3rp+ra)/5 β 8rp/5+ra/5 ν -3rp/5-ra/5 * q rq+3rp/5+ra/5 ρ -18rp/5-6ra/5 * B r B +3rp/5+ra/5

(k a k q aq M a M A a A β a β ρ aρ * ν aν * ) -rp (B.9b)

One can then express the dimensionless scaling of the SOL width as:

λ q a ∝ k k q q M M A A β β ρ ρ * ν ν * B B (B.9c)
where the exponents are given by: Using the coefficients as given in the scaling one finds again that the scaling law is not homogeneous leading ot B = 0.16. The only correction on a single exponent that stands within the error bars of the scaling law is a decrease of the exponent of B from -0.8 to -0.956. In that case one obtains: λ q a ∝ q 2.15 ρ 0.85 * A 0.03 β 0.46 ν -0.17 *

(B.11)

This result departs rather strongly from that assumed to hold qρ * and used to determine the parameter V of the κ-ε model (22d). Note that regarding the dependence on the aspect ratio A, we have used the expression of ν * without the trapped particle effect, hence for s ν = 0. For s ν = 3 2 , the exponent for the aspect ratio scaling is found to be -0.01.

Appendix B.3. Adimensional scaling of the κ-ε SOL width

The parameter scans performed with the 1-D κ-ε model indicate that:

λ q ∝ B v B T B -vq pol P vp R v R (B.12a)
v q v B v p v R 0.844 0.277 0.12 -1.06 Table B3. Exponents of the engineering scaling law obtained with the 1-D κ-ε model.

The coefficients v of the scaling law obtained are given in Table B3. With B pol = B T /(Aq) one can rewrite this dependence as:

λ q a ∝ q vq B v B -vq+1-v R T A vq+v R P vp (aB) v R -1 (B.12b)
We the introduce r B = v B -v q + 1 -v R , r q = v q , r p = v p , r a = v R -1, r A = v q + v R and therefore:

λ q a ∝ k k q q M M A A β β ρ ρ * ν ν * B B (B. 13
)
where the exponents are given by (B.10). Using the coefficients as given in Table B3, one finds again that the scaling law is not homogeneous leading to B = 1.15. However, the variation of the density during the various scans, governed by the chosen feedback scheme, and that does not exhibit a power law behavior, can explain this issue.

λ q a ∝ q 0.95 ρ 2.26 * A -0.48 β -0.05 ν -0.05 * B 1.15

T (B.14)

One finds therefore that the result exhibits a dependence on B T ρ * , hence typically inverse to the size a, that appears to be spurious. Such a size effect can be related to a relative change of the particle source and can have therefore have an impact on the turbulent transport. One can also observe that the κ-ε 1D model appears to exhibit a too strong dependence on R, see Table B4, which could be the drive for the difference between the ρ * 2.26 of the 1D model SOL width scaling and the input ∝ ρ * used to determine V . This difference in scaling law can also be related to the particle source localization and consequently the particular feedback scheme used in the present simulations. For a complete comparison further 1D simulations using a proper database of relevant experiments are required, with the same control scheme of the particle source as in the experiments. Furthermore, uncertainty propagation will be needed to identify the parameters that have a strong impact on the observed difference in the ρ * dependence, and then address the physics and measurements that are crucial. Finally, this work could also lead us to modify the input assumption used to determine V and its ρ * dependence to obtain a simulation output matching the experimental database. Assuming that the SOL width λ q /a is of the form: One can also use q = aB T /(AI p ) to obtain:

λ q a ∝ q q ρ ρ * A A (B.
λ q a ∝ aB T AI p q P τ E 1/2 na 3 Ak 1/2 B T a ρ A A (B.19)
Given the expression of τ E (B.1), one then has:

λ q a ∝ I α I p B (α B -2) T
n (αn-1) P (1-αp) M α M a (α R -5) A (α A +α R -1) k (α K -1) ρ/2 aB T AI p q A A (B.20a) which then yields the scaling law of the SOL width given in terms of engineering parameters:

λ q ∝ I α I ρ/2-q p B (α B -2) ρ/2+ q T P (1-αp) ρ/2 n (αn-1) ρ/2 a 1+(α R -5) ρ/2+ q A A -q +(α A +α R -1) ρ/2 M α M ρ/2 k (α K -1) ρ/2 (B.20b)

Let us consider the case ρ = q = 1, A = 0 that corresponds to the qρ * SOL width scaling:

λ q ∝ I α I /2-1 p B (α B -2)/2+1 T P (1-αp)/2 n (αn-1)/2 a 1+(α R -5)/2+1

A -1+(α A +α R -1)/2 M α M /2 k (α K -1)/2 (B.21a)

For the L-mode confinement scaling, Table B1, this leads to: λ q ∝ B -0.52 pol B 0.015 T P 0.135 R -0.555 n -0.3 a 0.45 M 0.1 k -0.18 (B.21c)

λ q ∝ I -0.
These results can be compared to those obtained with the scans performed with the 1-D κ-ε model, see Table B4. The agreement is rather poor, but at least there is no disagreement regarding the sign of the exponents. As discussed previously, the feedback process, leading to a change of the core density at fixed separatrix density, is an issue since the density n, assumed to be an independent engineering parameter changes during each of the scans that have been performed. exponents v q v B v p v R κ-model 0.844 0.277 0.12 -1.06 expected scaling 0.52 0.015 0.135 -0.555 Table B4. Exponents of the engineering scaling law obtained with the 1-D κ-ε model compared to its expected scaling.

Appendix C. The SolEdge2D model

The model [START_REF] Bufferand | Numerical modelling for divertor design of the WEST device with a focus on plasma-wall interactions[END_REF] is typically a system of Braginskii drift-reduced fluid equations [START_REF] Braginskii | Transport processes in a plasma[END_REF] that govern the evolution of the plasma density n, the parallel momentum and the total energy temperature for both electrons and ions assuming quasi-neutrality n e = n i and ambipolarity v ,e = v ,i . The latter constraint can be relaxed but the charge balance equation must then be solved to determine the electric field and consequently take into account the large scale E × B convective transport. In the simplified case with ambipolarity, the equations write as follows: . The latter is used in SolEdge2D as immersed boundary condition in the penalization technique. The parallel heat flux impinging on the solid target q ,BC will read q ,BC = (γnT v ) BC , where γ is the total sheath transmission coefficient.

∂n ∂t + ∇ • nv b + nv ⊥ = ∇ • D n ∇ ⊥ n + S n (C.1a) ∂(nv ) ∂t + ∇ • nv v b + v ⊥ = -∇ nT i m i + q i nE m i + R ei m i + ∇ • (νn∇ ⊥ v ) + S nv (C.1b) ∂ ∂t 3 2 nT i + 1 2 m i nv 2 + ∇ • 5 2 nT i + 1 2 m i nv 2 v b + v ⊥ = ∇ • κ i ∇ T i b + χ i n∇ ⊥ T i + νn∇ ⊥ 1 2 m i v 2 + q i nv E + R ei v + Q ei + S Ei (C.

Figure 1 .

 1 Figure 1. Variation of κ (left hand side) and ε ( right hand side) with the control parameter γ, for different values of the exponent η: η = 1 black line head down closed black triangles, η = 3/2 black line closed black circles, η = 1.25 blue line closed head up blue triangles, η = 2 blue line closed blue triangles and η = 0.5 dashed blue line open blue circles. Reference values η = 1, RANS κ -model [12] and η = 3/2 Kolmogorov theory [54].

Figure 2 .

 2 Figure 2. Variation of κ blue line closed circles left hand side axis, and ε black line open circles, right hand side axis, with the control parameter γ κ τ ql .

Figure 3 .

 3 Figure 3. Sketch of the 2D transport pattern in TCV in the large X-point height configuration. Pattern of the particle flux in the poloidal plane, simulations with SolEdge2D-EIRENE but without self-consistent turbulent transport, sketch from Ref.[75]. Panel a: standard transport modeling in 2D plasma-wall interaction. Panel b: modified choice of transport modeling enforcing ballooned transport. Strong outward particle flux is labeled "High > 0", moderate outward particle flux is "Medium > 0", and moderate reversed particle flux is "Medium < 0".

Figure 4 .

 4 Figure 4. Sketch of the 2D transport pattern in TCV in the large X-point height configuration. Panel a: effective particle diffusivity given by TOKAM3X, sketch of the results published in Ref.[75]. Panel b: Particle turbulent diffusion coefficient D n selfconsistently determined in SolEdge2D-EIRENE simulation using the κ-only model, sketch of the results published in Ref.[START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. Regions with largest particle diffusion coefficient are labeled "high", those with still significant particle diffusion are labeled with "medium".

Figure 5 .

 5 Figure 5. Radial profiles of plasma quantities when increasing the input power with feedback control on the separatrix density. (a) density. (b) Electron temperature. (c) Effective diffusion coefficient estimated from κ and ε. (d) Electric potential.

Figure 6 .

 6 Figure 6. Left hand side: Dependence of the core density n 0 and of the core pressure p 0 on the input power P in . Right hand side: scaling of the core T 0 and separatrix T a thermal energy in terms of the input power. The values in bracket are the scaling exponents.

Figure 7 .

 7 Figure 7. At the separatrix, power law dependence on the input power P in , Left hand side: particle diffusion coefficient D n and density gradient length L n = (∇ r n/n) -1 , Right hand side: neutral particle flux Φ n .

Figure 8 .

 8 Figure 8. Power law dependence of the plasma pressure e-folding length, used as a proxy for the SOL width λ SOL black open squares, approximate power law following these results dashed blue line, and power law governed by the experimental scaling law, black dashed-dot line. Left hand side scaling with P in ∈ [0.1 M W, 6 M W ]. Right hand side scaling with B pol ∈ [0.1 T, 0.25 T ]. The gray-shaded area stands for the 95% confidence interval of the regression used to determine λ p .

Figure 9 .

 9 Figure 9. Power law dependence of the plasma pressure e-folding length, assumed as a proxy for the SOL width λ SOL black open squares, approximate power law following these results dashed blue line, and power law governed by the experimental scaling law, black dashed-dot line. Left hand side scaling with R ∈ [1 m, 4 m]. Right hand side scaling with B T ∈ [2 T, 3 T ]. The gray-shaded area stands for the 95% confidence interval of the regression used to determine λ p .

  [s], for I p , the plasma current in [M A], B T the toroidal magnetic field in [T ], P the loss power in [M W ], n the line-averaged density in [10 19 m -3 ], M the average ion mass in [AM U ], R the major radius in [m] and finally A and k the aspect ratio and elongation, respectively.

Figure 10 .

 10 Figure 10. Comparison of the power law dependence of the thermal energy confinement time τ E on the injected power P in , Left hand side, and on the plasma current I p given B pol ∝ I p , Right hand side. The 1D simulation data, open black squares is fitted with the power law plotted with blue short-dash lines. The ITER L-mode scaling law (33a) yields the black dashed line. The κ-ε scaling stemming from the empirical scaling law of the SOL width (32b) is plotted with dash-dot black lines while the theory based κ-ε models (34a) are plotted with blue long-dash lines for q = 2 and blue dotted line for q = 1. Note that for the injected power scaling law, these two cases yield the same behavior.

Figure 11 .

 11 Figure 11. Stabilization by E × B shear. Left hand side: Radial profiles of the ion thermal energy T i when ranging the input power P in from 0.5 M W to 8 M W , closed square P in = 0.5 M W , close circle P in = 1 M W , closed head down triangle P in = 2 M W , open head up triangle P in = 2.5 M W , ope circle P in = 4 M W and open square P in = 8 M W . Right hand side: variation of plasma diffusivity D n with input power P in at different radial positions. Within the separatrix: open square ρ = 0.96; open circle ρ = 0.98, head up open triangle ρ = 0.999, closed head down triangle open square ρ = 0.998, and outside the separatrix: closed circle open square ρ = 1.0, open square ρ = 1.2.

Table 4 .

 4 Control parameters for the SolEdge2D-Eirene simulation of the shot W EST -55049, density at the inner boundary n OM P e , heating power into the simulation domain P in , poloidal field at the outer midplane B OM P pol . n OM P e [m -3 ] P in [M W ] B OM P pol governed by parallel advection from the low field towards the high field side. The SolEdge2D-Eirene radial profiles for the density n e and the electron thermal energy

Figure 13 .

 13 Figure 13. SolEdge2D-Eirene simulation of the WEST shot 55049. Left hand side: variation in the poloidal plane of κ normalized by its maximum value with logarithmic scale. Right-hand side: poloidal plane variation of D, units m 2 s -1 , with logarithmic scale.

Figure 14 .

 14 Figure 14. Divertor profiles and radial profiles at the outer midplane of the electron thermal energy, Left hand side, and density, Right hand side, obtained with the SolEdge2D-Eirene simulation of the shot WEST-55049. The divertor profiles are remapped at the midplane. Midplane profile, blue curve closed circles, divertor profile black curve open circles.

Figure 15 .

 15 Figure 15. Radial profiles at the outer midplane for the shot WEST-55049. SolEdge2D-Eirene (black curve closed circles), probe data (blue dashed line and open triangles). The experimental data measured with the reciprocating Langmuir probes (RLP) located at the top of the device with vertical plunge are mapped at the midplane.The radial extent is restricted to the region where the experimental data is reliable, in the far SOL R -R sep ≥ 0.05 m, the density is small and the interpretation of the probe characteristics is difficult, when approaching the separatrix R -R sep ≤ 0.01 m the probe head interacting with a hot plasma starts behaving as a secondary limiter. Left hand side: electron thermal energy. Right hand side: density.

  the private flux. The e-folding lengths in the SOL are λ W EST SOL = 9.1 mm and λ SE2D SOL = 7.3 mm, with a difference of order 20 %, and the decay length towards the private flux are λ W EST pf = 1.1 mm, smaller than estimated from the simulation λ SE2D pf = 1.5 mm.

Figure 16 .

 16 Figure 16. Experimental and simulation divertor profiles remapped at the outer midplane for the shot WEST-55049. Experimental data from Langmuir probes, labeled WEST, embedded in the divertor floor, closed blue circles, SolEdge2D-Eirene plain black line, labeled SE2D. Left and right logarithmic scales have the same ratio between max and min, but are shifted to ensure overlap of the peak values. Left hand side: ion saturation current. Right hand side: parallel heat flux.

Figure 17 .

 17 Figure 17. Experimental and simulation divertor profiles remapped at the outer midplane for the shot WEST-55049. Experimental data from Langmuir probes, labeled WEST, embedded in the divertor floor, closed blue circles, SolEdge2D-Eirene plain black line, labeled SE2D. Linear right and left scales are adjusted to ensure the overlap of the maxima. Left hand side: plasma density. Right hand side: electron thermal energy.

Figure A1 .

 A1 Figure A1. For WEST like parameters P in = 1M W , R = 2.5m, a = 0.5m, B T = 3.7T , Left hand side axis: Radial profile of the neutral density N 0 blue line open circles, Right hand side axis: radial profile of the ionization source term S n black line closed circles.

  a β = (α I + 3α n -8α p + α R ) 5(1 -α p ) (B.7b) a ν = (-α I + 2α n + 3α p -α R ) 5(1 -α p ) (B.7c) a ρ = (-6α I + 2α n + 18α p -6α R ) 5(1 -α p ) (B.7d)while the magnetic equilibrium properties are characterized by:a k = (α I + α K -α p ) 1 -α p (B.7e) a q = (α R -4α I -2α n -3α p ) 5(1 -α p ) (B.7f) a A = (5α A + 6α R -4α I -2α n -8α p ) 5(1 -α p ) -s ν a ν (B.7g)

+ s ν 3r p + r a 5 -

 5 B = 5r B + 3r p + r a 5 -r p a B (B.10a) β = 8r p + r a 5 -r p a β (B.10b) ν = -3r p -r a 5 -r p a ν (B.10c) ρ = -18r p -6r a 5 -r p a ρ (B.10d) k = -r p a k (B.10e) q =5r q + 3r p + r a 5 -r p a q (B.10f)A = 5r A + 8r p + r a5r p a A (B.10g) M = -r p a M (B.10h)

Appendix B. 4 .

 4 Dimensional scaling of the SOL width qρ *

52 pB 0. 015 T

 52015 P 0.135 n -0.3 a 0.415 A -0.555 M 0.1 k -0.18 (B.21b)Since B pol ∝ I p /a and A = R/a one can rewrite this result to obtain an expression comparable to (B.12a).

  b + v ⊥ = ∇ • κ e ∇ T e b + χ e n∇ ⊥ T e -env E -R ei v -Q ei + S Ee (C.1d)where ∇ = b • ∇ and ∇ ⊥ = ∇ -b • ∇ define gradients along the parallel and perpendicular direction, b = B/B being the unit vector along the magnetic field defining the parallel direction. v ⊥ = v -v b is the velocity component perpendicular to the magnetic flux surfaces due to the large scale flows when determined (for the present simulations v ⊥ = 0), E is the parallel component of the electric field, R ei is the parallel electron-ion friction force, S n , S nv , S Ei and S Ee are the particle, momentum and energy sources respectively.The boundary conditions are derived from the physics of the sheath: at the sheath entrance, the absolute value of plasma parallel velocity along the magnetic field lines reaches minimum sound speed. This is the Bohm condition for the plasma parallel velocity, expressed by |v | ≥ c s = e(Te+T i ) m i

  One also notices that the profiles shapes do not compare well. Differences in the density n e profiles, Figure

	8	7.3	7.3	9.1	| 4.1 8.2 2.9	18	0.3
	reasonably close.

  [START_REF] Schneider | Plasma edge physics with B2-Eirene[END_REF] one can step backwards to determine the expected scaling law in terms of the engineering parameters. Given ρ *

		ρ * ∝	T 1/2 B T a	(B.16)
	and:
		T ∝	P τ E nV	(B.17)
	so that:
	ρ * ∝	P τ E na 3 Ak	1/2 1/2 B T a	(B.18)

Acknowledgements

The authors thanks Prof. B.E. Launder for fruitful discussions on the original kepsilon model which inspired this work. This work has also been inspired by numerous discussions with Prof. P. Diamond, G. Dif Pradalier, X. Garbet and Y. Sarazin. The authors are also indebted to the Festival de Théorie, a meeting that enables in depth scientific exchange. The authors are very thankful to the referees who have contributed quite significantly to improving the paper. This work has been supported by the French National Research Agency grant SISTEM (ANR-19-CE46-0005-03). The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University-A*MIDEX, a French "Investissements d'Avenir programme". This work was granted access to the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01). This work has received funding from the Euratom research and training program 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A. The 1D transport model

Appendix A. [START_REF]ITER Physics Expert Group on Divertor Modelling and Database and ITER Physics Basis Editors. Chapter 4: Power and particle control[END_REF]

. The transport equations

The 1D model is obtained by averaging the transport equations Appendix C in the poloidal and toroidal directions. Transport is taken into account by diffusion that replace all convective transport processes at the microscopic level. The diffusion coefficients are chosen as discussed in Section 2.1 equation [START_REF] Rognlien | A fully implicit, time dependent 2-D fluid code for modeling tokamak edge plasmas[END_REF]. For the closed magnetic surfaces, this average must be considered as a flux surface average while for the open field lines a standard averaging on the parallel direction is used so that the parallel divergence terms yield contributions that are identified as the outflux onto the wall components. Since the parallel dynamics are ignored, the transport equations are then reduced to the density n and thermal energy evolution E i = 3 2 nT i for the ion and E e = 3 2 nT e for the electron. This simplified model assumes quasineutrality and a low Mach number regime. The Mach number of the ion flow u i /c s -where c s is the sound velocity-is thus assumed to be small on average, hence

and is only taken into account via the boundary condition where it is assumed to be finite and typically of order one. While the Bohm condition enforces |M | ≥ 1 with respect to the local boundary value of c s , it is not the case on average m i u 2 i T i . The brackets of the average are dropped in the following to simplify the notations. The plasma transport equations solved in the model are then:

In these equations, the function H(r -a) is the Heaviside step function used as a mask that defines the SOL region, H(r -a) = 0 for r < a and H(r -a) = 1 for r ≥ a so that the parallel loss terms apply. For these averaged equations, the convective loss term appear to be governed by the effective SOL confinement time for the particles, τ , typically of order L /c s . The coefficients γ e and γ i are reminiscent of that computed for the kinetic sheath transmission, but also take into account effects governed by the relationship between the sheath values of the density and thermal energy and their parallel average. The same remark holds for the transport terms since the average of the product of the local values of χ, n T and ∇T is not equal to the product of the average of these fields. The transport coefficients must therefore be considered as effective. For the heat conductivity, the transport coefficient χ eff to be used is such that χ eff < n >< ∇T >=< χn∇T >. This issue is not specific of the present model. The plasma equations are completed by a vorticity equation for Ω

where the electrostatic potential is φ, e the electron charge, and p i the ion pressure. This equation is derived from the charge balance equation For a typical plasma volume V = 2πRπk a 2 , where R is the major radius, a the minor radius and k the plasma elongation, the chosen heat source profile S E , Figure A2, is defined as a function of the normalized radius ρ = r/a:

where the constant p 0 = 7.425 10 -2 ensures that the cylindrical integral of VS E is P in , the injected power.

Appendix B. Adimensional scaling laws

Appendix B.1. Adimensional scaling of the confinement time

The generic form of the scaling laws that are considered for the energy confinement time is given in terms of the so-called engineering parameters as follows:

where τ E is obtained in [s], for I p , the plasma current in [M A], B T the toroidal magnetic field in [T ], P the loss power in [M W ], n the line-averaged density in [10 19 m -3 ], M the average ion mass in [AM U ], R the major radius in [m] and finally A and k the aspect ratio and elongation, respectively. It must be underlined that in fact this kind of expression is a mix including dimensional parameters I p , B T , n, P and R and dimensionless parameters M , k and A. Some parameters govern the magnetic equilibrium including its geometry, I p , B T , R, A, k, while others are less directly controlled such as the density n and the average mass ratio M that are really governed by particle confinement and sources, and the heating power coupled to the plasma P .

One can understand such a scaling law either as a fit that holds over a limited range of the parameter values, hence a local tangential fit, or as an effective scaling holding for all values of the parameters. Implicitly, the latter understanding of the underlying similarity appears to be favored since the former tangential form would generically introduce offsets, typically (1 + G/G 0 ) α for a given parameter G. For G G 0 this does not make any difference but the asymptotic behavior for G → 0 for α > 0, respectively G → ∞ for α < 0 is different, and, more importantly has different implications in terms of the underlying physics. Regarding the latter, a more relevant formulation is that given in terms of dimensionless parameters. For fusion plasmas these are the safety factor q, ν * the normalized collisionality, β the normalized plasma pressure and ρ * the normalized ion Larmor radius. If only four dimensionless control parameters exist, then the normalized confinement time Ωτ E must depend on these parameters, hence the 5 dimensional parameters must be expressed as functions of 4 dimensionless parameters.