One-dimensional modelling of limit-cycle oscillation and H-mode power scaling

To understand the connection between the dynamics of microscopic turbulence and the macroscale power scaling in the L–I–H transition in magnetically conﬁned plasmas, a new time-dependent, one-dimensional (in radius) model has been developed. The model investigates the radial force balance equation at the edge region of the plasma and applies the quenching effect of turbulence via the E × B ﬂow shear rate exceeding the shear suppression threshold. By slightly ramping up the heating power, the spatio-temporal evolution of turbulence intensity, density and pressure proﬁles, poloidal ﬂow and E × B ﬂow self-consistently displays the L–H transition with an intermediate phase (I-phase) characterized by limit-cycle oscillations. Since the poloidal ﬂow is partially damped to the neoclassical ﬂow in the edge region, the numerical results reveal two different oscillation relationships between the E × B ﬂow and the turbulence intensity depending on which oscillation of the diamagnetic ﬂow or poloidal ﬂow is dominant. Speciﬁcally, by including the effects of boundary conditions of density and temperature, the model results in a linear dependence of the H-mode access power on the density and magnetic ﬁeld. These results imply that the microscopic turbulence dynamics and the macroscale power scaling for the L–H transition are strongly connected.


Introduction
It is well known that turbulent transport on the edge of the plasma column strongly influences the confinement of particles and energy of a fusion device [1,2].It is also well known that the typical transition of the low confinement mode (L-mode) to the high confinement mode (H-mode) is associated with the formation of an edge transport barrier (ETB) [3,4].After decades of H-mode studies [5], there are two important phenomena that should be pointed out from the experiments.Firstly, that the L-H transition requires sufficient heating power, known as the L-H transition power threshold, and secondly, if the input heating power is close to the transition threshold, an intermediate phase, the so-called I-phase or limit-cycle oscillation (LCO), appears prior to the final transition into the H-mode.The former is related to the macroscale transport physics for the H-mode power scaling and the latter corresponds to the microscopic turbulence dynamics physics for the triggering of the L-H transition.Based on the experimental observations, it is challenging and necessary for a dynamical model to capture the basic behaviour of the L-I-H transition and predict the threshold power.
Based on experimental I-phase studies one can generally define two types of LCOs.The temporal dynamics of the turbulence-flow interaction was studied in several devices such as NSTX [6], TJ-II [7,8], AUG [9] and EAST [10], where experimental evidence was found for supporting the predatorprey relationship [11] between turbulence and flows.However, more recently, probe measurements in HL-2A [12] show the existence of another kind of LCO dominated by pressuregradient-induced diamagnetic drift, which shows opposite time sequences with respect to the predator-prey model.In addition, it is noted that the time sequences in the DIII-D experiments [13] are consistent with this opposite LCO, where one first observes a growth in the radial electric field causing a reduction in fluctuations.In parallel, significant effort has been devoted to the study of the condition and parameter dependence for the L-H transition power threshold [14,15].In general, this power threshold tends to have different behaviour in lower and higher density regions [16], but appears to increase with higher density, magnetic field and the size of the tokamak.Therefore, it is aimed at bridging the gap between microscopic turbulence dynamics physics and the macroscale power scaling in the L-I-H transition in close compliance with these relevant experiment results.
In the experiments described above, it is confirmed that the appearance of the ETB of the H-mode is both robust and ubiquitous, which naturally indicates there could be a universal mechanism leading to a strong reduction in turbulence and the formation of a transport barrier.To understand the mechanism, bifurcation theory models based on the suppression of turbulent transport by sheared flows successfully describe the time evolution of key variables that characterizes the L-H transition in different regimes [11,17,18], even including the radial profile evolution [19,20].Particularly, in order to understand the connection between microscopic turbulence dynamics physics and the macroscale power scaling in the L-I-H transition, the radial force balance equation is specially investigated to apply the decorrelation theory of E × B flow shear suppressing turbulence transport [21,22] driven by the pressure gradient.In order to make the connection between the flow shear intensity and turbulence intensity self-consistent and the whole dynamic system self-contained, the new feature of this one-dimensional model integrates the generally accepted understandings.(1) The heat transport equation is derived from the Braginskii energy conservation equation coupled with the particle conservation equation to determine the pressure gradient profile where the turbulence intensity is primary to determine the transport, which dominates the depth of the radial electric field well at the edge.(2) The poloidal momentum balance equation consists of the fluctuationdriven Reynolds stress, the collisional bulk plasma diffusion, neoclassical poloidal flow damping and the ion-loss cone near the separatrix, which is related to provide a negative feedback control for LCO during the I-phase.(3) The Reynolds stress is enhanced by increasing turbulence assisting the diamagnetic flow to lead the E × B flow shear intensity to prematurely achieve the critical point of the microscopic turbulent quenching, thus triggering the L-I transition.Finally, a novel time-dependent one-dimensional (in radius) dynamical model is composed of these five coupled macro fields to investigate the L-I-H transition and H-mode power scaling.
This paper is organized as follows.In section 2, the reduced five-field macroscale model related physics is introduced and the boundary conditions are discussed.In section 3, the numerical investigations of the L-H transition and the H-mode power threshold are presented.Section 4 contains our conclusions and remaining issues are summarized.

