S Baschetti

H Bufferand

G Ciraolo

Ph Ghendrih

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

P Tamain

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

Keywords: prey models, 1, 3 Quasilinear transport models

The 2D mean-field plasma edge transport description of plasma wall interaction is completed by a κ-ε model as in Reynolds Average Navier Stokes simulations for neutral fluids. The local evolution of the turbulent kinetic energy κ and its dissipation rate ε are revisited and slightly modified. It is shown that the κ-ε extends the quasilinear approach by self-consistently determining κ and the relevant time κ/ε leading to the diffusion coefficient κ 2 /ε. The κ-ε evolution is also shown to be equivalent to both coupled Ginzburg-Landau amplitude equations and predator-prey systems where κ is the prey and ε the predator. The dissipation process ε we enforce describes the small scale dissipation of Kolmogorov cascades. It depends on a free parameter akin to a velocity V . The chosen closure relates V to the the parallel connection time and the normalized Scrape-Off layer width qρ * , q is the safety factor and ρ * the ratio of the characteristic Larmor radius and plasma minor radius. A 1D model with κ-ε self-organized transport is used for comparison to empirical scaling laws of SOL width and energy confinement time in L-mode plasmas. Shortfalls of the scaling laws are analyzed. Possible changes of the closure for V are discussed. The 1D model is also used to test the transport response to a dependence of V to large scale velocity shear. Spontaneous confinement improvement when increasing the heating power 1 is observed with the development of an interface barrier at the separatrix. Plasma-wall interaction simulations for TCV and WEST are analyzed. A single scalar free parameter tunes the cross field transport. Experimental midplane and divertor profiles are compared to the simulations. Remarkable agreement is observed. The SOL width determined by the simulation for WEST is close to the experimental value with less than 20 % difference, while the scaling law for L-mode is off my more than a factor 3. The turbulent transport described by the κ-ε is not homogeneous, but ballooned as experimentally observed together with cross-field transport in the divertor SOL on the low field side, and nearly none in the private flux region. Contents 1 κ-ε and predator-prey models of turbulent transport 1.1 Bridging turbulence and transport models . . . . . . . . . . . 1.2 Predator-

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. ITER divertor operation in high performance scenarios is already expected to face such an issue [START_REF] Tamain | on Divertor. 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 chain 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, full-f 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 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 [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]. In the chain of models, between that required for optimization and first principle turbulence simulations, there remains a need for a model suitable for designing experiments, from adjusting plasma start-up, edge plasma conditions, and burn window, with short return times simulations, and generating information suitable for comparison to experiments as achieved with transport codes. However, these require stepping towards predictive transport modeling, freezing on a physics basis all the free parameters that account for cross-field plasma transport [START_REF] Baelmans | Efficient parameter estimation in 2d transport models based on an adjoint formalism[END_REF].

Despite the constant growth of computational power, engineering simulations for routine use in ITER size machines and in ITER relevant parameters will rely on the so-called transport codes providing only mean flow solutions. This is similar to Reynolds Averaged Navier-Stokes (RANS) models 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. Such a reduction on the degree of freedoms allows one taking into account additional equations describing surface physics processes and atomic physics in realistic tokamak geometries. Of particular interest are the various ionization stages and breakdown processes of D 2 molecules that govern the particle fueling process and the ionization-recombination stages of various impurities, both intrinsic and extrinsic, that must be taken into account. 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 montecarlo model of neutral-particle transport in diverted plasmas[END_REF] for neutral particles), SONIC [START_REF] Kawashima | Simulation of radiative divertor plasmas by ar seeding with the full wwall in jt-60sa[END_REF], EDGE2D [START_REF] Simonini | Models and numerics in the multi-fluid 2-d edge plasma code edge2d/u[END_REF] (coupled to the NIMBUS [START_REF] Cupini | Nimbus-monte carlo simulation of neutral particle transport in fusion devices[END_REF] for neutral particles), 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, developed in the team for more than 10 years, 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. The code SolEdge2D-EIRENE is operated with configurations of most existing tokamaks in Europe.

A key challenge for this class of codes is the description of cross field transport induced by the drifts. The small scale fluctuations of the latter are taken into account by ad'hoc diffusive processes, while large scale flows are found to also play a role and require solving the vorcity equation derived from charge balance. The predictive capability of such codes, as needed for ITER operation, depends to a large extent on the description of the small scale drift velocity fluctuations and their impact on cross-field transport. The effective particle diffusion 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. In present transport simulations, these coefficients 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. It is 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 scrapeoff-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. Another issue of interest is the ability of such transport codes to model transients. For most problems a steady state solution is computed and effective diffusion can model properly such steady states. However, 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] Fedorczak | Flow generation and intermittent transport in the scrapeoff-layer of the tore supra tokamak[END_REF], one expects the reduction to diffusive transport to become inappropriate when addressing time dependent processes. Finally, 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 targets plates, and consequently the operation constraints governed by the exhaust limitations. As a reference approximation of this SOL width, called simple SOL, we consider the standard e-folding lengths for cross field particle transport, typically λ 2 n = D n τ with τ = L /c s , L being the characteristic length of the open field lines and c s the plasma sound velocity.

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 neutral fluid turbulence [START_REF] Launder | The numerical computation of turbulent flows[END_REF] 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 coefficients mentioned above. A similar approach was recently proposed by [START_REF] Coosemans | A new mean-field plasma edge transport model based on turbulent kinetic energy and enstrophy[END_REF], but in a simpler configuration for isothermal plasma and 2D closed field lines domain. In this work, the perpendicular transport coefficients are determined from the turbulent kinetic energy κ and the enstrophy which is an invariant in 2D turbulence. 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].

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. The evolution equations for κ and ε are also compared to coupled Ginzburg-Landau amplitude equations. 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 L-mode plasma discharges and two different tokamaks, namely TCV and WEST. Finally, a discussion and conclusion Section closes the paper.

1 κ-ε and predator-prey models of turbulent transport 1.1 Bridging turbulence and transport models

The inherent complexity of turbulence, and consequently of the transport that it governs, is not readily modeled. When aiming at first principle simulations, significant computing resources are required. For plasmas, in view of ITER operation or reactor designs, reduced models must be considered to provide relevant information on appropriate return times. Predictive modeling capabilities are needed for ITER, to prepare, operate, possibly adapt the experiment scenarios on the fly and analyze the results. Coarse graining and building reduced model with predictive capability must then be regarded as crucial and fundamental research. Preserving locality, hence the dependence of transport on local properties, is a straightforward step to preserve simplicity of the numerical schemes. 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). These quantities are usually ad'hoc coefficients tuned to match available empirical evidence, which considerably impedes the predictive capability. A major mathematical change in the structure of the transport equation is the step from microscopic convective turbulent transport to a larger scale diffusive transport. 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 if determines in fact a projection operator and consequently the physics that one retains and that which is orthogonal. 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 times 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].

As recalled in the introduction, the standard approach of transport models is the calculation of the fluxes together with a fitting procedure of experimental profiles to determine the gradients and consequently, for each measurement point, the effective transport coefficient is defined as the flux divided by this gradient. For steady state conditions the calculation of the diffusion coefficients can therefore be seen as an alternative representation of the data. The use of these coefficients to analyze transients and to investigate transport at other locations than that used to determine them is a key step in the modeling effort. This is especially true for plasma-wall interaction where the profiles are typically measured in the midplane and where poloidal symmetry does not hold. Furthermore, in the divertor volume the plasma can enter very different parameter regimes. This situation can therefore be seen as reminiscent of that encountered in neutral fluids where turbulent properties can depart significantly from that observed in regions where measurements are more readily available. As final remark, one must keep in mind that the apparent universality of diffusive transport models is connected to the fact that it is in practice the simplest local transport model, together with ballistic transport, which can be implemented in equations. While ballistic transport questions the connection to the departure from thermodynamic equilibrium [START_REF] Ghendrih | Thermodynamical and microscopic properties of turbulent transport in the edge plasma[END_REF], non-linear dependencies of the diffusion coefficient, as with the κ-ε model, can drive transport properties that depart significantly from that of constant diffusion processes. As stated in the introduction, one can then expect that the use of a κ-ε model for plasma-wall interaction can improve the model consistency and predictive ability. We follow the standard point of view where κ is the kinetic energy per unit mass of the fluctuating transverse velocity, typically ≈ v 2 E and ε is a damping process acting on κ. In the two following Sections we discuss the connection of the κ-ε approach to models presently used to investigate plasma turbulence, both predator-prey and quasilinear transport models.

Predator-prey models

In neutral fluids the generic form given to the κ -ε model is:

∂ t κ + ∇•K = S κ -ε (2a) ∂ t ε + ∇•E = S ε (2b)
In Eq.( 2a), κ stands for the kinetic energy of the velocity fluctuations per unit mass, ∇•K stands for the divergence of the flux of κ, S κ is the source term generating the velocity fluctuations, and ε is the term driving the damping of the fluctuations. The evolution equation for the latter Eq.( 2b) is built to be similar to that of κ, with a term ∇•E that accounts for transport, a source term S ε . Regarding dimensionality and role, the field ε is the rate of energy dissipation. This quantity is also introduced to determine the universal power law dependence of the energy spectrum function, se so-called energy cascade, yielding the well known law 5/3 ε -1/4 where is the characteristic scale, hence the inverse of the wave vector. In such a framework, one obtains κ ∝ 2/3 ε 2/3 . It is to be noted that when considering the evolution equation of the fluctuating velocity field, typically starting from the Navier-Stokes equation, one finds that it couples the average of the velocity fluctuations to some power m to that at power m + 1. The problem is akin to that of the fluid hierarchy coupling moment m to moment m + 1 and, consequently, leads to the same closure issue. One can then read the κ-ε equations as a particular closure.

Ignoring transport, and therefore the coupling to neighboring positions in space, the system Eq.( 2) takes the following form:

∂ t κ = γ κ κ -β κ κ 2 -ε (3a) ∂ t ε = γ ε ε -β ε ε 2 (3b)
The source term for a field k, k standing for either κ or ε, appears therefore to be expanded in terms of k up to order 2,

