Validation of edge turbulence codes against the TCV-X21 diverted L-mode reference case

Self-consistent full-size turbulent-transport simulations of the divertor and SOL of existing tokamaks have recently become feasible. This enables the direct comparison of turbulence simulations against experimental measurements. In this work, we perform a series of diverted Ohmic L-mode discharges on the TCV tokamak, building a first-of-a-kind dataset for the validation of edge turbulence models. This dataset, referred to as TCV-X21, contains measurements from 5 diagnostic systems -- giving a total of 45 1- and 2-D comparison observables in two toroidal magnetic field directions. The dataset is used to validate three flux-driven 3D fluid-turbulence models: GBS, GRILLIX and TOKAM3X. With each model, we perform simulations of the TCV-X21 scenario, tuning the particle and power source rates to achieve a reasonable match of the upstream separatrix value of density and electron temperature. We find that the simulations match the experimental profiles for most observables at the OMP -- both in terms of profile shape and absolute magnitude -- while a poorer agreement is found towards the divertor targets. The match between simulation and experiment is seen to be sensitive to the value of the resistivity, the heat conductivities, the power injection rate and the choice of sheath boundary conditions. Additionally, despite targeting a sheath-limited regime, the discrepancy between simulations and experiment also suggests that the neutral dynamics should be included. The results of this validation show that turbulence models are able to perform simulations of existing devices and achieve reasonable agreement with experimental measurements. Where disagreement is found, the validation helps to identify how the models can be improved. By publicly releasing the experimental dataset, this work should help to guide and accelerate the development of predictive turbulence simulations of the edge and SOL.

A summary of nomenclature is provided in Appendix B.

Introduction
Magnetic-confinement-fusion devices cannot provide perfect confinement of the plasma in the core. Turbulent fluctuations and collisions lead to heat and particles being transported across field-lines, eventually reaching the solid walls of the device. The peak heat flux reaching plasma-facing components must be kept below engineering limits to prevent damage to the vessel. One key element to address the power exhaust problem is the diverted magnetic field geometry. In this configuration, an X-point is introduced in the tokamak boundary via shaping coils, diverting the boundary plasma to heat-resistant targets. Divertor geometries provide several benefits in comparison to simpler limited geometries. By increasing the separation of the confined region and the plasma-facing components, the divertor geometry improves the screening of impurities generated through plasma-wall interactions or intentionally injected to radiatively cool the plasma [1]. They also increase the volume for scrape-off layer (SOL) plasma cooling and the connection length over which cross-field transport can broaden the heat flux channel [1], as well as helping to reach improved-confinement and detached regimes [2]- [4].
Divertor plasmas are, however, challenging to model and predict, due to the interplay of turbulence, drifts, plasma gradients, coherent filaments and interactions with neutrals and the walls [1], [5]. Most divertor modelling is performed with transport codes such as SOLPS-ITER [6] or SOLEDGE2D [7], which treat the cross-field transport as an effective heat and particle diffusion, rather than directly modelling the small-scale convection due to turbulence. This reduces the computational cost, which in turn allows transport models to be run for large machines and over long time-scales. However, one limitation of transport modelling is that the diffusive transport coefficients are not self-consistently determined, and instead must be either heuristically fitted to experimental data, computed via reduced models [8] or determined via a coupled turbulence code [9], [10]. This can provide a reasonable match to the mean plasma profiles of existing devices, but cannot describe the timedependent behaviour of the plasma. Furthermore, particularly in the SOL, the plasma can form coherent structures called filaments, or 'blobs', which transport heat and particles ballistically rather than due to local gradients [11], breaking the assumption of diffusive transport.
Turbulent self-organisation also leads to non-linear behaviour, which complicates direct extrapolations from current to future devices. Therefore, to model the time-dependent dynamics and make predictive simulations of the divertor and SOL, it is necessary to simulate the turbulent nature of the transport * .
Direct numerical turbulence simulations require more sophisticated physical models than transport codes, and orders-of-magnitude more computational resources due to the 3D multi-scale nature of turbulence.
Full-size turbulence simulations of existing machines allow us to validate our turbulence codes, which is an important step towards the development of predictive simulations for future devices such as ITER. Validation (in combination with verification) is a common technique in software testing. In the fusion community, a set of best practices for model validation was proposed by Terry et al, 2008 [22] and Greenwald, 2010 [23] -outlining a rigorous validation methodology based on the quantitative comparison of multiple measurements at different 'primacy hierarchies'. Importantly, validation here is not a binary result, but rather a tool for checking the fidelity of the simulations and guiding targeted development of the models -with repeated validation suggested as part of a model testing and development cycle [23]. This methodology has already been used to test boundary turbulence simulations in basic plasma physics devices [24]- [26] and limited tokamak plasmas [17].
In this paper, we extend these previous works to the validation of edge turbulence codes in diverted tokamak geometry, with the goal of guiding the development of the models and assessing how close simulations are to reproducing realistic plasma behaviour. For this purpose, a diverted, Ohmic Lmode scenario has been developed on the Tokamak a Configuration Variable (TCV) [27], performed in both toroidal magnetic field directions. Thanks to the large suite of edge/SOL diagnostics available on TCV, an extensive experimental dataset has been collected, allowing for a stringent assessment of the simulationexperiment agreement. We refer to this scenario and dataset as the TCV-X21 validation reference case. This is used to test three 3D boundary turbulence codes -namely GBS, developed at the Swiss Plasma Center at EPFL, Lausanne [28], [29], GRILLIX [12], [21], developed at the Max Planck Institute for Plasma Physics, Garching and TOKAM3X [30], [31], developed at the CEA, Cadarache in collaboration with Aix-Marseille University. The codes solve subsets of the drift-reduced Braginskii fluid equations [32], which require sufficient plasma collisionality such that each plasma species is close to a local thermodynamic equilibrium [32], [33]. As such, the codes are not suitable for modelling the reactor core, and we focus our validation on the edge and open field-line region only.
By validating several codes against a common reference case, we can investigate how differences between the codes affect the results of divertor modelling, assess the importance of physical processes, and -using the results of the validation -guide the development of the codes. Validation against the TCV-X21 case could also benefit other boundary turbulence codes -such as XGC [18], COGENT [34], GENE-X [35], Gkeyll [36], ORB5/PICLS [37], GYSELA [38], FEL-TOR [39], BOUT++ [40], STORM [16], Hermes [41], GDB [42], HESEL [43] and SOLEDGE3X [44] -and could eventually be used as a common divertor reference case, similarly to the CYCLONE base case used for core modelling [45]. To enable future testing against the TCV-X21 case, we provide the experimental and simulation results in a Findable, Accessible, Interoperable, Reuseable (FAIR) data repository, along with additional documentation and data (such as the magnetic equilibrium) to help set up and post-process future validations. The data is available both as NetCDF files, and (for the experimental data only) in ITER Integrated Modelling and Analysis Suite (IMAS) format [46]. A dynamic repository is provided at gitlab.mpcdf.mpg.de/tcv-x21/ tcv-x21, which will be updated with the results of future comparisons against the reference case (this is encouraged, to provide an evolving picture of the stateof-the-art of divertor modelling). We additionally provide a static repository for the version used in this paper at zenodo.org, and a web-interface to the processing routines at mybinder.org. Throughout this paper, wherever extended analyses are available through the repository, we indicate this via a file-path relative to the root of the TCV-X21 repository.
The remainder of this paper is organised as follows: we first outline our validation methodology in Sec.2, following the methodology in Ricci et al., 2015 [25]. We then discuss the development of the experimental reference scenario and the collection of the experiment dataset in Sec.3. Next, we briefly introduce the three participating codes and discuss how the experimental reference scenario was simulated in Sec. 4. The results from the experiment and the simulations are compared both graphically and via a validation metric in Sec. 5. We then discuss the results of the overall validation and the physics observed in Sec.6.

Validation methodology
We start our validation by outlining the methodology. In this study, we perform both qualitative and quantitative validations in Sec.5, and consider the overall result in Sec. 6. For our qualitative validation, we simply mean graphically comparing the simulated results and the experiment. This is helpful for evaluating the ability of the codes to make predictions of the dominant physical processes, of the shape and magnitude of the profiles, and for building an understanding of why the simulations agree or disagree. Qualitative validation can, however, be imprecise or subjective -which is why we also perform a quantitative validation.
The goal of quantitative validation is to provide a single numerical value of the level of agreement between simulation and experiment. This is performed using a validation metric, which is a type of summary statistic like the average-absolute-difference or Pearson correlation coefficient. Since a validation is more meaningful if more observables are compared, a complementary quality metric is typically also given, which provides a measure for the number and precision of the observables used. For this study, we use the validation methodology presented in Ricci et al., 2015 [25], which is based on Terry et al., 2008 [22]. We briefly review here the concepts and terms.
To quantify the level of agreement between simulation and experiment in a validation involving several observables, a 'composite metric' is useful. In Ricci et al., 2015, a composite metric χ is computed from the individual levels-of-agreement for each observable j, denoted here R(d j ), which are combined via a weighted average. Each observable is weighted according to its 'primacy hierarchy' H j and its 'sensitivity' S j .
The individual level of agreement R(d j ) is computed from the root-mean-square of the errornormalised experiment-simulation difference (roughly equivalent to the RMS Z-score) for some observable denoted 'j' where the experimental measurement has estimated values e j,i and uncertainties ∆e j,i defined at some set of discrete data points i = {1, 2, ..., N j }. The simulation result is assumed to be continuous and so is interpolated to the experimental measurement positions, giving computed values s j,i and uncertainties ∆s j,i . The uncertainties of the experimental results is evaluated for each diagnostic in Sec.3. A rigorous estimate of the simulation uncertainty is difficult to determine, so we simply set the simulation uncertainty to zero. We note here that this has the effect of increasing our d j values, which we discuss in Sec.6. For a discussion of sources of simulation uncertainty, see Ricci et al., 2015 [25] and Chapter 5 of Computer Simulation Validation [47]. The level of agreement R(d j ) is computed from d j via a smooth-step function for d 0 the 'agreement threshold' and λ the 'transition sharpness', which are usually set to d 0 = 1 and λ = 0.5 [25]. The shape of this function and the effect of the d 0 and λ parameters are shown in Fig.1.
The level of agreement is combined with the primacy hierarchy and the sensitivity of the observable. The H j hierarchy weighting is computed as where h Exp and h Sim are the primacy hierarchies of the observable for the experiment and simulation, and h Comp is the combined hierarchy of the comparison. Higher values of the primacy hierarchy indicate observables which require stronger assumptions, or calculation from a model combining multiple measurements (for h Exp ) or combinations of multiple directly simulated quantities (for h Sim ). An extended discussion of the primacy hierarchy can be found in Ref. [48]. By using the inverse of h Comp in Eq.3, we weight directly available observables more than indirect observables. The observables used in this validation together with their primacy hierarchies are given in Tab.1.
The S j sensitivity weighting is computed as using the same notation as in Eq.1. The sensitivity is a measure of the relative total uncertainty of the observable. It approaches 1 for observables with very high precision, and 0 for observables which have very high uncertainties.
Finally, we compute the composite metric via the weighted average which gives values between 0 (quantitative agreement) and 1 (quantitative disagreement). This is combined with the overall 'quality' which gives higher values for validations considering more directly-computed, high-precision measurements.

