Validation of a full-plasma integrated modeling approach on ASDEX Upgrade

In this work we present the extensive validation of a refined version of the integrated model based on engineering parameters (IMEP) introduced in reference (Luda et al 2020 Nucl. Fusion 60 036023). The modeling workflow is now fully automated, computationally faster thanks to the reduced radial resolution of the TGLF calculation, and it includes the modeling of the toroidal rotation, which was still taken from experimental measurements in our previous work. The updated model maintains the same accuracy as its previous version when tested on the cases presented in the initial publication. The confined plasma, from the magnetic axis to the separatrix, is simulated without using any experimental information from profiles measurements, and the inputs of IMEP are the same engineering parameters used when programming a plasma discharge. The model validation database consists of 50 ASDEX Upgrade (AUG) stationary (over a few energy confinement time) H-mode phases, which largely cover the entire AUG operational domain. The prediction of IMEP is compared with experimental measurements and with scaling laws, such as the IPB98(y,2), the ITPA20-IL, and AUG specific regressions. This modeling framework has proven to be very accurate over the entire set of 50 cases, with a significantly lower mean relative error with respect to each of the scaling laws considered, accurately reproducing the change in pedestal and core confinement caused by a change in plasma current, heating power, fueling rate, triangularity, magnetic field, NBI voltage (i.e. the effect of a change in the core particle source), and heating mix (e.g. correctly predicting the effect on confinement caused by a change in T e/T i). Plasma confinement is correctly described by IMEP also for two particular operating regimes, such as the ITER baseline scenario, and the QCE regime (quasi continuous exhaust, also referred as type-II and small ELMs). This work clearly demonstrates the power of this approach in pulling out physics mechanisms to interpret subtle interdependencies and that a 1D integrated model can reproduce experimental results over very large parameter variations with a higher accuracy than any statistical regression. This approach has therefore the potential to improve the prediction of the fusion performance in future tokamak reactors.


Introduction
The prediction of the plasma global energy confinement plays an essential role in the design of future fusion devices and in the optimization of their operational scenarios. So far, the main tools for this type of predictions have been provided by empirical scaling laws, which have been obtained through multivariate regressions on large databases of observations from multiple tokamaks [1][2][3]. Recent progress in the understanding of the transport properties [4][5][6][7] and MHD limits [8][9][10][11] determining the confinement, as well as in the integrated modeling activity [12][13][14][15][16][17][18][19][20][21], have opened the possibility to obtain predictions of the plasma stored energy which use a limited number of data from the experimental measurements. The integrated modeling approach cannot only be made progressively independent of experimental measurements, but also increasingly based on theoretical models for the prediction of the turbulent transport and the MHD stability. Moreover, the integrated modeling approach with transport codes provides predictions of radial profiles of electron and ion temperatures and densities separately, a much more valuable information in the design of the operational point of a reactor with respect to a 0D approach, which only provides the confinement time or the total stored energy [22]. However, in the integrated modeling approach, a still critical and challenging step is that of obtaining predictions of the plasma kinetic profiles over the entire minor radius without using any information from actual measurements, that is, only using those (mainly global) parameters that enter the pulse program for the plasma discharge, such as magnetic field, plasma current, plasma shape, auxiliary heating powers and gas puff rates. Furthermore, a systematic validation of the predictions obtained from this approach over extended databases, covering wide variations of multiple parameters, have been rarely performed so far, even from one single device.
The prediction of the pedestal performance for ITER is based on models, such as EPED [23], which still contain assumptions on the input parameters and the boundary conditions. As an example the pedestal top density must be given as input, and the temperature and density values at the separatrix are given by simple assumptions (e.g. n e,sep = n e,ped /4). This poses two limitations: first, different assumptions on the pedestal density profile and its gradient can have a strong impact on the result of the pedestal stability calculation (as observed in ASDEX Upgrade (AUG), JET-ILW and TCV [24][25][26]), and second, it is important to take into account the effects of the scrape-off layer (SOL) and divertor conditions on the densities and temperatures at the separatrix [20]. This is because in tokamak operation, and particularly in a reactor, SOL parameters have to respect the constraints imposed by power exhaust. It is therefore important to take into account self-consistently the impact that different regimes of operation can have on the SOL and on the pedestal, and ultimately on global confinement. The model we developed allows the description of the coupling between SOL, pedestal and core, with a self-consistent treatment of the boundary conditions. The only assumed parameter is the effective charge Z eff , which is set as a fixed value (Z eff = 1.2) for most of the cases, or taken from experimental measurements when available. The model has proven to accurately capture the effects of some critical parameters on confinement, with fully predictive capabilities. It is important to extend the validation of the model on a large database to assess the reliability of such approach to correctly describe the coupling between the different plasma regions and their role in defining plasma confinement. This could help to understand if this approach has the potential to provide more reliable predictions for ITER and future tokamaks.
In this paper we present further developments of the integrated model based on engineering parameters (IMEP) introduced in reference [27]. The first test of the integrated model, as discussed in the previous publication, established its reliability to capture some fundamental dependencies of plasma confinement, specifically, that on the current, heating power, and fueling rate, and therefore also on the density. The positive outcome of this first validation suggested that no major ingredient of the model required a significant revision. Therefore, it has been possible to concentrate on a small refinement of the physics description included in IMEP, and an improvement of the modeling workflow execution in order to automatically run it over a large experimental database in an efficient manner. Such increased efficiency allowed a validation of the model on a database of 50 AUG experimental cases, spanning a wide range of AUG operating parameters. This application allowed us to better address the reliability of this approach in correctly describing the coupling between the different plasma regions and their role in defining plasma confinement.
This work describes the extensive validation performed, which also demonstrates the power of this approach in the identification of physics mechanisms which can explain observed complex inter-dependencies, particularly connecting different plasma regions like the SOL and the pedestal. We also describe the improvements and the new features that we have introduced in IMEP, in order to obtain an automatic modeling procedure which does not use any experimental measurement of the plasma densities and temperatures for the prediction of their complete radial profile, from the magnetic axis to the last closed flux surface.
For convenience of the reader, we briefly recall what has been done previously. As presented in our previous work, one of the main novelties of our integrated modeling workflow consists in a new pedestal transport model (implemented in the ASTRA transport code [28]), which is based on the experimental observation found by analyzing a pedestals database of AUG, DIII-D, and JET. The pedestals from these different tokamaks all showed to have a common feature: a constant ratio between the averaged pedestal electron temperature gradient (ETG) (in real space units) and the pedestal top temperature ∇T e /T e, top = −0.5 (cm −1 ), as shown in figure 2(a) of reference [29]. Since the gradient is computed with a non-normalized radius, this quantity is not dimensionless. This observation can be qualitatively related to a critical pedestal averaged electron temperature scale length, which can be interpreted as an R/L T e drive for turbulent transport, and therefore can be associated to ETG modes or micro tearing modes. The inclusion of this constraint into the ASTRA transport code allows us, for a given pedestal width, to find the electron heat conductivity (χ e,ped ) that fulfills such condition. We assume that the turbulent part of the ion heat conductivity and the particle diffusivity are proportional to the found χ e,ped which is consistently determined in the ASTRA simulation to fulfill that condition. The transport simulation (for a given value of the pedestal width) evolves the kinetic profiles until stationary conditions are reached. We launch many ASTRA simulations in parallel, each with a different value of the pedestal width. Figure 1 shows the resulting pedestal top pressure (a), and the profiles of pressure (b) and electron temperature (e) for different values of the pedestal width. By increasing the pedestal width, the pedestal top electron temperature increases, together with the density and the ion temperature, and therefore also the total pressure increases. Since the pedestal temperature and its average gradient are linked by the constraint ∇T e /T e, top = −0.5 (cm −1 ), also the temperature and the pressure gradients increase with increasing width, as illustrated in figures 1(c) and (f ). Thereby, with increasing pedestal width, a reduction of χ e,ped , illustrated in figure 1(d), allows the development of stronger temperature gradients. This is an important feature for instance in the understanding of the mechanism which allows the model to capture the change in confinement caused by different values of the NBI voltage, as presented in section 5.3. The core transport coefficients are calculated using the electromagnetic quasilinear turbulent transport model TGLF [5] and the neoclassical transport code NCLASS [30]. The latter is also used to compute the poloidal rotation, and the current conductivity and the bootstrap current for the simulation of the current density profile. The heat and particle sources, and the non-inductive current drive are computed by the TORBEAM code [31], and the NBI [32] and neutrals modules included in ASTRA. The TORBEAM and NBI modules require as input all the engineering parameters of the heating systems (e.g. power, angles, voltage of the beams, frequency of the gyrotrons, etc) and can predict the power distribution (and the particle source given by the NBI). The ASTRA neutrals module calculates the radial profile of the neutrals density taking as input the energy and the density of the neutrals at the separatrix. To account for sawteeth we consider the transport they cause in a 'sawtooth time averaged' way. We increase by a fixed amount the heat and particle transport coefficients (χ e,sr = 0.2 m 2 s −1 , χ i,sr = 1 m 2 s −1 , D n,sr = 0.5 m 2 s −1 ) in the plasma region inside of the inversion radius, so that the resulting kinetic profiles are a time averaged profile over the sawtooth period. We also use a model for the sawtooth crash [33] to take into account the effect of internal kink modes on the current density profile. The plasma equilibrium is calculated by the SPIDER [34] equilibrium code coupled to ASTRA [15]. Another important part of the model consists in a simple yet effective SOL description, composed by a set of analytical formulas, which set the boundary conditions (at the separatrix) of the electron and ion temperature and of the electron density. The MISHKA code [35] is then used to calculate the stability of peeling-ballooning modes (PBM) of the pedestals corresponding to the different values of the pedestal width. The final pressure profile, providing the result of the modeling workflow, is the one with the highest stable pedestal pressure, corresponding to pre-ELM conditions.
The modeling workflow is now computationally faster and more efficient, and it includes the modeling of the toroidal rotation, which was still taken from experimental measurements in our previous work. Experimental scans in operational parameters of critical importance have been included in this analysis, such as the plasma triangularity, the magnetic field, and the NBI voltage. The triangularity and the magnetic field scans address the robustness of the pedestal stability calculation (performed by MISHKA). As it will be shown, the NBI voltage scan provides a stringent test on the consistency among the different components of the integrated model, addressing in particular the reliability of the transport assumptions adopted to describe the pedestal, and their consistency with the MHD stability. We also explore the validity of the model on the QCE regime (quasi continuous exhaust, also referred as type-II and small ELMs) [36][37][38][39] and ITER baseline scenarios [40]. It is important to understand whether the model can reproduce the experimental behavior of these operational parameters, since these are potentially relevant for future fusion reactors. We analyze the accuracy of the model in reproducing the experimental pedestal and core confinement. Although the core and pedestal regions of the plasma are connected, and their interaction is fully taken into account in the modeling approach, we present detailed comparisons of the quality of the predictions for both core and pedestal separately. By this, we show that the modeling approach applied in this work in general allows predictions of the pedestal which are closer to the experimental observations than predictions for the core, somewhat reversing the common view that core transport is better predicted than the pedestal pressure. Section 2 describes the new features introduced in IMEP, explaining how it has been possible to increase the computational speed, and which changes have been made in the assumptions and physics description, as compared to the version of the model introduced in reference [27]. Section 3 presents the AUG experimental database considered for the extensive validation of IMEP, describing the variations of plasma parameters which have been considered. Section 4 presents the results obtained from this application, with a comparison of the prediction of the model with scaling laws and experimental data. Section 5 analyzes in details the results Figure 1. The new pedestal model included in ASTRA gives a transport constraint which determines the pressure at the pedestal top for a given pedestal width (a). Multiple parallel ASTRA simulations calculate the kinetic profiles, we show here the resulting total pressure (b) and electron temperature (e), from the magnetic axis to the separatrix (we show here the pedestal region only), for different values of the pedestal width. The corresponding pedestal gradients of pressure (c) and electron temperature ( f ) increase with increasing pedestal width, as the pressure increases more than linearly with the width, as shown in (a). This is achieved by imposing a value of the pedestal electron heat conductivity which decreases with increasing pedestal width (d).
obtained for each of the experimental scans studied. Summary and conclusions are presented in section 6.