S k = γ k k -β k k 2
, where γ k is a growth rate and where -β k k 2 is the first non-linear term. Although, from the point of view of an expansion of S k there is no constraint on the sign of the parameters γ k and β k , we consider the case such that both are positive. The non-linear term then ensures that the point k → +∞ is not a stable point. The dynamical system Eq.( 3) is then a rather generic local evolution system for two fields κ and ε in R + with comparable internal dynamics. While the chosen form for the source S k can be argued to be generic, the coupling between the two fields is specific. Indeed, the coupling only appears in Eq.( 3a) with the term -ε, therefore a damping of κ by ε.

But for this particular coupling, the system Eq.( 3) is quite generic. Setting ε ∝ A 2 in Eq.( 3b) one recovers the wave-amplitude equation for A known as the Ginzburg-Landau equation (with the cubic non-linear term) [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]. Alternatively, one can consider 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 predator-prey dynamics [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF] have been introduced in plasma-physics to investigate bifurcation-like phenomena of transport properties [START_REF] Diamond | Selfregulating shear flow turbulence: A paradigm for the l to h transition[END_REF][START_REF] Floriani | Selfregulation of turbulence bursts and transport barriers[END_REF]. We thus use this framework to analyze the κ-ε system Eq.( 3). The field ε is the predator and κ the prey. The predator appears to evolve independently from the prey with fixed points ε = 0 and ε = γ ε /β ε . Following the simplest case of the Ginzburg-Landau system for β ε > 0, one thus finds that for γ ε < 0, ε = 0 is the stable fixed point and when γ ε > 0 the stable fixed point switches to ε = γ ε /β ε . When considering ε as the rate of energy dissipation and assuming that there is a loss mechanisms analogous to viscosity at the smallest scales one must then relate the fixed point ε = γ ε /β ε to κ according to κ ∝ 2/3 ε 2/3 , consequently:

γ ε /β ε ∝ κ 3/2 , therefore γ ε /β ε = γκ * (κ/κ * ) 3/2 .
Here γ > 0 is a proportionality factor, dimensionaly the inverse of a time, and κ * ∈ R + is a convenient normalization of κ. The fixed point for the latter is then determined by:

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

κ κ * + κ κ * 1/2 - γ κ γ γ β κ κ * = 0 (4b)
Let us step towards a more general relationship by letting ε = γκ * (κ/κ * ) η .

In the standard neutral fluid η = 1 has been considered, while argument used when analyzing the scaling law of turbulent energy leads to η = 3/2 as discussed above. We then define κ * = γ * /β κ > 0 so that the fixed point is either κ/κ * = 0 or the solution of:

κ κ * + γ γ * κ κ * (η-1) - γ κ γ * = 0 (4c)
Depending on the problem of interest one can either define γ * = γ or γ * = γ κ . Assuming the former, hence for given γ > 0, one then finds, as with the standard Ginzburg-Landau wave amplitude equation, that for γ κ < 0, the only possible stable fixed point is κ = 0, while for γ κ ≥ 0 the stable fixed point is κ = 0 or the solution given by Eq.( 4c). For the two reference values of η, η = 1 and η = 3/2, analytical solutions are obtained: Let us now assume that γ κ > 0 is given, we then set γ * = γ κ , varying the coupling between the two fields κ and ε via the proportionality factor γ, Figure 1. For γ = 0, ε = 0 and the two fields are decoupled, κ/κ * = 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.

η =1 ; κ = γ κ -γ β κ (5a) η = 3 2 ; κ 1/2 = - γ 2(γ κ β κ ) 1/2 + γ 2(γ κ β κ ) 1/2 2 + γ κ β κ 1/2 (5b)
Compared to the case without coupling, γ = 0, adding the damping by ε introduces some form of self regulation as well as a possible delay in time because it is governed by an independent equation with different time scales. The latter effect will not occur if the time scales that govern the evolution of ε are small, hence γ → +∞ so that ε exhibits an adiabatic response. The regime with γ/γ κ ≤ 1 can thus be expected to exhibit such time delays and their impact of the dynamics of the system. It is to be underlined that in the standard κ-ε the non-linear saturation contribution in the ε equation is ∝ ε 2 /κ rather than ∝ ε 2 /κ 3/2 as retained here. This value is a marginal value, splitting η ≥ cases from that with η < 1. As shown on Figure 1, this governs a slight modification of the behavior, in particular instead of obtaining an asymptotic convergence towards zero when the control parame-ter γ is increased, both κ and ε switch to zero above the critical value γ = γ κ .

Adding transport to the evolution equation is important in the framework of deriving a transport model with as few as possible free parameters, however the key feature is the local dynamics that we have analyzed above. 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. Furthermore, the loss term ε, here governed by the losses at the Kolmogorov dissipation scale via the energy cascade, 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. Indeed, it is understood that either the nonlinear coupling of the source term to the free energy, hence the dependence of γ k 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 | Selfregulation 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 | Selfregulation 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]. In this paper we show that including the shearing effect, without complex and highly nonlinear dynamics, allows one generating an interface barrier self-consistently, Section 4.5.

Quasilinear transport models

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 (6)
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 ( 6) 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 (7a) ∂ t f 1 + L 0 f 1 + L 1 f 0 + L 1 f 1 -L 1 f 1 = 0 (7b)
The quasilinear theory simplifies Eq.( 7) by dropping the non-linear fluctuation terms in Eq.( 7b) hence removing L 1 f 1 -L 1 f 1 while retaining L 1 f 1 in Eq.( 7a). 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 Eq.( 7b) while retaining an order 2 term in Eq.( 7a) 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 Eq.( 7b). The latter can then be rewritten as:

∂ t f 1 + L * 0 f 1 + L 1 f 0 = 0 (8a) L * 0 f 1 = L 0 f 1 + L 1 f 1 -L 1 f 1 (8b) f 1 = - 1 ∂ t + L * 0 L 1 f 0 (8c)
The formal solution Eq.( 8c) 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 (9)
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 Eq.( 9), 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 (10a) D QL = 1 2 v 2 1 τ QL = κ τ QL (10b) τ QL = 2 ∂ t + v ∇ (10c)
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 Eq.( 10b). 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 ε ( 11
)
where C QL is a proportionality factor. It is also interesting to note that the time κ/ε also appears in the steady state analysis of the local evolution of the Predator-Prey system Eq.( 4c). The relations we have obtained are rather complicated because we have chosen the proportionality factor γ as free parameter. When considering as free parameter the characteristic time τ = κ/ε, one obtains:

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

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

D D * = κ 2 ε * εκ 2 * = γ κ τ -1 (13c)
The form Eq.( 11) 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 ε.

The SolEdge2D-EIRENE suite of codes

The 2D transport equations and boundary conditions implemented in the code SolEdge2D-EIRENE have been derived in [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF], and are recalled in appendix A. The model 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 condition can be relaxed then requiring the charge balance equation determining the electric potential to be solved. The geometry of the magnetic field plays a crucial role owing to the large anisotropy between the transport in the directions parallel and transverse to the magnetic field B = B b, b = B/B is the unit vector along the magnetic field defining the parallel direction. Indeed, the drift velocities of the charged particles are induced by the magnetic field, and of order ρ * = ρ 0 /a, the particle velocity in the parallel direction is of order 1 and transport in that direction is classical weakly collisional transport. The dimensionless parameter ρ * is a measure of the importance of the strength of the magnetic field accounted for by the typical Larmor radius ρ 0 compared to the tokamak minor radius a: in a strongly magnetized plasma ρ 0 a. The ρ * dimensionless parameter also characterizes the number of degrees of freedom of the system typically 1/ρ 2 * . For the purpose of using reduced models it is important to stress that the reference Larmor radius also depends on the typical thermal velocity. For ITER, one can expect that it is about 100 times smaller in the divertor region than in the core plasma, implying an increase or the number of degrees of freedom of the order 10 4 .

The structured mesh is based on a grid aligned onto the magnetic flux surfaces for numerical efficiency. An explicit domain decomposition technique allows one handling any complex magnetic geometry. An example of such decomposition is shown on Fig. 3 for a WEST plasma configuration [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]. 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] to embed any realistic tokamak wall geometry within the computational domain. An implicit/explicit Eulerian temporal scheme is used for the time integration. The code is parallelized using openMP and MPI libraries. Magnetic measurements at different locations surrounding the vacuum vessel are used as real-time inputs to achieve the numerical reconstruction of the plasma current density and the magnetic equilibrium, see for example [START_REF] Blum | Reconstruction of the equilibrium of the plasma in a tokamak and identification of the current density profile in real time[END_REF]. 3: 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.

In all equations of A, perpendicular turbulent fluxes have been modeled by diffusive and conductive perpendicular fluxes. Their coefficients D n , χ e , χ i and ν have to be tuned to provide the right balance between perpendicular and parallel transport in the model. These coefficients are defined by a matching procedure using experimental midplane profiles as input to the code. In the following, Section 3, a model is proposed and implemented in SolEdge2D-EIRENE to self-consistently estimate these coefficients in the edge and SOL non-isothermal diverted plasmas in steady-state.

3 The κ-ε model for edge and SOL plasma

The κ-ε model implemented in SolEdge2D

In the present work, we assume as in many references in neutral fluids that the turbulent Schmidt number Sc t = ν t /D n is of order unity [START_REF] Gualtieri | On the values for the turbulent schmidt number in environmental flows[END_REF], therefore D n = ν t /Sc t ≈ ν t . Furthermore, the empirical observation of the midplane profiles of density and temperature are matched by steady state transport simulations with the constraint χ e,i = 2 × D n = 2ν t /Sc t [START_REF] Chankin | SOLPS modelling of asdex upgrade h-mode plasma[END_REF]. The turbulent Prandtl number P r t = ν t /χ is therefore P r t = 0.5 × Sc t ≈ 0.5. As discussed in Section 1, the transport coefficients are determined in terms of κ and ε Eq.( 11) and set ν t accordingly:

ν t = C ν κ 2 ε ; D n = ν t Sc t ≈ ν t ; χ e = χ i = ν t P r t ≈ 2ν t (14)
with C ν as a constant parameter to be determined with an appropriate constraint, typically a fitting procedure of either empirical or simulation data. This approach allows one taking into account the scales of both turbulent production and dissipation independently. Indeed, the time evolution of κ and ε is governed by two characteristic transport equations Eq.( 2). This general form is given by:

d(•) dt = Production -Saturation + Diffusion (15)
It is to be noted that the diffusive transport Diffusion will depend on the two fields according to Eq.( 14). This adds a further non-linearity to the system that is not addressed in Section 1. The terms Production and Saturation determine the local behavior of the two fields Eq.( 3). The choice is made to base the production on the physics of the linear interchange instability, which is understood to be one of the the main drives of turbulence in the edge and Scrape-Off Layer (SOL). In particular, its linear growth rate γ leads to the poloidal asymmetry observed in the turbulent transport between the LFS and the HFS, the so-called ballooning of the turbulence which is not consistently taken into account in current approaches used in 2D transport codes. As already mentioned, the physics attached to the energy dissipation being complex in plasmas and still poorly understood, a quadratic saturation term the dissipation analogous to that used in CFD is chosen here. The local behavior is thus identical to that considered on a generic footing in Section 1.2, Eq.( 2). This leads to the following 2D transport equations for κ and ε:

∂ t κ + ∇ (κu b) -∇ ⊥ • (D κ ∇κ) = γ κ κ - 1 D ω κ 2 -ε (16a) ∂ t ε + ∇ • (εu b) -∇ ⊥ • (D ε ∇ε) = γ ε ε - V κ 3/2 ε 2 (16b)
The left hand side of these two coupled equations is the evolution induced by the convective and diffusive transport while the right hand side governs the local dynamics, the effective source and sink terms. In this system as well as in Eq.( 3), the coupling term is defined with no control parameter. The relative magnitude of the two fields is therefore fixed by this form of this coupling term. Furthermore, the generic terms driving the quadratic saturation in Eq.( 2), both β κ and β ε are replaced respectively by 1/D ω and V . These notations are consistent with the dimensionality of these control parameters, V being a velocity and D ω a diffusion coefficient.

The second term on the left hand side of these equations is the parallel component of the advection term. Consistently with the other equations of the model and quasilinear theory Eq. [START_REF] Tamain | The tokam3x code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF], and in the absence of large scale drifts, the cross-field transport of κ and ε is taken into account by a diffusive process, governed by the diffusion coefficients D κ and D ε , third term on the left hand side of Eq.( 16). The first term on the right-hand side of both equations 16a and 16b is the growth rate of the leading instability.

Model closure: defining the growth rates

When closing the system by choosing the values of the free parameters, we shall consider the simplification γ κ ≈ γ ε = γ I , and we shall consider γ I to be given by 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].

γ I = c s R R 2 ∇p i •∇B T p i B T ( 17
)
with R is the tokamak major radius and B T the toroidal magnetic field. The proportionality factor α 0 is the control parameter that determines the characteristic growth time of κ-ε given the functional dependence on the other fields, here p i playing a key role given that for an axisymmetric magnetic geometry RB T is constant. Although one can consider that using ( 17) is only consistent with the choice we make that plasma turbulence is driven by the interchange instability, one must stress that while Eq.( 17) is a global feature, determined for a flux surface and a particular geometry of a global eigen-mode (usually approximated by a Fourier mode), the control parameters γ κ and γ ε are local. Consequently they determine the local growth rate of the fields κ-ε. The implicit rule in Eq.( 17), namely that γ I = 0 when (∇p i /p i ) • (∇B T /B T ) ≤ 0 governs a ballooning of the local drive. As will be shown in Section 3.4, a vanishing drive combined to the other transport properties will govern an effective damping of the turbulent energy. The result 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.

The second term on the right-hand side of Eq.( 18a), κ 2 /D ω , stands for a self-saturation of the turbulence energy. It keeps the system stable, by preventing κ from growing to infinity. The parameter D ω is a free parameter, with dimension m 2 s -1 is a diffusion parameter, which can also be written as D ω = κ ω /∆ω. These quantities can be understood as the spectral width ∆ω in terms of frequencies of the turbulent energy, and a characteristic value of the turbulent energy. This single tuning parameter can thus be split into two components to provide the connection to published predator prey models [START_REF] Miki | Spatio-temporal evolution of the l → i → h transition[END_REF]. As discussed in Section 1, the right hand side term ε is a sink term given that one assumes ε to be positive. Finally the non linear term (V /κ 3/2 )ε 2 is the sink term in the ε evolution equation. The latter is controlled by the parameter V and has been designed such that in steady state conditions of the local evolution, Section 1, one gets ε ∝ κ 3/2 . Assuming that the velocity V scales like √ κ, one can show that the control parameter does not depend on the chosen normalization of κ.

With this choice of closure parameters, the local evolution equation for κ and ε Eq.( 3) takes the form:

∂ t κ = γ κ κ - 1 D ω κ 2 -ε (18a) ∂ t ε = γ ε ε - V κ 3/2 ε 2 (18b)
It is to be noted that a more complex form of the growth rate, including in particular a threshold in pressure gradient as well as other physics such as a Doppler shift governed by parallel transport is 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 features of parallel convection and diffusive damping appear independently in the evolution equations of κ-ε Eq.( 16), the present approach consists of using the growth rate proportional to γ I Eq.( 17) and let the physics included in the κ-ε system Eq.( 16) modify this drive. This aspect will be further discussed in Section 3.4.

Closure constraints for the parameters D ω and V

The control parameter of the quadratic term of the equation for κ, D ω in Eq.( 18a), 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. For the present case, it ensures the existence of a fixed point even when ε = 0, which is a steady state solution of Eq.( 18). As shown in Section 1, the steady state solution is obtained as the solution of a second order equation. This equation can be recast in an equation for the diffusion coefficient D = κ 2 /ε:

D 2 + DD ω -D V D ω = 0 (19a) D V = γ κ V 2 γ 2 ε (19b)
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 ω (19c)
Depending on V , one finds therefore two regimes, the weak turbulence regime for D V D ω , such that D ≈ D V ans 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 ω is the Kubo number that governs this transition. As a possible closure of the system, we relate the diffusion coefficient D in the weak turbulence regime to the SOL width determined by an empirical scaling law proposed in [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF]. In such a framework, we consider that the transition to the Bohm regime for D occurs for much larger values of D than required to match the experiments. The parameter D ω must therefore be chosen such that D ω D n , in the simulations D ω is defined in terms of κ max = 10 10 , D ω = κ max /γ κ . The exact value of this parameter will then have little impact on the behavior of the κ-ε system.

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. We therefore enforce the relationship:

D = C ν D V ≈ λ 2 SOL τ (20a)
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]. Accordingly, we assume λ SOL /a = α s qρ * , α s being the proportionality factor. Note that the chosen data base used for the empirical analysis does not properly account for the aspect ratio A variation. However, the dependence on such a parameter would be useful to discriminate various models. Alternative closures are discussed together with other theoretical considerations in Section 4.4. Given the closure chosen here, we obtain D V :

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

V = c s ρ * v 0 τ γ 2 ε γ κ (20c)
To obtain this expression of V we have considered L = qR and defined τ = L /c s . 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.

As discussed for the growth rate, it is important to underline that the SOL width is a global parameter since it describes the large scale balance between parallel and transverse transport. However, we introduce the constraint determining V Eq.( 20c) as a local property. Consequently, although we use the global interchange growth rate and the empirical scaling of the SOL width to constrain the local free parameters of the κ-ε model, this does not enforce that the global properties of the simulation will exhibit the empirical trends and in particular that the local property will govern the appropriate scaling of the global SOL width.

Linear analysis of the κ-ε model

Let us consider the linearized κ-ε model in Fourier space:

∂ t κ + ik u κ + k 2 ⊥ D κ * κ = γ κ κ - 2 D ω κ * κ -ε (21a) ∂ t ε + ik u ε + k 2 ⊥ D ε * ε = γ ε ε - 2V κ 3/2 * ε * ε + 3V κ 5/2 * ε 2 * κ (21b)
where κ * , ε * are the steady state solution homogeneous in space, hence solutions of the local system Eq.( 18) 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 Eq.( 21) describes the evolution of the system disturbed from the steady state solution by a perturbation varying in space. The growth rate γ of this perturbation is then determined by the non trivial solution of the following linear system:

γ + ik u + k 2 ⊥ D κ * -γ κ + 2 D ω κ * κ + ε = 0 (22a) - 3V κ 5/2 * ε 2 * κ + γ + ik u + k 2 ⊥ D ε * -γ ε + 2V κ 3/2 * ε * ε = 0 (22b)
The dispersion relation is then given by setting the determinant of this system to zero and takes the generic form:

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

A = γ κ -k 2 ⊥ D κ * - 2 D ω κ * -ik u (23b) B = γ ε -k 2 ⊥ D ε * - 2V κ 3/2 * ε * -ik u (23c) C = 3V κ 5/2 * ε 2 * (23d)
One can then observe the Doppler shift yielding the frequency k u as well as the stabilizing effect governed by the diffusion coefficients. This term governs a damping of the perturbation of the form k 2 ⊥ D * . The contribution of the transport terms thus generate damping processes that inhibit the drive terms and consequently govern the occurrence of threshold effects.

As a final remark regarding the dynamics of the system, it is important to note that 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 ∇ ⊥ κ κ - ∇ ⊥ ε ε ( 24
)
The divergence is simplified by considering here a slab geometry. Such a structure is comparable to that induced by a convective motion, hence 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 ∇ 2 ⊥ F + V pinch ∇ ⊥ F .
We have assumed here for the sake of simplicity that D pinch and V pinch are constant. Identifying the latter form of the flux divergence with expression Eq.( 24), one obtains

V pinch = D(∇ ⊥ Log(ε) -2∇ ⊥ Log(κ)
). The fact that κ and ε controlling the diffusion coefficient can vary in space thus generates transport features that are reminiscent of convection although the 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. It is then possible that the reduced cross-field transport model implemented in transport codes such as SolEdge2D-EIRENE 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 behaviors have been reported [START_REF] Loarte | Self-sustained divertor plasma oscillations in the jet tokamak[END_REF].