Experimental Scenario
In this work, we developed an experimental scenario in TCV for the validation of boundary turbulence codes. TCV is a medium size tokamak (R axis = 0.88m) with nominal vacuum toroidal magnetic field of B φ,axis 1.45T equipped with 16 independentlypowered shaping coils providing extreme plasma shaping capabilities [27]. The experimental scenario, referred to as the TCV-X21 reference scenario, is a Lower-Single-Null L-mode Ohmic plasma. A poloidal cross-section showing the magnetic flux surfaces of this scenario, obtained from the LIUQE magnetic reconstruction code [49], can be seen in Fig.2 together with the set of diagnostics used to collect the experimental dataset.
The discharges were performed in Deuterium and at a reduced toroidal field of B φ,axis 0.95T.
This has the benefit of increasing the characteristic perpendicular scale length of the turbulence, which is given in terms of the sound drift scale where m i is the mass of the Deuterium ions. To locally resolve the turbulence drive due to ballooning, driftwave and ITG instabilities, a numerical grid resolution in the direction perpendicular to the magnetic field on the order of the local drift scale is required (the exact resolution requirement depends on the instability and the numerical scheme). By reducing the toroidal magnetic field by a factor of 1.5, we can resolve the drift scale (or some multiple of it) with a 1.5× lower poloidal and radial resolution, reducing the number of simulation grid points and therefore the computational cost of the simulations. To avoid MHD modes and Ohmic H-mode transitions in the forward-field case, we also reduced the plasma current to I p 165kA, giving an edge safety factor of q 95 ≈ 3.2.
Since neutrals were not included in the simulations, the effect of their dynamics in the divertor volume was minimised in the TCV-X21 scenario by using a low electron line-average density (determined from the FIR chord shown in Fig.2) of n e ∼ 2.5 × 10 19 m −3 , corresponding to a Greenwald fraction of ≈ 0.25. The discharges were fuelled with D 2 from the top valve, indicated in Fig.2. Fig.3 shows no significant difference between the SOL profiles of electron temperature and density in the divertor entrance and the LFS target profiles, suggesting a sheath-limited regime and, therefore, negligible temperature parallel gradients due to recycling in the divertor region [50]. The substantial difference of the electron temperature near the separatrix in Fig.3 is because the divertor entrance is connected to the hot confined plasma, while the LFS target is disconnected from the confined plasma.
To investigate the effect of drifts, experiments were performed in both the "forward" and "reversed" toroidal field directions. In the forward ('favourable') field case, the ion ∇B drift points downwards from the plasma core towards the X-point, whereas in the reversed ('unfavourable') field case it points upwards, away from the X-point.

Diagnostics & Observables
In the following section, we introduce the diagnostics and the basic analyses used to compute the experimental profiles and uncertainties. Tab.1 shows all observables, divided between the different regions of the SOL and the respective diagnostic used to determine them. For all profiles, we use as the radial coordinate R u − R u sep , which is the distance between the measurement location and the separatrix, mapped along flux Poloidal cross section showing the magnetic surfaces (dark blue lines) reconstructed by LIUQE and the diagnostics used to gather the dataset. We indicate the location of the wall-embedded Langmuir probes at the high-field-side and low-fieldside divertor targets (HFS-LP and LFS-LP -red and blue circles), the Reciprocating Divertor Probe Array (RDPA -black L-shaped structure) and its swept area (black watermark box), the Thomson Scattering system (TS -orange squares), the Fast Horizontal Reciprocating Probe (FHRP -solid and watermark boxes in purple), the far infrared interferometer (used to estimate the line-average density, FIR -cyan), the field-of-view of the Infrared system (IR -green), and the position of the top valve used to fuel the plasma (magenta). at the low-field-side and high-field-side targets Figure 3: Comparison of upstream and target profiles in reversed-field. No significant reduction of T e in the LFS target (right) is observed. n e shows the expected factor 2 drop, characteristic of the sheathlimited regime.
surfaces to the outboard midplane. The flux surfaces in the private-flux region (PFR) are not connected to the outboard midplane (OMP), and so the upstream mapping is carried out using the corresponding surface with the same poloidal flux ψ in the confined region * . * A more involved method for computing the flux-surface label in the PFR is presented in Ref. [51]. Our method is simpler The separatrix, flux surfaces and poloidal flux are obtained from the LIUQE magnetic reconstruction [49]. This approach of using R u − R u sep removes small differences in the plasma positioning and magnetic geometry when combining data from repeated discharges, between different time-stamps, or in the comparison with the simulation profiles.
The magnetic reconstruction reveals a rapid oscillation of the strike-point, with a period of ∼ 35 ms. This movement has a peak-to-peak value of ∆(R u − R u sep ) <∼ 2.5mm in forward-field and ∆(R u − R u sep ) <∼ 1.5mm in reversed-field, affecting spatially fixed measurements which average over a long time interval. This is the case for the parallel heat flux estimated by the infrared vertical camera (Sec.3.3) and density, electron temperature and plasma potential estimated by the wall-embedded Langmuir Probes in swept-bias mode (Sec.3.2). In these cases, the rapid oscillation adds a broadening of the order of the peakto-peak movement to the profiles.
For the experimental measurements, small variations of the separatrix position contribute, to some extent, to the reproducibility uncertainty of the experimental profiles (introduced in the next paragraph). For the to compute, while the method in Ref. [51] may be preferable for detailed analysis of the PFR. simulations, quantifying the effect of the uncertainty in the magnetic geometry would require a sensitivity scan. This analysis is not performed in this work due to computational cost, and as such for the simulations we neglect the uncertainty in the magnetic reconstruction.
In general, we categorise the experimental uncertainty into three sources. The first is ∆e f it , the uncertainty related to fitting experimental data to a model. The second is ∆e dia , due to inherent characteristics of the diagnostics, e.g. uncertainties in the effective ion collection surface of Langmuir probes. Finally, the third is ∆e rep , the uncertainty related to the reproducibility of the observable assessed by comparing repeated discharges, typically performed in separate experimental sessions. Then, the total experimental uncertainty is evaluated as ∆e tot = ∆e 2 f it + ∆e 2 dia + ∆e 2 rep . Depending on the diagnostic and operation mode, not all the sources of uncertainty defined above are present.

Wall-embedded Langmuir probes
In TCV, both targets are covered by wall-embedded, dome-shaped Langmuir probes (LP) as shown in Fig.2. The probes are operated in four different modes: swept-bias, ion saturation current, floating potential, and ground current mode. The details of the basic probe analysis can be consulted in Refs. [54], [55]. Mean profiles of the electron density n e , the electron temperature T e , the floating potential V f l , the plasma potential (obtained as V pl = V f l + 3T e ), the ion saturation current density J sat parallel to the magnetic field, and the parallel current density J || are obtained from the swept-bias mode. Since we additionally have separate discharges where we use the Langmuir probes to measure V f l , J sat , and J || as a function of time, the time-averaged profiles for these quantities obtained from swept-bias mode are only used to determine the uncertainty ∆e rep associated with the reproducibility. The quantities obtained from the nonswept-bias modes are evaluated over time windows of 1ms. In ion saturation current mode, a constant negative bias of −100V is applied to the probes, resulting in a direct measurement of J sat , which is used to estimate the mean, J sat , the fluctuations (standard deviation), σ(J sat ), the skewness, skew(J sat ), and the Pearson kurtosis, kurt(J sat ). Direct measurements of the V f l and J || time histories are performed using, respectively, the floating potential mode (measuring the potential of the probe when floating with respect to the plasma) and the ground current measurement mode (current measured when the probe is biased to the vessel potential), allowing the estimate of mean (V f l and J || ) and fluctuation (σ(V f l ) and σ(J || )) profiles. For an improved spatial resolution, both strike-point positions were swept during the discharges.
For most of the quantities, ∆e rep is estimated using different shots, the only exceptions being σ(V f l ) and σ(J || ), where the single discharge available is divided into 100ms intervals and ∆e rep is estimated using the profiles resulting from the comparison of these sub-intervals. ∆e dia is estimated by assuming an uncertainty of ±0.1mm in the height of the probes (the wall LPs have 4mm of diameter and the domed-shape head protrudes from the tile shadow by 1mm [54]). This source of error affects the quantities n e , J sat , σ(J sat ), J || , and σ(J || ). The last source of uncertainty is ∆e f it , which affects only the swept-bias mode and is estimated as the 95% confidence interval of the IV four-parameter fit.

Infrared Cameras
The vertical Infrared thermography system (IR) covering the TCV outer target (see Fig.2) uses a camera operating with a frame rate of 160Hz and a spatial resolution is 2.5mm [51]. We estimate the heat flux at the LFS target for every camera frame and then average the results to obtain the averaged parallel heat flux q || as a function of R u − R u sep . We also use the standard parametrisation of the heat flux profiles [56] to determine the SOL power fall-off length λ q and spreading factor S.
The only source of uncertainty for q || is ∆e rep which is estimated comparing profiles from different time frames. This accounts for uncertainties related to the strike point oscillation mentioned in Sec.3.1 and the spacial calibration of the IR.