The criterion for the L-H transition
In order to investigate the dynamics of the L-H transition, we assume that if the E × B flow shear rate γ E = (r/q)∂(qV E /r)/∂r exceeds the shear suppression threshold it will trigger the quenching of turbulence [22], which leads to the transport bifurcation.Here it is assumed that the turbulence is driven by ion pressure gradient.Therefore, the shear suppression threshold would be near the maximum of turbulence growth rate, which is proportional to the pressure gradient and scaling as γ E,c ≈ γ max ∝ |∇p| α , where 0 < α < 1.This will be expected for ion temperature gradient driven mode [23,24], and specifically α = 0.5 is adopted in the model agreeing with the interchange mode.
In order to evaluate the E × B flow shear rate, the radial force balance equation is adopted: which indicates the connection between the radial electric field E r , the pressure gradient ∇p, particle density n, magnetic field B and the perpendicular flow V ⊥ .Particularly, this model focuses on the plasma edge physics of the L-H transition, which implies that the ∇p and poloidal flow V θ dominate E r [25].Then the approximation of the E × B flow is given by equation ( 2), where the sign convention of the model is that the positive value is in the ion diamagnetic drift direction and the negative value is in the electron diamagnetic drift direction.
It is obvious that equation ( 2) requires access to the cross-field heat, particle transport and cross-field poloidal momentum transport.For the transport coefficients of heat and particle, the model self-consistently presents the evolution of turbulence intensity I to characterize the dominating turbulence transport, Equation ( 3) directly integrates the properties of the reduced model [1,19], which describes the nonlinear damping by self-saturation of turbulence due to local broadening of the spectrum given by and the nonlocal turbulence spreading via the turbulent diffusivity κ [26].
Furthermore, this model emphasizes that the pressure gradient driving term with an equivalent growth rateγ , and the response of γ to the E × B flow shear rate, leads to turbulence bifurcation corresponding to the quenching effect of turbulence.
where γ E and γ E,c indicate averaging over the turbulence suppression region, and γ 0 is a certain constant to adjust the driving term γ ∇p to the order of turbulent decorrelation time scale [23].In the model, the width of the suppression region is determined by the radial correlation length of a turbulent eddy in the edge region (shear layer) as some multiple of the mesoscale given by L corr ⊥ ∝ L p ρ i , and L p = |∇p/p| −1 is the mean pressure gradient scale length in the edge region, ρ i is the ion gyroradius evaluated at the plasma thermal velocity [27].

The heat and particle radial transport
The following task is to obtain the diamagnetic flow ∇p/ZeBn from the global profiles of p = nT and n, where T is the ion temperature.Generally, the evolutions of the pressure and density profiles are given by the following one-dimensional reduced transport equations with external sources: where the external heat source profile with the deposition width L Q and the amplitude Q 0 are located in the core and given by a Gaussian function The particle source flux profile peaks at the edge and decreases exponentially with the ionization depth L s , given by Here a is the minor radius of the device, 2T /m i is the ion thermal velocity, and σ v ion is the ionization rate (the cross section averaged over a Maxwellian distribution) which is a strong function of the temperature T , given by [28,29] where φ T is the ionization potential (in deuterium φ T is 13.6 eV).Finally, the ionization depth is estimated as The heat flux p and particle flux n are expressed as functions of the ion temperature gradient and density gradient, respectively, and also the particle convection (pinch) is considered: where χ 0 I and D 0 I are proportional to the turbulence intensity corresponding to turbulent thermal and particle diffusivities, respectively.χ neo and D neo correspond to neoclassical transport coefficients which depend on different ion-ion collision frequency regimes in a tokamak [30]: Here ω t = V th /(qR) is the transit frequency, and is the inverse aspect ratio.The radial profile of safety factor q usually has its minimum value q 0 at, or close to, the magnetic axis and increases outwards to the edge value q a .A simple cylindrical model is adopted here to give a typical profile q = q a (r/a) 2 /(1 − (1 − (r/a) 2 ) qa /q0 ), and the ion-ion collision frequency is given by [31] Finally, in order to support the transport at the L-mode dominated by turbulent processes, reasonable ratios between the coefficients are to be compared with typical neoclassical values [31].
Furthermore, to consider the peaking of the density profile, the particle convection velocity is given by [32,33] Here the particle convection is mediated by the anomalous pinch including two mechanisms, where the first term indicates a velocity proportional to the curvature of the magnetic field with the turbulent equipartition (TEP) coefficient C q [34,35], the second term indicates a velocity proportional to the inverse scale length of temperature gradient with the thermo-diffusion coefficient C T [36,37].In this model, the safety factor q profile is fixed but the temperature gradient is self-consistent with the heat transport.In order to get a reasonable profile of particle pinch velocity, it is suggested to calibrate the ratios between these two coefficients by comparing with the observation in [38].