Refined integrated modeling workflow
In comparison to the previous work [27], we introduced a few important improvements in the modeling workflow. Each ASTRA simulation is run on a single CPU with 16 cores, over which TGLF is parallelized. We reduced the radial gridpoints in the TGLF calculation from 80 to 48, allowing us to move from 5 to 3 serial calls, with a reduction of the computational time of one ASTRA simulation from about 3 h to 1 h and 50 min. We tested on many different experimental cases that the reduced radial resolution, employed for the calculation of the core turbulent transport fluxes, does not affect the final result of the simulation, leading to identical final kinetic profiles. For clarity, we recall that the ASTRA simulations are performed using a linear radial grid consisting of 601 points and a time-step of 2 ms.
We also switched to the more recent version of TGLF (saturation rule 1 geo with CGYRO units), with a new calibration of the saturation rule on CGYRO nonlinear simulations. The main difference we noticed by using the new version (for the range of parameters considered in this work) consists in the higher values of the ion heat flux at the pedestal top region. We reported in our previous work how TGLF was systematically underpredicting the ion heat flux in this region, which forced us to artificially increase the ion heat diffusivity at the pedestal top by a predetermined amount, and prevented us from fixing the boundary of the TGLF calculation further out (ρ tor,bc = 0.78, where ρ tor = √ Φ N is the normalized toroidal flux), since the further we went the lower was the predicted ion heat flux. The new version predicts a larger and more realistic value of the ion heat flux in the vicinity of this location, allowing us to set the TGLF boundary at ρ tor,bc = 0.85, thus decreasing the radial extent of the transition region (TR) between the pedestal and the core (TGLF domain), in which the diffusion coefficients are predetermined in order to get continuously regular profiles from the pedestal to the core region (as in the experimental profiles). By this the radial extent of the region with our ad hoc diffusivities is reduced, and the theory-based TGLF calculation has more authority on the final prediction of the kinetic profiles. The definition of the value of the diffusion coefficients in the TR is described in reference [27] (section 2.1.4).
We also included the effect of the dilution due to fast ions in TGLF, by simply removing the fast ion density from the thermal main ion density, although not treating the fast ions as a separated species in TGLF [41] (we include 3 species in TGLF: electrons, ions, and one impurity). This implies that quasi-neutrality is not respected in TGLF, although it is in ASTRA. We do this to at least take into account the effect of the fast ions on the turbulent fluxes via dilution [42]. This effect has shown to give negligible changes in most of our cases, probably because of the usually large thermal density values of AUG operation, which leads to a generally low ratio of fast ions energy to thermal energy (W fast /W th ). However, we observed that for cases with central deposition of the NBI power and with W fast /W th > 0.45 the effect can be considerable, causing a difference in the central ion temperature up to 1 keV and in the core thermal energy up to 15%, consistent with other recent applications of TGLF to AUG plasmas [43].
Another novelty introduced in this work is the modeling of the toroidal rotation. By doing this we calculate all the components of the radial electric field core profile (the poloidal rotation is calculated by NCLASS), which is used by TGLF to calculate the effect of the E × B shear on turbulent transport. The TGLF prediction for AUG plasmas is not very sensitive to the E × B shear, but by doing this we self consistently include the effect of the toroidal rotation in our simulations, without the need of using any experimental information from profile measurements. Since we only need the core component of the toroidal rotation, and the modeling of this quantity in the pedestal region is not trivial, we simply simulate it from the magnetic axis to the pedestal top. ASTRA allows us to set the boundary condition of the toroidal rotation directly at the pedestal top, while keeping the boundary condition of the other kinetic profiles at the separatrix. The toroidal rotation from the pedestal top to the separatrix is simply ignored, while its value at the pedestal top is not taken from experimental measurements, but is calculated using the formula from references [44][45][46][47][48], where f pass = 1 − 1.46 √ + 0.46 √ is the passing ion fraction given by equation (15) in reference [48], where = r/R 0 is the inverse aspect ratio. d c is a term of the diffusion coefficient in the ballooning form D(θ) = D 0 (1 + d c cos θ) and is assumed constant and equal to d c = 1.0,R X is the normalized major radius of the X-point, which varies from −0.46 to −0.32 across the cases considered in this work. The other terms in the formula are the effective safety factor q eff = 5 B t a 2 √ (1+κ 2 )/2 R 0 I p [49], where κ is the plasma elongation, the ion temperature at the pedestal top T i,top (keV), the ion charge state Z, the turbulence inhomogeneity length L φ (cm), which we assume proportional to the electron temperature scale length L φ = 0.5L T e . μ i is the ion atomic mass number (2 for deuterium), R 0 (m) the major radius, Q i (MW) the ion heat flux through the pedestal, and τ (Nm) the applied torque, by the NBI in this case. The first term in the formula is derived from a theoretical description of the intrinsic rotation based on a transport model that incorporates orbit-loss, transport-driven SOL flows, and turbulence intensity gradient. Such a model is able to capture the basic T i /I p scaling and X-point-position dependence of experimental observations. The second term takes into account the effect of applied torque to the edge rotation. We then simulate the core rotation by assuming an effective Prandtl number P r = 1 (cf references [50][51][52]), which means that we assume a momentum diffusivity χ φ equal to the ion heat diffusivity associated to turbulent transport χ i,turb , calculated by TGLF. The torque provided by the NBI is calculated by the ASTRA NBI module. All of the AUG cases considered in this work are heated by means of the NBI system. Therefore, we never face the problem of predicting an intrinsic core rotation profile. Thanks to the fact that an external torque is always present, the assumption P r = 1 provides a satisfactory description of the experimentally measured core toroidal rotation. We also changed the estimation of the radiated power, since the one previously applied in reference [27] has proven to be not reliable on the extended validation database, especially at low plasma current, and therefore low densities. This is because a fixed value of the tungsten concentration c W (the main source of radiation in AUG confined plasmas for cases without seeding of impurities) was assumed to evaluate the radiated power P rad = c W n 2 e L Z , where n e is the electron density, and L Z the tungsten cooling factor. In low density cases this frequently led to underpredicted radiated power when comparing with the bolometer measurements. Instead, we noticed that the measured radiated power in the database of experimental cases considered in this work, which does not include impurity seeded discharges, nor any case which exhibited impurity accumulation, robustly followed a linear dependence on the heating power (P heat ), with no dependence on the density. Physically, this could be connected to the increase of the divertor temperature with increasing heating power and the consequent increased tungsten sputtering yield. Therefore we adopted the simplified solution of capturing this behavior by assuming a radially constant tungsten concentration, and scale its value to obtain a total radiated power of P rad = 0.25(P heat − 2) (MW), obtained from a regression on experimental data. It is important to point out that this relation has been derived on stationary (over a few energy confinement time) H-mode phases of AUG discharges, therefore it might not be applicable to other plasma regimes or other devices. This also allows us to be independent from assuming the tungsten concentration. This assumption can be replaced in the future with theory based impurity transport modeling.
Regarding pedestal transport, we modified our description of pedestal particle transport, and now we make use of the electron neoclassical particle diffusivity and not of that of the ions (calculated by NCLASS). We recall that we assume the particle diffusion coefficient to be equal to D n,ped = c D/χ χ e,ped + D n,neo , where the term c D/χ χ e,ped represents the turbulent component of D n,ped , being proportional to χ e,ped through c D/χ . The value of the coefficient c D/χ = 0.03 was obtained through an optimization procedure trying to match a set of different experimental pedestal density profiles, as explained in reference [27] (section 2.1.2). Since in the presence of impurities D n,neo,e is smaller than D n,neo,i , we increased c D/χ = 0.06 (instead of c D/χ = 0.03) to compensate for this, in order to obtain a final D n,ped of similar magnitude. This difference has an almost negligible impact on the final results, but gives a more consistent description of the electron particle transport in the pedestal region, which, in this form, is also less affected by the assumed impurity content.
Finally, we implemented an automatic solver for the pedestal width, in order to efficiently find the highest stable pedestal width, without the need of running a fine scan on an extended range of pedestal width values. This consists in running first a rough scan with only 8 pedestal width values (Δρ tor = 0.097, 0.093, 0.089, 0.085, 0.081, 0.077, 0.073, 0.069). This results in finding a transition point from stable to unstable conditions between two different consecutive pedestal widths. At this point a finer scan in pedestal width is run with the highest accuracy (Δρ tor = 0.001). If 4 consecutive pedestal widths (each with a difference of Δρ tor = 0.001 from each other) are found such that the 2 lower ones are both stable and the 2 higher ones are both unstable, the iterative procedure has converged. The highest stable pedestal width is the final one. An example of such procedure is illustrated in figure 2, which depicts the normalized ideal growth rates of the MHD instabilities (γ MHD /γ Alfven ) and their spectra (for a set of toroidal mode numbers) as a function of the pedestal width (in Δρ tor ) and pressure. The circles (red) represent the first iteration, which finds a transition from stable to unstable, with a quite large gap in pedestal pressure (from ∼14 kPa to ∼18 kPa). This gap is filled by the second iteration, represented by the crosses (blue), which finds 3 more stable pedestal widths, reducing the gap and finding a pedestal width right below the stability limit (γ MHD = 0.04γ Alfven ). Since there are now 4 consecutive pedestal widths, 2 stable and 2 unstable, the procedure has converged. Figure 2 (right) shows how the most unstable toroidal mode number decreases in value with increasing pedestal width (and pressure).
This procedure has proven to be robust when applied on the database of the 50 cases presented in this work, allowing to efficiently run the whole modeling workflow and to make the most out of the computing resources available and without requiring human intervention.
If the highest stable pedestal pressure is found within the second iteration of the pedestal width solver, the whole workflow requires ∼ 6 h to simulate the plasma profiles corresponding to a time point of a discharge. This is achieved by running in parallel all the different ASTRA simulations, corresponding to different values of the pedestal width, on different nodes (CPUs) of the computing cluster, where every ASTRA simulation requires 1 h and 50 min on 16 cores of a single CPU (Intel Xeon Gold 6130). The computation time could be reduced by parallelizing the TGLF calculation (the slowest part of the ASTRA simulations) over a larger number of cores (e.g. 1 h and 15 min on 24 cores, 40 min on 48 cores), but we are currently limited by the 16 cores available on each CPU. Once the ASTRA calculations for a certain iteration of the pedestal width solver are finished, all the MISHKA calculations are run simultaneously, where not only the calculations corresponding to different pedestal widths are run in parallel, but a further parallelization is performed over the toroidal mode numbers. This allows us to run the whole MHD stability calculation in ∼1 h, for a given iteration of the pedestal width solver.
A further reduction of the total computational time required for the workflow execution can be achieved by accelerating the calculation of the core turbulent transport fluxes, which is at present the slowest process in the ASTRA simulations. For this purpose the TGLF model can be replaced by neural network surrogates [18,19,53]. Although the validity of these approaches on AUG parameters has not been addressed yet, a significant computational speedup is expected, which could bring the workflow execution time from ∼6 h to ∼2 h.

