Incorporating nonlocal parallel thermal transport in 1D ITER SOL modelling

Accurate modelling of the thermal transport in the ‘scrape-off-layer’ (SOL) is of great importance for assessing the divertor exhaust power handling in future high-power tokamak devices. In conditions of low collisionality and/or steep temperature gradients that will be characteristic of such devices, classical local diffusive transport theory breaks down, and the thermal transport becomes nonlocal, depending on conditions in distant regions of the plasma. An advanced nonlocal thermal transport model is implemented into a 1D SOL code ‘SD1D’ to create ‘SD1D-nonlocal’, for the study of nonlocal transport in tokamak SOL plasmas. The code is applied to study typical ITER steady-state conditions, to assess the relevance of nonlocality for ITER operating scenarios. Results suggest that nonlocal effects will be present in the ITER SOL, with strong sensitivity in simulation outputs observed for small changes in upstream density conditions, and drastically different temperature profiles predicted using local/nonlocal transport models in some cases. Global flux limiters are shown to be inadequate to capture the spatially and temporally changing SOL conditions. Introducing impurity seeding, under conditions where detached divertor operation is achieved using the flux-limited Spitzer-Härm models used in standard SOL codes, simulations using the nonlocal thermal transport model under equivalent conditions were found to not reach detachment. An analysis of the connection between SOL collisionality and nonlocality suggests that nonlocal effects will be significant for future devices such as DEMO as well. The results motivate further work using nonlocal transport models to study disruption events and low collisionality regimes for ITER, to further improve accuracy of the nonlocal models employed in comparison to kinetic codes, and to identify more appropriate boundary conditions for a nonlocal SOL model.


Introduction
An important limiting factor for the design of future tokamak fusion devices is the thermal heat flux onto the exhaust target plates. Modelling of the tokamak Scrape-Off Layer (SOL) Original Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
for the next major fusion experiment, the ITER tokamak, predicts that the unmitigated steady-state heat load on the divertor target plates will be of the order of 40 MWm −2 [1], which will have to be mitigated to 10 MWm −2 -the accepted upper limit that the materials can withstand [2]. This heat flux challenge will be even greater for future pilot power plants such as DEMO [3].
Accurate modelling of the thermal transport in the SOL is therefore of vital importance, both in accurately predicting the heat flux through the SOL and onto the divertor targets, and in assessing the design and effectiveness of heat flux mitigation methods such as detachment. However, it has long been known that in the presence of steep temperature gradients, classical local transport theory breaks down [4,5]. The thermal transport becomes 'nonlocal', depending on conditions in distant regions of the plasma, where this nonlocal kinetic effect on thermal transport has significant impacts on the temperature and heat flux profiles [4,6]. The necessary conditions can be found in SOL plasmas (as well as in laser heated plasmas relevant to inertial confinement fusion (ICF) where this is a long standing problem [4][5][6][7]), where temperature gradients can become steep approaching the divertor target plates [8]. Nonlocal effects could therefore potentially have significant impact on plasma conditions across the SOL and the divertor target, with implications for the target plate conditions and access to detached divertor operation. Accurately capturing nonlocality in thermal transport models is therefore a major challenge to be addressed in predictive modelling of the SOL.
Many methods have been proposed to capture nonlocality without resorting to full kinetic modelling, which is not feasible for realistic simulations required to assess the performance of future devices (although some progress is being made in this area [9][10][11]). The most rudimentary and commonplace is the use of flux limiters on the electron heat flux in plasma fluid models -tuned limitation of the heat flux when the local Spitzer-Härm model predictions become unphysically large. However, these do not account for all the kinetic effects of nonlocal transport such as 'pre-heat' (where long mean free-path electrons stream ahead of the main heat fron and heat distant regions of the plasma), and the accuracy and predictive capability of flux limiters is questionable. The use of a flux limiter calibrated to smaller-scale ICF experiments was one reason for the discrepancy observed in the thermal transport in the laserheated plasma on the National Ignition Facility from that predicted by simulation codes [7], and the models continue to be revised [12]. Several more advanced nonlocal transport models have been developed to attempt to address such discrepancies [13][14][15][16]. These have had varying success at reproducing results from kinetic simulations [17][18][19], but have not yet been implemented into large scale complex SOL simulation codes.
One such nonlocal thermal transport model, first developed by Ji, Held and Sovinec [16] and then further developed by Omotani [20], takes an analytical approach to solve the electron drift kinetic equation to determine a model for the heat flux [16]. This model uses a moment expansion of the kinetic equation in a polynomial basis, truncated to a specified number of moments. The resulting set of equations are solved as an eigenvector problem to produce moment equation solutions, of which the 3rd-order corresponds to the thermal transport. This model has been tested and compared with both kinetic and local flux-limited thermal transport models [17], and whilst the model has some limitations -for example being unable to predict preheat effects -it has been shown to have the advantage of reproducing flux-limited cases of varying magnitude flux-limiter in time-dependent simulations over a range of plasma conditions [17] without requiring hand-tuned fitting parameters.
In this paper we describe the 1D SOL code, 'SD1D', that has previously been implemented by Dudson et al [21] using the BOUT++ framework and our incorporation of the Ji-Held nonlocal model for approximating the nonlocal heat flux. This version of the code, referred to here as 'SD1D-nonlocal', has been used to study the impact of including nonlocal thermal transport on the output of 1D SOL simulations for the ITER tokamak. Simulations employing the nonlocal model were then compared to those using local Spitzer-Härm and fluxlimited (FL) heat flux models, to analyse the implications of any differences in SOL and reactor target plate conditions. This paper is structured as follows: section 2 describes the nonlocal Ji-Held model, and the previous work that has been performed with this model. Section 3 describes the combined SD1D-nonlocal model, section 4 discusses some simple model benchmarking work performed and section 5 then analyses results of this code for a 1D ITER-like simulation case, comparing results against local and flux-limited models. A density scan is then performed as a method of altering the collisionality of the plasma being studied, with results shown in section 6. Section 7 discusses an analysis of the nonlocality of a tokamak SOL in terms of the connection length and electron mean free path, and in section 8 a fixed-fraction impurity model is introduced to both the local and nonlocal model simulations. Discussion of the results and further work is given in section 9, and conclusions are outlined in section 10.