4 Physics of κ-ε transport in a 1D model

The 1-D κ-ε transport model

The 1-D model is obtained by averaging the SolEdge2D transport equations in the poloidal and toroidal directions. 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

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

                 ∂ t n - 1 r ∇ r rD n ∇ r n = S n -H(r -a) n τ ∂ t E e - 1 r ∇ r rD n T e ∇ r n + χ e nr∇ r T e = S E -H(r -a) nT e γ e τ ∂ t E i - 1 r ∇ r rD n T i ∇ r n + χ i nr∇ r T i = S E -H(r -a) nT i γ i τ (25a)
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 . However, 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. It will be be further discussed when addressing the 2D simulations, Section 5, where the ballooning nature of turbulent transport is an important property recovered with the κ-ε model, yielding in particular χ =< χ >.

The plasma equations are completed by a vorticity equation

Ω ≡ m i ∇ ⊥ • en∇ ⊥ φ/B 2 + ∇ ⊥ p i /B 2 .
It will be used when addressing the impact of velocity shear on the transport properties. In this definition of the vorticity Ω, the electrostatic potential is φ, e the electron charge, and p i the ion pressure. This equation is derived from the charge balance equation including the polarization drift current for the ions.

∂ t Ω - 1 r ∇ r rν∇ r Ω = H(r -a) 1 τ eφ T e -Λ (25b)
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. Neutrals are also considered since they govern the particle source by ionization S n . Only the particle balance equation is used for the neutrals with density n 0 , diffusive transport with constant diffusion coefficient D n0 , a source term Φ n0 and the ionization sink S n .

∂ t n 0 - 1 r ∇ r rD n0 ∇ r n 0 = Φ n0 H(r -a) -S n (25c)
The transport model for the neutrals can readily be replaced by convection.

The source term Φ n0 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 Φ n0 constant throughout the SOL region, hence of the form Φ n0 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 (25d)
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 in fact 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 1-D model Eq.( 25d) must therefore be understood as effective parameters, as discussed for plasma transport Eq.( 25a) and Eq.( 25b). The values of the control parameters in Eq.( 25d) are thus modified to account for such loss terms as well as for corrections stemming from the approximation made in averaging process.

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. On the 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 Φ n0 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 4 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 source term, Figure 4 right 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 5, is defined as a function of the normalized radius ρ = r/a:

VS E = P in p 0 exp - 1 1 -ρ 2 (26a) p 0 = 1 0 ρdρ exp - 1 1 -ρ 2 = 7.425 10 -2 (26b)
where the constant p 0 = 7.425 10 -2 ensures that the cylindrical integral of VS E is P in , the injected power.

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 (27a)
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 (27b)
Within the error bars of Eq.( 27a), 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 Eq.( 27b) 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]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]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 3.3 Eq.( 20c) 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 and the effective 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 1-D 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 γ e and χ e /D n , Eq.( 14). 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 of particularly suitable to investigate the 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 Eq.( 27a). 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 6. 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 when increasing P in at constant separatrix density, Figure 7 left hand side. However, one can observe that the density decrease is less pronounced for P in ge1M W compared to the very low power cases. Conversely, the thermal energies, Figure 6 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 , 7 right hand side. The overall effect on the plasma pressure p = nT is an increase. In the core, Figure 7 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 that D n exhibits a power law dependence on the injected power, D n ∝ P 0.644 in , Figure 8 left hand side. Furthermore, from the density profiles Figure 6 topleft, 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 8 left hand side. Assuming that this scaling is consistent with the choice made to determine the parameter V Eq.( 20c), 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 Eq.( 27a). 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 Eq.( 27a). There 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 the 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 8 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 6 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 efolding length λ p and compare it to the experimental scaling law Eq.( 27a). On Table 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 dimensionaly correct experimental scaling law.

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 scaling law since all the power law exponents have the same sign. However, the magnitude of these exponents observed with in the κ-ε simulations are systematically smaller. We first discuss the results for the power scan, Figure 9 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 Eq.( 27a) 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 trend is recovered with the κ-ε simulations, Figure 9 right hand side, but the exponents is nearly a factor 1.5 smaller than that in Eq.( 27a), and closer to ≈ -0.5 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]. The depen- dence on R, keeping all other engineering parameters constant, hence minor radius, poloidal and toroidal magnetic field, density and input power, Figure 10 left hand side, is mainly governed in Eq.( 27a) 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 10 left hand side, but the power law behavior is clearly different. Finally, when varying the toroidal magnetic field alone, Figure 10 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 Eq.( 27a) several issues are raised which can readily explain the observed differences. First, as we have just recalled, the reference expression is not dimensionaly 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

In the Section, we discuss the theoretical background of the proposed κ-ε model. The starting point of the analysis of turbulent transport common to neutral fluids and plasmas are the velocity fluctuations. 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. One recovers therefore that V Eq.( 20c) is the relevant velocity driving turbulent transport1 . To obtain the expression of V , we have used the characteristic time scale τ = qAa/c s . In fact two time scales play a role, τ and the linear time scale associated to the growth rate of the interchange instability, namely 1/γ I Eq. [START_REF] Reiter | The eirene and b2-eirene codes[END_REF].

γ I τ = α 0 q R 2 ∇p i •∇B T p i B T ∝ qA 1/2 (28)
This ratio of relevant time scales is thus linear in the safety factor which can prove important when stepping from the Larmor radius ρ 0 as standard turbulent scale to qρ 0 that appears to be the relevant SOL width scaling.

To obtain this simplification, we have assumed here that ∇B ≈ B/R while that ∇p i ≈ p i /a, this scaling property, commonly considered can however be questioned. Given a characteristic time τ and the velocity amplitude V one can define the characteristic ballistic length scale λ = V τ , while a characteristic diffusive length scale is √ D V τ , provided the diffusion coefficient is D V . In Section 3.3, the time scale τ is chosen to be τ . However, when retaining τ = 1/γ I the ballistic and diffusive length scales defined above are identical. With the present notations, setting γ κ = γ ε = γ I , one obtains consequently:

V = c s ρ * q √ A (29a) κ = V 2 = c 2 s ρ 2 * q 2 A (29b) ε = γ I κ = c 2 s τ ρ 2 * q 2 q √ A (29c)
We have assumed here that ε scales according to the steady state solution of the local drive for κ Eq.( 18) with D ω D n , see Section 3.3. These expressions are then used to define the diffusion coefficient D V Eq.( 19b) and the SOL width:

D V = κ 2 ε = ac s q 2 ρ 2 * 1 √ A (30a) λ = a qρ * (30b)
Compared to the expression of D V implemented in the 1-D code, Section 3.3, one finds different exponents for both safety factor q and aspect ratio A. The choice of the characteristic time scale has an impact of the final result and appears to be the most difficult to define and determine without ambiguity on a physics base. It is to be noted that the choice of the closure based on the SOL width also involves both the choice of the characteristic time and that of the dominant transport process, either ballistic with velocity V or diffusive with diffusion coefficient D V .

To account for the fact that defining the appropriate time scale is not straightforward, we introduce three undetermined exponents ρ , q and A in the definition of D V .

D V = ac s ρ 2 * ρ q q A A (31)
In view of the discussion of the characteristic time scale, q = 2, A = -0.5 stands for the choice made above, and q = 1, A = -1 for the choice made in Section 3.3. In the chosen framework, the exponent of ρ 2 * allows one taking into account a ρ * 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 Eq.( 31) will be used to determine the energy confinement time and SOL width scaling laws from the κ-ε model. Prior to that step, let us first revisit the connection between the energy confinement time and the SOL width.

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 (32a)
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 ( 32b
)
This expression is obtained with the approximation λ q ≈ qρ 0 with ρ 0 = c s /Ω 0 , Ω 0 being the reference ion Larmor gyration pulsation, leading to τ ≈ q 2 R/(Ω 0 λ q ). Such a connection between the global energy confinement time and the width of the boundary layer has been investigated experimentally [START_REF] Brunner | Highresolution heat flux width measurements at reactor-level magnetic fields and observation of a unified width scaling across confinement regimes in the alcator c-mod tokamak[END_REF]. Using the scaling of λ q Eq.( 27a) 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 (33a)
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 (33b)
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: 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 (34b)

τ L E,
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 appear to be have a small impact but the elongation k at given safety factor. Plasma shaping appears to enhance confinement Eq.( 34a) which consequently drives a reduction of the SOL width Eq.( 34b).

In the framework of the κ-ε approach, one can also estimate the energy confinement time given the diffusion coefficient D V Eq.( 30a) 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 Eq.( 31), one can estimate the confinement time τ

(κ-) E considering that τ (κ-) E ∝ a 2 /D V : Ω 0 τ (κ-) E ∝ ρ -1-2 ρ * q -q A -A (35a)
The Larmor pulsation Ω 0 = eB/m is used to enforce the standard time normalization. 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 (35b)
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]Chapter 2: Plasma confinement and transport[END_REF]. Comparing Eq.( 35a) and Eq.( 35b) 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 3.3. 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 Eq.( 35a) 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 (35c)
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 closer to the values obtained with τ ∝ 1/γ I , q = 2, A = -0.5, than with τ ∝ τ , q = 1, A = -1. 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 1-D model allows one comparing the predicted energy confinement time to that achieved in the steady state 1-D κ-ε simulations. The appropriate trend of the dependence on the input power P in is recovered, Figure 11 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 Eq.( 34a), dashed black line. One finds that the simulation data also compares well with the power law Eq.( 33b) 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 Eq.( 35a), yields P -0.6 in blue long-dash line. The scan in B pol , hence in I p appears to qualitatively right, Figure 11 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 Eq.( 34a). 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 Eq.( 33b) black dash-dot line. The theory based κ-ε models Eq.( 35a) 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 33b) is plotted with dash-dot black lines while the theory based κ-ε models Eq.( 35a) 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. 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.

Changes in the closure constraint based on zonal flow control of the turbulent energy

