Drift wave theory for transport in tokamaks

We revisit the theory of tokamak transport due to drift waves. We recall the whole development, starting with simple theories from the 1950s till today advanced theories using advanced fluid models taking full account of kinetic effects in the frequency regime of drift waves, also including the effects of zonal flows, and also fully nonlinear kinetic theory itself. Traditionally drift waves have been described by either kinetic theory or expanded fluid theories. The expansions have usually been made in the ratio of magnetic drift frequency and frequency. As was just shown toroidal drift waves are mainly driven by the magnetic drift resonance so such an expansion is usually not allowed. It is, in fact, natural that toroidal effects play an important role since they originate from bending the system to a torus, thus eliminating the third direction in which the magnetic field does not confine plasma. Due to an exact fluid closure, we are now able to use fluid theory completely without expansion, thus maintaining the fluid resonances due to magnetic drifts in the denominators. This gives us a new normalization of drift wave equations and this enables us to recover nonlinear (Dimits) upshifts, spinup of poloidal rotation in internal transport barriers and the L–H transition. The principle of our reactive closure is that we include all moments with sources in the experiment. Here, it is usually enough to use the diamagnetic heat flow as closure term but more general cases are, of course, possible and remain parts of our general closure. The fact that the model is self-consistent then also leads to the overall experimental power scaling e ∼ P . Using a full transport matrix we are then also able to obtain adequate particle pinches and apparently nonlocal phenomena such as the heat pinch on DIII-D. Recently an extension of the derivation of the fluid closure also showed that quasilinear theory works well for the real part of the eigenfrequency in our fluid description of driftwaves. A fluid description using an exact closure with a fully valid quasilinear approach also lets us cover cases with nonlinear thresholds for gradients then including also sand pile thresholds. Finally, we use a correlation length at the maximum of the E × B drift in k-space. Then we can go beyond the gyro-Bohm scaling if we allow the parameter k to vary. Thus, we are, in practice, able to relax all limitations of traditional drift wave theories. This means that we can recover all aspects of low-frequency tokamak Extended author information available on the last page of the article Reviews of Modern Plasma Physics (2019) 3:8

Although neoclassical effects can still be important in some regions e.g. close to the magnetic axis where turbulence usually is weak, or in transport barriers, turbulent ion energy transport is typically a factor 3-10 stronger than neoclassical transport. Turbulent transport is usually strongly coupled to drift motions set up by the confining gradients. Drift velocities associated with such gradients are naturally associated with frequencies which are typical of the turbulent motion. Since turbulent motion in tokamaks is usually dominated by E × B motion, the turbulent electron energy transport is usually comparable to the turbulent ion energy transport. However, since neoclassical transport is proportional to the square of the gyroradius, electron neoclassical transport can almost always be neglected.
Then, having established the overall picture of transport we also have to look at the problem of global stability. Here the most dangerous instabilities are global kink modes and interchange modes. As it turns out the worst of these instabilities can be avoided using a suitable rotational transform usually expressed in terms of the safety factor q = r B t /RB p where t and p stand for toroidal and poloidal, respectively. Now kink modes are driven by the current which generates B p . Thus, kink mode stability requires a large enough q. Interchange modes, on the other hand, are driven by curvature. However, in a tokamak the most dangerous interchange modes localize in regions with unfavourable curvature. These are called ballooning modes. The length of the regions with unfavorable curvature increase with q. Ballooning modes are because of this more stable for small q. Fortunately, we can usually find a regime with 1 < q < 6 where both ballooning modes and kink modes are reasonably stable. There is also a stability boundary against external kink modes which is q > 2 at the outer boundary. However, this typically limits the plasma beta to below around 3% including also neoclassical effects. This concerned global magnetohydrodynamic (MHD) stability where the rotational transform q was the most basic stability parameter. Concerning transport, a suitable q turned out not to be enough. As it turned out, multipole confinement systems had a very strong transport, Bohm transport, scaling as T/B. This was for a while considered as a real threat to the possibility to achieve nuclear fusion. This very strong transport turned out to be due to large electrostatic vortices called Convective cells (Taylor and McNamara 1971;Okuda and Dawson 1973;Sagdeev et al. 1978;Hassam and Kulsrud 1979;Weiland 1980). The scaling T/B was first derived by Taylor and McNamara (1971) and was later described in more detail by Okuda and Dawson (1973). Basic fundamental work on convective cells was also made by Sagdeev et al. (1978). Fortunately, it turned out that convective cells could be short circuited by magnetic shear. This opened up for new optimism about the achievement of nuclear fusion. Another improvement of confinement could be obtained due to favourable average curvature of the magnetic field lines (Hassam and Kulsrud 1979). However, this effect could be partly overcome by higher order nonlinear effects (Weiland 1980).