Experimental results considered for the extended validation of the model
We extended the validation of IMEP on a database consisting of 50 experimental cases, 10 from the previous work, which have been recomputed with the improved model described in section 2, plus 40 new ones, which include: a triangularity scan, a magnetic field scan, a plasma current scan at similar line averaged density, a NBI voltage scan (obtaining the same total heating power but with different particle source), scans in the heating systems (on and off axis heating, ratio of ion to electron heating), ITER baseline scenarios (2 discharges), and QCE scenarios (2 discharges). In comparison with our previous work, the extended database allows us to test the model on the variation of the magnetic field, which is another important parameter included in the IPB98(y,2) scaling law which is possible to test on a single device such as AUG. The entries in the database have been chosen in order to obtain the largest possible variation in the main engineering parameters, and in order to include some particularly relevant cases. This results in a variation across the database in plasma current I p = 0.6-1.2 MA, magnetic field B t = 1.5-2.85 T, safety factor q 95 = 3-8, heating power P heat = 2-14 MW, triangularity We selected a triangularity scan to address the robustness of the MHD stability calculation, since the change in pedestal pressure with different plasma shapes is caused by a change in the peeling-ballooning stability boundary. These effects have been shown to be captured by linear ideal MHD stability codes [54][55][56][57], like MISHKA, so our model should be able to capture this effect. Another important test is provided by the scan in the NBI voltage, which causes a change in plasma parameters in every plasma region: it changes the core particle source provided by the NBI, causing a change in particle flux crossing the separatrix which affects the amount of recycled neutrals in the SOL, ultimately driving a change in pedestal pressure. This NBI voltage scan provides a stringent test on the consistency among the different components of the integrated model, addressing in particular the reliability of the transport assumptions adopted to describe the pedestal, and their consistency with the MHD stability.
The values for the effective charge Z eff were taken from the integrated data analysis [58] when available (if previously calculated), however for most of the cases this was not available so for simplicity we assumed a fixed effective charge equal to Z eff = 1.2 (typical AUG value for experiments without seeding of impurities, in the intermediate to high density range). The result of the model is anyway very little sensitive to the effective charge in this range of values (Z eff = 1.1-1.5) relatively close to Z eff = 1.