The proposed κ-ε model for plasma turbulent transport is closely linked to that used in neutral fluids. In particular, the field ε, which governs the dissipation of the turbulent energy, is consistent with the damping process of the turbulent energy spectra to the small scales. However, a key physics mechanism that is now understood to control the turbulent transport in plasmas is that of turbulent eddy shearing by large scale flows, typically the zonal flows [START_REF] Kim | Zonal flows and transient dynamics of the l → h transition[END_REF]. In particular, the regimes of improved confinement observed is 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 | Selfregulating shear flow turbulence: A paradigm for the l to h transition[END_REF][START_REF] Kim | Zonal flows and transient dynamics of the l → h transition[END_REF][START_REF] Miki | Spatio-temporal evolution of the l → i → h transition[END_REF][START_REF] Floriani | Selfregulation 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 the ballistic transport at the smallest scales. As discussed, this velocity is found to be of the order of the E × B drift velocity, thus proportional to c s ρ * . The discussion on the relevant time scale τ , and consequently on the length scale V τ used for the closure, typically V τ = qρ * a, indicates that a further dependence of the form q ηq /A η A can be expected. 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 for 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 ρ * 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 q ηq A η A (36a)
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 fixed point for κ-ε is such that κ = V 2 and ε = γ I κ so that D V = κ/γ I . Given 1/γ I = a √ A/c s , the modified diffusion coefficient is then:

D V = V 2 a √ A c s = c 2 s ρ 2 * a √ A c s 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 2 q 2ηq A 2η A = ac s ρ 2 * 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 2 q 2ηq A 2η A -0.5 (36b) λ = a qρ * 1 + (ω z τ z ) 2 1 + r(ω z τ z ) 2 q ηq-1 A η A -0.5 (36c)
In the 1-D simulation case the parameter ω z is defined given the electric potential profile solution of the vorticity equation Eq.( 25b). Two free parameters the 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, 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 12. The effect of shear appears to be localized just inside the separatrix, close to the inversion point of the electric field Figure 12 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 12 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 a relatively simple way of generating a feedback on turbulent transport 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 mechanism for the flow shear impact on transport is devised to achieve a transition. Alternative models governing the local behavior of κ ε driven by physics insight have been proposed to investigate the transition to transport barriers [START_REF] Diamond | Selfregulating 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]. One issue that is readily raised by the model used here is the stabilization of the microscopic turbulent velocity V Eq.( 36a). Indeed, at given amplitude of the fluctuating electric potential, the reduction by shearing of the eddy sizes drives an increase of the velocity amplitude. A reduction of turbulent transport must then be driven by a reduction of the characteristic microscopic time scale and consequently of the typical length scale governing turbulent transport. The importance of such issues to step towards a predictive capability of plasma transport models can be investigated with 1D-models as done here and then further developed with 2D modeling.

Confrontation to experiments: SolEdge2D-EIRENE modeling

We present here the results of 2D modeling with the transport code SolEdge2D-EIRENE assuming pure deuterium plasma, hence no impurity physics, and neglecting the mean field velocity drifts. The κ-ε equations Eq.( 16) are solved coupled to the plasma transport equation, see Appendix A. The parallel advection of κ appearing in the 2D equations Eq.( 16) contribute to smoothing the solution. The simulations presented in this Section illustrate the novel aspects of 2D transport modeling when taking into account the transport generation by the κ-ε fields. For the sake of simplicity, the simulations are run without solving the vorticity equation, therefore without the large scale drift terms and consequently without transport regulation by the shear of the large scale flows. Typically 10 6 iterations are needed to reach steady state conditions with, as well as without, the κ-ε model.

Table 2: TCV parameters used for the SolEdge2D-EIRENE simulations.

I p [M A] B T [T ] B pol [T ] R[m] a[m] P in [M W ] n OM P sep [m -3 ] k 0.

Confrontation to TCV data

We are interested here in lower single null (LSN) Ohmic L-mode discharges in TCV -Tokamak à configuration Variable with good coverage of plasmawall diagnostics [START_REF] Coda | Overview of the tcv tokamak program: scientific progress and facility upgrades[END_REF][START_REF] Maurizio | Divertor power load studies for attached l-mode single-null plasmas in tcv[END_REF][START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF]. Our main interest is a set of three discharges with three magnetic configurations achieved by varying the length of the outer divertor leg h X , namely the vertical distance from the X-point to the divertor strike point: the short-leg h X = 0.21 m shot #51262, the medium-leg h X = 0.36 m shot #51333, and the long-leg h X = 0.64 m shot #51325. The other plasma parameters remain nearly unchanged: major radius R = 0.89 m, minor radius a = 0.22 m, elongation k = 1.4, plasma current I p = 210 kA. For the present purposes, we refer to the TCV data published in 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] and Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. In the SolEdge2D-EIRENE simulations, the gas-puff is adjusted to obtain the same density at the separatrix outer midplane (OMP) for the three cases, n OM P sep = 0.7 10 19 m -3 and the input power crossing the separatrix into the SOL P in is set at 100 kW , see Table 2. The flexibility of the penalization technique to implement the boundary condition at the wall allows one running realistic simulations of the complete volume, thus including any effect governed by the relative positions of the wall and magnetic surfaces. We also use here patterns of transport properties from TOKAM3X simulations with the TCV geometry and published in 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]. This set of TCV experiments are particularly interesting because they allow one varying the relative role of SOL cross-field transport in the vicinity of the core plasma and that along the divertor leg from the X-point to the strike point. The relative change can be characterized by the ratio h X /(ka) equal to 0.68, 1.17 and 2.89 for the three chosen experiments. It is to be underlined that for homogeneous transport properties, the weight of the divertor cross-field transport should scale from 0.68 to 1.17 and 2.89 compared to the "main chamber" SOL transport and that this effect should be identical in the private flux and divertor SOL region. The simulation results will therefore be sensitive to the explicit or implicit choice made regarding the extrapolation of the diffusion coefficient to the whole poloidal plane. The κ-ε model is particular in that respect since the transport coefficients throughout the poloidal plane is generated consistently by the transport, source and sink of κ and ε. This process evolves in time with the plasma pressure gradients that control the growth rates of κ and ε. For these first tests with κ-ε transport self-organizations, the optimization of the free parameters was not addressed. Only tuning the magnitude of a single scalar, the proportionality factor C ν in Eq.( 14), for all transport coefficients and for the whole simulation domain was performed.

Midplane and divertor profiles

The confrontation of the simulations and the experimental profiles remapped to the outer midplane has been presented in Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. In the outer mid-plane plasma one finds in both the experiments and the simulations that the profiles do not exhibit any significant variation with the X-point height h X . This appears to indicate that when maintaining the same plasma control parameters for the core, hence same input power, plasma current, gas puff, toroidal magnetic field, there is no dramatic change in the cross-field properties. In Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF], the midplane radial profiles from the simulations and from the experimental data for the different lengths of the divertor leg are observed to be in good agreement for the SOL density and the electron temperature profiles. Some difference, typically in the range of 15 %, can however be observed regarding the density profiles in the edge plasma region, from the separatrix to the inner boundary of the simulation domain. Given the fact that a single free parameter has been tuned for the three configurations and the whole simulation domain, this level of agreement is already a demonstration of the capability of the κ-ε model to self-consistently determine the transport properties in SolEdge2D-EIRENE.

Let us now compare the heat flux profiles measured by the Infra-Red system to the output of the SolEdge2D-EIRENE simulations, Figure 13, for the three configurations, from left to right panel, small, medium and large X-point heights. Here, the profiles are normalized to the peak value to highlight the comparison of the shape of the deposition footprint. Also, the background heat flux to the divertor floor in the private flux region and far SOL is subtracted from the experimental data for the sake of this comparison since the simulations are performed without impurity radiation. Note that the data points are not actual TCV data but a sketch of the measurements published in [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF][START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. The overall agreement between model predictions and experimental data is fair, the peak values agree within 10 % for the short leg, the error is larger for the long leg case, reaching 30 %. For the medium case, the data at the peak value is missing making the comparison more difficult. Regarding the shape of the divertor heat flux deposition, Figure 13, Figure 13: Numerical and sketch of the experimental profiles published in [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF][START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] of the perpendicular heat flux remapped at the outer midplane for the three TCV magnetic configurations. The sketched Infra red data from TCV is indicated by the open squares, SolEdge2D-EIRENE output plain line. When needed, the latter is slightly shifted radially using the upper axis. For the medium leg case, the experimental profile in the peak region is extrapolated, dashed line.

the simulations tend to underestimate the fall-off length towards the private flux region, also the SOL fall-off exhibits a close to exponential decay while the experimental data suggests a variation of the decay lengths between the near and far SOL. The case of the medium divertor height, center panel, is more difficult to analyze because of the missing data. To characterize the results we consider the estimate of the SOL heat flux channel using the now standard fitting procedure [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF], Figure 17. The corresponding values are labeled λ SOL for both the simulation and experimental results.

Transport properties in the poloidal plane

In this Section we revisit the transport analysis that has been performed for this set of TCV experiments. In 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] SolEdge2D-EIRENE simulations without the κ-ε model are presented together with micro-turbulence isothermal simulations with the code TOKAM3X [START_REF] Tamain | The tokam3x code for edge turbulence fluid simulations of tokamak plasmas in versatile magnetic geometries[END_REF]. For the purpose of the present discussion we reproduce in this paper the salient features of the transport patterns in the poloidal plane. Comparable SolEdge2D-EIRENE simulations with the κ-ε model are found in Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF].

Let us first analyze the SolEdge2D-EIRENE simulations 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 14. 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. To highlight the key properties of the simulation outputs, three levels of the particle flux have been identified, strong outward particle flux labeled "High > 0", moderate outward particle flux "Medium > 0", and moderate reversed particle flux labeled "Medium < 0". The latter highlights transport into the private flux region and is located in the vicinity of both divertor legs. 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 scrapeoff-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 14, panel (b). A modest improvement when Figure 15: Sketch of the 2D transport pattern in TCV in the large Xpoint height configuration. Panel a: effective particle diffusivity given by TOKAM3X, sketch of the results published in 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]. Panel b: Particle turbulent diffusion coefficient D n self-consistently determined in SolEdge2D-EIRENE simulation using the κ-ε 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". 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], but still the key experimental trend when varying the X-point height is not recovered.