Ji-Held nonlocal thermal transport model
An outline of the derivation for the Ji-Held nonlocal thermal transport model is given here; the full derivation can be found in References [16,20]. The electron distribution function f e is expressed as the sum of a Maxwellian part f (0) e and non-Maxwellian part δf e , such that f e = f (0) e + δf e . δf e can be thought of as the nonlocal contribution -the distortion of the distribution away from Maxwellian caused by nonlocal electron transport. The drift kinetic equation for a collisional plasma can then be written as: where ⟨·⟩ denotes the gyroaverage, v || is the component of the velocity parallel to the magnetic field, l is the distance along the field line and C(·,·) is the linearized Fokker-Planck collision operator. This equation is expanded in fluid moments on the basis ) are tensor Harmonic polynomials and L

Te
) are associated Laguerre polynomials. The polynomial basis is typically truncated to a high order (e.g. 20 × 20), reducing the kinetic equation to a set of first-order ODEs for non-Maxwellian fluid moments: where z is a dimensionless length defined by ∂z ∂l = λ ee −1 (where λ ee is the electron-electron collision length), n A are the parallel fluid moments of δf e , Ψ A B is the matrix coming from v || , and g A is the drive term from the gradient of the Maxwellian part.
Equation (2) is then transformed to the eigenvector basis wheren A andĝ A are the components of n A and g A on the eigenvector basis and ζ(A) are the corresponding eigenvalues, and the equations are solved as an eigenvalue problem. Transforming back to the original basis and selecting the (1,1) moment (giving the heat flux) gives the nonlocal thermal transport model: where W A B is the matrix formed from the eigenvectors W (B) (whose eigenvalues are ζ (B) ) connecting equations (2) and (4). This model has been previously studied in a 1D SOL model, for a test case representative of a edge-localised mode (ELM) crash in the JET tokamak SOL [20]. The model was shown to be able to self-consistently calculate the heat flux reduction based on plasma parameters, and respond to changes in the plasma conditions over time. It was observed that over the course of an ELM-crash the degree of flux limitation varied by up to two orders of magnitude. The flux limitation calculated by the model has also been shown to be in reasonable agreement with results from a kinetic code across a range of collisionality regimes [17]. However, these studies did not capture any of the important atomic and neutral physics that would likely have a significant impact on the realised heat flux, which become increasingly important when modelling high-power tokamaks such as ITER. In the following section, we shall introduce the SD1D code, which is capable of taking these effects into account.