The poloidal flow equation
The other key factor to change the E ×B flow in equation ( 2) is the poloidal flow component, which is governed by the poloidal momentum balance equation coupling various processes [23,39].In the model, it is focused on terms for ion and neoclassical poloidal flow damping as follows: where bV is from the collisional bulk plasma viscosity with the diffusivity [40] proportional to the ion Larmor radius and collision frequency D p ∝ ρ 2 i ν ii ∝ T ν ii /B 2 .And V ∇V originates from the turbulent Reynolds stress [41] with an equivalent negative viscosity [42] proportional to the turbulence intensity D rr ∝ I/B 2 , which brings an enhancement of momentum flux across the edge during the increase in turbulence intensity.
Meanwhile lc is from the ion-loss cone where the edge region is smaller than the poloidal gyroradius attributed to the direct ion loss.In response, the poloidal rotation will be driven by the torque associated with the ion-orbit loss [43,44] where 2qV E /εV th .Specifically considering the ion-loss cone located very near the edge, the ion loss flux has a deposition width λ b ∝ qε −0.5 ρ i proportional to the width of the banana orbit of trapped particles given by lc = lc 0 exp The last term of equation ( 17) is the contribution of the neoclassical poloidal flow damping, which means that the induced poloidal rotation will decay towards the neoclassical value within a decay time following the general expression [19,45] neo where the neoclassical poloidal velocity V neo θ = −K p ∂ r T /ZeB has been derived from neoclassical theory [46,47].However, the exact values of damping rate γ damp and the coefficient K p depend on the collisional transport regime.In the model an approximate formula [48] is adopted to evaluate K p , which is positive in the rare collisional regime and negative in both the plateau and collisional regimes.On the other hand, the damping rate γ damp is a piecewise function of the ion collision frequency, which is proportional to ν ii in the rare collisional regime and is controlled by the transit frequency ω t in the plateau regime, and finally is inversely proportional to ν ii in the collisional regime, given by [1]

Boundary conditions and parameters
So far, the four coupled fields (including turbulence intensity I , density n, temperature T and poloidal velocity V θ ) are obtained to identify the critical point of the microscopic turbulent quenching in which the E × B flow shear is governed by the radial force balance equation.In order to numerically solve this system of equations, the model is supposed to achieve first the L-mode equilibrium from arbitrary initial conditions; therefore, the boundary conditions become crucial to affect the spatio-temporal evolution of the pedestal through the L-H transition and the scaling of the H-mode access power with magnetic field and plasma parameters.The simulation region is a one-dimensional space normalized by the minor radius a, so r = 0 and r = a correspond to the centre and the edge of the plasma.Due to the cylindrical geometry, it is natural that the particle and the heat fluxes are zero at the centre of the plasma where the poloidal velocity will be zero as well.For the turbulence intensity, the natural boundary conditions are applied both at the core and the edge.These boundary conditions are given as follows: However, boundary conditions at the edge are more complex especially when the plasma parameters depend on the edge values.In order to be consistent with experimental results, scaling of edge density and temperature based on local and multi-machine databases [49,50] are referenced in the model instead of explicit modelling of the scrape-off layer (SOL) region.It is indicated that the separatrix density n sep is a fraction of the line-averaged density n and is also related to the toroidal magnetic field B, so the edge density condition has the form: where B 0 = 2 T is the reference value, and the scale factor c of magnetic field will affect the dependence of the H-mode threshold power on the magnetic field.We should point to, without directly dealing with SOL region physics, the fact that there is also a so-called two-point model of SOL transport to estimate the upstream SOL value as the edge temperature [27,51].This edge boundary condition gives an explicit connection between the edge temperature and input power P in as well as the separatrix density given by where A 1 and A 2 are related to the connection length in the near-separatrix SOL and the parallel energy flow (crosssectional surface) area [27], which means that these parameters may be constants for a certain tokamak plasma configuration.Finally, it is important to note that the boundary condition for the poloidal flow equation at the edge is constrained by the radial force balance equation in the model.The radial electric field, E r , at the last closed flux surface (LCFS) is mostly constant during the L-H transition, see, e.g., the review [5]; therefore, we assume, as a simple case, that the fixed value zero of the E × B flow at the edge constrains the poloidal velocity to vary with the diamagnetic flow, Finally, the boundary conditions are integral to the model.But the edge values of plasma parameters as well as the constant radial electric field at the LCFS are just logical compromising applications to simplify the physics.Actually, the recent DIII-D and EAST tokamak experiments reveal flow and turbulence amplitude oscillations outside the separatrix [13,52], which suggest that open field line physics may play an important role in the oscillation of L-H transition behaviour.Therefore, a self-consistent separatrix SOL-edge model is required to match the edge dynamics; however, this is beyond the capacity of the present model and will be addressed in future work.