The effective particle cross-field diffusion coefficient as been estimated from micro-turbulence simulations in the TCV geometry running the TOKAM3X code [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF]. The characteristic transport pattern is sketched on Figure 15 panel a. The same kind of leveling for the amplitude of the diffusion coefficient has been chosen, "High" for the regions with large values and "Medium" for the regions with moderate values. On the same Figure 15 panel b, the 2D pattern of a SolEdge2D-EIRENE simulation with the same TCV geometry and using the κ-ε turbulent transport model [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] is presented. The self-consistent calculation of D n in the latter simulation readily exhibits turbulent transport ballooning. As reported for all choices of the X-point height, these TCV simulations are characterized by a peaking of D n in the plasma low field side equatorial plane [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF]. The first observation that can be made on Figure 15 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 coefficient is observed to increase from the in-Figure 16: Sketch of the 2D transport pattern in TCV in the large X-point height configuration. Panel a: particle flux from TOKAM3X, published in 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]. Panel b: Particle flux determined with SolEdge2D-EIRENE simulation, 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]. ner 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 edge-core interplay in fusion plasmas[END_REF]. However, the particle diffusion coefficient determined by the κ-ε model appears to be maximum at the separatrix and then decays into the SOL. This feature can be observed in Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] is not reported on Figure 15 panel b. Such a roll-over appears to be observed in some experiments [START_REF] Prisiazhniuk | Density fluctuation correlation measurements in asdex upgrade using poloidal and radial correlation reflectometry[END_REF]. Regarding the ballooning aspect, one finds that the poloidal extension 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]. 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. This feature is specific of the SOL divertor region on the low field side and does not appear in the private flux region. Some transport enhancement, localized on the high field side divertor separatrix is also reported in Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] and not shown here. Comparing the micro-turbulence and κ-ε transport pattern indicates that qualitative agreement is achieved. A possible route for optimization of κ-ε transport in fusion plasmas could be achieved by stepping to quantitative agreement based on a set of simulations for different tokamak geometries and control parameters.

We now compare the particle flux determined by a TOKAM3X simulation Figure 16 panel a, 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], and the SolEdge2D-EIRENE simulations from Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF], Figure 16 panel b. One recovers more or less the pattern observed for the magnitude of the turbulent particle diffusion coefficient but with further peaking in the regions with larger density gradient, hence on the low field side separatrix and extending more or less along the separatrix on the low field side divertor leg. A region with large negative cross field particle flux is localized on the high field side divertor leg, in the very vicinity of the separatrix, can be observed in Ref. [START_REF] Baschetti | Study of the role of the magnetic configuration in a k-model for anomalous transport in tokamaks[END_REF] but is not reported on Figure 16. The overall turbulent transport features agree with experimental observations such as the ballooned feature [START_REF] Labombard | Transport-driven scrapeoff-layer flows and the boundary conditions imposed at the magnetic separatrix in a tokamak plasma[END_REF], the localization along the low field side divertor leg, as 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]. Some differences are observed between the κε transport modeling and the 3D turbulence simulation with TOKAM3X. The latter for TCV-like geometry [START_REF] Gallo | Impact of the plasma geometry on divertor power exhaust: experimental evidence from TCV and simulations with SolEdge2D and TOKAM3X[END_REF] as well as that for Compass-like geometry [START_REF] Galassi | Drive of parallel flows by turbulence and large-scale e×b transverse transport in divertor geometry[END_REF] are performed with isothermal conditions in include the large scale electric drift flows which are not computed in the presently available SolEdge-EIRENE simulations with κ-ε turbulent transport. In particular, the 3D micro-turbulence simulations exhibit large transport in the regions with large flux expansion in the X-point region as well as in the top region towards the secondary X-point. Such features could be a signature of poloidal spreading from the large turbulence regions into stable regions. This important facet of self-organized turbulent transport is accounted for explicitly in the chosen evolution equations of the fields κ and ε where convective transport is restricted to parallel transport. The latter appears to exhibit some features of spreading by smoothing the pattern of turbulent transport. Improving the large scale convective transport in the κ-ε turbulent transport model could drive turbulence spreading and consequently improve the predictive capability of the simulations.

Regarding the particle transport in SolEdge2D-EIRENE simulation, Figure 16 panel b, one can notice that the large amplitude region of the cross-field particle flux extends in the plasma midplane up to the main chamber wall. In the simulation, the latter therefore acts as a large surface secondary limiter which scrapes-off the parallel flux into the divertor volume. The shadow of this secondary limiter thus tends to narrow the particle deposition pattern on the divertor floor, this shadowing effect being all the more important that the connection length is reduced. Another issue is the impact on the core high performance scenarios of such spurious particle recycling in the main chamber volume. On short pulses the response of the wall, which depends considerably on the wall material, will play a major role. When considering long pulses, the effective recycling coefficient in that area will tend towards unity, Table 3: TCV parameters and SOL width 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 into account the error bars. inevitably modifying the overall particle recycling pattern. This physics can be addressed in the SolEdge2D-EIRENE transport simulations because the plasma volume is meshed up to the wall and because the latter is not aligned on a magnetic surface. This geometry cannot be handled with the present TOKAM3X setting.

A q 95 B T /B pol λ H [mm] λ L [mm] λ q [mm] λ max q [mm] λ min q [mm]

SOL width modeling for TCV

In this Section we compare the predictions of the SOL width based on SolEdge2D-EIRENE simulations with κ-ε turbulent transport to the experimental analysis performed with the TCV data. For the latter, we reproduce here the results reported in 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], including the error bars, closed circles on Figure 17. The uncertainty in the simulation data, open circles, is smaller except for the medium divertor height. One can notice that the simulation data does not capture the increasing SOL width with divertor height. The SOL width also appears to be overestimated. A similar behavior with no clear dependence on divertor height is found for the simulation output when analyzing the decay length towards the private flux region λ pf . The magnitude of this decay length is always smaller in the simulations, open squares, than in the experiment, closed head down triangles. However, regarding the latter, one finds a clear trend indicating a reduction of the decay length with increasing divertor height. The closest TCV configuration to that considered for ITER is the so-called short leg configuration. For the latter the scaling laws found in the literature [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 | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] should apply, see Table 3. The safety factor, estimated as q 95 ≈ 6 leads to a comparable SOL width for twice the scaling law [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] in H-mode and for the scaling law in L-mode [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF]. These values match for q 95 ≈ 6.6. The SolEdge2D simulations yield a value in the range 7.8 mm 8.7 mm, therefore in reasonable agreement. However the TCV data for the same configuration yields λ SOL ≈ 4.1 mm ± 1.1 mm, closer to the values expected in H-mode for TCV according to [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF]. Although we have found that using the κ-ε model allows one recovering turbulent transport in the divertor SOL region, one finds that the effect in terms of transport does Also indicated are the predicted SOL width given by the scaling in Ref. [START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] for the H-mode dashed line labeled by H, and twice that value as expected for the L-mode, dotted line labeled by L. Extrapolation to zero X-point height is also indicated, dashed-dot line, for both private flux and SOL characteristic scales.

not govern the experimentally observed increase of the SOL width. Also, one finds that the drop of the heat deposition towards the private flux is typically a factor 2 smaller than the decay reported from the Infra-Red data.

One could argue that a drive for specific turbulence transport in the divertor region is missing in the κ-ε. It would be particularly important in terms of modeling if its growth rate does not exhibit the ballooning character of interchange Eq.( 17). However, the fact that the experimental decay length towards the private flux region is found to decrease with X-point height does not precisely fit with such a possibility. When extrapolating the experimental data to a flush divertor configuration, with the X-point on the divertor floor, one finds that the SOL width and the decay length towards the private flux region are approximately the same 2.75 mm. While the SOL width behavior is found to scale more or less according to the geometrical feature of the various configurations, the result for the private flux region is altogether unexpected, in particular one would expect zero transport for the flush divertor where the private flux region disappears altogether. It appears therefore that this behavior does not support the idea that a transport model without the ballooning character of interchange must be included. An alternative would be that the large scale flows in the vicinity of the X-point, reported in numerical simulations of micro-turbulence [START_REF] Galassi | Drive of parallel flows by turbulence and large-scale e×b transverse transport in divertor geometry[END_REF] vary with X-point height. Such physics can be addressed with SolEdge2D-EIRENE [START_REF] Bufferand | Implementation of drift velocities and currents in soledge2d-eirene[END_REF]. It is to be underlined that a few free control parameters exist in the κ-ε model. All but C ν Eq.( 14) have been set to one in the present work. They offer a possibility of finer tuning according to experimental evidence as has been done in neutral fluids. This optimization step could modify the relative importance of SOL turbulent transport between main chamber divertor regions. It is to be reminded that in the spirit of the neutral fluid use of the κ-ε model, the tuning of the free parameters is performed with reference experiments and then kept as universal input for all simulations. We aim at the same philosophy for plasmas and step out of tuning transport coefficients at need. The considerable uncertainty of the scaling law for L-mode SOL width [START_REF] Scarabosio | Outer target heat fluxes and power decay length scaling in l-mode plasmas at jet and aug[END_REF], leads to λ max q /λ q = 6.1 and λ q /λ min q = 7.4 , see Table 3. One thus finds that the κ-ε modeling implemented in SolEdge2D-EIRENE has a better prediction capability with 0.96 ≤ λ SE2D SOL /λ T CV SOL ≤ 1.7, the superscript T CV standing for the experimental results and the superscript SE2D for the SolEdge2D simulation results.

Confrontation to WEST experiment