Comparison with experimental measurements and IPB98(y,2) scaling law
In figure 3 we compare the measured and the predicted thermal energy, as well as its pedestal and core components. The experimental thermal stored energies have been calculated as the average between the values obtained using two different techniques. These are: (1) the integral of the experimental kinetic profiles obtained in the stationary phase before an ELM crash, and (2) W MHD − W fast , where W MHD is the plasma energy obtained from the equilibrium reconstruction, averaged over the multiple stationary phases before the ELM crash in a time interval of 0.1 s, and W fast is the energy of the fast ions, calculated by ASTRA. The error-bars represent the maximum and minimum values between the integral of the experimental kinetic profiles and all the data-points of W MHD − W fast within the time interval considered (retaining only the stationary phases before an ELM crash). This implies that, in the figures, the value of the experimental thermal energy in not necessarily centered at the middle of the error-bar. We observed that over the database the experimental values of W MHD agree well with W th + W fast , with W th calculated as the integral of the experimental kinetic profiles. This is also reflected by the rather small error-bars of the thermal energy. The core and pedestal experimental thermal energies are calculated as integrals of the pressure profile. Specifically, the pedestal stored energy W th,ped is the energy in the entire plasma volume, with pressure constantly equal to the pedestal top pressure from the pedestal top inward. The core stored energy W th,core is the total stored energy minus the pedestal energy W th,core = W th,tot − W th,ped . The same procedure is applied for the modeled and the measured profiles, assuming for both the same value of Z eff for the calculation of the ion density. The experimental pressure profile is constructed by fitting the measurements of the electron and ion temperatures, and of the electron density. The fits are obtained combining several diagnostics, in particular we used the two Thomson scattering systems (one for the core and one for the edge of the plasma) [59] and the electron cyclotron emission (ECE) system [60] for the electron temperature, the lithium beam emission spectroscopy diagnostic [61,62] and the deuterium-cyanide-nitrogen laser interferometer [63] for the electron density, and the charge exchange recombination spectroscopy systems [64] for the ion temperature. The edge Thomson scattering data is used to fit both electron temperature and density, together with the ECE and lithium beam data, so the profiles in every time window are automatically aligned to each other. The electron temperature and density are shifted together such that the separatrix electron temperature is ∼100 eV, as is typical of AUG [65]. The profiles are fitted using a modified hyperbolic tangent for the temperatures, and the two-line fit [66] for the density. We found the two-line fitting tool to be more robust for evaluating the pedestal top density, while the modified hyperbolic tangent can sometimes underestimate it due to the usually large scatter of the experimental data in this region of the plasma.
The accuracy of the prediction is evaluated using the mean relative error (MRE) where the index i represents a case in the database, and n is the total number of cases. The error is calculated between the predicted quantity y i,pred and the experimental measurement y i,exp . Figure 3(a) shows a comparison of the measured thermal energy with the prediction of IMEP (purple squares), and IPB98(y,2) (gray diamonds). In the experimental cases considered in this work, the time averaged plasma energy content calculated over multiple ELMs cycles differs from that just before the ELMs crash by less than 2%. This difference is negligible with respect to the error-bars of the measurements and with respect to the differences among the various predictions. It is therefore acceptable to compare the time averaged thermal energy predicted by the IPB98(y,2) with the thermal energy just before the ELMs onset predicted by the model and calculated from experimental data. The MRE for the model MRE = 5.78% suggests that the prediction is accurate in reproducing the change in energy confinement caused by the different parameters of plasma operation, and is much more accurate with respect to the IPB98(y,2) scaling law which has a MRE = 22.25% and which exhibits the tendency to overestimate the thermal energy. The upper and lower dashed lines in figure 3(a) represent a value of the relative error of RE = |y pred − y exp |/y exp = ±20% and help to better visualize the distribution of the relative error across the database. We underline that no data is taken from the measurements of the kinetic profiles in the ASTRA simulations. The model and the scaling law use similar input parameters, with the exception of the density, which is required by the IPB98(y,2), while it is predicted by IMEP, using the fueling rate as additional information. However, the integrated model makes use of empirical elements, in particular for the description of the pedestal and SOL regions, that could be AUG specific, and that the scaling law cannot benefit of. Figure 3(b) shows a comparison between the predicted and measured values of the thermal energy separated in the pedestal and core components. It can be noticed that the prediction of the pedestal pressure is very accurate and robust (MRE = 5.25%), while the core thermal energy exhibits more scatter (MRE = 12.16%). This means that the core energy is the term responsible for the largest portion of the error on the total thermal energy, somewhat reversing the common view that core transport is better predicted than the pedestal pressure. The core thermal energy tends to be overestimated at low energies, and underestimated at high energies. This could be due to a slightly inaccurate stiffness of TGLF [67]. Figure 3(c) shows a comparison between the predicted and measured values of n e at the pedestal top (blue triangles) and n e at the center of the plasma (red triangles), at ρ tor = 0.1. We chose this radial location instead of ρ tor = 0 since at the magnetic axis some of the experimental measurements are missing and fitted profiles result from extrapolations. An important aspect is that IMEP can correctly predict the pedestal top electron density. This is an improvement over the EPED model, which lacks the capability of predicting the pedestal top electron density, as this must be given as input. Figure 3(d) shows a comparison between the predicted and measured values of T e (orange diamonds) and T i (green squares) at the pedestal top, both of which are well reproduced by IMEP. Figure 4 shows again a comparison of the measured thermal energy with the prediction of IMEP (squares), and the IPB98(y,2) scaling law (crosses), but with color-maps which show the change in the main plasma parameters across the data set. It is clearly visible that the plasma current is the plasma parameter most strongly correlating with the thermal energy. This is because energy confinement has a strong dependence on this parameter (it scales almost linearly with I p ). The effect of plasma current is well captured by both the model and the scaling law. As one can see, the points for which the IPB98(y,2) is less accurate correspond to high fueling rate (Γ D ) levels, which can cause strong confinement degradation. This effect is not captured by the scaling law. Our model is instead capable of capturing this phenomenon, as reported in our previous work [27]. We observed that in our simulations a different fueling rate causes the peeling-ballooning stability boundary to be crossed at different locations (different toroidal mode numbers). The toroidal mode number (n) of the most unstable modes calculated by MISHKA, for the pedestals predicted by the model, becomes larger at higher fueling level. The transition from PBM stable to unstable conditions is located at high toroidal mode numbers (n 20) for high fueling cases (Γ D 2 × 10 22 e/s), therefore at the ballooning-dominated (low current and high pressure gradient) part of the stability diagram, illustrated in red in figure 5. The pedestals sitting in this region are limited to a lower pressure, since they are closer to the ballooning limit. At low fueling rate instead the stability boundary is crossed at the 'nose' of the diagram (corresponding to lower toroidal mode numbers and illustrated in green in figure 5), which is a more optimal region where higher pressures are allowed. This is an indication that the model is correctly capturing the effect of fueling on confinement, which usually causes a reduction in pedestal pressure, and that  the transport description is consistent with the MHD stability calculation, well reproducing the physics of this dependence.
We also noticed that the cases for which the peelingballooning stability boundary is crossed at the peelingdominated (high current and low pressure gradient) part of the stability diagram, corresponding to low toroidal mode numbers (n 3), are cases at low plasma current I p = 0.6 MA (q 95 ∼ 7), low fueling rate Γ D < 0.6 × 10 22 e/s, and high heating power P heat 8 MW. This combination of parameters produces pedestals with low density and high temperature, and therefore low collisionality, which maximizes the bootstrap current, and therefore causes the pedestal to be closer to the peeling limit, illustrated in blue in figure 5. The model accurately reproduces the experimental measurements also for the 5 cases in these conditions.
It is also visible that most of the points in the database have a value of the magnetic field B t = −2.5 T. This is because the AUG ECRH system is optimized for central heating (the ECRH power is deposited between 0 < ρ tor < 0.3 for most of the cases in the database) at such magnetic field value, which is therefore the one used for most of the experiments. It is therefore difficult to obtain a scan in B t , since this requires different settings of the ECRH system in order to achieve similar profiles of the power deposition at the two different B t values. Nevertheless, we found and included into the database two similar discharges providing a scan in B t , where no ECRH power was used, without being affected by tungsten accumulation thanks to the low value of the plasma current. Figure 6 shows that the thermal energy prediction of IMEP is even accurate for the two ITER baseline scenarios (cyan filled crosses) and the two QCE cases (yellow filled stars) considered. The pedestal pressure is also very well reproduced, with a relative error lower than RE < 10% for each case. This means that the model is robust also when applied on ITER relevant range of parameters (q 95 = 3) and that (at least for these cases) a linear ideal MHD stability code like MISHKA is capable of reproducing the correct pedestal pressure for the cases in the QCE regime. In our selection we chose a case with a high  The model proved to be capable of capturing the effects of these different operating parameters on pedestal (and global) confinement also for the QCE cases. The model is conceived to describe conditions close to the PB stability boundary, as usually expected for type-I ELMs. However, any regime that is not characterized by type-I ELMs, but which reaches pedestal values close to the PB stability boundary, can be in principle modeled. The fact that this approach correctly predicts plasmas in the QCE regime suggests that these cases are not far from the PB stability boundary. Also, these results show that the negative impact of high fueling can be consistently reproduced by IMEP in both the type-I and QCE regimes. From the modeling standpoint, this can be understood through the fact that, in both conditions, this effect is mainly connected with the radially outward shift of the peak pressure gradient, in regions where the safety factor is higher and thereby leading to limiting values of the α max,MHD parameter at lower maximum pressure gradients.
Another important experimental case included in this analysis is the stationary phase at 1.9 s of the discharge #35681, with B t = −2.5 T, I p = 0.6 MA, P heat = 7.5 MW, δ = 0.24, Γ D = 0.14 × 10 22 e/s. This experiment aimed at studying the effect of ELMs on the fast-ion densities at the plasma edge using a new edge FIDA spectrometer installed at AUG [68]. From a modeling perspective it is particularly interesting for two reasons: first, in the experiment tangential NBI sources were used, obtaining a strongly off-axis power deposition, and second, the edge FIDA spectrometer is capable of measuring the intensity of the edge D-alpha line simultaneously with the FIDA radiation, which allows the determination of the background density of neutrals. The first point provides a test of IMEP when using off-axis heating, which can have an impact on core confinement. The second point allows us to compare the model prediction of the density of neutrals with the one inferred from the experiment. The discharge has been performed with a considerable amount of ECRH power (P ECRH = 2.5 MW), and strongly off-axis NBI heating (P NBI = 5 MW), at relatively low electron density (n e = 4 × 10 19 m −3 ). This results in a high electron to ion temperature ratio (T e /T i ∼ 2) at the center of the plasma, which has a strong impact on core turbulent transport, destabilizing ion temperature gradient (ITG) modes [69][70][71][72]. As a result, this causes a decrease in core confinement. The IPB98(y,2) scaling law has been created including cases with dominant on-axis heating and without ECRH, and therefore it cannot be expected to correctly reproduce such effects on confinement. The relative error is very large for this case RE W th ,IPB98(y,2) ∼ 50%. IMEP instead can reproduce this effect, as this physics is captured by TGLF, well reproducing the electron to ion temperature ratio (RE T e /T i ,IMEP < 10%) and proving to be very accurate also for this case, with RE W th ,IMEP < 5%.
The experimental density of neutrals has been obtained combining the calculation of the KN1D code [73] and the FIDA measurements, as described in reference [74]. This sophisticated procedure yields to an accurate estimate of the density of the neutrals at the separatrix n 0,sep,exp ∼ 3.0 × 10 15 m −3 . IMEP predicts a value of n 0,sep,IMEP ∼ 2.2 × 10 15 m −3 . Given the simplicity of the model describing the neutrals in the SOL, this result is quite encouraging. This case has a very low fueling rate, meaning that recycling is probably the major component contributing to the density of neutrals at the separatrix. In this way we are directly testing the recycled neutrals term included in our model, without the need of disentangling the contributions of two different components (the one coming from recycling and that coming from the fueling valves) on the density of neutrals. Even though n 0,sep is slightly underpredicted by the model, the relative error on the pedestal top density prediction is very low RE n e,top ,IMEP ∼ 10% and within the uncertainties of the experimental measurements. Overall, these results are providing further confidence on both the core transport modeling and the description of the neutrals in the SOL, although the latter is approximated with quite simple assumptions.