Basic properties of physical quantities in L-mode and H-mode
With the structure described above, the model takes EAST tokamak plasma parameters [53] as reference, which is helpful in calibrating the free parameters in the model.Therefore, the initial L-mode equilibrium is reached as the standard case characterized by the line-averaged density of 2 × 10 19 m −3 , the ion temperature of about 1.3 keV in the core and 0.08 keV at the edge with roughly 0.8 MW critical heating power and major radius R 0 of about 1.9 m, minor radius a of about 0.45 m.The safety factor, q, profile is fixed in the model with q increasing from the minimum value q 0 = 1 in the core to q a = 4 at the edge.The corresponding profiles of heat source with an amplitude Q 0 of about 0.89 MW m −3 located in the core and particle source with an amplitude S 0 of about 3.7 × 10 22 m −3 s −1 localized at the edge are shown in figure1, where the power deposition width L Q = 0.2a is fixed, but the ionization depth L s varies with the temperature and the density given by equation (10) with a value of about 0.03 m in this standard case.
Then the heating power starts to ramp up linearly with a small increment and finally the model locks in to the H-mode state identified by the quench of turbulence and formation of the ETB.In the model, the heat and particle diffusivities are both dominated by the turbulence intensity with their corresponding proportion.In order to visualize the difference of transport properties figure 2 shows that in the edge region the particle diffusivity D = D 0 I + D neo is governed by turbulent transport in L-mode but drops a lot in H-mode due to the suppression of turbulence, where D 0 = 2.5 m 2 s −1 is the optimized value and the other ratios of thermal diffusivity and particle convection velocity are given by equations ( 15) and (16).
As a consequence of transport coefficient structures in different confinement regimes, pedestal and steep profile gradients for density, pressure and temperature in the edge region are formed in the H-mode, as shown in figure 3. Also, it is clear that the profiles of density and pressure maintain stiffness in the radial direction except for the edge region, where the turbulence is suppressed.But in contrast, the core temperature becomes lower in H-mode because of the significant increase in thermal conductivity (χ 0 I + χ neo )n in equation (11) for the contribution of density while the heating power only has a small increment during the study.We should emphasize that these profiles are theoretically self-consistent in the model in which the pedestal width is determined by the turbulence suppression region.We point out that the model is intended to study the physics of the pedestal, which involves the variation of width and height of the pedestal with varying plasma parameters.In the model, the pedestal depends on how the turbulence is suppressed.Figure 4 shows profiles of turbulence intensity and the inverse pressure gradient scale length L −1 p = |∇p/p| at L-mode and H-mode, respectively.At L-mode equilibrium, the turbulence intensity has a maximum value at the edge region, where the pressure gradient is significant to support the turbulence driving term γ 0 ∇p.In the model, it is assumed that the width of the suppression region is determined by the radial correlation length of turbulent eddies in the edge region, which is proportional to the mean pressure gradient scale length.And the numerical calculation of theL −1 p -profile at L-mode suggests that the radial width of 0.94 < r/a < 1 is a reliable edge region.This corresponds to about 2.7 cm close to the corresponding values in the DIII-D [13] and JFT-2M [54] tokamak experiments.
When the suppression region is empirically specified, the model always traces the average of the E×B flow shear |∂ r V E | and |∇p| α over this edge region.If the ratio of |∂ r V E |/|∇p| α exceeds a certain constant corresponding to the criterion of γ E > γ E,c , then equation (4) would switch γ to zero just in this region, which is the direct application of triggering the turbulence quench.Naturally without the drive from the pressure gradient at the edge region, the turbulence intensity would be strongly damped by the nonlinear self-saturation term I with a time scale of 0.2 ms in this case; also, the nonlocal turbulence diffusivity κI is important for controlling the degree of turbulence suppression, and it is expected that the minimum value of the total heat diffusivity χ = χ 0 I + χ neo in equation (11) would approach the neoclassical level in the turbulence suppression region.Furthermore, besides the profiles of density and pressure, the poloidal flow is another key component to affect the E ×B flow shear.Governed by equation ( 17), the profiles of the poloidal flow and neoclassical poloidal flow at L-mode and Hmode are shown in figure 5, where the poloidal flow at the edge shifts the velocity from about 3 km s −1 in the L-mode to about 8 km s −1 in the H-mode.The significant effect here is that the poloidal flow is highly damped towards the neoclassical flow almost all over the radial region except is the region constrained to vary with the diamagnetic flow at the edge layer.It is obvious that there is a sharp poloidal flow shear near the edge since there is a big gap between the neoclassical flow and the diamagnetic flow.Above all, as the comprehensive result of the evolution of the radial force balance equation, it is natural from the model to visualize the profile of the E × B flow and its two important components from equation (2); this is shown in figure 6, which indicates that the negative radial electric field well, E r , is formed right inside the edge, associated with the negative contribution of pressure gradient to determine the depth of the well and the positive contribution of poloidal rotation to be consistent with the fixed value of E r at the edge.In the Lmode, the E r -well is relatively low and flat due to the shallow pressure gradient.However, the poloidal flow provides a sharp flow shear component near the edge, which implies that before the positive feedback of the diamagnetic flow becomes dominating, the negative feedback between turbulence and poloidal flow will sustain the oscillatory behaviour of the whole system for a period of time.
As a complementary explanation, based on the plasma parameter profiles above, the model is able to self-consistently identify the different neoclassical transport regimes according to equation (13).Taking the L-mode in the standard case as an example, figure 7 shows that all over the radial region both the ion-ion collision frequency ν ii and ε 3/2 ω t are less than the minimum value of the transit frequency ω t implying that the whole model does not evolve in the collisional regime, but in