Reciprocating Divertor Probe Array
The Reciprocating Divertor Probe Array (RDPA) installed at the bottom of TCV (see Fig.2) provides 2D measurements of a variety of quantities by combining a fast, vertical linear motion and a radial array of 12 rooftop Mach probes [57]. A typical RPDA plunge takes approximately 350ms and its Mach probes are operated in three different modes: sweep-bias, J sat , and V f l . The quantities obtained in this way are timeaverages of n e , T e , V f l , V pl , M || and time histories of J sat and V f l . The time histories are used to estimate mean and fluctuation profiles of J sat and V f l , as well as the skewness and kurtosis of J sat .
The vertical positions Z of the RDPA dataset are translated to the coordinate system Z − Z X , where Z X is the vertical X-point position determined by LIUQE. Analogously to the radial profiles, this approach removes small differences of the plasma positioning when combining data from repeated discharges or when comparing with the simulation data.
The two main sources of uncertainty in sweptbias and J sat mode are ∆e dia , which is estimated considering an uncertainty of 10% in the probe area, and ∆e rep , which is estimated comparing different shots. For V f l mode, the only source of uncertainty is ∆e rep , since the value of V f l does not depend on the probe area. This is also the case for T e and V pl estimated from swept-bias mode.

Thomson Scattering
The Thomson scattering system (TS) installed on TCV consists of 109 observation positions (chords) covering the region between Z = −69cm to Z = +55cm at a radial location of R = 0.9m (see Fig.2). This TS system can provide spatial profiles of the n e and T e covering a range of T e from 6eV to 20keV [58] (and even down to 1.4eV in the divertor [59]). In our analysis, we use TS data measured near the separatrix in the divertor entrance to produce divertor-entrance profiles, seen in Fig.3. The uncertainty sources are ∆e f it , estimated from the analysis procedure used to obtain n e and T e , and ∆e rep , which is estimated from the comparison of the profiles obtained in different shots.

Fast Horizontal Reciprocating Probe
The horizontal reciprocating probe mounted at the outboard midplane (FHRP, at Z = 0), shown in Fig.2, consists of a probe head with ten electrodes, which are used in different configurations and operation modes to provide measurements of time-averaged and fluctuation quantities [60]. The double probe configuration, operated in swept-bias mode, is used to determine T e and n e . Similarly to the wall LPs, direct, time-resolved measurements of J sat and V f l , performed at 2.5 − 5MHz, are used to estimate the mean and the fluctuations for both quantities, and the skewness and kurtosis for J sat . Two electrodes are used to determine the parallel Mach number M || . The sign convention of M for the FHRP is such that positive values refer to flows in the counterclockwise direction if the torus is seen from top. For the magnetic helicity used in the present experiments (standard helicity), this corresponds to a parallel flow directed towards the LFS target. For T e and n e , all three sources of uncertainty are present and estimated following the same procedure employed for the wall LPs. For the quantities related to J sat , the main sources of uncertainty are ∆e dia , estimated considering an uncertainty of 10% in the probe area, and ∆e rep , which is estimated by comparing different discharges. For the fluctuations and mean value of V f l , only ∆e rep is relevant since this quantity is independent of the probe area.
The ion collection area of the FHRP electrodes is calculated here using the total probe surface area (and accounting for sheath expansion as detailed in [60]) rather than its projection along the magnetic field. This weak magnetic field assumption is made because the FHRP electrode dimensions (cylinders of length 1.5mm and radius 0.75mm) are comparable to the ion Larmor radius at the outboard midplane (ρ s = c s /Ω i ≈ 1.2mm for T e = T i = 20eV and B φ = 0.76T at R = 1.1m). This is not the case for the wall LP and RDPA electrodes, which are larger than the local ion Larmor radius and therefore are treated as strongly magnetised.

Simulation codes
The TCV-X21 scenario was simulated with the GBS, GRILLIX and TOKAM3X 3D two-fluid drift-reduced Braginskii turbulence codes.
In this section we first provide a brief introduction of the models in each of the codes, and then discuss the choice of sources and parameters to model the TCV-X21 scenario. For the complete physical models (excluding boundary conditions) of the respective codes, the reader is directed to Giacomin and Ricci, 2020 [61] for GBS, appendix A of Zholobenko et al, 2021 [21] for GRILLIX, and Tatali et al, 2021 [13] for TOKAM3X. For a discussion of the sheath boundary conditions used in this study, see Sec.4.1 and appendix A. The codes have all previously been verified via the Method of Manufactured Solutions [12], [30], [62], to ensure that the model equations have been correctly numerically implemented. TOKAM3X has additionally been verified via the a posteriori iPOPE method [63].
The codes all solve versions of the drift-reduced Braginskii fluid equations [32] -giving the evolution of a plasma density n (under the assumption of quasineutrality), separate electron and ion temperatures T e and T i , the parallel ion velocity u , the parallel electron velocity v or parallel current density J = en(u − v ) and the electrostatic potential V pl . Additionally, GRILLIX evolves the parallel component of the electromagnetic vector potential A . In this work, the codes neglect the neutral dynamics. To approximate the particle source due to neutral ionisation, a simple confined-region source (discussed in Sec.4.3) is used.
Rigorously, the fluid theory does not allow for modelling low-collisionality plasma regions [32], [33]. The collisionality in the plasma core is too low for fluid models to be formally valid, and as such we do not expect to have a good agreement with the experiment in this region. Nevertheless, both GBS and GRILLIX include the plasma core in their simulations (see Fig.4). This circumvents the need to apply boundary conditions at the core, which do not have a clear physical analogue [61]. Additionally, the ion drift approximation breaks down at the entrance to the magnetic presheath [64] and as such the codes aim to mimic the effects of the plasma sheath via 'sheath boundary conditions', rather than directly modelling the sheath.

Contrasting the codes
There are a number of significant differences between the codes, despite the fact that they all are based on drift-reduced Braginskii models. We consider a few of the most significant differences here, to help interpret the differences found between the simulated results. Considering first the models, the codes apply different assumptions to simplify the numerical implementation of the model. In this work, both GBS and TOKAM3X use the Boussinesq approximation (although with different flavors, for details see Ref. [65] for GBS and Ref. [13] for TOKAM3X) and treat the electrostatic limit of the dynamics, while GRILLIX does not apply the Boussinesq approximation and includes the effects of electromagnetic induction. Additionally, GRILLIX and TOKAM3X include terms for electron-ion heat exchange, in contrast to GBS. Since the start of this project, new versions of the codes with extended models have been developed -these are discussed in Sec.6.5.
The codes employ different sheath boundary conditions at the magnetic presheath entrance near the divertor targets. The details of the boundary conditions used for each code are given in Appendix A, with the key differences summarised here. For GBS, the parallel ion velocity u is set to c s . In GRILLIX and TOKAM3X, corrections for the E × B drift (both codes) and curvature drift (TOKAM3X only) are included in the u boundary condition, and the drift-corrected ion velocity is set to ≥ c s to allow supersonic transients to be freely advected across the boundaries (see Appendix.A and Ref. [7], [66]). For the electron and ion temperatures, GBS enforces ∇ T = 0, while in GRILLIX and TOKAM3X sheath heat transmission factors are used. GBS and TOKAM3X determine coupled expressions for the current and plasma potential in the electrostatic limit, whereas GRILLIX assumes that internally-generated currents freely flow across the boundaries (via a ∇ J = 0 or 'free-flowing' boundary condition) and sets the plasma potential such that V f l → 0. The effect of the boundary conditions is discussed in Sec.6.3.
The codes use different methods to discretise their model equations -from a fourth-order non-fieldaligned scheme in GBS [67], to a domain-decomposed flux-aligned scheme in TOKAM3X [13], to a locallyfield-aligned scheme in GRILLIX [12]. In GBS and TOKAM3X, the boundary conditions are enforced at boundary grid points, while in GRILLIX an immersed boundary method is used [12]. Since all the codes have been verified, the choice of discretisation will have no impact on the solution which the codes will converge to at arbitrarily high grid resolution. For a given grid resolution however, the discretisation error can vary depending on the choice of discretisation. Furthermore, the choice of discretisation affects the geometrical flexibility, computational efficiency and scalability of the codes.

Physical parameters
The physical and numerical parameters of the simulations can be varied to permit coarser spatial resolutions and a larger time-step, which reduces the computational cost of the simulations. The simulations set their resistivity in terms of the Braginskii value [68] where m e , e and τ ee are the electron mass, elementary charge and electron collision time, and we have taken the Coulomb logarithm to be equal to a constant value of 13 (the weak parametric dependence of the Coulomb logarithm is dropped). GRILLIX used the value of η as defined in Eq.8, while GBS and TOKAM3X increased η by factors of 3 and 1.8 respectively, to permit the use of a larger time-step (see Ref. [69]) and avoid numerical instabilities. The codes also set their electron and ion heat conductivities in terms of the Braginskii values [68]. The electron heat conductivity is 16 nT e τ ee m e = 23.9 T e 40eV 5/2 MW/(eV m) (9) and the ion heat conductivity is which we have evaluated for Deuterium ions. GRILLIX used the heat conductivities as defined in Eqs.9-10, with a periodicity limiter (Eq.B.63 from the SOLPS-ITER manual [6]) to limit the core heat flux. Conversely, both GBS and TOKAM3X used reduced heat conductivities, to avoid time-step limitations. TOKAM3X reduced the heat conductivities given by Eqs.9-10 by a factor of 1.8. GBS used an effective heat conductivity of χ ,e = 1.29(n/n ref )MW/(eV m) for the electrons and χ ,i = 0.037(n/n ref )MW/(eV m) for the ions, with n ref = 6×10 18 m −3 (corresponding to Eqs.9-10 evaluated at n = n ref and T e = 41.3eV and then reduced by a factor of 20). Using the experimental values, we see that this gives a heat conductivity reduced by a factor of 20 at the OMP, and by a factor of 4.8 near the targets * . The effect of using relaxed parameters is discussed in Sec.6.2.