Comparison with recent scaling laws
A new scaling law for the prediction of ITER plasma performance has been recently created, the ITPA20-IL (ITER-like) [3]: The new ITPA20-IL scaling law has been obtained from the ITER-like subset of the updated ITPA global H-mode confinement database, including plasma triangularity as variable.
The updated ITPA H-mode global confinement database contains new data with metallic wall from JET-ILW and AUG W (tungsten) wall and an increased number of observations which better envelope the ITER operational conditions. Important differences between the IPB98(y,2) and the ITPA20-IL are the strongly reduced dependence on the line averaged density, the decreased dependence on the major radius, the increased I p exponent and the slightly negative B T exponent, with the sum of these two exponents remaining fairly close around ∼1.1.  Figure 7 shows a comparison between the measured thermal energy with the prediction of IMEP (purple squares), IPB98(y,2) (gray diamonds), and ITPA20-IL (orange circles) scaling laws. The new ITPA20-IL scaling law is more accurate than the IPB98(y,2) on the selection of AUG experimental cases considered in this work, reducing the MRE from MRE IPB98(y,2) = 22.25% to MRE ITPA20−IL = 14.40%. In particular, it can be seen that the ITPA20-IL is very accurate at lower energies, while at higher energies its error also tends to be larger. A more detailed analysis hints that the error becomes larger with larger values of Γ D , as shown by figure 8 (left). This is consistent with the observations from figure 4, but this feature is even more evident for the new scaling law since the ITPA20-IL is more accurate than the IPB98(y,2) at low fueling rate, while the error of the two scalings is comparably large at high Γ D . This could be interpreted as a limitation of the ITPA20-IL in capturing the confinement degradation caused by an increase in Γ D , which cannot be directly described as a dependence on the plasma density, figure 8 (right).
Because the ITPA20-IL becomes less accurate at higher energies, and the thermal energy is strongly correlated with I p , one could argue that the reduced accuracy of the scaling law could correlate not only with Γ D but also with I p . However, figure 9 shows that the ratio between the measured and predicted thermal energy (H 98(y,2) and H 20−IL ) is more strongly correlated with Γ D for both the IPB98(y,2) and ITPA20-IL scaling laws. This analysis has been carried out on a subset of the database, selecting only cases with a triangularity lower than δ < 0.24, in order to decouple the effect of δ from Γ D and I p , which can have a strong effect on the H-factor. In figures 9(b), (c), (e) and ( f ) it is clearly visible that (especially for values of the plasma current larger than I p > 0.9 MA) at fixed I p the H-factor (H 98(y,2) and H 20−IL ) decreases with increasing Γ D . Instead, figures 9(a) and (d) show that at fixed Γ D the H-factor does not change with I p . Figure 7 also shows a comparison with the prediction of the AUG-2W scaling law [75], which is obtained performing a regression on the W wall subset of the AUG confinement database The AUG-2W scaling law gives very similar results to the ITPA20-IL, with the tendency to predict lower values of the thermal energy, and therefore resulting in a lower MRE MRE AUG−2W = 12.19%. The AUG-5C + W is another AUG only scaling law proposed in [75], containing also data from AUG operation with C (carbon) wall, and using as variables also the plasma density and triangularity. The predictions of the AUG-5C + W are very similar to the ITPA20-IL (MRE AUG−5C+W = 14.38%), and thereby it has not been included in the figure. Remarkably, IMEP, which is capable of capturing the effect on confinement caused by different fueling rate levels, results more accurate than all of the scaling laws considered, even those specifically developed on the AUG confinement database only. This provides further demonstration that an integrated model has the capability to capture hidden dependencies beyond the possibilities of scaling laws.

Accuracy of kinetic profiles prediction
Since a comparison of integral quantities, such as the stored energy, and local plasma quantities, such as the pedestal top density, does not give information on the quality of the prediction of the plasma profiles, two figures of merit are also adopted to describe how well the predicted profiles reproduce the fits of the experimental measurements. The two figures of merit are defined as the root mean square relative error (RMSRE) and the MRE. The RMSRE is calculated as where the index i represents a given radial location, and n is the total number of radial points, equidistantly distributed on the radial coordinate ρ tor . The error is calculated between the plasma quantities predicted by the model y i,pred and estimated by the fits of the experimental measurements y i,exp . The MRE is calculated as Figure 10 shows histograms of the RMSRE (top) and MRE (bottom) calculated for the T e (left), T i (middle), and n e (right) profiles corresponding to the 50 different experimental cases considered. The errors are calculated for the core (black) 0.2 < ρ tor < 0.85, which is the TGLF computational domain, avoiding the sawteeth region, and the edge (orange) 0.85 < ρ tor < 1 regions of the plasma. The core profiles of T e , T i , and n e exhibit similar low values of the RMSRE, which indicates the good  quality of the prediction, especially for n e . The edge profiles exhibit a higher value of the RMSRE. The large scatter of the experimental measurements and the steep gradients typical of the pedestal region cause difficulties in obtaining fits with a high level of accuracy. In particular, the evaluation of the pedestal width contains very large uncertainties. Therefore, the fits of the experimental measurements and the profiles predicted by the model can be radially misaligned to each other, and because of the steep gradients, this can result in large relative errors. This is the reason why the RMSRE is larger in the edge than in the core of the plasma, which does not necessarily mean that the predicted pedestal profiles are in bad agreement with the experimental measurements. The experimental pedestal top plasma quantities are instead very robustly estimated by the fits, due to the very low gradients typical of this region. This allows an accurate and reliable comparison between the predicted and measured pedestal top values. When analyzing the accuracy of the pedestal prediction it is therefore more important and reliable to compare the pedestal top quantities (T e , T i , and n e ), which are very accurately reproduced by the model, rather then the RMSRE between the predicted profiles and the fits of experimental measurements. However, this comparison on the edge part of the profiles is also shown for completeness. Finally, the pedestal thermal The MRE gives information if the model tends to overestimate (MRE > 0) or underestimate (MRE < 0) the plasma profiles. It can be noticed that all edge quantities are quite well centered around MRE = 0. For the core T i and T e profiles, the majority of cases exhibit a MRE > 0, which suggest that TGLF tends to overestimate the core temperatures profiles, especially T i . The core density profile instead, tends to be slightly underestimated, as already deduced from figure 3(c), partly compensating the overprediction of the temperature profiles, and producing therefore a MRE on the core thermal energy which is quite well centered around zero.
The quality of the toroidal rotation prediction is illustrated in figure 11, which shows a comparison between the measured and predicted rotation values at three different radial locations: pedestal top (blue), mid-radius (green), and center of the plasma (red). The prediction of v tor exhibits a larger scatter with respect to the other quantities previously analyzed, maintaining nevertheless a reasonable MRE. It is important to point out that the evaluation of the experimental value of the pedestal top rotation can contain large uncertainties, which can be comparable to the error of the predicted values. The scatter does not increase significantly towards the center of the plasma, meaning that the estimate of core momentum transport is quite accurate and robust. The prediction of the toroidal rotation profile has the only purpose of calculating the radial electric field in the core, without using any experimental information from profile measurements. This is needed to calculate the E × B, which is an input parameter required by TGLF, and can have an impact on the final prediction of the kinetic profiles. However, since TGLF is weakly dependent to the E × B in the parameter domain of these experiments, it is not important that the prediction of the toroidal rotation is very accurate. Therefore we can accept a relatively larger error on the predicted toroidal rotation with respect to the predictions of density and temperatures.