SD1D-nonlocal
The SD1D model [21] in BOUT++ evolves standard 1dimensional Braginskii plasma fluid equations: where n is the plasma electron density, p is the total plasma pressure, V || is the plasma flow velocity, T e /T i are the electron/ion temperatures respectively, q is the heat flux (with former and latter terms in equation (6f ) being convective (q conv e|| ) and conductive (q cond e|| ) terms respectively), κ ||e is the electron heat conduction coefficient, j || is the plasma current density, andb gives the unit vector for the direction of the magnetic field. The model assumes equal electron and ion temperatures, such that only a single set of fluid equations are required to model the plasma species. Source terms in the equations are given by S n and S p for the particle and pressure sources respectively representing cross-field transport from the core into the flux tube of interest. A fluid-diffusive neutral model is employed, with neutrals being evolved with their own equivalent set of fluid-diffusive equations similar to equations (6a)-(6c) [21]: where n g , p g , T g and V g are the neutral gas density, pressure, temperature and velocity respectively, D g is the neutral diffusion coefficient given by and κ g is the neutral gas heat conduction coefficient, defined by where v g,th = √ eT g /m i is the neutral thermal velocity, ν cx is the charge-exhange frequency and ν gg is the neutral-neutral collision frequency. The V g is an effective parallel velocity, given by the sum of a parallel flow and parallel projection of a perpendicular diffusion [21]: Transfer channels in the equations are given by S, R, E and F terms, where S is the net recombination (i.e. plasma sink/neutral source), R is the radiation energy losses, E is the energy transfer to neutrals (including recombination, ionisation, charge exchange and elastic collision terms), and F is the ion friction losses from charge exhange and recombintation. The inclusion of transfer channels that capture energy/particle transfer between the plasma, neutrals, impurities, and encompassing processes such as recombination, ionisation, charge exchange, recycling etc, make this a complex SOL model that goes beyond that in which the Ji-Held nonlocal transport model has previously been studied in [20].
SD1D models half the tokamak SOL, from stagnation point (plasma flow = 0, assumed to be at the outboard midplane) to the divertor target. A source region is defined between the upstream stagnation boundary and a specified X-point location, defining the spatial extent of source terms S n and S p . The S p term is fixed to a constant across the source region such that the integrated power flux equals a specified parallel power flux density at the X-point. The particle source S n can either be set to a specified constant input particle flux into the domain, or alternatively can use a PI feedback controller to set the upstream density at s || = 0 to a specified value.
A symmetry (zero flow) boundary condition is applied at the upstream stagnation boundary, corresponding to: At the downstream boundary a standard plasma sheath edge boundary is employed: q || = γ SH T e n i V || (13) where γ is the ratio of specific heats and γ SH is the sheath heat transmission factor (assumed to be 5/3 and 7.8 respectively in this study) [22]. For our study, we compare the outputs of SD1D simulations using three thermal transport models: Spitzer-Härm, 'flux-limited' Spitzer-Härm, and the Ji-Held nonlocal model.

Spitzer-Härm
The heat flux calculation in the standard SD1D model is calculated using equation (6f ), where the latter term is the local Spitzer-Härm model for diffusive conduction [23], with κ ||e given by where ln(Λ) is the Coulomb logarithm.

Flux-limited Spitzer-Härm
When temperature gradients become large, the Spitzer-Härm predictions for the heat flux can start to greatly exceed the 'free-streaming flux' q fs = n i T e V Te (where V Te is the electron thermal velocity), seen as an approximate physical limit. A simple correction that is widely employed is the 'flux-limited' model, where the maximum parallel heat flux predicted by Spitzer-Härm is limited to a specified fraction α of q fs using the equation The flux-limited model originates in laser-plasma and inertial confinement fusion studies [7], but have also been applied to tokamak SOLs in modelling low collisionality regimes [24,25]. To implement a flux-limited model into the SD1D code, the thermal conductivity in the heat flux calculation is adjusted to an 'effective thermal conductivity' term, which is mathematically equivalent to equation (15).

Ji-Held nonlocal model
To create a new model with nonlocal thermal transport, the heat conduction term (second term in equation (6f )) was replaced with the nonlocal heat flux given by equation (5), to create SD1D-nonlocal. Boundary conditions in the nonlocal model are adapted to match those in Equations (12) and (13), which impose this condition on the boundary moments. A notable difference in how the code fills the guard cells at the sheath boundary in this case is that the local and FL models employ a zero-gradient Neumann condition for T e across the boundary, whereas the SD1D-nonlocal code employs a constant-gradient Neumann condition which is instead applied to log(T e ), preventing negative guard cell T e in the case of low target temperatures. In both cases the Neumann conditions are 'free-floating' conditions that do not impact on the boundary/domain values of T e .