Equilibrium, resolution and sources
We select a single 'reference equilibrium' -TCV shot 65402 at time t = 1.0s -which is representative of the experimental discharges. The magnetic field structure of the reference equilibrium is computed by LIUQE and is provided as an input to the codes for the simulations in both toroidal field directions. By using the magnetic field from LIUQE and the electron temperature from Thomson scattering, we can approximately determine the drift scale (Eq.7) as a function of position. We find that the confined region should have ρ s ≥ 1mm, while the open field-line region has ρ s as small as 0.3mm. Therefore, a perpendicular resolution of the order of a few mm should resolve most of the 'primary' turbulence drive in the confined region, which ballistically drives SOL turbulence, and partially resolve 'secondary' instabilities, which locally drive turbulence in the open field-line region. In this work, GBS used a perpendicular resolution of ≈ 2mm, GRILLIX used a perpendicular resolution of ≈ 1mm and TOKAM3X used a resolution approximately equivalent to 1mm radially and 4mm poloidally at the outboard midplane. In the toroidal direction, GBS used 128 planes for 2π of toroidal angle, GRILLIX used 16 planes for 2π of toroidal angle and TOKAM3X used 32 planes for π of toroidal angle (half-torus). Due to the different discretisation methods the resolution requirements may vary dramatically between the codes, although this is difficult to quantify without resolution scans. Due to computational cost, different resolution were tested with GRILLIX only, with the results discussed in Sec.6.1.
The simulations are flux-driven, with freely evolving profiles determined by the balance of * Here, we see that including the strong temperature dependence reduces the heat conductivity in the SOL (which makes the simulations less expensive). In this work, it was not included in GBS due to the divergence of the Braginskii heat flux in the core, while a new version of GBS includes the temperature dependence and a heat flux limiter. source terms, transport mechanisms and sinks at the device walls. The sources for this study are selected to provide simple approximations of Ohmic heating and neutral ionisation. The temperature source is selected to be close to the magnetic axis, since this is the position where Ohmic heating is primarily expected [70]. The density source is placed just inside the confined region, as shown by the green shaded region in Fig.4A. Ionisation in the divertor is not taken into account. For GBS and GRILLIX, the source positions are indicated in Fig.4 and given in TCVX21/grillix_post/components/ sources_m.py. For TOKAM3X, pressure sources (i.e. combined density and temperature sources) are located in the vicinity of the core limiting flux surface, which corresponds approximately to the same position as the GBS and GRILLIX density source. Treating the neutral dynamics via a simple confined-region source is a strong approximation in this work -this is discussed in Sec.6.5.
In addition to the confined-region sources, small additional sources are added in the open field-line region to prevent numerical instabilities which occur at very low temperatures or densities. For GBS, additional particle sources are added in the private-flux region and at the inner boundary where flux surfaces become tangent to the wall. These sources are intended to prevent the density from dropping below 10 16 m −3 . For GRILLIX, point-wise adaptive sources are used to prevent the density from dropping below 5 × 10 17 m −3 and the electron and ion temperatures from dropping below 2eV. For TOKAM3X, additional sources were not required in this work.

Constraints and free parameters
Since the simulations self-consistently evolve the plasma profiles as well as the turbulence, they have only a few free parameters which can be tuned to match the experiment. The most significant free parameters are the power and particle source rates. In all codes, the density source rate is adjusted such that the separatrix value of the simulated density profile approximately matches the separatrix value measured by Thomson Scattering. In GBS and GRILLIX, both density and temperature sources act as sources for power (since adding particles at non-zero temperature requires energy). The total power added to the plasma is where n, T e , T i , S n , S Te and S Ti are all functions of R, Z and φ, and the ionisation energy, stored in the plasma as potential energy, is not included here. In GBS, the electron temperature source rate is adjusted such that the separatrix T e value at the outboard midplane matches the value measured by Thomson Scattering, and the ion temperature source rate is set to 25% of S Te . In GRILLIX, S Ti = 0 and S Te = 1 n (S P − (T e + T i )S n ), which simplifies Eq.11 to P = 3 2 S P d 3 V . Therefore, a negative T e source is applied at the n source position to maintain a constant power. The S P can be considered a power source for electrons, which is adjusted to match the Ohmic power, and the ions are heated via the equipartition term. In TOKAM3X, sources are used for the pressure instead of for the temperatures, such that the total power is P = 3 2 (S pe + S pi )d 3 V . Equal pressure sources are used for the electrons and ions, and the electron pressure source is adjusted such that the separatrix T e value matches the value measured by Thomson Scattering.
Therefore, in all simulations, the density value at the separatrix should approximately match the Thomson scattering separatrix value as a result of tuning, while the rest of the profile is free to vary. Additionally, in GBS and TOKAM3X, the electron temperature value at the separatrix should match the Thomson scattering separatrix value (while the rest of the T e profile is free), while in GRILLIX the total injected power is set to approximately match the experimental Ohmic-heating power and the whole T e profile is free.
The order-of-magnitude of the resulting source rates can be compared to the experiment.
For the power injection, the simulation source rate was equivalent to 170kW for GBS, 150kW for GRILLIX and 30kW for TOKAM3X -compared to a total Ohmic-heating power of 150kW, of which approximately 120kW crossed the separatrix (estimated from a tomographic reconstruction of the radiated power measured with bolometry).
For the particle source, the simulation source rate was equivalent to ∼ 2 × 10 21 s −1 for GBS and 1.85 × 10 21 s −1 for GRILLIX and TOKAM3Xcompared to ≈ 3 × 10 21 s −1 inferred from the total out-flux to the Langmuir probes, assuming perfect recycling. Therefore, the simulated and the expected experimental source rates come out at similar ordersof-magnitude. However, the power varied by more than a factor of 5 between the simulations, despite each simulation achieving upstream separatrix T e values which are similar to the experiment. From a simple two-point model, we expect that the upstream separatrix T e will be weakly dependent on the input power (T e,upstream ∝ P 2/7 sep , Eq.5.7 in Ref. [50]) -and as such it is possible to achieve roughly the same upstream T e value with a wide range of input powers. However, for the target T e value, a stronger P 10/7 sep dependence (Eq.5.10 in Ref. [50]) is expected.

Simulations and post-processing
Each of the three codes performed simulations of the TCV-X21 scenario for a physical time of at least 2ms, allowing the sources, cross-field turbulent transport, plasma profiles and sinks to approach a dynamic equilibrium (or saturated ) state. Statistical moments were calculated over the last 1ms of each simulation, sampled at approximately 1µs intervals. In Sec.5 we compare results from both field directions for GBS and GRILLIX, while TOKAM3X performed a simulation only in the forward-field direction. During the setup of the simulation, an issue in the aspect ratio resulted in the R coordinate of the TOKAM3X being effectively shifted inwards by −25cm, which was corrected for in post-processing. In addition, the normalisation parameters of the TOKAM3X simulations were adjusted in post-processing to improve the match of the outboard midplane separatrix values of n and T e . This renormalisation can be performed consistently since the equations are implemented in a dimensionless form which retains the parametric dependencies * .
Renormalisation also changes the effective value of other physical parameters such as the resistivity, heat conductivity and source rates (the values given in Sec.4.2 and Sec.4.4 are computed after renormalisation).
A worked example showing how the observables in Tab.1 are calculated is given in TCV-X21/notebooks/ simulation_postprocessing.ipynb.

Validation
In this section, we present the overall result of the validation and show individual profiles from the experiment and simulations. We start by giving the overall quantitative result, to quickly indicate which observables agree particularly well or poorly. We then show the profiles obtained at the outboard midplane and divertor entrance (Sec.5.1), which are found to give reasonable agreement. This is contrasted to the divertor target profiles (Sec.5.2), where a reduced level of agreement is found. Finally, we show the divertor volume profiles from the RDPA (Sec.5.3). Note that due to space limitations it is not possible to show figures for all observables. Figures for all observables (and additionally the simulated ion temperature) can be found in TCV-X21/3.results.
We limit our validation analysis to the range R u − R u sep < 2.5cm for all diagnostics, and to R u − R u sep > −0.9cm for the wall Langmuir probes. This removes regions where the signal acquired by the probes is very low, which can prevent the correct fitting of * Setting the Coulomb logarithm equal to a constant is one exception. the IV curve by the four-parameter model. Additionally, this range ensures that the comparison points are on the grid of all simulations (since the codes use different radial grid extents, indicated in Fig.4). We also note that the RDPA 2D profiles are limited to Z − Z X > −0.32m, with points close to the targets cropped to avoid possible effects of the pre-sheath entrance. The data over an extended range is available in TCV-X21/1.experimental_data.
For each simulation and each observable, the normalised distance d j (Eq.1) and sensitivity S j (Eq.4) is computed. Points with a very low experimental uncertainty ∆e j /e j < 10 −3 (typically due to a lack of repeat discharges to estimate the reproducibility error) are removed from the calculation of d j . The values of d j and S j are used, together with the primacy hierarchies H j given in Tab.1, to compute the overall composite metric χ (Eq.5) and quality Q (Eq.6) for each simulation and including all observables. In addition, the effective composite metric and quality values χ diag and Q diag , are computed taking observables from a single diagnostic. The result is given in Tab.2. We find that the level of agreement computed from individual diagnostics varies significantly. Both the reciprocating midplane probe (FHRP) and divertor-entrance Thomson scattering (TS) show appreciable quantitative agreement, while the reciprocating divertor probe array (RDPA) and the divertor target profiles (HFS-LP/LFS-LP/LFS-IR) show poor agreement. This suggests that the change in agreement is due to the measurement location rather than the diagnostic itself: better agreement is found for diagnostics which are close to the confined region (where the TS values at the separatrix were used to tune the sources) than for diagnostics in the divertor volume or at the targets. To understand the quantitative result, we show comparisons of several observables, grouping the results by location.  Table 2: Quantitative validation result for each observable. For each code and field direction (indicated by (+) for forward field and (−) for reversed field), the d j ('normalised distance', Eq.1) and S j ('sensitivity', Eq.4) terms are given. The normalised distance d j gives the root-mean-square Z-score of the difference between the experiment and simulation, with green cells indicating good agreement (d j → 0) and red cells indicating poor agreement (d j → ∞, with the colour scale limited to d j ≤ 5.0). The sensitivity S j indicates the precision of each observable, with S j → 0 for observables with high relative uncertainty and S j → 1 for observables with low relative uncertainty. The combined level-of-agreement χ (Eq.5 for d 0 = 1.0 and λ = 0.5) and the comparison quality Q (Eq.6) are given for each diagnostic individually as well as for the overall validation.  Table 3: Near-SOL decay lengths (in cm), for the mean density (n) and electron temperature (T e ) measured at the outboard midplane (OMP) and divertor entrance (DE), in forward (+) and reversed (−) toroidal field direction. Profiles are fitted in the range [0cm, 1.5cm]. The observable uncertainty is included in the fitting uncertainty.