Ion
Dia.  the rare collisional or plateau regime; especially near the edge due to the relatively low temperature, ν ii is even larger than ε 3/2 ω t indicating that the related region in radius belongs to the plateau regime.

The LCO during the L-H transition
In addition to the explicit equilibrium of L-mode and H-mode as shown above, the model is aimed at displaying the spatiotemporal dynamics of the L-H transition.With the parameters calibrated by the L-mode equilibrium standard case, at an arbitrary time of L-mode equilibrium we start to slowly ramp up the heating power linearly with 1% increment in 200 ms, which directly increases the pressure gradient especially at the edge.Meanwhile, the turbulence intensity keeps growing to enhance the momentum flux across the edge by the Reynolds stress term in the poloidal flow equation.This will assist the diamagnetic flow to lead the E × B flow shear intensity prematurely achieve the critical point of the microscopic turbulent quenching, thus triggering the L-I transition.During the I-phase the physical quantities of the system show quasiperiodic behaviours with a certain phase relation, i.e. the LCO.When the heating power keeps ramping up, the period of turbulence suppression increases and the bursting time decreases gradually.The system evolves into the final stage where the pressure gradient is high enough to dominate the flow shear intensity, which is identified as an H-mode with strong shear suppression of turbulence and the ETB.For instance, these continuous evolutions of the L-I-H transition are characterized by the turbulence intensity and the E × B flow, shown in figure 8.
Particularly during the I-phase, there are two different oscillation relationships presented in the model corresponding to whether the poloidal flow is damped to the neoclassical flow or not.As the E × B flow profile shows, in figure 6, the E rwell bottom is more close to the edge at L-mode equilibrium.Therefore, radial positions r/a = 0.97 located inside of the corresponding radial electric field increase first, leading to the decrease of the turbulence.In figure 9(b) the limit cycle rotates clockwise, which is similar to the predator-prey model, where the turbulence increases first leading to increase in the turbulent Reynolds stress.
In order to understand the details of the LCO, it is necessary to consider the competition between the poloidal flow and the diamagnetic flow.Figure 10 shows the time evolution of the turbulence intensity, E × B flow, diamagnetic flow and poloidal flow at the same radial locations as in figure 9.It is observed that the oscillation frequency is about 1 kHz and, specifically, that the E × B flow has a phase difference of almostπ between these two radial locations in figure 10(b).In the model, the E × B flow is only affected by the diamagnetic flow and the poloidal flow.From figures 10(b)-(d), it is clear that at the inner position (such as r/a = 0.97) the oscillation amplitude of the poloidal flow is relatively small due to the neoclassical damping, but there is no constraint to the evolution of pressure gradient, then the diamagnetic flow could oscillate freely.Therefore, the oscillation of the E×B flow is dominated by the diamagnetic flow and increases in magnitude when the turbulence intensity gets suppressed.In contrast, at the outside of the E r -well very close to the edge, where the poloidal flow is free from the neoclassical damping and strong in oscillation amplitude, the oscillation of the E × B flow is dominated by the poloidal flow.These observations explain the two different oscillation relationships between the E × B flow and turbulence intensity.Also, it is important that the negative feedback between the poloidal flow and the turbulence intensity originates from an equivalent negative viscosity provided by the Reynolds stress (compare figure 10(a) with figure 10(d)).This negative feedback enables the model to maintain the oscillatory behaviour until the poloidal flow oscillation is not able to reduce the E × B flow shearing rate below the shear suppression threshold.Furthermore, it is a positive feedback between the diamagnetic flow and the turbulence intensity, which implies that the edge transport decreases, the pressure gradient sharpens, enhancing the equilibrium shear and finally the model would access the H-mode with ramping up of the heating power.During the whole I-phase the plasma pressure and its gradient are not only taken as equilibrium parameters but are also oscillating through transport interaction, according to the radial force balance equation, the local phase of E × B flow depends on which of the oscillations diamagnetic flow and poloidal flow is dominant.Figure 11 shows the spatiotemporal evolution of the diamagnetic flow and poloidal flow, which confirms that only very near the edge the oscillation of poloidal flow is significant, otherwise the diamagnetic flow leads the oscillation of the E × B flow.