Benchmarking SD1D-nonlocal
Convergence testing and comparison with kinetic code results for the Ji-Held model and Omotani's implementation of it has been performed in previous studies with the model [16,17,20]. As a check on the successful integration of Omotani's Ji-Held model implementation into SD1D, benchmarking simulations were performed for both high and low collisional SOL scenarios, comparing the outputs using the Ji-Held heat flux model with the standard Spitzer-Härm. SD1D-nonlocal simulations were run on a 1D grid with 320 evenly-spaced gridcells for a domain with connection length L = 100 m from stagnation point (s || = 0) to divertor target, with X-point location at s || = 80 m. The Ji-Held simulations were performed using 20 × 20 moments in the Legendre-Laguerre basis, which has been shown to be sufficient for convergence of the nonlocal model output [17,20]. Simulation parameters were set such that the particle source produced an upstream density of~3.7 × 10 19 m −3 , and power source such that the parallel power flux density at the X-point was 0.24 GW m −2 . High recycling fraction of 0.95 was applied at the target boundary,  and the sheath transmission factor set to 7.8. Such conditions are similar to those of 2D simulations of the ITER SOL for the 3rd/4th grid cell ring outside the separatrix for 80-100 MW exhaust power, albeit at a higher upstream density and longer connection length to enforce high collisionality.
Simulations were run to steady-state, and the resulting temperature/heat flux profiles produced using the different heat flux models are shown in figure 1. These simulation parameters produced a SOL collisionality ν * SOL = L/λ e of~170, indicating a SOL that is highly collisional. Under such conditions the results from the Ji-Held model should reduce to be the same as for Spitzer-Härm. This is indeed observed, with the simulation outputs converging on the same result for both models.
Decreasing the input particle source so that the upstream density decreased to~2.1 × 10 19 m −3 , the Ji-Held and Spitzer model outputs start to significantly diverge (figure 2). The conductive heat flux remains approximately the same in both simulations, since the total power being transported through the 1D domain is fixed by the specified input power. However the calculated temperature profiles required to carry this heat flux are dramatically different between the Spitzer-Härm and Ji-Held simulation results, with an elevated T e profile and steeper temperature gradients for the results using the Ji-Held model. Given the only difference between the simulations here is the heat flux model employed, the descrepancies observed can be attributed to nonlocal transport effects predicted by the Ji-Held model. The change in conditions reduced the collisionality of the simulations to ν * SOL ∼ 55-70. For lower collisionality conditions such as this, nonlocal effects would be expected to become more relevant, and so the diverging outputs observed is aligned with our expectations.
One issue that arises for the simulations using the Ji-Held model is the appearance of a peak in the T e profile at s || ∼ 10-20 m. The calculated heat flux from the Ji-Held model remains positive throughout this region (starting at q || = 0 at s || = 0), despite the positive temperature gradient. The cause of this is currently unclear, but could be attributed to difficulties of the model handling a stagnation-to-target SOL domain. The initial implementation of the Ji-Held model in BOUT++ used a target-to-target simulation domain [20], and the conditions on the boundary moments are optimised for this setup. A stagnation-to-target domain would require a different mathematical treatment of the boundary moments, which may be introducing errors that cause this observed T e peak. However, we believe the error introduced here to be a small factor -the temperature profile is almost flat over several collision lengths at the midplane end of the simulation domain in both cases, therefore the influence of the midplane boundary condition should not be too significant.

1D ITER-like tokamak
Simulations were run with relevant parameters for 1D ITER steady-state conditions: upstream density set to 4.0 × 10 19 m −3 using the density feedback controller in SD1D; parallel power flux density at the X-point set to 0.8 GW m −2 ; connection length L = 70 m from stagnation point (s || = 0) to divertor target, with X-point location at s || = 42 m. These parameters are consistent with 2D modelling studies of ITER with a SOL exhaust power of 100 MW and SOL e-folding length λ q = 3.5 mm, for the first grid cell ring outside the separatrix [26,27]. High recycling fraction of 0.95 is applied at the target boundary, with a neutral diffusion factor set to account for a~2 • field line angle with the divertor target.
No impurity seeding is applied to the simulations at this point. Simulations for this ITER study were again run on a 1D grid with 320 evenly-spaced grid-cells, and the Ji-Held model simulations used 20 × 20 moments in the Legendre-Laguerre basis. Spitzer-Härm, flux-limited (FL) and nonlocal SD1D heat flux models were run to stationary steady-state solutions, with the FL runs performed with an α value of 0.2 (typical value based on previous limited kinetic simulation results [28,29]). Particular focus in these studies will be given to comparison of the nonlocal model with the FL α = 0.2 case, as this is the typical model applied for simulations of the ITER SOL in large scale fluid codes. The SD1D-nonlocal code does take significantly longer to run than the standard local SD1D -running on 16 cores the SD1D model took 8-10 hours to run a real-time simulation of 40 ms for the ITER conditions used in this study, whereas SD1D-nonlocal needed 5-7 days to complete the same. This in part motivated the choice to limit the simulations to 20 × 20 moments, as the simulation time scales significantly with the number of moments employed.
Resulting SOL temperature profiles for our ITER 'basecase' steady-state conditions are shown in figure 3 for each heat flux model. The electron temperature (T e ) profile from the Ji-Held nonlocal model is notably hotter across the entire domain than the FL α = 0.2 model output. This difference is at a minimum at the domain boundaries, with a 5% discrepency (~10 eV) for upstream temperatures, and target boundary temperatures within 5 eV for the two models. For the region immediately in front of the divertor target (s || = 60-68 m), the difference increases to >20 eV (>10%). This is an important region in the SOL for neutral-interactions/radiative energy losses necessary for achieving detachment, and therefore for the T e profile to be notably hotter here using the nonlocal model could be a concern.  However towards the target the flux-limitation increases significantly, with the q cond e|| profiles diverging and the equivalent α dropping to <0.1, demonstrating the spatially-dependent nature of the flux-limitation with the nonlocal model.