The level of transport
The level of transport in tokamaks is, of course, strongly dependent on the level of the turbulence.
The usual mixing length estimate of the saturation level is to use the condition that the nonlinear convection equals the linear (Kadomtsev 1965). This is possible because the perturbation varies much more rapidly in space than the background. Thus, introducing the cross field drift v E = E × B/B 2 we get: where L n is the scale length of the background variation and k r is the inverse radial correlation length.
We are now interested in determining the maximum displacement corresponding to the turbulence level given by (2). The condition (1) leads to the usual mixing length estimate for transport (Kadomtsev 1965;Dupree 1967) where γ is a typical linear growth rate at the correlation length. We will, however, here use a somewhat more general approach. We then balance the linear growth rate with the convective nonlinearity (Weiland and Nordman 1988). Fortunately, this gives the same result for the continuity and energy equations: We now see these as equations in real space, replacing gradients by 1/L where L is the correlation length. We also consider the mode number corresponding to the inverse correlation length. Then cancelling the perturbations (4a, 4b) leads to the level of excitation (Weiland and Nordman 1988;Weiland et al. 1989;Connor and Pogutse 2001): where φ is the electrostatic potential, ρ s is the ion gyroradius at the electron temperature and c s is the ion acoustic velocity. Equation (5a) corresponds to the E × B displacement: Here the displacement is reduced by the real eigenfrequency and turns into the wave breaking limit when γ ≫ ω r . Thus (5a, 5b), which is non-Markovian, gives a more accurate result than (2). Here a finite ω r means that the displacement is reduced due to rocking of the eddies. Since the drift frequency and the growth rate typically are of the same order for drift waves the level (5a) is usually comparable to the level (3). Both are thus comparable to the experimental result (Liewer 1985): We note that calculations are usually made to the cubic power of the small parameter ε corresponding to the nonlinear Schrödinger case. However, then higher order resonant terms would be of order ε 5 which would be ignorable.
As shown by (4a, 4b) the continuity and energy equations lead to the same level of turbulence, (5a, 5b). Then a formal calculation of the transport level, just using the relation between E × B displacement and drift, gives us The kernel (7) was first obtained as a part of the effective thermal diffusivity of ITG modes (Weiland and Nordman 1988) but the density part was later rederived for resistive ballooning modes (Connor and Pogutse 2001) using a reactive fluid closure at zero ion temperature. Using Boltzmann electrons, the thermal conductivity we derived in Weiland and Nordman (1988) was: but has later been recovered in several other derivations (Weiland et al. 1989;Bateman et al. 1998). Equation (8a) was derived from the quasilinear temperature perturbation: by assuming Boltzmann electrons. In this case, the pinch terms cannot give a net pinch. However, in more general cases, including electron trapping or electromagnetic electrons we can have pinches of all involved free energies. The only restriction is that the total energy moves outward. Note that we have here kept the magnetic drift resonance in the denominator. This is possible using our exact fluid closure derived in . Expansion in the magnetic drift, which has been customary, is usually not allowed. This can be seen from the real part of the ITG mode frequency: where ε n = 2L n /R and flr indicates gyroradius corrections. We note that the fact that the ITG mode propagates in the ion direction (ω r < 0) is a toroidal effect which dominates in almost all the discharge Another feature of (8a) is that it involves only one mode number at the inverse correlation length. Of course, the transport coefficient will first be expressed as a convolution sum in k-space. In Weiland and Nordman (1988) we evaluated the convolution sum numerically and found that this sum could be replaced by one term at the correlation length. The correlation length was there found numerically to be the inverse mode number of the linearly fastest growing mode as normalized by the drift frequency. The detailed analytical interpretation of this came as late as 2016 in Weiland (2016) (Eq. 15) and concludes that the correlation length is at the maximum of the E × B drift in k-space. Then, using this mode number, Eq. (8a) is in full agreement with all our previous equations. However, we have here used the assumption that E × B convection is entirely stabilizing. This is only true at the correlation length. Thus, we must take care to use (8a) only at the correlation length.
We here used the quasilinear approximation. However, as has been found later (Weiland 2018; this is a very good approximation (of order ε 2 ). This means that from the transport point of view our quasilinear treatment is equivalent to that of a fully nonlinear and nonlocal turbulence code [agreement with a local turbulence code (Rogers et al. 1998) was found in Weiland and Nordman (1988), Weiland et al. (1989) and Weiland (2014)]. The kernel (7) was also part of the diffusivity in the Hasegawa-Wakatani system (Wakatani and Hasegawa 1984;Connor and Pogutse 2001). Also in this case the quasilinear approximation was found to be very good (Connor and Pogutse 2001). In (8b) the Doppler shift in the denominator is due to the magnetic drift of ions which is the average velocity the ions are moving with. Thus, this is a signature of the magnetic drift resonance. When electron trapping is added (Liu 1969;Weiland et al. 1989) we find the corresponding drift resonance for electrons. These different resonances mean that the ITG mode mainly drives ion transport and the Trapped Electron (TE) mode mainly drives electron transport (Weiland et al. 1989). Equation (8a) shows the effective transport coefficient for the simplest ITG mode with Boltzmann electrons. The factor in front includes off-diagonal contributions. Here, the pinch term prop to ε n = 2L n /R means that transport is reduced towards the axis. Actually also the growth rate contributes to this. Now Eq. (8a) gives us the overall level of transport in the simplest possible case for illustration. More general cases calculated along the same principles have turned out to usually be in good agreement with both experiments (Bateman et al. 1998 A very important aspect is here that we have not expanded in the magnetic drift (Weiland 2014;. Actually, the pinch term in (8a) is due to quadratic effects of the magnetic drift. As is evident this pinch will reduce the heat flux towards the axis and tends to give a transport coefficient that grows with radius. This had previously been one of the main problems for drift wave models as also shown in Scott et al. (1990). We studied this case in Weiland (2012) with much better result than obtained in Scott et al. (1990). The result in Weiland (2012) was obtained from an interpretative code. Predictive codes (Weiland and Nordman 1991) usually give a more clear radial growth of transport. The fact that the radial growth of transport coefficients is a toroidal effect was also found in kinetic simulations by LeBrun et al. (1993).

Role of the fluid resonance
As has been emphasized in many of our papers we usually need to keep the fluid magnetic drift resonance. This is also true for the kinetic magnetic drift as pointed out in Weiland (2014) and .
The best overview of the need to keep the magnetic drift resonance is probably Weiland (2014). We start by observing that Eq. (8b) holds both in the adiabatic and isothermal limits. Thus, it includes the drift resonance. That this is exactly correct follows from . However, this has previously been anticipated on good grounds. Originally almost all drift wave models expanded in the adiabatic limit. This gave fairly good agreement with the physics in the outer third of the tokamak cross section but did not work at all closer to the axis. Our model was actually the first to solve this problem (Weiland et al. 1989;Bateman et al. 1998). This also meant that it worked for both L and H mode equilibria (Weiland and Nordman 1991;Bateman et al. 1998). Of course, this has also been necessary for obtaining the full L-H mode transition in the present code (Weiland 2014). Our first transport code was published in 1991 (Weiland and Nordman 1991). It was a 1.5 D transport code as mainly used in the ITER Expert group (ITER Physics Basis Editors, ITER Physics Expert Groups 1999) and ITPA (ITER Physics Basis Editors 2007) tests of transport models. It uses a flat boundary near the axis and fixed edge boundaries at the separatrix. The independent coordinate is the radial flux coordinate and usage of local elongation and flux transformation provides the additional half dimension. We use an implicit method to integrate profiles of temperature and density in time until we reach a steady state. Sources of temperature and density are numbers taken from a database. We obtained an L-H transition already with our first transport code (Weiland and Nordman 1991) but in these simulations the only stabilizing mechanism was finite larmor radius (FLR). There was no parallel ion motion and no resistive ballooning modes so it was the drift waves that were stabilized. In this situation, the transport barrier was not resolved but the increase in confinement time was a factor 2.5. Thus the interior had the character of a real H-mode equilibrium. Then we also found the experimental scaling of the confinement time with power

Fluid closure
A very important property of (4a, 4b) is that the nonlinearity has been assumed to be entirely stabilizing. This means that the turbulent waves always propagate out from the source at the correlation length. To avoid reflections, we then have to require the inverse cascade towards longer wavelengths to be absorbed. The condition for this was recently found to be fulfilled due to zonal flows (Weiland and Zagorodny 2016) (Fig. 1).
The strong flowshear at marginal stability is due to the exact fluid resonance in (8b). We realize that the waves with the wavelength approaching system size will always hit marginal stability somewhere, thus becoming strongly damped. Thus, the form (7) of the transport kernel strongly depends on the fluid closure. Since dissipation due to wave particle resonances strongly reduce the zonal flows, we need a reactive fluid closure to have transport kernels of the type (7). In Wakatani and Hasegawa (1984) and Connor and Pogutse (2001) the reactive closure was due to zero ion temperature but in our case it is due to detuning of the wave particle resonance by nonlinear frequency shifts .
The basic coherent papers which our fluid closure is based on are those by Mattor and Parker (1997) and Holod et al. (2002). These were generalized in . We point out here that it is trivial to include the magnetic drift in Mattor and Parker (1997) and Holod et al. (2002) since it just introduces a phase shift. We will here look at a multiple three-wave system as described in .
u 3l sin( jkl + Δ jkl t). Such systems have two attractors which are the random phase levels and the center of force for trapped particles. As was found by Mattor and Parker it is really the nonlinear frequency shifts that cause the phase mixing. Since we can here go to high order systems we can reach the situation where we can derive a Fokker-Planck equation. Also in this limit, it is the nonlinear frequency shifts which cause the phase mixing.
Here, the friction and diffusion are due to turbulence. For constant coefficients (12a, 12b,c, 12d) has an analytic solution given by Chandrasekhar (1943). Here, we have added a source S H which is supposed to represent heating or current drive. This source would correspond to a continuous modification of the slope of the velocity distribution in such a way that the nonlinear oscillation of the phase velocity of turbulent waves would not cancel the wave-particle energy exchange. We can regard this as a source in velocity space. Thus, in this case we would need to describe transport in velocity space and that, in fact, requires a kinetic description. However, the drift waves that cause bulk transport in configuration space have frequencies that are about two orders of magnitude below that of the heat sources so at these frequencies we will have We also recall that the analytical solution by Chandrasekhar is only valid for constant coefficients in the Fokker-Planck equation. We have, however, tested more general cases (Holod et al. 2005) including also the magnetic drift frequency numerically and again recovered an asymptotic situation without average energy exchange between resonant particles and waves. Since we here approach a stationary asymptotic state the Chandrasekhar solution will be valid asymptotically. Equation (12d) is thus the only limitation of our closure which then is exact for bulk transport.
We can also include the Chandrasekhar solution into an orbit integration (Weiland 2018;: Page 10 of 21 This means that u is the fluid velocity and w is the kinetic velocity. Now (14a) will contain both fluid velocity part and kinetic velocity parts. We denote these by subindices f and k as and analogously for the nonlinear part. Then taking the ensemble average in the saturated state in Fig. 8 we know that the only possible correlation is due to waveparticle interactions so we must have: Thus, we have and we have only fluid parts left. Now we know that Accordingly, we can conclude that Now this is mainly due to ω = ω r + iɣ. Then, if we are close to marginal stability (17) would not necessarily apply to the imaginary part. This is e.g. the case when an instability is stabilized by nonlinear flowshear (Waltz et al. 1995;Zagorodny and Weiland 2009) However, the real part would still fulfil (17). Thus we can conclude that the linear real eigenfrequency is generally a very good approximation for the real part of ω in fluid theory as also found by Connor and Pogutse (2001).
We also have interesting support for our nonlinear upshift from the fully nonlinear kinetic simulations in Dimits et al. (2000). This is from Hatzky et al. (2002) showing that the resolution in number of particles is particularly sensitive close to wave-particle resonances, as is the case in the nonlinear upshift region and from the fluctuation-dissipation theorem, Sitenko (1982). When the number of particles is insufficient the fluctuations are increased since particles can not effectively screen each others space charge. For larger fluctuations, we also have larger dissipation. Thus, we effectively introduce dissipation into the nonlinear particle simulations. Thus, the width of the upshift was reduced in the Parker simulations in a similar way to that from the IFS-PPPL model (Dimits et al. 2000).

Particle and heat pinches
As mentioned above our fluid model includes a full transport matrix. For drift waves, this basically means that we include ion and electron heat pinches and electron and (15a) (18) → − E×B .
impurity particle pinches. This property was originally tested on the electron heat pinch in DIII-D (Luce et al. 1992;Weiland and Nordman 1993).
Here, we used a rather early version of our model which was electrostatic and did not include zonal flows explicitly (Weiland et al. 1989). It did, however, of course use the normalization (7) which assumes the implicit presence of zonal flows. We here studied a DIII-D L mode which had only ECH heating (Luce et al. 1992;Weiland and Nordman 1993). Thus, there was no externally generated ion rotation and we had fairly low beta. Thus, the electrostatic version of our model without explicit zonal flows (Weiland et al. 1989) should apply. However, the real part would still fulfil (17). Thus, we can conclude that the linear real eigenfrequency is generally a very good approximation for the real part of ω in fluid theory as also found by Connor and Pogutse (2001). Thus, this was a sensitive test of the most fundamental properties of our model including the fluid closure and the normalization (7). Our transport equations in matrix form were: The experimental setup was an ECH source at half radius and gas puffing at the edge. The result is shown in Fig. 2.
In Fig. 2a we show the electron temperature profile. It has a very small wiggle at the source. There is also a wiggle in the experimental points The very small trace of the experimental point was, by several observers, seen as a sign of a nonlocal mechanism. However, we are here using a local drift wave model. The sensitivity is due to the balance between outgoing and incoming flux components as seen in Fig. 2b. This model is still the limiting case of our present model and thus includes exact reactive fluid drift resonances. The heat pinch requires a density gradient. The density was not shown in the original PRL by Luce, Petty and Haas which was our only background material at the simulation. Thus, we entered a particle source at the edge which we kept as a free parameter. However, since there is no particle source in the interior the total particle flux has to vanish in the interior. This means that incoming and outgoing flux components have to vanish also for the particle flux. This completely defines the density profile for a given temperature profile. Thus, in this simulation the electron temperature and density profiles build up self consistently and in the end agree quite well with the experiment. The expression for the particle flux is comparatively simple and proportional to n (20). We can here see that for moderately large ε n it is the ITG mode that will tend to drive a particle pinch.
It follows from the second part of Eq. (20) that an inward particle flux (Γ < 0) will be driven by an electron temperature gradient for an ITG mode (ω r < 0). We note, however, that also trapped electron modes (ω r > 0) can drive a particle pinch for sufficiently 1 3 flat density. However, there is a built in tendency for the system to move towards a state where length scales of the dynamic gradients do not differ too much. Thus, the second term in (20) usually dominates. Nevertheless, we should not ignore the third term. The simulated particle profile for our DIII-D shot is shown in Fig. 2c. The importance of the density gradient was later verified by Austin et al. (2009) redoing this experiment in H mode with flat density. In this case, there was no heat pinch.
The fact that we here get a strong particle pinch is also important. The equations used here were from Weiland et al. (1989). This was the first model to achieve a strong particle pinch of experimental order. This problem was expressed by Wagner and Stroth (1993). Fig. 2 a Experimental points and inferred profile of electron temperature and simulated T e and T i profiles. b Electron and ion simulated heat fluxes. c Simulated density which was obtained completely without knowing the experimental result. Reproduced from Weiland and Nordman (1993) with the permission of AIP publishing

Parallel ion motion
Our first papers on tokamak transport (Weiland 1980;Weiland and Nordman 1988;Weiland et al. 1989;Nordman 1991, 1993) all ignored parallel ion motion. Nevertheless, there was often very good agreement with experiment, Nordman 1991, 1993;Bateman et al. 1998;Tardini et al. 2002;Weiland 2012). Parallel ion motion is usually only important near stability boundaries but then, of course, plays a major role in connection with transport barriers. The model for parallel ion motion presently used in our model, with general ballooning, was introduced in Weiland (2004). The eigenvalue equation for the pure ITG, using the ballooning formalism (Connor et al. 1978) was written: where η is the generalized poloidal angle (Connor et al. 1978). Our eigenvalue equation then takes the form: With asymptotic solution where Thus we get our final solution (Weiland 2004): where ⟨⟩ means average over the asymptotic eigenfunction given by (22b). Here, we can see the reason for the formation of transport barriers for small magnetic shear. This is the fact that our eigenfunction becomes flat for small shear which leads to a (21a) 1 3 small average of the harmonic function g. The usual technique to combine asymptotic and interior solutions is to use a constant φ approximation in the interior. Here, since our eigenfunction is varying slowly in the interior we can use this eigenfunction everywhere. In fact, the solution (23) turned out to be in excellent agreement with shooting code results (Weiland 2004). Small shear was also important in obtaining internal transport barriers on JET (Weiland et al. 2011). A more recent result is the simulation of the particle pinch on Tore-Supra by the Cadarache and SWIP groups (Zhong et al. 2013). We recovered very similar results but also could identify the mode profile   (Fig. 3). There we got a region of strong ballooning, enhancing the ITG mode at the experimental pinch.
As mentioned above the direction of the particle flux usually depends on which mode is driving the flux. In the experiments on Tore-Supra the particle pinch appeared in the region where ITG was dominant. This was found both by the Cadarache and SWIP groups (Zhong et al. 2013) using QualiKiz (Bourdelle et al. 2007) and our group . In our case, the dominance of ITG was shown to be enhanced by strong ballooning. There are actually two trapped electron modes which could play a role here. These are the Ubiquitous mode (Coppi and Pegoraro 1977) driven by charge separation and the mode driven by compressibility and electron temperature gradient (Liu 1969;Weiland et al. 1989). The ITG mode (Coppi and Pegoraro 1977;Horton et al. 1981) is always driven by compressibility and is symmetric to the TE mode driven by compressibility. Nice overviews of these modes were shown in Garbet et al. (2004), Weiland (2012), , Kikuchi and Azumi (2015) and Weiland (2018) where several considered the Tore-Supra experimental regime.

Momentum transport
We have already indicated that an implicit effect of zonal flows is behind the transport kernel (7). When zonal flows are important at the correlation length we need to calculate them explicitly. The most important effect is on poloidal rotation. It can be described through the transport equation: where the flux p is generated nonlinearly by the Reynolds stress: The toroidal momentum flux has, is our simulations, mainly been important for the spinup of poloidal momentum in internal transport barriers. The generating equation is ): Fig. 3 a Effective diffusivity as a function of temp. gradient for a Tore-Supra shot in.
b The mode width parameter α for the same scan as in a. Reproduced from  with the permission of AIP Publishing where we here added electromagnetic effects. We calculate the combined influence of toroidal and poloidal momentum through their generation of a radial electric field (Staebler et al. 1997)

Transient transport
To test different transport models, it has become usual to use transient mechanisms such as oscillatory strength of heating. The first test of our model in this respect was made by the ASDEX group, using modulated ECH heating on ASDEX upgrade (Tardini et al. 2002). As a background, we here repeat the way we calculate the trapped electron transport (Weiland et al. 1989). The total electron density response is calculated as: Here, subindex t stands for trapped and subindex f for free. For f t we use the deeply trapped particles i.e. f t = (2ε/(1 + ε)) 1/2 . The free particles are considered to be completely free. Thus, we ignore barely trapped and barely free particles. These populations have opposite effects on the total response and, in fact, this representation has usually turned out to work well. One example is the calculation of the heat pinch on DIII-D above. The results by the ASDEX group (Tardini et al. 2002) are shown in Fig. 4.
As is evident from Fig. 4 the result is almost perfect. This confirms the validity of our treatment of trapped electrons.
Later tests for ion modulation have not been quite as good. However, for ion modulation we have to remember the importance of the kernel (7) which depends sensitively on the absorption of the inverse turbulent cascade by zonal flows according to Fig. 1. This mechanism could be disturbed by a modulation in the ion heating in the experiment since it is the ions which generate the zonal flows.

The L-H transition
The zonal flow generation in our code has turned out to, in addition to the spinup of poloidal rotation in internal transport barriers and the Dimits shift, to also give us an experimental. L-H transition (Weiland 2014).
This was first demonstrated on EAST (Weiland 2014) but has later also been tested on JET and ITER (Weiland 2018).
A problem has been that we need to reduce the fixed edge boundary conditions to get nice pedestals. However, the boundary conditions usually have quite small temperature and density so measurement errors can easily give large relative errors. Then, for numerical reasons, we can not go closer to the separatrix than to r/a = 0.97. This will have a tendency to enlarge the difference. In our original EAST simulation, the temperatures were 30% of the experimentally inferred values at r/a = 1. However, the density was about 0.1 times the experimental. An interesting aspect is here that the density gets much flatter in H-mode. This has also been seen experimentally. Another interesting result is that alpha particle heating reduces the slope of the H-mode pedestal (Weiland 2016). This will be favourable by reducing the wall load. We have updated the calculations from Weiland (2014) by including also ion viscosity. This has been done in the same way as was introduced for trapped electrons in Nilsson and Weiland (1994). The result is shown in Fig. 5. In Fig. 5, in particular the density is interesting. There is a particle pinch on the edge barrier which is the location of the only particle source (gas puffing). Thus, the pinch carries particles inward. However, this is the dominant pinch in the system. Some pinch fluxes are also present locally in the interior and the central density is actually higher than the density on top of the pedestal. We note the difference to Fig. 2c. The density in H-mode is much flatter than that in L-mode although there are no internal particle sources in either case. The main difference is due to the very steep pedestal in H-mode. In H-mode there is no real steady state but the interior density builds up very slowly.

Discussion
We have in this work given an overview of the development of the theory for turbulent transport in magnetically confined fusion plasmas. We emphasize the importance of the fluid closure and the magnetic drift resonance. In fact, all problems with early drift wave models are solved by introducing a model which does not expand in the magnetic drift frequency. We here focus on achievements made by a two fluid model with an exact reactive fluid closure. This fluid closure enables us to include unexpanded denominators with magnetic drift resonances included. In combination with a full transport matrix this means that all limitations on simpler drift wave theories are eliminated. The full transport matrix is often needed and it is not advisable to try to calculate particle transport without including a possible drive from the electron temperature and vice versa. Off diagonal elements are usually due to toroidal effects which enforces our previous advice to keep the magnetic drift resonance at full strength. At the same time, this is a fluid theory so numerical integration over local particle positions and velocities are eliminated. Our only limitation is that we still need kinetic theory for heating and current drive. It is also interesting to see that the importance of the magnetic drift resonance is similar in kinetic and fluid theory ). One of our most recent results is that we, by combining Dupree's method of integration along nonlinear orbits with Chandrasekhars analytical result (which we can generalize numerically) on the phase mixing of wave particle interaction due to nonlinear frequency shifts, we can integrate also the effects of the real frequency in the phase mixed case. Note that this is still strongly nonlinear so it goes beyond quasilinear theory.