We present now a confrontation to the 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 differences with the Table 4: 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 LH, 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

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]
3.4 1.25 5.0 0.132 TCV case are additional heating, higher toroidal magnetic field and larger plasma size. The aim of this comparison is to complete the predictive capability for the divertor deposition pattern in a different operating regime in a configuration with reduced X-point height. The plasma parameters for this WEST discharge are presented in Table 4. In order to set the input power for the SolEdge2D-EIRENE simulation we must determine the power P in used as boundary condition for the simulation. With P LH = 4 M W of Lower Hybrid power injected by two antennas in a low density plasma, the ohmic heating becomes small enough to be neglected. A fraction of the injected power by the LH system is not coupled to the plasma due to various loss mechanisms, such as ripples losses. The core radiation measured by the WEST bolometer system amounts to P core rad = 2 M W , while the edge radiation including the outer plasma and divertor region amount to P edge rad = 0.6 M W . As input for the simulation we consider P in = 1 M W having in mind that the present set-up to test the κ-ε model does not handle impurities and their radiation. The density value at the inner boundary of the SolEdge2D-EIRENE simulation domain is set at, n BC = 2.5 10 19 [m -3 ], in agreement with the expected density value in the WEST edge plasma for this kind of discharge and consistent with the density that can be extrapolated from the reciprocating probe. The profiles for the density n e and the electron thermal energy T e obtained

Midplane profiles

The density and electron thermal energy profiles obtained with the SolEdge2D-EIRENE simulation are compared to the experimental profiles measured with the reciprocating Langmuir probes (RLP) located at the top of the device with vertical plunge. These profiles are mapped at the outer midplane for comparison, Figure 19. 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. 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 cross-field transport in the whole simulation domain the overall agreement is remarkable.

Divertor profiles

As analyzed with Figure 19, the SolEdge2D-EIRENE simulation of the shot WEST-55049 is in the hot divertor regime with a small thermal energy variation between the midplane and the divertor, and similarly for the density. However, in the experiment 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 a change of behavior in the parallel transport. 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 that the actual values. On Figure 20 and Figure 21 are plotted 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 20, 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 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. 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 21 right hand side. In practice one finds a ratio of the order of 2.2, therefore reasonably close. One also notices that the profiles shapes do not compare well. Differences in the density n e profiles, Table 6: 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]. 20 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 6 yield quite comparable values 8 ± 0.8 [mm] in the SolEdge2D-EIRENE simulations and 8.2 ± 0.9 [mm]. 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]. As for TCV, 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, both TCV and WEST experiments analyzed here are 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 4.4, 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

The plethora of plasma parameters and physical phenomena as well as the possibility of various geometries of the magnetic equilibrium and consequently its relative position to the wall will require a chain of models to prepare control and analyze the ITER experiments. On the one hand, physics based reduced models will be used to explore routinely the plasma scenarios at modest computational cost and on the other hand first principle simulations based on a couple of fundamental equations will model at large computational cost some facets of the intricate non-linear physics at play. The 2D transport codes based on averaged fluid equations and used to investigate plasma-wall interaction fulfill requirements of the reduced models. However, their predictive capability is still acknowledged by the community to be insufficient and hindered by the large number of free parameters compared to the relatively reduced set of measurements. In that respect a major challenge for such models is improving the way turbulent transport of heat and particles is accounted for. The latter governs transport across the magnetic flux surfaces, and, consequently, from the core plasma into the edge region, generating the plasma boundary layer in interaction with the wall components in both main chamber and divertor.

We present here an approach to self-consistently estimate the cross-field fluxes in the edge and SOL regions of diverted plasmas with a model that can be handled in a transport code such as the SolEdge2D-EIRENE suite of codes. This approach is inspired by the family of two-equations models developed in Reynolds Averaged Navier-Stokes (RANS) for neutral fluids, in particular the so-called κ-ε model. In such models turbulent diffusive transport is assumed and the transport coefficients are computed at each time and each position, according to the values of the two additional fields. In the reference neutral fluid models, a couple of free parameters have been tuned to ensure the best match with experimental evidence. They are then considered as fixed once for all and kept identical in all simulations. Our starting point in this paper has been a slightly generalized κ-ε framework.

Consequently, 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, the κ-ε concept is to express the effective particle diffusion, viscosity and heat conductivity as functions of the two fields κ-ε governed by two additional transport equations algebraically derived. The turbulent kinetic energy κ has been chosen completed by the field ε. These equations are quite generic and are comparable to key theoretical models such as the Ginzburg-Landau amplitude equations [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 plasma transport reservoir models [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] and predator-prey models [START_REF] Murray | Mathematical Biology: I. An Introduction[END_REF][START_REF] Diamond | Selfregulating shear flow turbulence: A paradigm for the l to h transition[END_REF][START_REF] Floriani | Selfregulation of turbulence bursts and transport barriers[END_REF]. Both predator-prey 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 are shown to provide a theoretical background to the κ-ε description of turbulent transport. For these aspects we concentrate on the local evolution. We then find that the turbulent kinetic energy is a prey and ε the predator. As in the neutral fluid model, ε is the turbulence dissipation rate. As a consequence κ and ε are both defined to be positive. One can then define a diffusion coefficient D = κ 2 /ε and a characteristic time τ = κ/ε. We have bridged the κ-ε model to the quasilinear theory where the diffusion coefficient is shown to depend on κ the kinetic energy of the fluctuating velocity and a characteristic time usually associated to the spectral frequency width. While in the quasilinear theory the system is closed by assuming κ and τ to be defined by other means, the κ-ε model provides equations for a self consistent evolution of these fields as with predator-prey models. 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 positive growth rates, therefore a drive by a given instability. In the present work we have considered 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. In the standard neutral fluid κ-ε model this quadratic saturation term is of the form ∝ ε 2 /κ. We show that any saturation term of the form ∝ ε 2 /κ d , with an exponent d ≥ 1, leads to comparable behavior. Our choice allows us identifying ε as the dissipation channel via the Kolmogorov energy cascade to small scales, thus providing a clear picture for the meaning of the ε field. This choice departs from that usually made for fusion plasma predator-prey systems where the main predator is associated to zonal flows [START_REF] Busse | Geofluiddynamical Wave Mathematics: Research Contributions[END_REF][START_REF] Diamond | Zonal flows in plasma-a review[END_REF], hence eddy shearing by large scale flows. We also show that the two parameters D ω and V together with the chosen growth rate define a Kubo number V 2 /(γ I D ω ) governing a transition from weak

D ≈ D V = V 2 /γ I to strong turbulence D ≈ √ D V D ω .
Given the growth rate γ I and setting D ω 1, the remaining free parameter of the model is the velocity V and describes transport in the weak turbulence limit. Our analysis shows that with very few changes, one can propose a theoretical background to the κ-ε model. We have also underlined that introducing the ε field is in fact a way to define a characteristic time. From our analysis we can conjecture that the key point for the success of the κ-ε model is that the ε field is very similar to the κ field, so that one can build the ε evolution equation according to that for κ but for the local dynamics.

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]. Normalizing the SOL width to the plasma minor radius a, we have therefore assumed the scaling relationship to be λ SOL /a ∝ qρ * , q being the relevant safety factor and ρ * the ratio of a characteristic Larmor radius and plasma minor radius. Given the SOL characteristic time, namely the parallel connection time τ ∝ qAa/c s we have then assumed λ 2 SOL = D V τ , which defines V = λ SOL /τ √ γ I τ . In these expressions A is the aspect ratio and c s is the sound velocity, namely the SOL characteristic velocity governed by the Bohm condition at the sheath [START_REF] Stangeby | The plasma boundary of magnetic fusion devices / Peter C. Stangeby[END_REF]. We have then obtained V /c s = ρ * √ γ I τ . We have found therefore that the velocity V scales like a drift velocity V ∝ ρ * c s with extra dependence on A and q. This result therefore appears to be rather generic of plasma turbulence. It is used in the simulations presented in this paper, both for 1D and SolEdge2D-EIRENE modeling.

Results of the 1D model are confronted to experimental evidence by analyzing the computed SOL width and comparing the result 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]. Given the chosen closure for V , one could expect to achieve perfect agreement. This is not the case. This has led us to analyze the scaling law and several shortfalls have been identified, the most striking being that its dimensionless form departs significantly from the reference qρ * dependence and that in fact ii 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. We have then found a less intuitive picture of V that departs in particular from the straightforward scaling properties of the drift velocity. The latter features are reminiscent of the results for the SOL width scaling determined from a data base analysis of micro-turbulence simulations [START_REF] Fedorczak | Width of turbulent sol in circular plasmas: A theoretical model validated on experiments in tore supra tokamak[END_REF][START_REF] Fedorczak | On the dynamics of blobs in scrape-off layer plasma: Model validation from two-dimensional simulations and experiments in tore supra[END_REF]. Since λ SOL /a ∝ qρ * appears to be a robust result, the change in V must therefore be associated to a change in the characteristic time, the latter then appearing as a key issue in understanding and modeling turbulent transport. 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 the V exhibits a monotonic decrease from the low shear limit to the large shear limit, the ratio between these values being ≈ √ 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 | Selfregulating shear flow turbulence: A paradigm for the l to h transition[END_REF][START_REF] Kim | Zonal flows and transient dynamics of the l → h transition[END_REF][START_REF] Miki | Spatio-temporal evolution of the l → i → h transition[END_REF][START_REF] Floriani | Selfregulation of turbulence bursts and transport barriers[END_REF]. Furthermore, the model has been devised to exhibit the appropriate trend, namely a decrease 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 self-consistent manner.

Stepping to the 2D plasma-wall simulations in L-mode, we find that the κ-ε 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 both TCV and WEST are 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 scrapeoff-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 L-mode simulations of TCV and WEST. Further adjustment of the model can be achieved using experiments and/or micro-turbulence simulations. It is to be noted also that the model can be used to model the core plasma as well as self consistent transitions to improved confinement regimes. As also analyzed in the paper, 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. The model can be completed using other instability growth rates and one could also split the model according to fluctuations on either ion or electron scales.

In the present paper we have thus developed the theoretical background of the κ-ε model. Using this reduced model to investigate turbulent transport we have shown its relevance and the possibility for complete transport modeling, including the core plasma and boundary layer up to the wall, as well as the possibility to handle non-linear self-organization of plasma transport. It thus appears that the κ-ε framework increases quite significantly the predictive capability of the class of RANS-like transport models for plasmas. spired 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. 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 programme 2019-2020 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

A The SolEdge2D model