Varied ITER SOL collisionality
The upstream density conditions were varied to investigate different collisionality regimes with the nonlocal model, by increasing and decreasing the upstream density controller to values of 3 × 10 19 m −3 and 5 × 10 19 m −3 . These values are consistent with the range of upstream densities that are typically investigated for 2D ITER simulations [26]. Increasing n sep to 5 × 10 19 m −3 , a 25% increase on our ITER base scenario, the simulation output becomes less nonlocal (figure 5) as a result, with greater agreement between the nonlocal and FL α = 0.2 models observed upstream, but again diverging towards the divertor target, with significantly higher target plate temperatures that have better agreement with the Spitzer-Härm output. Flux limiter simulations for various α values are performed again to identify if an equivalent global flux limiter can reproduce the nonlocal model results under these conditions. Upstream temperatures were best reproduced using  Reducing the upstream density to n sep = 3 × 10 19 m −3 (a 25% decrease on the ITER base scenario), however, resulted in a drastically different temperature profile using the nonlocal model than for both the Spitzer and FL models. Upstream, the Ji-Held model produced a hotter T e profile (>30 eV) across most of domain (figure 7), before signigicantly dropping to lower T e at the target boundary (~30 eV), with much steeper dT e /ds gradients. This behaviour is in contrast to Ji-Held model outputs for the base ITER conditions and raised upstream density cases. Given ITER will not operate at a single set of plasma parameters, but rather these will vary both in time and across flux surfaces, this result indicates that nonlocality could have importance in regions/regimes for the ITER SOL where the upstream density and/or SOL collisionality are low.

Predicting SOL nonlocality
It would be desirable to be able to determine if a tokamak SOL will exhibit nonlocal effects without having to run a full nonlocal SOL code first, since these simulations are computationally demanding and time-consuming. In this section, we explore the potential of typical metrics for their ability to predict/assess the nonlocality of a tokamak SOL.

Assessing nonlocality using local temperature scalelength
One obvious candidate is the L T /λ e metric (the 'inverse Knudsen number' K −1 n ), which gives the ratio of the temperature gradient scalelength L T and the electron mean free path λ e .  This is the typical metric used in nonlocality studies, for which L T /λ e < 100 is usually taken as indication that nonlocal kinetic effects would become present. In the context of the tokamak SOL, this metric will typically take its minimum value close to the divertor target, where temperature gradients are steepest. For our study, we take λ e ≈ λee √ 2 , assuming λ ee ≈ λ ei for T e ≈ T i (where λ ee and λ ei are the electron-electron and electron-ion collision lengths). Calculating the value of L T /λ e across the domain for the Ji-Held model solutions over the density scan (shown in figure 9(left)), the minimum value of this metric notably decreases as the upstream density decreases, with minimum values of~30,~20 and~10 obtained for the higher density, base-case and lower density ITER scenarios respectively. This trend reflects the level of importance of nonlocality that was observed in the results in sections 4 and 5, with increasing significance of nonlocal effects on the temperature/heat flux profiles as the SOL upstream density (and therefore overall collisionality) was decreased. In all cases, the L T /λ e minimum is low enough to suggest that significant nonlocal effects should be present, as was observed in the simulation results in the discrepencies of the nonlocal model with the local Spitzer-Härm and FL models.
As mentioned, the aim is to be able to predict SOL nonlocality without having to run the full nonlocal code, so applying the L T /λ e metric to the Ji-Held model outputs is not satisfactory to this end. Instead, it would be ideal if applying the metric to the outputs from codes using local thermal transport models could be used for this predictive capability. However, we find that for SOL codes with local Spitzer-Härm or FL thermal transport models and typical sheath boundary conditions, the L T /λ e metric is unable to predict the level of nonlocality of the SOL that would be observed in a nonlocal simulation. Calculating L T /λ e for the Spitzer-Härm (figure 9(middle)) and α = 0.2 FL steady state ITER solutions (figure 9(right)), the minimum value does not show the same trend of decreasing in magnitude with decreasing collisionality, and shows no obvious trend across the density cases. The lowest density case of 3 × 10 19 m −3 appears to have the highest minimum L T /λ e value, and does not drop below 100 significantly at all, in contrast to the trend and high impact of nonlocality observed in the simulations. It is believed that the explanation to this lies in the interaction between the sheath boundary conditions and the thermal conduction models, and its impact on the value of L T /λ e in the simulation domain (an example and discussion of this is given in appendix A). This undermines the use of K −1 n for evaluating the potential nonlocality of a system from standard local heat flux model simulations with typical sheath boundary conditions.