Outboard midplane and divertor entrance profiles
Experimentally, we see from Tab.3 that the OMP (FHRP) and divertor entrance (TS) give similar λ n and λ Te fall-off lengths in forward-field, while in the reversed-field the fall-off lengths at the OMP are 1.5 − 2× larger than at the divertor entrance, although with a much higher fit uncertainty. For the forward-field λ n , TOKAM3X predicts broader n profiles, while GBS and GRILLIX match the experimental fall-off length within the uncertainty. All simulations predict too broad forward-field T e profiles, with GRILLIX matching the closest, then GBS, and then TOKAM3X. In the reversed-field case, GBS reproduces the narrowing of the n and T e profiles between the OMP and divertor entrance (although to a smaller extent than in experiments), predicting λ n within the uncertainties. λ Te is slightly too large, but still within uncertainty at the OMP. Conversely, GRILLIX predicts similar fall-off lengths as in the forward-field case. It does not show a narrowing of λ n,OM P between the OMP and the divertor-entrance, but still matches all falloff-lengths except λ − Te,DE within uncertainty. We see that the higher fit uncertainty for the TCV reversedfield λ n is because the experimental n profile (shown in Fig.5.B) is not a simple exponential decay in the range R u − R u sep ∈ [0cm, 1.5cm]. Instead, a flat region around around R u − R u sep = 0.25cm is seen in both the FHRP and OMP measurements. This is not, however, reproduced in the simulations.
The use of relaxed parameters is likely to be part of the cause of the broadened profiles for TOKAM3X and GBS. However, if relaxed parameters were the only cause for broadening, TOKAM3X (which uses values closer to the Braginskii values than GBS) should predict narrower T e profiles than GBS, while the opposite is found. Additionally, GRILLIX (which uses the Braginskii parameters directly) predicts broadened forward-field T e profiles -indicating that relaxed parameters alone cannot explain the behaviour. It is likely that the lack of neutral dynamics is contributing to the broadened T e profiles, since neutral ionisation would add an additional sink of energy in the open field-line region.
Another possible cause for the different profile widths is the different energy source rates (given in Sec.4.4).
For the V pl profiles (Figs.5.E-F), we see in the experiment that the profiles are monotonically decreasing with a steeper slope in forward-field than in reversed-field. For all codes, V pl is positive in the SOL and of a similar amplitude as in the experiments, with a very good match in forward-field. In GBS, V pl peaks near the separatrix in forward-field case, and further into the SOL in the reversed-field case. This is within uncertainty for the forward-field case, while in the reversed-field case a disagreement is found close to the separatrix. In GRILLIX, the peak of V pl is at the separatrix in both field directions, such that the agreement is within the error bars in forward-field, but not in reverse-field, where V pl is decreases too steeply into the near-SOL. In TOKAM3X (forwardfield only), V pl is very flat, but still within uncertainty except in the vicinity of the separatrix. The reason for the worse agreement found for the reversed-field V pl profiles (compared to the reasonably good agreement in forward-field) is difficult to determine, since the potentials can be affected by both the confined region dynamics and the sheath boundary conditions [21], [71].
For the M profiles (Figs.5.G-H), we see in the experiment that the parallel flow changes direction

Low-and High-field-side target profiles
Overall, a worse match between simulation and experiment is found for the target profiles compared to the midplane profiles. Generally, for most observables, the codes capture the correct peak order-of-magnitude and the features visible in the experiment are (roughly) reproduced. However, the majority of observables are not matched within experimental uncertainty, and the broadness of the profiles is seen to vary significantly amongst the simulations. All simulations fail to accurately predict the V f l target profile, and under-predict the σ(V f l ) and σ(J sat ) by factors of 2 or more. Furthermore, experimentally, a significant effect of the toroidal field reversal is seen for the n, T e and V f l and J profiles at the targets; at the HFS target the plasma is colder and denser in forward-field than in reversed-field, and on the LFS target a prominent private-flux peak in n and J sat is observed in forward-field. The simulations are not able to reproduce the double peak n profile at the LFS target in forward-field, missing the primary peak lying in the PFR. Generally, the codes provide very different predictions of the divertor target profiles. As such, we present the results from each code separately, highlighting general trends as well as observables with particularly good or poor agreement.
Starting with GBS, the broadness and peak values of the profiles vary depending on the toroidal field direction and between the targets. A significant effect of the toroidal field reversal is seen in the simulated n, J sat and J profiles at LFS and HFS targets. At the LFS target, the SOL n peak value approximately matches the experiment in the reversed-field case, but the profile is broadened with respect to the experiment. Conversely, in the forward-field case, the width in the SOL matches more closely, but the PFR peak is missed. At the HFS target, the forward-field n peak value is underestimated and shifted towards the SOL, while the reversed-field profile is again broadened with respect to the experiment. Conversely, the T e and V pl profiles are much broader than the experimental profiles. This is likely due to the reduced heat conductivity, which changes the ratio between the parallel and perpendicular heat fluxes. Since the plasma potential is related to the electron temperature, it will also be broadened. Of the three codes, GBS is the only code which predicts a significantly-nonzero V f l , but the predicted profiles do not match the experiment. At least for the reversed-field case, this could be a consequence of the profile broadening.
For GRILLIX, the density at the targets is too low, while the shape is loosely recovered, except for the PFR peak observed at the LFS in forward-field. The SOL T e profile is matched well at both targets and in both field directions, while in the private-flux region the T e profile drops off too sharply. Due to the V f l = 0 ⇒ V pl = ΛT e Figure 7: Comparison of averaged profiles at the low-field-side divertor target. Rows from top: Profiles of the mean plasma density (n e ), electron temperature (T e ), plasma potential (V pl ) and the parallel current density (j ) in forward (left column) and reversed field (right column). The experimental data from the low-field-side Langmuir probe array (LFS-LP) is indicated by the blue line, with its uncertainty indicated by the blue shaded region. The other lines give the corresponding simulated profiles from the three codes. The corresponding d j values (Eq. 1) are given in the legend for each code. Profiles of the mean plasma density (n e ), electron temperature (T e ), plasma potential (V pl ) and the parallel current density (j ) in forward (left column) and reversed field (right column). The experimental data from the high-field-side Langmuir probe array (HFS-LP) is indicated by the blue line, with its uncertainty indicated by the blue shaded region. The other lines give the corresponding simulated profiles from the three codes. The corresponding d j values (Eq. 1) are given in the legend for each code. Figure 9: Comparison of direct experimental measurements and their standard deviations at the low-field-side divertor target. Rows from top: Profiles of the mean ion saturation current (J sat ) and its standard deviation (σ(J sat )), mean floating potential (V f l ) and its standard deviation (σ(V f l )) in forward (left column) and reversed field (right column). The experimental data from the low-field-side Langmuir probe array (LFS-LP) is indicated by the blue line, with its uncertainty indicated by the blue shaded region. The other lines give the corresponding simulated profiles from the three codes. The corresponding d j values (Eq. 1) are given in the legend for each code. boundary condition, the sharp drop towards the PFR in T e leads to a corresponding drop in V pl . This will in turn give a strong electric field across the separatrix. The V f l = 0 boundary condition is clearly incorrect, which causes the reversed-field SOL V pl to disagree regardless of the good match in T e . Despite applying a potential boundary condition which corresponds to an insulating (J = 0) sheath, the J profile (which is set equal to the internal currents) is both significantly nonzero and also, surprisingly, a reasonably good match to the experiment (except for the HFS profile in forwardfield). The current boundary conditions and the effect of the electric field along the target on u is discussed further in Sec.6.3.
For TOKAM3X, the HFS target is substantially better matched than the LFS target. Good matches are found for the HFS n, T e and V pl , although the far-SOL profile is too broad. For the LFS target, it is seen that the T e and V pl profiles are increasing into the far-SOL. Since there are no mechanisms to heat the electrons in the far-SOL in the TOKAM3X simulation, this appears to indicate a numerical issue such as pollution due to a buffer zone applied at the limiting flux surface. The LFS n and J sat have the correct amplitude in the SOL, but both miss a prominent PFR peak. Both V f l and J at both targets are significantly lower than observed.
We finally consider the heat flux and decay length of all codes. The parallel heat flux profiles measured at the LFS target by the infrared camera are shown in Fig.10. Since the simulations do not directly evolve the parallel heat flux, its value is calculated from the sum of the electron and ion contributions. In GBS, the parallel heat flux is computed from the parallel-convective heat flux 5 2 n v T e + u T i + 1 2 m i u 3 , where v and u are the electron and ion velocities respectively, and the last term accounts for the kinetic energy associated with the bulk ion flow. In GRILLIX and TOKAM3X the heat flux to the sheath entrance is computed from the sum of a conductive heat flux −χ ,e T 5/2 e ∇ T e − χ ,i T 5/2 i ∇ T i and the total convective heat flux 5 2 n (vT e + uT i ) ·b, for v and u the total electron and ion velocities (the vector-sum of the parallel, E × B and -for TOKAM3X only -the ∇B components of the velocity). Note that the conductive component for GRILLIX and TOKAM3X is intended to mimic the electron cooling by the sheath, and that the χ values used are the Braginskii values, which do not account for the modifications of the distribution function due to the wall [74]. The calculation of the heat flux for GRILLIX is given in TCV-X21/tcvx21/ grillix_post/observables/heat_flux_m.py.
We see that GRILLIX predict the peak heat flux in forward-field reasonably well while GBS under-predicts Figure 10: Parallel heat flux profiles and Eich profile fits at the low-field-side target. The experimental parallel heat flux data from the infrared camera for the low-field-side target is indicated by the blue line, with its uncertainty indicated by the blue shaded region. The other solid lines give the corresponding simulated profiles from the three codes, and the corresponding d j values (Eq.1) is given in the legend for each code. For each profile shown, the corresponding Eich-type fit [73] is indicated by a dashed line of the same colour. The fitted parameters are given in Tab.4. Profiles are fitted over the entire range of data.
the peak heat flux. In reversed-field both GBS and GRILLIX predict a peak heat flux which is less than half the measured value. By comparison, in forwardfield TOKAM3X predicts an extremely broad, lowamplitude heat flux profile, with a peak heat flux ∼ 25% of the measured peak heat flux. We use a Levenberg-Marquardt algorithm [75] to fit Eich-type profiles of the form [73] where r = R u − R u sep and λ q , S, q 0 , q BG and r 0 are fitted parameters corresponding to the (upstream mapped) heat flux decay length, the power spreading factor, the peak heat flux, the background heat flux, and a free radial shift. The parameter r 0 was introduced to account the uncertainty in the spatial calibration of the IR, which can change the strike-point position, affecting the fit of the Eichprofiles. We found that in forward and in reversedfield, respectively, r 0 = −2.1mm and r 0 = −0.6mm. The fitted profiles are indicated by dashed lines in Fig.10, and the fitted parameters for λ q and S are given in Tab.4. During fitting, the experimental uncertainty is neglected since this is found to improve the match between the raw and fitted profiles. The fitting procedure is provided in TCV-X21/tcvx21/analysis/ fit_eich_profile_m.py.  Table 4: Eich profile fitted parameters (in mm) for the mean parallel heat flux profile, for the heat flux decay length (λ q ) and power spreading factor (S), in forward (+) or reversed (−) toroidal field direction.
For the simulations, the simulated double-peak seen by GRILLIX for the forward-field clearly cannot be described by Eich-type fits. In the forward-field direction, GBS reproduces the experimental spreading factor S, but predicts a larger heat-flux decay length λ q , while GRILLIX predicts smaller values for both parameters. For TOKAM3X, the fitted parameters have very high uncertainty, indicating a low quality fit. In the reversed-field direction, GBS again reproduces the experimental spreading factor and predicts a larger λ q , while GRILLIX predicts a smaller spreading factor and a slightly larger λ q . The low-amplitude TOKAM3X profile may be related to the very low injected power in that simulation (see Sec.4.4). The broadened GBS heat flux profiles are consistent with the profiles at the LFS target. For GRILLIX, the smaller S in both field directions is consistent with the interpretation of narrowed profiles due the fast parallel advection, discussed in Sec.6.3.