Analysis of specific dependencies
The following subsections describe in greater details the effects of changes in plasma triangularity, magnetic field, NBI voltage, and plasma current with similar line averaged density.

Triangularity scan
Among the cases included in the database, we focus on the two stationary phases of the AUG discharges #33195 and #33194, respectively at 4.2 s and 5.1 s, providing a triangularity scan (δ = 0.22-0.39) at B t = −2.5 T, I p = 0.8 MA, P heat = 11.4 MW, Γ D ∼ 0.9 × 10 22 e/s. An increase in plasma triangularity usually causes an increase in confinement which stems from the pedestal, due to a stabilization of the MHD modes [54][55][56][57]. In particular, the higher triangularity alters the shape of the peeling-ballooning boundary, allowing the pedestal to reach a higher pressure gradient and current density.
In our model the transport in the pedestal is not affected by triangularity, as it does not have any impact on the pedestal transport coefficients (apart from neoclassical transport which has a very small contribution). Therefore, without considering the MHD stability calculation yet, when comparing the ASTRA simulations results obtained at different triangularity for the same value of the pedestal width, we get the same pedestal pressure (also the TGLF results show a weak dependence on triangularity). This means that the relation between the pedestal width and the pedestal pressure does not change when changing triangularity. The change in confinement has then to be captured by the MHD stability calculation. Figure 12 shows the normalized spectra of the ideal growth rates (γ MHD /γ Alfven ) calculated by MISHKA for a set of toroidal mode numbers (1 n 40), corresponding to low (left) and high (right) plasma triangularity. The different colors correspond to different values of the pedestal width. The pedestal pressure and its gradient increase with the pedestal width, as explained in section 1, thereby also the growth rates increase. We selected the same range of pedestal widths values for this scan, in order to perform a better comparison. As one can see, for the same pedestal width the growth rates are strongly reduced at higher triangularity, especially for the higher toroidal mode numbers, which correspond to ballooning modes. This results in a higher value of the highest stable pedestal width at higher triangularity, and therefore a higher pedestal pressure. Figure 13 (right) shows that the predicted pedestal pressure is in good agreement with the measurements both at high and low triangularity, while the core thermal energy is overpredicted at high triangularity. As a result the change in total thermal energy predicted by the model is slightly overestimated, as it can be seen in figure 13 (left). The reason for this can be attributed to an underestimate of turbulent transport by TGLF for the case at high triangularity. The IPB98(y,2) captures the change in confinement due to the positive dependence on the density, which increases with triangularity [56]. The ITPA20-IL scaling law (not included in the figure) has a weaker dependence on density, but an explicit dependence on triangularity, and gives a prediction very similar to that of the IPB98(y,2), featuring the same trend but reducing the error from 27% to 21%. Figure 14 shows a comparison between experimental measurements (dots) and predicted profiles (lines) of electron and ion temperature and electron density at high (green) and low (purple) triangularity. As one can see, the values at the pedestal top are well reproduced, capturing the increase of both temperature and density with triangularity. The core gradients are very well reproduced at low triangularity (with the exception of T i inside ρ pol < 0.4, where ρ pol = √ Ψ N is the normalized poloidal flux), while at high triangularity the core gradients are overestimated by TGLF, especially for the electron and ion temperature. As a result TGLF does not predict a change in core profiles gradients as strong as the one observed from the measurements.
Overall, IMEP well captures the change in confinement with triangularity, especially for pedestal confinement.

Magnetic field scan
A scan in the toroidal magnetic field is provided by the two stationary phases of the discharges #31555 and #34955. In these two similar discharges with I p = 0.6 MA, P heat = 5 MW, δ = 0.22, the magnetic field B t is varied between 1.5 T and 2.8 T. The fueling rate is slightly different, being higher at lower magnetic field Γ D,1.5(T) ∼ 0.97 × 10 22 e/s, Γ D,2.8(T) ∼ 0.70 × 10 22 e/s. A scan in the magnetic field is another variation included in the IPB98(y,2) that is possible to test on AUG. The two stationary time windows selected represent the most similar pair of AUG discharges with the largest variation in B t we found. At these values of the magnetic field, it is difficult to get a similar distribution of the ECRH power, so we selected two cases without ECRH power, therefore with NBI heating only.
The measured thermal energy, showed in figure 15 (left), exhibits a negative dependence on the magnetic field W th ∝ B −0.45 t , that is not captured by the IPB98(y,2) scaling law, which instead predicts a slightly higher energy with increasing magnetic field. The IPB98(y,2) has a small positive exponent of B t , α B = 0.15. The experimental dependence on B t is consistent in sign with the ITPA-20IL scaling law (α B = −0.13), and is similar in value with that of the AUG-2W scaling law (α B = −0.344). IMEP predicts the correct trend, but does not predict a change in thermal energy as strong as observed experimentally (α B = −0.135). This is because the model, while it correctly captures the change in pedestal thermal energy (higher at lower B t ), predicts the same or even decreasing core thermal energy when going to lower magnetic field values, as it can be seen in figure 15 (middle), while the measured core energy is higher at lower B t values. Figure 15 (right) shows the predicted (solid) and measured (dashed) pressure profiles for the different magnetic field values, where it can be seen that the core gradients are underpredicted at low magnetic field (green lines), while the pedestal pressure is accurately matched. In particular, the prediction starts to deviate significantly from the measurements at ρ pol ∼ 0.6. The reason why TGLF fails in capturing this change in core confinement is most likely associated to the effect of β on turbulent transport. Since β has a quadratic dependence on B t , by reducing the magnetic field at constant pressure by almost a factor of 2, β increases by a factor of ∼4. In addition to this effect, the pressure decreases with increasing field, leading to an even stronger increase in β by reducing the magnetic field (β 1.5(T) = 2.1% ∼ 5 × β 2.8(T) = 0.42%).
In non-linear gyrokinetic simulations turbulent transport is strongly reduced by an increase in β, whereas in linear simulations this effect is weaker [76,77]. A comparison between linear gyrokinetic simulations, performed with the GENE code [78,79], and TGLF at mid radius (ρ pol ∼ 0.6), illustrated in figure 16 shows that TGLF underestimates the effect of β even against linear GENE. While the spectra of TGLF and GENE exhibit similar growth rates at high B t (low β), a large difference between the growth rates predicted by the two models is noticeable at low B t (high β). In particular, by increasing β, GENE (linear) predicts a reduction of the maximum growth rates by almost a factor 3 (γ 1.5(T) ∼ 0.37 × γ 2.8(T) ), while TGLF predicts a change in the maximum growth rate by only a factor ∼1/3 (γ 1.5(T) ∼ 0.66 × γ 2.8(T) ). The instabilities which are found by both codes at these wave numbers (k y ρ s ) are ITG modes. As already mentioned, the non linear effect of β on the growth rates (in particular for ITGs) is stronger than the linear one, and since in TGLF this effect is even weaker than in linear gyrokinetics, we can conclude that with our model we strongly underestimate the turbulence stabilization with β, and therefore the change in core confinement with B t is underestimated. Another element in support of this hypothesis is represented by figure 17, which shows the measured and predicted electron and ion temperatures and the electron density at the different magnetic field values. As one can see, the larger disagreement between the model prediction and the measurements is in the ion temperature for the low field (high β) case. Also the predicted core density profile is flatter than the measured one. Turbulent transport driven by ITG modes usually affects both temperature and density profiles, limiting their gradients.
The core density profile is not very sensitive to the particle source, and it is rather determined by turbulent transport [16]. These considerations suggest that TGLF is underpredicting the ITG stabilization caused by the increase in β. Thereby, the correct prediction of TGLF at high field shows that when the electrostatic limit is approached, the TGLF model is more reliable, whereas when β becomes large the electromagnetic stabilization is underpredicted.
Overall IMEP has shown to capture the effect of the magnetic field on confinement more accurately than the IPB98(y,2) scaling law, and consistently with the ITPA20-IL scaling law. The experimental dependence on B t is more consistent with the AUG-2W scaling law. This confirms the validity of this analysis. However, in this pair of discharges both the magnetic field and the fueling rate are varied. Since the fueling rate can have a strong impact on the pedestal pressure, we will present in a future work an experimental scan of B t at fixed fueling rate, to more completely assess the impact of B t on confinement.
The pedestal prediction has proven to be robust also with respect to a variation in B t , while the accuracy of the core thermal energy prediction is affected by a too weak stabilization of ITG modes by β in TGLF, which is responsible for a larger error on the total thermal energy evaluation at low field (high β).