Assessing nonlocality using connection length
An analysis is instead attempted for the SOL collisionality parameter ν * SOL = L/λ e . Figure 10 shows this value calculated across the domain under the ITER base and raised/lowered density conditions. We are again interested in the minimum value of L/λ e for our assessment of nonlocality. The results show a clear trend across all three thermal transport models; the raised density scenario has a minimum L/λ e value of~50, the ITER base scenario has a minimum of~35, and the lowered density case has a minimum of~20. This reflects the trend observed in the simulation results in sections 5 and 6, with the nonlocality becoming more important as the density decreases. Since these results have proven to be relatively consistent across all three models, it allows for local and flux-limited simulations to be useful in assessing the nonlocality of tokamak SOL conditions. Rather than running the full nonlocal model (which takes significantly longer to run than the local codes), running SOL simulations with the Spitzer-Härm or flux-limited models and then calculating the minimum L/λ e value can be used as a rough initial assessment for whether the SOL conditions for a particular tokamak would be relevant for nonlocal transport.
A simple formula can be applied for this purpose also, if theoretical, experimental or simulation estimates of upstream temperature/density conditions are known. Inserting the definition of λ e = β Te 2 ne (where the constant β = 1.5×10 54 √ 2ln(Λ) and taking ln(Λ) = 15) and adjusting the units in the terms, a formula for ν * SOL can be written as where T u and n u are the upstream electron temperature/density respecitvely. Using parameters representative of the FL α = 0.2 model ITER cases studied in Sections 4 and 5, the L/λ e values estimated from this formula are 53.3, 34.9, and 19.2 for the raised density, ITER-base and lowered density scenarios respectively. These are consistent with figure 10 and the relevance of nonlocality to the three cases observed in this study previously. Figure 11. Plot of Lnu against Tu for DIII-D [30] and JET [31] H-mode inter-ELM data, the FL α = 0.2 ITER scenarios investigated in Sections 5 and 6, and predicted SOL conditions for a future DEMO [32]. Contours of ν * SOL = L/λe = 100, 50 and 30 are shown.
Equation (17) is applied to a typical DIII-D case [30] and JET shot data [31] for SOL parameters in H-mode (between ELMs), as well as to approximate SOL conditions predicted for DEMO [32]. Results are shown in table 1 and figure 11. Only one of the JET shots shown here has L/λ e that drops notably below 100 (to 71.8), but even this case does not reach a collisionality that is comparable to any of the ITER cases investigated here (though it is stated in reference [31] that some of the lowest collisionality cases were excluded from their analysis). This indicates why local thermal transport models have reasonable success with the modelling of these experiments. However, experiments to study nonlocal thermal transport on existing devices may be possible for lower collisionality conditions than in the JET data shown here. Some nonlocality has been observed in kinetic modelling using the KIPP code for some JET H-mode discharges [33].
The predicted DEMO conditions result in a low L/λ e of 18.3, a similar value to the ITER low density case we studied. This similarity occurs despite DEMO having a higher upstream temperature, due to a longer connection length of a larger device which compensates. This result suggests that if nonlocal effects will have some importance for ITER conditions, as our results suggest, they will pose at least a similar level of concern for future pilot-plant reactor relevant devices like DEMO, if not greater. Nonlocal transport will almost certainly have a huge impact for modelling of ELMs/disruption events, but the results here indicate that in DEMO and lowercollisionality cases in ITER it is unlikely that nonlocal effects could be ignored even in steady-state.

Impurity seeding
To mitigate the high anticipated heat loads on the divertor target plate, ITER will employ impurity seeding in the SOL to induce radiation energy losses in the plasma, to lower the target plate temperature and enable detached divertor operation. To investigate the impact of nonlocal transport on code predictions with impurity seeding, simulations are repeated with a 'fixed-fraction' carbon impurity model -where impurity concentration is set at a percentage of the plasma electron density throughout the domain -aiming to reach detached conditions. Unfortunately, whilst the SD1D code could successfully model detached conditions (target plate T e < 2 eV), the SD1Dnonlocal code using the Ji-Held model was unable to do so; in all cases under detached conditions the CVODE solver aborted code execution due to failure to reach convergence in simulation timesteps. It is therefore not possible to do a comparison under fully detached conditions of the Ji-Held model with the local transport models, and we are forced instead to investigate differences in model outputs as detached conditions are approached, and for changes in the detachment threshold. The impurity fraction f imp was steadily increased for the ITER base case scenario in section 4 to increase radiative energy losses in the SOL, once again for simulations using all three thermal transport models. Using the α = 0.2 FL model, the onset of detachment occured at an impurity fraction of~28%. This is an excessively high impurity fraction -it has been estimated that a impurity fraction of 10% in the divertor volume will be required in ITER to achieve divertor detachment [34] -but this is due to: 1) using a carbon impurity in SD1D-nonlocal, whereas the 10% estimate is for nitrogen; 2) no use of target gas-puffing in our simulations to reduce target temperatures; and 3) a large hysteresis effect using the impurity model -that once lower target temperatures are obtained (for which the impurities radiate more strongly), detachment is achieved and maintained when the impurity fraction is then decreased to <5%, more in-line with the 10% estimate. Given that SD1Dnonlocal could not reach detached conditions with the Ji-Held model, this hysteresis effect could not be taken advantage of for our study, so comparisons at the higher impurity fractions (f imp > 25%) under attached conditions had to be used instead. For these reasons, in this study we are not concerning ourselves with the absolute value of impurity fraction required to reach detachment, but instead simply using the impurity model as a tool to disproportionately remove energy from the divertor target region as a method to approach detached conditions and to study what happens when nonlocal heat flux and temperature-dependent impurity radiation effects are combined. Figure 12(upper) shows T e , n e and n g profiles for the three thermal conduction models for the ITER-base scenario with an impurity fraction f imp = 27%. This is immediately before the FL model detaches -further increasing f imp to 28% causes the target temperature to decrease below 50 eV in the FL model, where the impurities radiate much more effectively, and from there the target temperature cascade down to~0 eV and detachment is onset. Significant differences between the simulation outputs at f imp = 27% are apparent, in particular that the target plate T e in the FL model has decreased significantly from 130 to 60 eV between the zero impurity and f imp = 27% cases, but for the nonlocal model the target T e is still high at~120, only a 15 eV decrease from conditions without impurities ( figure 3). This is despite the fact that the two models produced target temperatures within 5 eV of each other without any impurity seeding. Increasing the impurity fraction to 30% ( figure 12 (lower)), the FL model becomes detached, with T e dropping to <2 eV and large decrease/increase in n e /n g respectively in the target region, but both the Ji-Held model and Spitzer-Härm remain attached with very high target temperatures over 100 eV. This demonstrates a clear change in the detachment threshold between using nonlocal and FL models for ITER conditions, which presents a potential concern for designs of divertor detachment systems on ITER.