Divertor volume profiles
In Fig.11 we show 2D profiles of the mean plasma density (n), plasma potential (V pl ), parallel Mach number (M ) and ion saturation current density (J sat ) fluctuations in the divertor volume, comparing the experiment and simulations.
The mean, we find that M is subsonic throughout the entire divertor leg, with M ∈ [−0.58, 0.67] in forwardfield and M ∈ [−0.78, 0.3] in reversed-field. This may indicate the presence of a significant neutral ionisation source, which will reduce parallel flows [77]. In forward-field, the flow is directed towards the target, on the order of ∼ 10 − 20% of c s in the SOL and ∼ 50% of c s in the PFR. In reversed-field, a flow away from the target is found for most of the divertor leg, which may indicate a returning flow for the strong poloidal E × B flow directed towards the target. Flow magnitudes ∼< 20% of the c s are seen in the SOL (some weak flow towards the target is also observed), while in the PFR the flow reaches 80% of c s directed away from the target near the X-point. .3G) the simulated parallel flows greatly exceed the local sound speed, while the direction of the parallel flow in the PFR appears to roughly correspond to the measurements. The fast parallel flows may be related to a feedback loop affecting the E × B-drift correction u boundary condition used in GRILLIX, which is discussed more in Sec.6.3. For TOKAM3X (Fig.11.3D), the effect of the E × B drift is weaker than in GRILLIX, consistent with the weak divertor V pl gradients (see Fig.11.2D), such that the values and direction of M are dominated by the drift-free component of the Bohm-Chodura condition. The measured J sat fluctuations are approximately constant along the divertor leg (Figs.11.4A and 4E), with a moderate increase observed towards the target in forward-field and near the X-point in reversed-field. Radially, the fluctuation intensity shows the same pattern as the other fields, with a broader profile in forward-field and a marked peak near the separatrix in the reversed-field case. The relative fluctuation level σ(J sat )/J sat displays a marked hollowed profile in both field directions (not shown here), with a typical σ(J sat )/J sat ∼ 10% around the separatrix, a σ(J sat )/J sat ∼ 50% in the far-SOL, and a σ(J sat )/J sat ∼ 100% deep in the PFR. We also remark that the region with low σ(J sat )/J sat is broad and extends from the SOL into the PFR in forwardfield, while in reversed-field, the region of low relative fluctuations is only in the SOL and sits in a narrow region around R u − R u sep ≈ 0.1cm. In contrast, the J sat fluctuations (the absolute fluctuations) in the simulations are significantly smaller, with the fluctuation profiles showing a strong drop when moving from the divertor entrance to the LFS target (see Figs11). We find that the lower simulated absolute fluctuation levels also translate to lower σ(J sat )/J sat levels (not shown here). In TOKAM3X, σ(J sat )/J sat < 10% throughout the divertor. In GBS, the simulations in both field directions recover the qualitative shape of the experimental σ(J sat )/J sat profiles in the SOL. The region of reduced relative fluctuations along the divertor leg is found with σ(J sat )/J sat ∼ 10% as well, as the increase in the radial direction. In the reversedfield simulation, this is up to ∼ 50% in the region at Z − Z X ≈ 0m and ∼ 15% at Z − Z X ≈ −0.35m. Conversely, the forward-field simulation shows a less steep poloidal gradient with σ(J sat )/J sat ∼ 15%−25%, where the maximum and minimum values are found, respectively, at Z − Z X ≈ 0m and Z − Z X ≈ −0.35m. In GRILLIX, in forward-field σ(J sat )/J sat ∼ 50% is seen localised near the separatrix in the PFR, while σ(J sat )/J sat < 10% is seen in the SOL. In reversedfield, a peaked relative fluctuation profile is observed with a σ(J sat )/J sat ∼ 25% along the separatrix, while the rest of the divertor displays σ(J sat )/J sat < 10%. We find that GBS simulations show a reasonably good quantitative agreement for the skewness and kurtosis of J sat , as shown in Tab.2, suggesting the blob dynamics are being reasonably well reproduced [78].

Discussion
We now consider the overall outcome from the quantitative and qualitative validations. The immediate result from both analyses is that outboard midplane and divertor entrance profiles are captured to a reasonable degree by all codes, whereas the divertor volume and target profiles show a worse agreement. The disagreement between simulation and experiment is particularly apparent for the parallel Mach number in the divertor volume (Figs.11.3), the standard deviation of J sat at the divertor targets (Figs.9.C and 9.D) and the floating potential at the divertor targets (Figs.9.E and 9.F). For observables where the simulations disagree with experiment, we find that the simulations also typically disagree with each other, which indicates that the simulations are sensitive to the differences between the models, model parameters, and numerical parameters.
In Sec.6.1 we discuss the result of the quantitative validation, and the observed sensitivity of the simulated profiles to the numerical resolution. Then, in Sec.6.2, we discuss how the choice of physical parameters was found to affect the results. In Sec.6.3, we discuss the sensitivity of the simulations to the choice of sheath boundary conditions -particularly for the parallel ion velocity, parallel current and floating potential. In Sec.6.4, we discuss the effect of toroidal field reversal, contrasting the results to previous transport modelling. We then close our discussion of the simulations in Sec.6.5 by assessing what additional physics might be required to match the experimental results, such as neutral dynamics, and to provide predictions beyond TCV-X21. Finally, in Sec.6.6 we suggest minor extensions to the experimental reference dataset, which could further constrain the models and test additional capabilities of the simulations. The results of the quantitative analysis correspond to a lower level of agreement and a significantly higher quality compared to previous works using the same methodology, such as Ricci et al., 2015 [25], which found a χ = 0.5 − 0.75 with a Q ≈ 4 for TORPEX plasmas and Riva et al., 2020 [17] which found a value of χ = 0.45 and Q ≈ 4 for limited TCV discharges. We note that a companion work to this paper, published as Galassi et al., 2021 [26] found χ = 0.85−1.0 and Q ≈ 4 in TORPEX plasmas with an internal X-point. The reduced level of agreement found in this work is partly due to our assumption of zero simulation uncertainty, in contrast to the previous studies which -except for Ref. [26] -estimated non-zero simulation uncertainties. We can explore the impact of non-zero ∆s j by setting the simulation uncertainty ∆s j to the experimental uncertainty ∆e j , which gives χ between 0.7−0.87. This is closer to the range of the previous studies, although still at a lower level of agreement. Beyond reducing our level of agreement, we see that our assumption of zero uncertainty also causes the metric to occasionally give non-intuitive results. This is particularly noticeable when the experimental uncertainty varies strongly across a profile, since low uncertainty points effectively dominate the quantitative level of agreement. To avoid this, we could introduce a finite ∆s j , or alternatively we could explore alternative metrics such as Gaussian process regression [79].
To estimate the simulation uncertainty, we could use error-estimation methods such as a sensitivity analysis [22] or Richardson extrapolation [47]. However, such techniques require repeat simulations, and due to computational cost of the simulations, a rigorous error analysis was not performed within this work. Instead, GRILLIX was used to explore the effect of varying the poloidal resolution. We performed a resolution scan, comparing simulations at 2mm, 1mm and 0.5mm perpendicular resolutions. At the OMP, similar mean profiles were found across the resolution scan, while a phase-shift analysis showed an increase of drift-wave turbulence relative to ballooning turbulence with increasing resolution. In the open field-line region, higher fluctuation levels (particularly in the PFR) and more transport into the far-SOL was found with increasing resolution. This suggests that a resolution of a few mm is sufficient to resolve most of the primary turbulencedrive by confined-region instabilities, while higher resolutions are required to capture more of the secondary turbulence driven locally in open field-line regions with low temperature (or alternatively, small drift scales). We also investigated the statistical uncertainty of the GRILLIX results by using the bootstrap method [80]. This indicated that higher-order statistical moments, such as skewness and kurtosis, had significant statistical uncertainty, while the mean and standard deviation had negligible statistical uncertainty.