The dependence of H-mode threshold power on the density and magnetic field
The other main purpose of the model is to investigate the dependence of H-mode threshold power on the density and magnetic field.This is also inspired by the radial force balance equation that makes a clear connection between the E × B flow and the pressure gradient, particle density and magnetic field.Just as the results shown in figure 6, the diamagnetic flow dominates the contribution to the depth of the E ×B flow-well near the edge, so essentially the shear suppression threshold is related to the diamagnetic term.Therefore to analyse the connection between microscopic E×B flow shear suppression threshold and the macroscale power scaling in L-H transition, the E × B flow in equation ( 2) is reduced to the critical approximation in equation ( 27) (below) to emphasize the pressure gradient as affected by heating power and neglect the negative feedback effect of poloidal flow, which qualitatively predicts the density and the magnetic field linear dependence of the transition power threshold.
Based on the parameters of the standard L-mode equilibrium case in section 3.1, where the reference toroidal magnetic field B 0 is about 2 T and the effective charge Z is fixed about 2.5.
During the study of H-mode power threshold, the boundary conditions given by equations ( 24) and ( 25) are applied in the model, which is to parametrize the SOL physics.Also the ionization depth L s for the particle source self-consistently depends on plasma parameters given by equation (10).First with the fixed magnetic field B/B 0 = 1, near the L-H transition, figure 12 shows that the ionization depth is almost inversely proportional to density with a value of about 3 cm in the standard case n 0 = 2 × 10 19 m −3 and edge temperature approaches a saturation at a higher line-averaged density with the fixed factor A 1 = 0.25 and A 2 calibrated by the standard case in the model.With these conditions the investigation of power threshold versus line-averaged density is shown in figure 13, which indicates that the modelling results are quite close to the linear scaling law P L−H ∝ n , especially in the higher density regime.However, the model is not able to achieve the lower density branch with minimum threshold power which is due to the missing physics in the model.We should note that recent investigations in ASDEX Upgrade demonstrated that the ion heat channel plays the key role in the L-H transition [16,55].In standard L-H transition experiments the electrons are dominantly heated during the power ramp-up.Thus, the resulting ramp-up of the ion energy flux across the LCFS will depend on the energy coupling between electrons and ions over the whole plasma volume and thus on the density.The heat transfer from electrons to the ions will decrease with decreasing density, which may explain the power minimum at a certain density.When the power input is directly through the ions as in the present model one should not expect to find a roll-over of the threshold power at a specific density.Towards this study, it is suggested to make the distinction between the electron and ion transport channels to describe the collisional decoupling of the electron and ion transport in the low density regime.
Actually near the L-H transition, the E × B flow profiles corresponding to the various line-averaged densities have their individual structures, which are resolved to show the details in the model.By way of illustration, figure 14(a) shows the variation of the E × B flow profiles with different density regimes by zooming in at the turbulence suppression region, which indicates that the E × B flow tends to deepen the well with the increase in line-averaged density.Meanwhile the statistics shown in figure 14 the E × B flow over the suppression region almost remains constant, besides the increase in the profile peak value.By checking all of the components in equation ( 2), it is suggested that in the higher density regime with higher collisionality the decrease in the factor K p in the neoclassical poloidal velocity, at the rare collisional region, leads the poloidal flow to decrease more with lower velocity in the ion diamagnetic direction, which makes the E × B flow-well deepened.
According to equation (27), it is helpful to understand the constant average of the E × B flow over the suppression region, which implies that T L −1 p would offset the influence of the poloidal flow on the fine structures of the E × B flow for different line-averaged densities.The comparison between the average inverse scale length of temperature and density gradient over the suppression region shown in figure 15 indicates that the temperature gradient scale length is dominant in this case, where even L −1 n increases due to the peaking of the particle source with shorter ionization depth in higher density regimes.But the increase in thermal conductivity would flatten the temperature profile to make L −1 T decrease strongly, which is consistent with higher heating power requirement for the higher density regime.Additionally, there is a density limit, when the heating power just exceeds the threshold power with a small fraction, which means in higher density regimes the finite heating power near the threshold is not enough to keep the E × B flow shear rate exceeding the shear suppression threshold due to the dramatic increase in density after the L-H transition flattening the temperature gradient further.
In order to investigate the L-H power threshold versus the magnetic field, the line-averaged density is specified as n / n 0 = 1, but the density boundary condition equation ( 24) is still applied as well as the same form of the temperature boundary condition.Just before the L-H transition, the ionization depth and edge temperature at different magnetic field scenarios are shown in figure 16, where the change of temperature is dominant in equation (10), due to the fixed lineaveraged density, which implies that the ionization depth is almost linearly depended on temperature.Also the change of heating power becomes dominant in equation ( 25) implying the edge temperature to vary linearly with the magnetic field variation without explicit saturation at higher magnetic fields.Particularly with the magnetic field exponent c = 0.3 in equation ( 24), the model succeeds in reproducing the linear dependence of the L-H transition threshold power on the magnetic field P L−H ∝ B/B 0 shown in figure 17.Besides the upgrade of the model for the physics of the lower density branch with minimum threshold power, the sensitivity of power threshold to the dependence of the edge density on the magnetic field also requires the model to account for more self-consistent boundary conditions by more accurate SOL physics.The response of E × B flow profiles to the different magnetic field values is different from the situation with varying density near the L-H transition, as shown in figure 18(a).The most remarkable effect of the increase in the magnetic field is that the E × B flow-well structure tends to move towards the edge with an almost constant peak value.To understand these fine structures of the E × B flow in the suppression region, it is necessary to consider the influence of the magnetic field on the poloidal flow.According to equations (18) and (19), the effective diffusivity is inversely proportional to B 2 , therefore even if the poloidal flow has a large value at the edge due to the boundary condition, the inward diffusion of poloidal momentum is comparatively weaker at higher magnetic fields.The consequence is that the minimum value of the poloidal flow profile has an outward shift, damped by the neoclassical velocity.
The results shown in figure 18(b) also confirm that the peak value of the E × B flow is independent of the magnetic field; however, the average of the magnitude of the E × B flow over the suppression region decreases with increasing magnetic field due to the decrease in the width of the E × B flow-well.In order to reveal further details, the variations of the average inverse scale length of pressure, temperature and density gradient over the suppression region with different magnetic fields are shown in figure 19, where the decrease in the inverse density gradient scale length becomes dominant due to the broadening of the particle source with longer ionization depth at a higher magnetic field.Since the line-averaged density has been fixed, it is reasonable that the temperature profile gets steepened with increasing heating power, which implies that the collisionality becomes lower and the factor K p increases in the rare collisional regime leading to a higher neoclassical poloidal velocity in the ion diamagnetic direction.In addition to the contribution of poloidal flow with higher velocity in the ion diamagnetic direction, the decrease in the diamagnetic flow due to the increase in the magnetic field and slight decrease in L −1 p , shown in figure 19, make the inside of the E × B flow decrease in the electron diamagnetic direction and the well structure gets narrower.
Generally, equation ( 27) also shows that the E × B flow is inversely proportional to the magnetic field B, which would then require a higher heating power at higher magnetic fields to drive the E × B flow and achieve the shear suppression threshold of the turbulence.Furthermore, it should be pointed out that the diamagnetic flow is the foundation to dominate the radial force balance equation in the model and it corresponds to a similar magnitude of the E×B flow under conditions slightly above the L-H transition.However, the fine structure of the E × B flow-well implies that in the regime of edge turbulence suppression the effect of the poloidal flow is strongly related to the H-mode power scaling, in addition to the negative feedback to the oscillatory behaviour during the L-I-H transition.