Discussion
The results in this paper highlight that nonlocal thermal transport may well be important for the ITER SOL and for ITER divertor designs, with the nonlocal model implemented showing notable differences in simulation outputs compared with local thermal transport models, and showing strong sensitivity to even small changes in parameters away from the ITER base scenario. The results also highlight the inadequacy of flux-limiters in being able to accurately capture the thermal transport in conditions where collisionality varies significantly across the simulation domain, which our results suggest is relevant to ITER in steady-state. Capturing changing collisionality regimes over time is also beyond the scope of flux-limiters, for events such as ELMs or transients, which are most certainly relevant for ITER.
Particularly concerning features are the hotter temperature profiles observed for the nonlocal model in the region between the X-point and divertor target in some scenarios (figures 3 and 5), in particular for the base ITER scenario where the change in detachment threshold with impurity seeding is observed between transport models. These factors are likely relatedwith the impurities radiating less effectively at higher temperatures, therefore removing less energy from the divertor region at the elevated temperatures for the Ji-Held model. In this case it could be expected that the Ji-Held model for the reduced density scenario (n sep = 3 × 10 −19 m −3 ), which reaches lower divertor temperatures than the other models (figure 7), may instead reach detachment earlier than the FL model. However, other explanations are also possible; reducing temperatures via impurity seeding and therefore making SOL conditions more collisional may result in the Ji-Held model predicting a lower level of flux limitation, with a resulting flatter temperature profile more similar to that predicted by Spitzer-Härm, which also did not detach at the 28% threshold of the FL α = 0.2 model. This is indeed observed in the impurity simulations, with the equivalent α increasing as f imp was increased ( figure 13). Indepth analysis of the Ji-Held model dependencies would be required to identify the predominant effect here. Regardless of the underlying cause, such a discrepency as seen in figure  12(right) could not be tolerated for ITER, and would cause significant damage to the divertor target if realised in practice.
It is also worth noting that the SOL conditions used in this study are representative of 2D ITER simulations with a SOL power width λ q of~3.5 mm. The expected value of λ q for ITER is not currently agreed amongst reseachers, and there is Figure 12. Comparison of Te, ne and ng profiles for SD1D and SD1D-nonlocal for the ITER-base conditions with fixed fraction carbon impurity f imp of (upper) 27% and (lower) 30%. All solutions shown are in steady-state besides the FL model solution for 30% impurity, for which the detachment front is moving upstream towards the X-point and at this impurity fraction eventually MARFEs. evidence from multi-machine scalings that λ q could be much smaller,~1 mm or even less [35,36]. For narrower λ q the parallel power flux density in the SOL would significantly increase, potentially as high as~5 GW m −2 . Nonlocal transport effects would become even more important for the ITER SOL in such a scenario.
Modelling these scenarios in 1D is very simplified compared to 2D models, and so it is important to highlight the need for full 2D modelling with nonlocal transport in order to properly assess the impact of these effects for ITER. The results obtained in this paper therefore motivate the study of a nonlocal transport model in 2D codes in future work: to determine if the discrepancies with Spitzer-Harm/FL models persists, or if the impact of cross field transport between flux tubes of varying collisionality reduces or removes these differences. If such discrepancies as observed in this work are also present in 2D simulations with nonlocal parallel transport, this would have impact on not only the separatrix temperature and target detachment physics, but also on cross field T e gradients and related ∇T e -driven turbulence and transport. In theory, the Ji-Held model could be incorporated into any 2D (e.g. SOLPS, UEDGE, EDGE2D) or 3D (e.g. GBS, HERMES, STORM, TOKAM3X) by replacing the relevant term for the parallel electron thermal conduction in the pressure/energy equation with the Ji-Held model (or other nonlocal thermal transport models). However, the factor >10 increase in simulation time to run the Ji-Held model over the Spitzer/flux-limiters found in this study limit the model's application to 2D/3D codes, which have already long computational time demands, and may cause problems in keeping simulations running in a reasonable timeframe.
These results should be viewed in light of the nonlocal model's limitations as well. The model has previously been compared to the kinetic code KIPP in a simplified problem of large temperature gradients, and was found to underestimate the level of flux limitation in comparison to the kinetic code results [17]. That said, the fact that the Ji-Held model underestimates the flux limitation suggests that it is underestimating the importance of nonlocality in our work. This would potentially mean the discrepancies observed in the T e profiles and detachment threshold may be further increased and be an even greater concern for ITER modelling than suggested here. This motivates further study of nonlocal models for the ITER SOL in low collisionality regimes and in the context of ELMs/disruptions. The nonlocal model also fails to reproduce any kinetic preheat effects, so only in part addresses issues for nonlocality in the SOL. The exact nature of the discrepancies observed in the simulation outputs from nonlocality may vary depending on the nonlocal thermal transport model employed, so other models that more accurately capture nonlocal flux-limitation and preheat effects should also be investigated for further SOL simulation research.
In addition, appropriate boundary conditions for nonlocal thermal transport models remains a concern that needs addressing. The nonlocal model was employed with a standard sheath boundary condition given by equation (13). This itself is a dubious assumption, since this condition employs local thermal transport assumptions. In collisionless plasma conditions, streaming thermal-energy-carrying electrons from regions further away in the plasma would impact the boundary, making this sheath condition invalid. More attention needs to be given to deriving appropriate boundary conditions that can be used with this model based on nonlocal transport theory.