NBI voltage scan
In our database we also included the stationary phases at 2.8 s of the discharges #35288 and #35289. In these two similar discharges with B t = −2.5 T, I p = 0.8 MA, P heat = 7.25 MW, Γ D ∼ 0.4 × 10 22 e/s, δ = 0.22, the NBI power is the same P NBI = 5 MW, but obtained with a different number of NBI sources at a different voltage (V NBI = 42 and 92 kV) [80]. The flux of particles across the LCFS for the two NBI configurations then differs by a factor of 2.5.
In these experiments the density remains almost the same, while one would have expected a change caused by operating the NBI heating system at different voltage, due to the resulting different particle source (as explained later in this section). We were then interested to see if IMEP can reproduce the experimental observation. These experimental cases provide an essential test of the model in validating the consistency between transport and pedestal stability and the connection between heat and particle transport. The prediction of IMEP can result accurate only if the change in the temperature and density profiles (and therefore also in confinement) caused by the different NBI voltage levels is well reproduced. We can test if the model contains the correct elements to achieve this. For this reason these cases have been chosen for a detailed analysis in this work. Another important point that motivated this analysis is that the IPB98(y,2) scaling law fails in reproducing the change in confinement observed experimentally, caused by the different voltage levels. It is therefore interesting to see if our model performs better than the scaling law by capturing these effects.
When operating at lower voltage the power provided by each beam is lower, so to deliver the same power, 8 beams are needed at V NBI = 42 kV, while only 3 beams are needed at  V NBI = 92 kV. Even though the power delivered by each beam changes with voltage, the amount of particles injected by each beam into the plasma remains unchanged. This means that for the lower NBI voltage case, in which more beams are used to achieve the same heating power, the particle source is significantly larger (by a factor ∼2.5). Therefore, one could expect that at lower voltage the density is higher, due to the higher particle source. Instead, as shown by the experimental data in figure 18(c), the density is similar for the two cases. An understanding of this experimental observation can be obtained with the integrated modeling. Going from high to low voltage, the experimentally measured pedestal pressure has decreased by ∼25%, due to the higher particle source, which causes an outward shift of the density profile. The electron temperature is lower at lower voltage, figure 18(b). The mechanism of the reduction in pedestal confinement is similar to that coming from an increase in the fueling rate, as presented in our previous work. The higher particle source from the NBI causes an increase in the particle flux leaving the plasma, leading to a higher density of neutrals via recycling. This causes a steepening of the density profile in the pedestal region, as observed from experimental measurements, and reproduced by the model. Figure 18(a) shows the scans in pedestal width, where the filled symbols correspond to PBM unstable conditions, while the open ones represent stable conditions. The largest open symbols identify the highest stable pedestal pressure, which correspond to the final result of the model for the pedestal prediction. As one can see, the predicted pedestal pressure is lower at lower NBI voltage. One can also see that the same pedestal pressure among the 2 different cases (low and high V NBI ) corresponds to different values of pedestal widths, and therefore the same pedestal width among the 2 different cases corresponds to different values of pedestal pressure. This is because, in the ASTRA prediction of the kinetic profiles, the higher particle source from the NBI (at lower voltage) causes an increase in the particle flux leaving the plasma, leading to a higher density of neutrals via recycling which causes a steepening (and therefore an increase) of the density profile. For the same value of the pedestal width, the electron temperature at the pedestal top is similar among the 2 different cases due to the constraint ∇T e /T e, top = −0.5 (cm −1 ), while the pedestal top density is higher with lower voltage (due to the higher particle source). Therefore, for the same value of the pedestal width, the pedestal top pressure is higher with lower voltage, or, for the same value of the pedestal pressure, the pedestal width is smaller with lower voltage. Since the peak of the pressure gradient is located approximately in the middle of the pedestal width, this moves outwards with lower voltage. Because the ballooning stability is sensitive to the location of this peak (the closer to the separatrix the more it is unstable due to the strong increase of the safety factor approaching the separatrix), the lower voltage case will be limited to a lower value of the pressure gradient, because of the lower ballooning stability limit. This is illustrated in figure 18(d), which shows the pressure gradient profiles of the highest stable pedestals for the 2 different cases. The legend of the figure also indicates the toroidal mode number (n) of the most unstable mode, which is higher at lower voltage (n = 40), meaning that it is more ballooning unstable.
In the pedestal width scan the pedestal pressure increases more than linearly with the width, in fact also the pressure gradient increases with the width. This comes from the imposition of the constraint ∇T e /T e, top = −0.5 (cm −1 ), which causes the pedestal electron heat diffusivity χ e,ped to decrease in value with increasing pedestal width, as already explained in section 1. Because of this, and because the highest stable pedestal width is lower at lower voltage, the value of the electron heat diffusivity χ e,ped is larger at lower voltage, and so is the particle diffusivity (since D n,ped ∝ χ e,ped ), as it can be seen in figures 18(e) and ( f ). The particle diffusivity becomes larger by a factor ∼2 at lower voltage, compensating for the almost double source of neutrals (the value of the separatrix neutral density is shown in the legend of figure 18(c)), and resulting in a similar pedestal top density.
Another effect to take into account is the different power deposition of the NBI with different voltage. The beams are absorbed more in the edge of the plasma at lower voltage, as it can be seen in figure 19 (left, solid lines). This effect is consistently included in the simulations. The profiles in the figure are calculated with the ASTRA NBI module. Figure 19 (left) also shows the electron (dotted) and ion (dashed) flux surface crossing powers (respectively Q e and Q i ). The different beam deposition could also contribute to the difference in Q e and diffusivity in the pedestal region at the different NBI voltage levels (higher at lower voltage). The difference in the flux surface crossing power profiles is instead small in the core. We have studied this effect directly by repeating the simulations and replacing the power density profiles from one case to the other. We have found that the different NBI power deposition alone does not have any impact on the pedestal pressure. Therefore, we conclude that this is not an important effect. Figure 19 (right) shows the radial profiles of the particle source provided by the NBI (solid) and by the SOL neutrals (dashed), produced mainly by recycling, for the two different NBI voltage levels. The NBI particle source changes significantly with the different voltage levels (by a factor ∼2.5). As a consequence, also the source provided by the SOL recycled neutrals changes (by a factor ∼2). It can be noticed that the differences in the particle source profiles, produced by the different NBI voltage levels, are significantly larger than the differences in the heating power densities (note the logarithmic scale). We have also studied in more detail this effect by repeating the simulations and replacing the particle source from one case to the other. We have found that the different particle source alone has a strong impact on the pedestal pressure. We can therefore conclude that the change in confinement is caused by the different particle source, which strongly impacts the pedestal stability. Figure 20 shows that the predicted pedestal energy matches the measurements, and IMEP well captures the change in pedestal confinement with the different NBI voltage levels. The change in core energy, and then in total thermal energy as well, is also well captured. The prediction of the IPB98(y,2) does not change with the different NBI voltage levels, since all its input parameters (including the line averaged density) do not change, strongly overestimating the plasma energy at low voltage.
A more complete and specific documentation of these simulations, with additional experimental cases, will be reported in a separate publication [80]. However, in the framework of the present paper, this case provides a demonstration of how important it is to take into account core, pedestal, and SOL effects self-consistently. The combined estimation of core particle transport and sources with the different NBI voltages affects the SOL neutral density via recycling, which (the neutral density) impacts the pedestal MHD stability. This allows us to capture the change in confinement caused by the different NBI voltage. This analysis is also a very stringent test for the assumption on the pedestal transport, and on the description of the neutrals in the SOL, which can only result accurate with a correct representation of the underlying physics. It would also not be possible to obtain such results without assuming that the pedestal particle diffusivity is proportional to the pedestal electron heat diffusivity (D n,ped ∝ χ e,ped ). This suggests that our assumptions on pedestal transport and its consistency with the pedestal stability give a correct representation of these phenomena, at least in the modeling of these experimental conditions. Therefore, the present modeling results can also be considered to provide hints on the properties that a theorybased pedestal transport model should have in order to also reproduce these experimental observations.