Conclusion
We have developed a time-dependent one-dimensional (in radial direction) dynamical model to describe self-consistently the evolution of the turbulence intensity, density and pressure profiles, poloidal flows and E × B flow in a magnetically confined plasma.The central point of the model is the application of the decorrelation theory of sheared E × B flow for suppressing turbulence transport for investigating the connection between microscopic turbulence dynamics and the macroscale power scaling for the L-I-H transition.The model focuses on the plasma edge physics, and the E × B flow is directly evaluated by the reduced radial force balance in which the perpendicular flow is dominated by the poloidal flow in the edge region and the diamagnetic flow is governed by the cross-field heat and particle transport.
The free parameters of the model are calibrated by the standard case of typical EAST tokamak plasma parameters.The boundary conditions at the edge are of key importance for H-mode power scaling and carefully considered to be consistent with experimental results.Instead of explicit modelling of the SOL physics, the scaling of the edge density by equation ( 24) is adopted, which is associated with the lineaveraged density as well as the toroidal magnetic field, and the expression for edge temperature from two-point models of SOL transport shown in equation ( 25) is applied.These are functions of heating power as well as the edge density.Another important issue is to obtain the edge value of the poloidal flow by considering the radial force balance equation as a constraint condition due to the assumption of the radial electric field with a fixed value at the edge which is determined by SOL physics.As a consequence, the diamagnetic flow provides a negative contribution in the electron diamagnetic drift direction and the poloidal flow is a positive contribution in the ion diamagnetic drift direction.
By slowly ramping up the heating power starting from typical L-mode conditions, the model achieved the I-phase characterized by LCOs during the L-H transition.The model also revealed the linear dependence of the H-mode threshold power on the density and magnetic field.The results in this paper are briefly summarized as follows: (a) Two different types of limit-cycle oscillation.The poloidal flow term provides a negative feedback control for limit-cycle oscillation during the I-phase, where there are two key points, the sharp poloidal flow shear near the edge due to the neoclassical flow damping and the effect of turbulent Reynolds stress on viscosity to modify the poloidal flow shear.The oscillations of the poloidal flow maintain the I-phase before the positive feedback of the diamagnetic flow dominates.However, because of the boundary constraint in the model the poloidal flow is partially damped towards the neoclassical flow in the edge region, which represents two different oscillation relationships between the E × B flow and the turbulence intensity.One is that the oscillation of the E × B flow is dominated by the diamagnetic flow and increases in magnitude when the turbulence intensity gets suppressed; the other is that only very near the edge does the turbulence intensity grow first followed by the increment of the localized flow due to the poloidal flow being free of the neoclassical damping and dominating the oscillation.
(b) The linear dependence of the transition power threshold on the density.The diamagnetic term in the radial force balance equation is found to dominate the threshold power for accessing the H-mode.But the pressure gradient at the edge region also depends on the edge value of density and temperature, therefore by adopting the appropriate boundary conditions the model achieves the linear dependence of the L-H transition threshold power on the line-averaged density with a fixed magnetic field.
Particularly, the E × B flow profile has individual fine structures in the turbulence suppression region due to the poloidal flow damped towards different neoclassical flows in different density regimes.This also results in that the average of the E × B flow over the suppression region is kept almost constant.But the peak value of the flow profile increases with density.Additionally, it should be noted that in the higher density regime, one encounters a density limit for the heating power near the threshold.This is because the finite heating power is not sufficient for the E × B flow shear rate to exceed the shear suppression threshold due to the dramatic increase of density after the L-H transition.
(c) The linear dependence of the transition power threshold on the magnetic field.In this case the dependence of edge density on the magnetic field is the key issue.The magnetic field exponent strongly affects the transition power threshold.Contrary to the situation of the transition power threshold variation with density, for the variation with the magnetic field the ionization depth is dominated by temperature rather than density.The average inverse scale length of pressure is dominated by the density gradient rather than the temperature gradient, and the fine structure of the E × B flow shows that the peak of the profile tends to move towards the edge with a constant value, but the average of the E × B flow over the suppression region decreases with the magnetic field.
Finally, there are still several issues remaining for further improving the model.These include the achievement of the lower density branch with threshold power increasing for decreasing density, which will require the model to include the physics for describing the collisional coupling of the electrons and ions for energy transfer from heated electrons to ions in low density regimes.In the present study, the heating is directly through the ions which imply a monotonic increase in the threshold power with density.The detailed dependence of the threshold power on the boundary conditions and the dynamics of the radial electric field at the LCFS require the model to apply more accurate SOL physics.These issues are the subject of further work.
The preliminary results of the model presented in this paper implied that the physics of the microscopic turbulence dynamics and the macroscale power scaling of the L-I-H transition are one thing with two aspects.The model provided convincing scaling relations and will be upgraded to study back-transitions as well as the effects of different fuelling and heating schemes with the final goal of providing ITER-relevant physics. bV