Conclusions
While much attention has been given to perpendicular transport in the SOL, relatively little has been done to improve the sophistication of parallel transport models. As demonstrated here, reduced-kinetic models of nonlocal parallel transport provide a potential solution to addressing this problem. Despite the limitations outlined in section 9, the results in this paper provide insight into how the incorporation of nonlocal thermal transport models could impact ITER SOL simulation outputs, in particular demonstrating the potential for the inclusion of nonlocal transport effects to change predictions for the ITER SOL and divertor target for low density scenarios, which are important to model accurately for divertor heat flux considerations. Simulations performed for the ITER steady-state conditions with varying upstream density suggest a significant relevance of nonlocality in the ITER SOL, with strong sensitivity observed for small changes in the upstream density conditions. Results show notable discrepencies in the temperature profiles predicted by the Ji-Held model against both the Spitzer-Härm and FL models, typically showing elevated T e profiles, as well as in the detachment thresholds between models when adding impurity seeding. Global flux limiters are shown to be inadequate to capture the changing SOL conditions/physics across the domain. An analysis of SOL collisionality and nonlocality suggests that nonlocal effects will also be very significant for future devices like DEMO, that can not be ignored even in steady-state. Further investigations into the applications of nonlocal models for time-dependent SOL simulations in ITER should be performed, in particular in application to full 2D modelling of the ITER SOL. Other nonlocal transport models that more accurately capture flux limitation as well as preheat effects should be investigated for further SOL simulation research. Work should also be carried out to determine appropriate divertor boundary conditions for nonlocal transport conditions. and λ e = β T e 2 n e (A2) where β = 1.5×10 54 √ 2ln(Λ) is a constant, and ln(Λ) is the Coulomb logarithm. In the local Spitzer-Härm model, the electron parallel heat flux q ||e is defined by the equation: where κ oe is the plasma electron thermal conductivity. K −1 n will be at a minimum where the temperature gradient is steepest, which occurs towards the divertor target. Rearranging equation (A3), inserting the value for dTe dx into equation (A1), and substituting values of T e and n e for their values at the target plate gives a new function for K −1 n : The heat flux at the target is subject to the boundary condition: where γ is the sheath heat transmission coefficient and the plasma sound speed c s is given by where µ is a constant. Substituting these equations into equation (A4), a value for the minimum is found to be: All terms in this resulting equation are constants. Therefore, for simulations run to steady-state with conduction only using the local Spitzer-Härm heat flux model and sheath boundary condition assumptions, the minimum value of K −1 n will simply be a constant, determined by the γ factor set for the target boundary. A constant value of the L T /λ e minimum would then not provide an indication of the degree of nonlocality of the SOL. In reality, the convective component of the heat flux and radiative losses would alter this result to not be an exact constant. But this analysis demostrates as an example how the sheath boundary can impact on the simulation results and the measured L T /λ e metric, irrespective of the collisionality of the system.