The model, which is presented in detail in [START_REF] Bufferand | Numerical modelling for divertor design of the west device with a focus on plasma-wall interactions[END_REF], includes the following equations describing the conservation of plasma density, momentum and energy: differential operators ∇ = b • ∇ and ∇ ⊥ = ∇ -b • ∇ to define gradients along the parallel and perpendicular direction . 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 + ∇ • nvb + nv ⊥ = S n (37a) ∂(nv) ∂t + ∇ • nv vb + v ⊥ = -∇ nT i m i + q i nE m i + R ei m i + ∇ • (νn∇ ⊥ v) + S nv (37b) ∂ ∂t 3 2 nT i + 1 2 m i nv 2 + ∇ • 5 2 nT i + 1 2 m i nv 2 vb + v ⊥ = ∇ • κ i ∇ T i b + χ i n∇ ⊥ T i + νn∇ ⊥ 1 2 m i v 2 + q i nvE + R ei v + Q ei + S Ei (37c) ∂ ∂t

B Adimensional scaling laws 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:

τ E ∝ I α I p B α B T n αn P -αp M α M R α R A α A k α K ( 38
)
where 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, 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 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. This generates a constraint, the so-called Kadomtsev similarity constraint [START_REF]Chapter 2: Plasma confinement and transport[END_REF].

β = nT B 2 β 0 ; β 0 = n 0 T 0 B 2 0 /(2µ 0 ) (39b)
ρ * = T 1/2 aB ρ * 0 ; ρ * 0 = m p T 0 eB 0 a 0 (39c) q = F q (k, δ)q cyl ; q cyl = aB T RB pol ; B pol = µ 0 2π

I p akF I (k, δ) (39d)
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 . The various constants that appear in Eq.( 41) 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 Eq.( 39d) 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 (42a)
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 (42b)
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 (43a)
Ωτ E ∝ k a k q aq M a M A a A β a β ρ aρ * ν aν * (43b)

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 ) (44a)
With the coefficients of the ITER1996th L-mode scaling, Table 7 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.

[?]. The three exponents a β , a ν and a ρ characterize the plasma:

a β = (α I + 3α n -8α p + α R ) 5(1 -α p ) (44b) a ν = (-α I + 2α n + 3α p -α R ) 5(1 -α p ) (44c)
a ρ = (-6α I + 2α n + 18α p -6α R ) 5(1 -α p ) (44d) while the magnetic equilibrium properties are characterized by:

a k = (α I + α K -α p ) 1 -α p (44e)
a q = (α R -4α I -2α n -3α p ) 5(1 -α p ) (44f)

a A = (5α A + 6α R -4α I -2α n -8α p ) 5(1 -α p ) -s ν a ν (44g)
and the main ion species effective mass dependence is:

a M = α M (1 -α p ) (44h)
In the ITER physics basis, 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.

B.2 Adimensional scaling of the empirical SOL width

Let us now consider the scaling law proposed for λ q :

λ q = λ * ,0 B u B T q uq P up R u R (45a)

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 the 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 (45b)

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

λ * ,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 8: 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

(46a)

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 (46b)

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

λ q a ∝ k k q q M M A A β β ρ ρ * ν ν * B B ( 46c
)
where the exponents are given by: 5r q + 3r p + r a 5 -r p a q (47f) A = 5r A + 8r p + r a 5 + s ν 3r p + r a 5 -r p a A (47g)

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 * [START_REF] Vedenov | Quasi-linear theory of a plasma[END_REF] This result departs rather strongly from that assumed to hold qρ * and used to determine the parameter V of the κ-varepsilon model Eq.( 20c). 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.

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 (49a)
The coefficients v of the scaling law obtained are given in Table ??. With

v q v B v p v R
0.844 0.277 0.12 -1.06

Table 9: Exponents of the engineering scaling law obtained with the 1-D κ-ε model.

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 (49b)
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 mq M M A A β β ρ ρ * ν ν * B B ( 50
)
where the exponents are given by Eq. [START_REF] Norscini | Self-organized turbulent transport in fusion plasmas[END_REF]. Using the coefficients as given in Table 9, one finds again that the scaling law is not homogeneous leading ot 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 ( 51
)
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 an impact on the turbulent transport. Further simulations are required, with a more appropriate particle source control, eventually addressing the particle confinement time, to have a complete analysis.

B.4 Dimensional scaling of the SOL width qρ * Assuming that the SOL width λ q /a is of the form:

λ q a ∝ q q ρ ρ * A A (52)
one can step backwards to determine the expected scaling law in terms of the engineering parameters. Given ρ *

ρ * ∝ T 1/2 B T a (53)
and:

T ∝ P τ E nV (54)
so that:

ρ * ∝ P τ E 1/2 na 3 Ak 1/2 B T a (55)
One can also use q = aB T /(AI p ) to obtain:

λ q a ∝ aB T AI p q P τ E 1/2 na 3 Ak 1/2 B T a ρ A A (56)
Given the expression of τ E Eq. [START_REF] Grisolia | Plasma wall particle balance in tore supra[END_REF], 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 (57a)

Figure 1 :

 1 Figure 1: Variation of κ and ε with the control parameter γ, for different values of the exponent η, the reference values η = 1 black line head down closed black triangles and η = 3/2 black line closed black circles, and for η = 1.25 blue line closed head up blue triangles, η = 2 blue line closed blue triangles and η = 0.5 dashed blue line open blue circles. Left hand side: variation of κ. Right hand side: variation of ε.

Figure 2 :

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

Figure

  Figure3: 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.

Figure 4 :

 4 Figure 4: 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. hand side, is peaked close to the separatrix and decays rapidly towards the plasma core. The ionization source in the SOL is localized in the vicinity of the separatrix.

Figure 5 :

 5 Figure 5: Radial distribution of the energy source VS E , where V is the plasma volume and S E the source of energy density in M W/m 3 .

Figure 6 :

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

Figure 7 :Figure 8 :

 78 Figure 7: 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.

  Figures 9 and 10 are plotted the variation of λ p from the present κ-ε 1-D model -black curves open squares-and power law variations: the experimental scaling law -dash dot black curve, and a fit of the results of the 1-D model, dashed blue curve. 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, 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 ].

Figure 10 :

 10 Figure 10: 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 ].

  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 (34a) 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.

Figure 11 :

 11 Figure 11: 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 Eq.( 34a) yields the black dashed line. The κ-ε scaling stemming from the empirical scaling law of the SOL width Eq.(33b) is plotted with dash-dot black lines while the theory based κ-ε models Eq.( 35a) 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 12 :

 12 Figure 12: 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.

Figure 14 :

 14 Figure 14: 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.[85]. Panel a: standard transport modeling in 2D plasma-wall interaction. Panel b: modified choice of transport modeling enforcing ballooned transport.

Figure 17 :

 17 Figure 17: Model scaling (left panel) and experimental scaling (right panel) of λ SOL against the outer divertor leg length L div normalized by the plasma vertical size ka where k and a are the plasma elongation and radius, TCV experimental data closed circles and error bars from [85], SolEdge2D-EIRENE simulations open circles. The private flux decay length are also shown: TCV data closed head down triangles, simulation results open squares.Also indicated are the predicted SOL width given by the scaling in Ref.[START_REF] Eich | Scaling of the tokamak near the scrape-off layer h-mode power width and implications for iter[END_REF] for the H-mode dashed line labeled by H, and twice that value as expected for the L-mode, dotted line labeled by L. Extrapolation to zero X-point height is also indicated, dashed-dot line, for both private flux and SOL characteristic scales.

Table 5 :Figure 18 :

 518 Figure 18: 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 19 :

 19 Figure 19: Radial profiles for the shot WEST-55049 at the outer midplane, SolEdge2D-EIRENE simulation black curve closed circles, reciprocating probe data blue dashed line and open triangles. Left hand side: electron thermal energy. Right hand side: density.

Figure 20 :

 20 Figure 20: 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. Let hand side: ion saturation current. Right hand side: parallel heat flux.

Figure 21 :

 21 Figure 21: 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.

=

  ∇ • κ e ∇ T e b + χ e n∇ ⊥ T e -envE -R ei v -Q ei + S Ee (37d) n, T and v are the moments of the averaged kinetic distribution function. b ≡ B/B is the versor aligned to the magnetic field, v ⊥ = v -vb is the velocity component perpendicular to the magnetic flux surfaces, 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 entrance of the pre-sheath, 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 se, ,i | ≥ c s = qe(Te+T i ) m i

  ( 40c), one thus obtains:n = n 0 B 8/5 ρ * ν * /A sν qA

B = 5r B + 3r p + r a 5 --3r p -r a 5 -r p a ν (47c) ρ = -18r p -6r a 5 -

 555 r p a B (47a)β = 8r p + r a 5 -r p a β (47b) ν = r p a ρ (47d) k = -r p a k (47e) q =

Table 7 :

 7 τ 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 Exponents of the engineering scaling lax of the energy confinement time of the IT ER1996 L-mode scaling.

	α I	α B	α n	α p	α M	α R	α A	α k
	0.96 0.03	0.4 0.73	0.2	1.83 0.06 0.64

In the framework of SOL physics, it is convenient to use the sound velocity c s as reference thermal velocity rather than a thermal velocity.

M = -r p a M (47h)

Acknowledgements

The authors thanks Prof. B.E. Launder for fruitfull discussions on the orignal k-epsilon model which inspired this work. This work has also been in-

which then yields the scaling law of the SOL width given in terms of engineering parameters:

Let us consider the case ρ = q = 1, A = 0:

For the L-mode confinement scaling, table 7, this leads to:

Since B pol ∝ I p /a and A = R/a one can rewrite this result to obtain an expression comparable to Eq.( 49a).

These results can be compared to those obtained with the scans performed with the 1-D κ-ε model, see Table 10. The agreement is rather poor, but at least there is no disagreement regarding the sign of the exponents. As

κ-model 0.844 0.277 0.12 -1.06 expected scaling 0.52 0.015 0.135 -0.555

Table 10: Exponents of the engineering scaling law obtained with the 1-D κ-ε model compared to its expected scaling.

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.

Bibliography