Figure 1 .Figure 2 .
Figure 1.The profiles of heat and particle sources in the reference case of L-mode equilibrium.

Figure 3 .
Figure 3. Profiles of (a) density, (b) pressure and (c) temperature at L-mode (blue line) and H-mode (red line), respectively.

Figure 4 .
Figure 4. Profiles of (a) turbulence intensity and (b) the inverse scale lengths of pressure gradient at L-mode (blue line) and H-mode (red line) respectively, where the vertical black line indicates the turbulence suppression region of 0.94 < r/a < 1.

Figure 5 .
Figure 5.The profiles of the poloidal flow and neoclassical poloidal flow at L-mode and H-mode.The ion (V θ > 0) diamagnetic direction is indicated.

Figure 6 .
Figure 6.The profiles of the E × B flow, the poloidal flow and the diamagnetic flow from the radial force balance equation at L-mode and H-mode.The ion (V θ > 0) and electron diamagnetic directions are indicated.

Figure 7 .
Figure 7.The profiles of the ion-ion collision frequency ν ii and ε 3/2 ω t in L-mode.

Figure 8 .
Figure 8. Spatio-temporal evolution of (a) E × B flow and (b) turbulence intensity during a power ramp from 4000 ms to 4200 ms with a 1% increment.

Figure 9 .
Figure 9. Plots of turbulence intensity I versus the absolute E × B flow |V E | at the radial location of r/a = 0.97 (a) and r/a = 0.996 (b).The filled symbols denote the final point.

Figure 10 .
Figure 10.Time evolution of turbulence intensity (a), E × B flow (b), diamagnetic flow (c) and poloidal flow (d) at various radial locations of r/a = 0.97 (blue line) and r/a = 0.996 (red line).

Figure 12 .Figure 13 .
Figure 12.With a fixed magnetic field B/B 0 = 1, the ionization depth and edge temperature versus line-averaged density.

Figure 14 .Figure 15 .
Figure 14.With the fixed magnetic field B/B 0 = 1, (a) the profiles of E × B flow in different line-averaged density situations at the turbulence suppression region.(b) The average and the minimum of the E × B flow over the suppression region vary with different line-averaged densities.

Figure 16 .Figure 17 .
Figure 16.With the fixed line-averaged density n / n 0 = 1, the ionization depth and edge temperature versus magnetic field.

Figure 18 .Figure 19 .
Figure 18.With the fixed line-averaged density n / n 0 = 1, (a) the profiles of the E × B flow in different magnetic field situations at the turbulence suppression region.(b) The average and the minimum of the E × B flow over the suppression region vary with different magnetic fields.