Impact of relaxed parameters
Whereas GRILLIX used the Braginskii values for resistivity and heat conductivities directly, relaxed parameters were used in GBS (η increased by factor 3 and χ ,e,i reduced by factor between 4.8 (at the targets) and 20 (at the OMP)) and TOKAM3X (η increased by factor 1.8 and χ ,e,i reduced by factor 1.8) in order to reduce computational costs and improve the numerical stability. Previous studies have indicated that artificially increased collisionality leads to broadened SOL profiles [13], [16], larger blobs and higher fluctuation levels [13]. Increased resistivity can directly affect the plasma potential (affecting the scale of potential fluctuations and the magnitude of the E × B transport) [13], extend the inertial regime of SOL blobs [81] and change the position where blobs are generated [82]. Reduced heat conductivity leads to a higher fraction of the heat flux being convected by parallel flows [83], broader profiles [84] and allows larger parallel temperature gradients to develop.
In the poloidal density profile (Fig.4), we see that the codes predict noticeably different blob sizes -large for GBS, intermediate for TOKAM3X, and small for GRILLIX -which is the ordering expected if physical parameters were the dominant factor determining blob size. Comparing the profiles, the difference between the codes at the OMP is small and within experimental uncertainty, while at the divertor targets the differences between the codes are more pronounced. In general, GBS approximately matches the width of the n, J sat , and J profiles in forward-field, while in reversed-field the profiles are slightly broadened. Conversely, the reversed-field simulations often predict the peak value more accurately than the forwardfield simulations. For GRILLIX, the target profiles are noticeably narrower, which is likely due to the fast parallel advection (see Sec.6.3). For TOKAM3X, the LFS target profiles returns a low quantitative and qualitative agreement, while the HFS target gives appreciable agreement -which may indicate that the effect of relaxed parameters is increased for longer leg lengths. However, we note that the electron temperature profile increasing into the far-SOL for the TOKAM3X low-field-side target profile cannot be explained by relaxed parameters, and may indicate numerical pollution via the buffer zone applied on the limiting flux surface. It is also unclear whether the use of relaxed parameters can explain the reduced fluctuations in the divertor leg, since in Ref. [13] fluctuations were found to increase with increasing collisionality. Here, it is likely that grid resolution (and the numerical discretisation scheme) is having an impact, but a combined sensitivity analysis and resolution scan was not performed in this work. It is also unclear which of the increased resistivity or the reduced heat conductivity is having a more significant effect on the target profiles. A parameter scan would be useful to determine the effect of relaxed parameters on the simulated target profiles, and reveal whether different parametric dependencies are found in limiter and divertor geometries.

Influence of sheath boundary conditions
The codes use different sets of sheath boundary conditions, detailed in Sec.4 and Appendix A. One effect of the different sheath boundary conditions can be directly seen in the simulated parallel Mach number (M ) in the divertor volume, shown in Fig.11.3. Here, GBS shows the expected flow profiles for a sheathlimited regime without drifts or ionisation in the divertor [77]. The Mach number is large relative to the experiment throughout the divertor, approaching M = 1 at the target, and similar in both forward and reversed field. This is consistent with the u = c s boundary condition, but it does not match the direction of the measured M profile, which is seen to be dependent on the toroidal field direction. The effect of the toroidal field reversal on M appears to be matched better by the E × B-drift corrected u boundary conditions used by GRILLIX, which are described in Ref. [66]. However, the GRILLIX simulations find excessively fast parallel flows with |M | 1). This might indicate a feedback mechanism, where the parallel flows affect the T e profile. Due to the V pl = ΛT e boundary condition, this in turn affects the poloidal E × B drift and the u boundary condition. It appears that this leads to self-steepening mechanism, eventually leading to very fast parallel ion velocities. The increased parallel flow rate in GRILLIX may explain part of the disagreement in other profiles. Faster parallel advection will lead to a lower target n for a given upstream n, and by modifying the balance of parallel and perpendicular transport it will also lead to narrower profiles (such as Fig.10).
The boundary conditions also set the potential and current at the sheath. For V pl particularly, since it is solved via an elliptic equation, the potential boundary condition can have a significant effect throughout the open-field line region [21], [64]. In GBS and TOKAM3X, the potential and current boundary conditions are coupled in the electrostatic limit, giving ) and a coupled V pl boundary condition. In GRILLIX, the current was allowed to freely flow and the potential was set to ΛT e . Considering the target V f l profiles, the GBS potential boundary condition is seen to give a non-negligible V f l (with a reasonable magnitude, albeit a non-matching shape), while in GRILLIX and TOKAM3X V f l ≈ 0. The GBS current profile approximately matches the shape (but not the amplitude) in the forward-field, while in reversed-field it is broadened and shifted with respect to the experiment (likely due to the used of relaxed parameters). In contrast, the freeflowing boundary condition used by GRILLIX gives a reasonable match to the experiment for J (except at the forward-field HFS), which is surprising since this boundary condition does not have a theoretical basis. To investigate this further, a GRILLIX simulation with an insulating current (J = 0) boundary condition was performed. It was found that the interior currents were similar to the free-flowing current simulation, but a very strong heating near the boundaries was caused by the compression required to force J → 0. Therefore, it appears that in GRILLIX the currents observed at the boundary are driven internally, independently from the boundary conditions. This result warrants further investigation, to determine whether the match to the experiment is fortuitous and -if not -what is driving the current internally. Finally, TOKAM3X predicts very low currents at the boundary. This gives a reasonable match at the HFS target (where the other profiles also agree well), but not at the LFS target, potentially due to the disagreement in the other profiles.

Toroidal field reversal and in-out asymmetry
The codes are seen to mostly under-predict the effect of the toroidal field reversal and the in-out asymmetry between the low-and high-field-side divertor targets. This could be due to an overall underestimation of large-scale plasma drifts, which reverse with toroidal field direction [5], [76], [85], or due to a lack of symmetry-breaking terms such as ion-orbit-loss [86]. Here, it is interesting to compare to previous modelling of TCV with the UEDGE transport code [76], which found that the HFS target profile was ∼ 3× colder and denser in forward-field than in reversed-field, primarily due to E × B drifts. This effect was also observed experimentally, although it appears that UEDGE was overestimating the effect of the background drifts, in contrast to this study where it appears that we are underestimating their effect. As such, it would be interesting to compare turbulence and transport simulations of the TCV-X21 scenario, to compare the E × B drifts, as well as to provide additional information about the neutral dynamics and power and particle sinks in the divertor.

Towards an improved match
In future works, we expect that the match to experiment should improve when GBS and TOKAM3X use more realistic values for the resistivity and heat conductivity. For all simulations (and for GRILLIX especially), the overall match should improve if the parallel ion velocity can be reduced towards experimental values, which might be achieved by using other boundary conditions for v and/or V pl . Theoretical work to identify a consistent coupled set of boundary conditions applicable to the full system of equations (including electromagnetic and ion thermal effects, potentially extending on Ref. [64]) could significantly improve the fidelity of the simulations in the divertor.
To continue improving the quantitative match between experiment and simulation will likely also require additional physics to be considered. Despite targeting a sheath-limited low-recycling regime in the TCV-X21 scenario, it is nevertheless likely that the most significant missing physics term is the neutral dynamics. For instance, the Mach number measurements in the divertor volume (Fig.11.3) suggest that there is a non-negligible neutral ionisation source in the divertor leg [77]. By adding neutral dynamics, we expect that the simulated parallel advection in the divertor will be reduced. The neutrals will also introduce plasma cooling in the vicinity of the X-point and an asymmetry in the density source, which can change the profiles and, thus, the turbulence drive [87].
As such, new versions of the codes are currently targeting simulations of the TCV-X21 scenario. For GBS, a new version has recently been developed. This version does not use the Boussinesq approximation, includes electromagnetic effects and couples to a kinetic neutrals model [67]. For GRILLIX, a refactored version will be used for testing boundary conditions, including neutrals and extending the resolution scan [87].
For TOKAM3X, or rather its successor SOLEDGE3X [44], the Boussinesq approximation will not be used, neutral and impurity physics will be provided via a coupling to EIRENE, and additional terms for the parallel viscosity and a complete Reynolds-stress tensor will be included.
Beyond TCV-X21, further extensions of the model beyond the drift-reduced fluid approach might be necessary, especially for future validations against larger experimental devices and more challenging plasma conditions. Kinetic corrections to the parallel heat flux [88], ion orbit loss effects [86] and finite Larmor radius corrections (in the spirit of gyrofluid or velocity-space-decomposed-gyrokinetic [89] models) may become important. In this context, comparisons to gyrofluid and gyrokinetic simulations of the TCV-X21 reference case and future divertor validation cases are of great interest.