Current scan at similar line averaged density
The large database considered for the extended validation of IMEP allowed us to also analyze a scan in plasma current at similar line averaged density, instead of fixed fueling rate as in the case presented in our previous work. The scan is represented by stationary phases of the discharges Figure 19. Radial profiles of the NBI power deposition (solid), and of the ion (dashed) and electron (dotted) flux surface crossing powers for the two different NBI voltage levels (left). Radial profiles of the particle source provided by the NBI (solid), and by the SOL neutrals (dashed), produced mainly by recycling, for the two different NBI voltage levels (right). #33195 and #30668 corresponding to a change in plasma current from I p = 0.8 MA to I p = 1.2 MA, respectively. The two discharges have otherwise the same parameters B t = −2.5 T, P heat = 12 MW, δ = 0.3, and similar line averaged density, being however higher at higher currentn e, 0.8 (MA) = 5.5 × 10 19 m −3 ,n e,1.2 (MA) = 6.6 × 10 19 m −3 . The corresponding Greenwald fractions f G =n e /n G , where n G = I p /(πa 2 ), are f G,0.8(MA) = 0.57 and f G,1.2(MA) = 0.45, respectively. As already discussed in our previous work, an increase in plasma current at fixed fueling rate produces an increase in the pedestal top density, due to the improved PBM stability. Therefore, to achieve a similar density, a higher fueling rate is required at lower current Γ D,0.8(MA) = 0.95 × 10 22 e/s, Γ D,1.2(MA) = 0.51 × 10 22 e/s. This scan allows a more appropriate comparison between the model and the scaling law on the accuracy in capturing the effect of the plasma current, while keeping all the other regression parameters of the IPB98(y,2) fixed, or within a small variation as for the line averaged density. However, as discussed in section 4, a change in fueling rate, such as the one required in these discharges to produce similar densities at different currents, is sufficient to produce a strong decrease in pedestal confinement. As a consequence, the decrease of global confinement in the plasma at low current is stronger than that predicted by the IPB98(y,2), which significantly overestimates the thermal energy at low I p (and high Γ D ), with H 98 = 0.7, as it is shown in figure 21. In contrast IMEP is capable of capturing the combined effect of the fueling rate and the plasma current on the pedestal pressure and on core confinement, leading to a prediction of the plasma stored energy with a significantly lower MRE MRE = 6.61% than the IPB98(y,2) scaling law MRE = 20.12%. The change in pedestal top density is also well reproduced by the model. These experimental results exhibit a more than linear dependence on I p . Assuming that the difference in the fueling rate is responsible for 15% of the change in the stored energy, as observed in fueling scans with the same variation in Γ D [27], one is left with a ratio of the thermal energy at high and low current of W th,1. These considerations in the comparison with recent scaling laws suggest the potential generality of this integrated modeling approach.

Summary and conclusions
In this work we have shown an extensive validation of IMEP, the integrated model presented in reference [27] which combines theory-based components and empirical elements to simulate the confined plasma, with some important updates which make it computationally faster and more accurate. The new improved version of the integrated modeling workflow does not require any experimental information from profiles measurements as input, as now also the toroidal rotation is selfconsistently predicted, while for the radiated power we use a simple but realistic general empirical scaling for stable Hmode phases at AUG. The only assumed parameter is Z eff , a quantity on which the result of the model is not very sensitive, at least in the range of relatively small values of Z eff obtained in a metallic wall device like AUG. In fact, Z eff has a limited variation in the AUG (with tungsten wall) database considered in this work. Thereby, similar results can be obtained assuming a fixed value of Z eff . This aspect is planned to be improved in the future by including a self-consistent modeling of the impurities, similarly as done in reference [19].
The validation of the model has been extended to a database of 50 experimental cases, which largely cover the AUG operational domain. The prediction of the model is compared with experimental measurements and with the scaling law of the ITER physics basis IPB98(y,2), regularly applied over the last two decades to qualify confinement in present experiments. Comparisons are also performed with the most recently derived scaling law ITPA20-IL from the international multidevice tokamak confinement database specific for ITER like plasmas, and with AUG specific regressions, made on the AUG confinement database only. It has been shown that the prediction of the model is very accurate and robust over the entire set of 50 cases, with a significantly lower MRE with respect to each of the scaling laws considered. This shows that indeed a sophisticated model which contains a description of the physics regulating plasma confinement can capture important dependencies beyond the possibilities of a log-linear regression of experimental data. In particular we have noticed that, for the AUG experimental cases considered, all of the scaling laws tend to be less accurate at higher values of the fueling rate (Γ D ), as they do not capture the confinement degradation caused by this parameter. The fact that IMEP is instead capable of capturing the effect of fueling correctly is a very important aspect for the study of the scenarios for ITER and future fusion reactors. Power exhaust sets constraints on the possible variation of the fueling rate, therefore it is important to take into account the effect that the operating conditions have on the fusion performance, by including it in the simulations.
Overall, the analysis carried out in this work shows that the model can accurately reproduce the change in confinement caused by a change in plasma current, heating power, fueling rate, triangularity, and magnetic field. Some important experimental cases have been included in the data set, which has been considered in this work to provide specific stringent tests for the model. In particular, the model has proven to reliably capture the effect of different heating mix (e.g. correctly predicting the effect on confinement caused by a change in T e /T i ), and of different NBI voltage (i.e. the effect of a change in the core particle source). This modeling framework has proven to also correctly describe plasma confinement for two particular operating regimes, such as the ITER baseline scenario, and the QCE scenario. This shows that plasmas in the QCE regime are not far from the peeling-ballooning limit provided by linear MHD stability codes, such as MISHKA, allowing accurate predictions also in this plasma conditions, at least for the representative cases considered here.
The pedestal transport model we presented can robustly capture the effect of the different plasma parameters on the pedestal pressure. It also brings important advantages, as it can accurately predict the pedestal top density with no need of experimental information, as opposed to the widely applied EPED model where this must be given as input. This increases the predictive capability of the entire integrated modeling approach. It also brings for the first time (to our knowledge) the capabilities of modeling separately the electron temperature and the ion temperature profiles in the pedestal region (without assuming T e,ped = T i,ped ). Finally, the inclusion of the formulas from references [44][45][46][47][48] into the model allow the prediction of the pedestal top toroidal rotation with reasonable accuracy.
Comparing the quality of the prediction on the pedestal and core thermal energy components, we found that the larger error is associated to the core. This could be surprising, as core transport is considered to be better understood and predicted than pedestal confinement. This does not mean that the description of core transport is performing worse than expected (the error of the predicted core energy is small), but rather that the empirical pedestal transport model included in this work provides a sufficiently realistic description and that the MHD stability provides a very robust and accurate constraint. The successful application of the empirical pedestal model could therefore give guidelines in the understanding of the actual physical properties of transport in the pedestal. We have seen that TGLF tends to overestimate transport for cases at high β. In particular, TGLF can underestimate the reduction of the ITGs growth rates when increasing β compared to linear gyrokinetic simulations, and it also does not include the nonlinear effect, which is usually considerably larger. The pedestal transport model is based on experimental observations that we expect to be applicable to other devices since it relies on a multi-machine analysis which has identified a common parameter. However, because such common parameter is dimensional, the validity of this assumption needs to be tested for machines with a different size than AUG. IMEP contains elements which could be specific and applicable only to AUG. The main aspects which produce a strong machine dependency are the geometry of the divertor and its baffles, the location of the gas valves, and the materials of the plasma facing components, which would probably make the estimation of the divertor neutral pressure not valid for different machines or divertors. A new scaling should then be derived for the tokamak of interest, using either experimental measurements if available, and/or synthetic data from simulations, particularly for non existing devices. Codes for the simulation of the SOL profiles are capable of predicting the divertor neutral pressure, therefore it would be interesting to study if a small database of such simulations would be sufficient to obtain a scaling for divertor neutral pressure p 0 , based on synthetic data rather than experimental measurements. If this approach results successful, leading to a similar accuracy of the prediction for AUG cases, the validation of the integrated model could be extended to other tokamaks. Testing the model on a larger device, such as JET, and on a smaller one, e.g. TCV, would allow us to also obtain a scan in the size of the device, which, if successful, would increase the confidence in the prediction for ITER, for which a large database of SOL simulations is already available, and therefore a dependence of p 0 on the main engineering parameters could be extracted.
The capability of simulating the kinetic profiles of the confined plasma increases the accuracy and the reliability in the prediction of energy confinement with respect to 0D scaling laws, as a 1D model can include the description of the physics phenomena which are strongly dependent on the gradients of the kinetic quantities, both in the core and in the pedestal regions. In particular we have shown that the density profile affects the shapes of the pressure gradients, which have a strong impact on the pedestal stability. Moreover, the prediction of 1D profiles is important for a reactor, since the shape of the kinetic profiles can influence the design and the prediction of the operation of a burning plasma [22], and therefore this is a great advantage that the integrated model provides over 0D approaches.
This work clearly demonstrates that it is possible to combine many different components for the description of plasma confinement into an integrated modeling workflow, which covers the entire confined plasma cross section within a computationally easily tractable framework and produces very realistic results. For the first time it has been proven that a 1D modeling approach can reproduce experimental results over very large parameter variations, as allowed by a single device, with a higher accuracy than any statistical regression, even those performed on the device itself only. The approach of integrated models, like the one presented in this work, has therefore the potential to improve the prediction of the fusion performance in future tokamak reactors.
Finally, the integration of the different modules can provide important insights to better understand interdependencies, particularly between different plasma regions, which are not possible to explore otherwise. Thereby, in addition to the increased predictive capabilities, which are promising also for applications to other devices, this approach has proven very helpful in the identification of hidden dependencies, specifically those resulting from effects that connect the different plasma regions (from SOL to core) which cannot be identified in 0D statistical studies nor in physics studies which focus on specific plasma regions only.