Summary and conclusions
The predictive capabilities of divertor turbulence simulations was rigorously assessed via a validation against a series of dedicated diverted TCV discharges in both forward and reversed field direction. The discharges were performed at a lower toroidal magnetic field than typical for TCV (0.95T vs 1.45T) to decrease the computational cost, thus allowing direct, full-size simulations of the experimental scenario. Moreover, the discharges were carried out in low density, sheath-limited conditions, designed to reduce the effect of the neutrals, which are currently not included in the simulations. The discharges were repeated several times to improve the statistics of the experimental measurements and to investigate the effect of toroidal field reversal. An extensive experimental dataset, which we refer to as TCV-X21, was collected for the purpose of validating the codes, which includes a broad range of 1D and 2D measurements of fluctuation and time-averaged quantities at the outboard midplane, at the divertor entrance, throughout the divertor volume, and at the target plates. The validation dataset is provided in an open repository, as a reference for future validations and to allow benchmarking of boundary turbulence simulations.
Full-size simulations of the TCV-X21 reference case, in both toroidal field directions, were performed by three separate edge turbulence modelling groups, using the GBS, GRILLIX and TOKAM3X codes. The simulations were flux-driven, meaning that the plasma profiles, turbulence and transport were evolved self-consistently. As such, the simulations were validated against the full set of experimental measurements, including both mean profiles and fluctuations. The only tunable physics parameters in these simulations were the position and strength of the density source and of the temperature or power source. The source positions were chosen to approximately match the positions of Ohmic power deposition and neutral ionisation, and the rates were adjusted to approximately match the experimental outboard midplane density and temperature (or power) at the position of the separatrix. To evaluate the quality of the match between the simulations and experiment, the different observables were compared graphically and via a quantitative validation metric.
An appreciable match between simulation and experiment was found, particularly at the outboard midplane. The simulations were able to match several of the outboard midplane profiles, including both mean profiles and fluctuations. By comparison, the simulation results lie mostly outside the experimental error bars at the divertor targets and in the divertor volume, while the order-of-magnitude and approximate shape of several profiles agree with the experiment. Additionally, although at the outboard midplane the three codes provided similar results, the difference between the simulations increased towards the targets. Since the simulations adjusted their sources for the outboard midplane separatrix values, the reduced agreement in the divertor was partially expected. However, it also appears that there are additional physical and numerical effects which are becoming important in the divertor.
By comparing the simulations against the experiment and each other, we find possible causes for the reduced agreement in the divertor, which may indicate how the simulations could be further improved. We see that simulations using relaxed parameters -artificially increased resistivity or reduced heat conductivity -found broadened divertor target profiles. The divertor profiles were also seen to be sensitive to the energy source rate, which is consistent with results from a twopoint model. Simulations which resolved the local drift scale found increased fluctuations away from the confined region. Since computational cost is the main limit on the choice of parameters and resolution, continued improvement of the numerical efficiency and scalability of the simulations should help to improve the match to the experiment.
The simulated divertor profiles are found to be sensitive to the choice of sheath boundary conditions applied at the divertor targets. The use of driftcorrections in the parallel ion velocity boundary condition was found to give much faster velocities than the experiment, while simulations which neglected drift corrections found a lower flow speed but did not reproduce the flow pattern observed in the divertor volume. Additionally, the codes differed in their description of the current crossing the sheath. Here, a simple 'free-flowing' boundary condition was found to give reasonable agreement when compared to the experiment or to simulations employing the usual Bohm-current boundary condition. In both cases, the tight coupling of the boundary conditions and profiles prevent generalisation of the results. Instead, we highlight that further work in investigating sheath effects and identifying an optimal, numerically-stable set of boundary conditions would be highly beneficial for divertor simulations.
Beyond this, it is expected that additional physics is required to achieve a quantitative match between simulation and the TCV-X21 reference. The addition of neutral physics is expected to be particularly important -since despite experimental efforts to reduce the divertor neutral pressure, our simple treatment of the neutrals as acting as only a confined-region density source neglects several important effects such as localised ionisation near the X-point and in the far-SOL. As such, the codes will target repeat simulations of the TCV-X21 scenario including neutral physics, more realistic physical parameters and improved boundary conditions. This validation shows that the TCV-X21 reference case developed in this work allows for the rigorous validation and bench-marking of turbulence simulations. The results of this first validation indicate that turbulence simulations are already providing a promising match to the experiment, which can be further improved by targeted development of the codes. This validation methodology can be extended to other codes since the TCV-X21 reference data has been publicly released. Additionally, similar validations could be performed at more challenging plasma parameters, with neutral physics included, in more complex magnetic geometries, at larger magnetic field strengths and on larger machines. Continued validation of turbulence simulations over a broad parameter space will accelerate the development of the codes, to improve their usefulness for interpreting results from existing machines, and eventually to enable predictive simulations of future fusion reactors. This work was supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project IDs s1126, s1129, s1028 and s882. The simulations in this work were performed within the 4th and 5th cycles of MARCONI-FUSION HPC. The required computational resources for GBS were allocated under the projects DIVturb and XTurb. We acknowledge PRACE for awarding us access to SuperMUC-NG at GCS@LRZ, Germany. The required computational resources for GRILLIX were allocated under the projects FUA34 EST3D, FUA35 EST3D and FUA35 TSVV3. TOKAM3X simulation were additionally performed on HPC systems provided by GENCI.

A. Boundary conditions
The model equations for each of the codes can be found in Giacomin and Ricci, 2020 [61] for GBS, in appendix A of Zholobenko et al., 2021 [21] for GRILLIX and in Tatali et al., 2021 [13] for TOKAM3X. The models equations require the use of boundary conditions at the edge of the computational domain. Of these, the choice of boundary conditions applied at divertor targets is particularly impactful, since these boundary conditions are used to include sheath effects in the simulations.
The sheath is known to have a significant effect on the plasma dynamics, affecting the plasma flows, currents and potential. However, the assumptions in the fluid approximation break down in the sheath (and magentic presheath [64]), and so the sheath cannot be self-consistently modelled by fluid codes (or by gyrokinetics -self-consistent sheath modelling requires fully-kinetic 6D simulations [94]). Instead, in this work, 'sheath boundary conditions' aim to mimic the effect of the sheath. The choice of boundary conditions for each code is given here.

A.1. GBS boundary conditions
GBS uses a set of magnetic pre-sheath boundary conditions for u , v , n, V pl , T e and for the vorticity ω = ∇ 2 ⊥ V pl [64], [67] that have been rigorously derived via a first-principles analysis of the sheath dynamics and which are in agreement with kinetic simulations of this region. In our simulations, we use the same boundary conditions as in Giacomin et al. [61], which neglect corrections due to the gradients of density and plasma potential along the wall.

A.2. GRILLIX boundary condition
For GRILLIX, the parallel velocity is set according to a flow-reversal-limited Bohm-Chodura boundary condition, written as where f + means the nearest parallel-neighbouring point in the direction towards the main plasma volume, b is the parallel unit vector andn is the wall-normal unit vector. This boundary condition modifies the parallel velocity to prevent inwards E × B-drifts across the boundary. The projection into the wall-normal directionn is to prevent cross-field flow oblique to the parallel direction, although we note that theb ·n projection of the parallel velocity can lead to very fast flows. We additionally set free-flowing boundary conditions for the parallel current density J = J + and density n = n + . The electrostatic potential is set to the floating potential V pl = ΛT e for Λ = 2.69 the sheath potential drop. No boundary condition is used for the parallel electromagnetic potential A . Finally, for the electron-and ion-temperatures we set sheathheat-transmission boundary conditions of the form ∇ log(T e ) = −γ e χ ,e nu (14) whereγ e = 2.5 is the anomalous electron sheath-heattransmission factor. An equation of the same form is used for the ions replacing T e with T i and using γ i = 0.1.

A.3. TOKAM3X boundary conditions
In TOKAM3X, the boundary conditions are developed in Tatali et al., 2021 [13]. Bohm-Chodura boundary conditions are enforced for the parallel velocity of ions: where the total first order ion drift velocities (ExB plus curvature) are taken into account, c s = Te+Ti mi is the local acoustic velocity andn is the normal to the surface in the outgoing direction. Here u + stands for the velocity at the nearest discretization point in the poloidal direction. The boundary condition on the parallel current J links it to the dimensionless plasma potential V pl and the floating sheath potential drop ΛT e = −0.5 ln 2π me mi 1 + Ti Te T e according to the following relation [50]: Finally, heat fluxes follow the standard sheath boundary conditions derived in [50], very similarly to what is used in GRILLIX: whereΓ i is the total particle flux andq e/i is the total heat flux for electrons (e) or ions (i). The sheath heat transmission factors were set respectively to γ e = 4.5 and γ i = 2.5 for these simulations.

B. Nomenclature, units and sign convection
The sign convection of the parallel Mach number M and the parallel current is defined to match the experimental measurements. For the parallel current, measured at the divertor targets, positive values indicate current flowing into the targets and negative values indicate values flowing from the targets. For the parallel Mach number, measured via immersed probes, the sign convention is such that positive M indicates flows towards the low-field-side divertor target, and negative M indicate flows towards the high-field-side divertor target.
Simulation codes -GBS: A fully-non-aligned 3D fluid turbulence code (developed at SPC, Lausanne) -GRILLIX: A locally-field-aligned 3D fluid turbulence code (developed at IPP, Garching) -TOKAM3X: A flux-aligned 3D fluid turbulence code (developed at CEA, Cadarche) -LIUQE: An equilibrium reconstruction code for TCV Diagnostics and locations -FHRP: Horizontally-reciprocating probe at the outboard midplane -RDPA: Reciprocating Divertor Probe array -LP: Langmuir Probe (in TCV, wall embedded Langmuir Probes) -TS: Thomson Scattering system -LFS-LP Low-field-side divertor target Langmuir probes -HFS-LP: High-field-side divertor target Langmuir probes -LFS-IR: Infrared camera measurements at the low-field-side divertor target -OMP: outboard mid-plane -DE: divertor entrance (TS measurement position) -SOL: scrape-off-layer -PFR: private flux region -LFS: low-field-side, or low-field-side divertor target -HFS: high-field-side, or high-field-side divertor target Plasma quantities (and base units) -R: radial displacement from the axis of rotational symmetry, in metres, or radial direction -Z: vertical displacement from magnetic axis, in metres, or vertical direction -φ: toroidal direction -: parallel-to-magnetic-field direction -R u − R u sep : flux surface label, giving the radial distance between the flux surface and the separatrix, in metres -Z − Z X : vertical coordinate relative to the Xpoint, in metres -B: magnetic field vector -B: magnetic field strength, in tesla -I p : plasma current, in amperes -ρ s : sound drift scale, in metres -β: the ratio of kinetic to magnetic pressure -η : parallel resistivity -χ ,e : parallel electron heat conductivity -χ ,i : parallel ion heat conductivity -n: plasma density for electrons (n e ) and ions (n i ), in particles per cubic metre -Λ: sheath potential drop, in units of per-coulomb -T e : electron temperature, in electron-volts -T i : ion temperature, in electron-volts -V pl : plasma electrostatic potential, in volts -J sat : ion saturation current density, in ampereper-square-metre -V f l : floating potential, in volts -J : parallel current density, in ampere-per-squaremetre -q : parallel heat flux, in watts-per-square-metre -λ q : parallel heat flux decay length, in metres -M : parallel advective velocity, normalised by the local sound speed -c s : sound speed, in metres-per-second -u : the parallel ion velocity, in metres-per-second -v : the parallel electron velocity, in metres-persecond -A : the parallel component of the electromagnetic vector potential, in tesla-metre -P : power, in watts. P sep is the power crossing the separatrix, in watts total uncertainty from the experiment and simulation respectively, for some observable j. ∆e j consists of the uncertainty related to fitting experimental data to a model (∆e f it ), inherent diagnostic uncertainty (∆e dia ) and the uncertainty related to the reproducibility across repeat discharges (∆e rep )