Validation of SOLPS-ITER simulations against the TCV-X21 reference case

This paper presents a quantitative validation of Scrape-Off Layer Plasma Simulation-ITER (SOLPS-ITER) simulations against the TCV-X21 reference case and provides insights into the neutral dynamics and ionization source distribution in this scenario. TCV-X21 is a well-diagnosed diverted L-mode sheath-limited plasma scenario in both toroidal field directions, designed specifically for the validation of turbulence codes (Oliveira, Body et al 2022 Nucl. Fusion 62 096001). Five new, neutrals-related observables are added here to the extensive, publicly available TCV-X21 dataset. These are three deuterium Balmer lines in the divertor and neutral pressure measurements in the common and private flux regions. The quantitative agreement metric used in the validation is combined with the conjugate gradient method to approach the SOLPS-ITER input parameters that return the best overall agreement with the experiment. A proof-of-principle test of this method results in a modest improvement in the level-of-agreement; the shortcomings impacting the result and how to improve the methodology are discussed. Alternatively, a scan of the particle and heat diffusion coefficients shows an improvement of 10.4% in the level-of-agreement, approximately twice as high as that achieved by the gradient method. This result is found for an increased transport coefficient compared to what is usually used for TCV L-mode plasmas. The simulations further indicate that ∼65% of the total ionization occurs in the SOL, from which ∼70% in the divertor regions, despite being a sheath-limited regime, motivating the inclusion of self-consistent neutral models in future turbulence simulations on the path towards improved agreement with the experiment.


I. INTRODUCTION
The power exhaust problem is one of the key challenges faced by the magnetic confinement fusion community.Significant progress has been achieved with the introduction of diverted magnetic geometries such as the single-null configuration planned for ITER 1 .In a single-null plasma, a poloidal magnetic null (X-point) is introduced in the plasma boundary, localizing the hot core plasma some distance away from the vacuum vessel wall and directing the "open" magnetic field lines and thus the majority of the heat flux to specially designed plates called divertor targets.However, it will still be challenging to keep the heat load deposited onto the target plates of future devices such as ITER 2 and SPARC 3 , and even more so of a DEMO 4 within technological limits.To address this issue, it is necessary to have reliable modeling of the transport and the underlying physical process governing the "open" field line region called the scrape-off layer (SOL).
The investigation of the plasma dynamics in the SOL is difficult due to the limited diagnostic access and the complexity of the most complete theoretical models.Validation of numerical simulations against experimental observables stands as a suitable methodology to test and improve the capabilities of the current models 5,6 .In this approach, the largest possible number of experimental observables are quantitatively a) yinghan.wang@epfl.chb) See Reimerdes et al. 2022 Nucl.Fusion 62 042018 compared to the simulation results, and the overall level of agreement is quantified by an agreement metric, i.e. a single numerical value denoted by χ.Such a validation provides a robust framework to assess the current agreement, effect of model improvements and, when agreement is judged satisfactory, provides confidence in the code outputs and their predictive capabilities.
The modeling of the tokamak boundary plasma can be carried out by edge plasma transport codes such as SOLPS-ITER 7 , SOLEDGE2D 8 , UEDGE 9 or edge turbulence codes, as in the recent multi-code validation involving GBS 10 , GRIL-LIX 11 , and TOKAM3X 12 , which presented the first full tokamak size edge turbulence simulations 13 of a diverted single-null plasma in the Tokamak à Configuration Variable (TCV) 14 .The scenario simulated by the turbulent edge codes in this study, referred to as TCV-X21 reference case, is a Lmode plasma in sheath-limited condition to minimize the effect of the neutral dynamics in the divertor volume, as these effects were not included in the turbulence simulations.Instead, the ionization source was prescribed as an input and assumed to be localized in the outer region of the confined plasma.However, a non-negligible effect of neutrals was suggested as a possible cause of the relatively low level of agreement in the divertor region and at the divertor targets in the simulation-experiment comparison 13 .
In this work, we validate the SOLPS-ITER code against the TCV-X21 scenario, as a first step to understand the role of the neutrals in this case, with kinetic neutrals simulated with EIRENE.For this purpose we extend the publicly avail-able 13 TCV-X21 dataset (Tab.I) with two new observables for the neutral dynamics: Balmer line intensities measured by the Divertor Spectroscopy System (DSS) 15 and neutral pressure measurements from the Baratron gauges 16 .We conduct a similar quantitative validation as in Ref. 13, using the methodology of Ref. 17 (briefly reviewed in Sec.IV).With this, we expect to help guide future edge turbulence simulations of the TCV-X21 reference case including neutral dynamics.We also provide proof-of-principle tests of using different approaches to determine the SOLPS-ITER input parameters that result in the best simulation-experimental agreement, in particular, the conjugate gradient method applied on the validation metric.
The paper is organized as follows: the TCV-X21 experiment reference case and the associated dataset are introduced and extended in Sec.II.The SOLPS-ITER simulations are described in Sec.III.Then in Sec.IV, we present the qualitative and quantitative validation results for simulations with the standard up-stream matching approach and the systematic methodology to determine the input parameters based on the quantitative validation result.In Sec.V, we analyze and discuss the results of the validation and the effect of neutrals in the TCV-X21 scenario.Finally, the conclusions are presented in Sec.VI.

II. TCV-X21 EXPERIMENTAL DATASET AND ITS EXTENSION
In this work we use the experimental dataset of the TCV-X21 reference case which is publicly available at https: //github.com/SPCData/TCV-X21.This scenario is a lower single null L-mode Ohmic plasma in the TCV tokamak 14 with a toroidal field of B φ ,axis ≃ 0.95T, a plasma current of I p ≃ 165kA and an electron line-average density ⟨n e ⟩ ∼ 2.5 × 10 19 m −3 measured by the far infrared camera (FIR, Fig. 1(a), vertical cyan line).TCV-X21 includes data in both toroidal field directions to study the effect of drifts, with the convention that "forward" (Forw) denotes the field direction where the ion ∇B drift (B × ∇B) points downwards, from the core towards the X-point and "reversed" (Rev) when it points upwards 13 .
Tab.I lists the diagnostics, observables and their respective hierarchies (used in the validation metric) of the TCV-X21 dataset.Fig. 1(a) shows the position of the listed diagnostics.The TCV-X21 dataset includes mean and fluctuation profiles of observables covering the divertor targets and volume, the divertor entrance, and the outboard midplane.In this work, we only consider the mean profiles because SOLPS-ITER does not predict fluctuation quantities.Detailed information about the observables and diagnostics can be found in Ref. 13.In this work we add new observables to the dataset, namely the divertor neutral pressure measurements and deuterium Balmer lines, which will be introduced in the following subsections.
Reciprocating divertor probe array (RDPA) for divertor volume  The baratron pressure gauges (BAR) considered here provide measurements of neutral pressure at the TCV floor (p div ) and in the turbo pump duct (p tmp -see Fig. 1(a)).The gauges are installed at the end of extension tubes to protect them from the magnetic field of the tokamak.Therefore, we need a model to relate the measurement to the in-vessel pressure 16 .The energetic atomic and molecular divertor neutrals flowing into the tube undergo thermalization through collisions with the walls and the atoms recombine to form molecules.At the end of the tube, after several bends, the pressure is determined solely by the molecular density at wall temperature, T wall = 300K.For the experiment-simulation comparision, we use the 0D model discussed in Ref. 19 and Ref. 20  to determine p div and p tmp from SOLPS-ITER outputs.To compensate for the pressure drop in p tmp induced by the turbo pump, the experimental data of p tmp is multiplied by a factor of 1.5 21 .The measured pressure is averaged over 0.8s of the flattop phase of several, repeat discharges of the TCV-X21 scenario.
The main source of uncertainty of the p div and p tmp measurements is ∆e rep , the uncertainty related to reproducibility 13 , estimated from repeat discharges.

B. Divertor Spectroscopy System
The Divertor Spectroscopy System (DSS) installed in TCV (Fig. 1(a)) provides measurements of the line-integrated visible radiation at different wavelengths, corresponding to different atomic processes 22 .The DSS system has 30 chords along which it is possible to measure the wavelength spectra with a spectral resolution of up to 0.02 nm 15 .The measurement is a line integral of the emission along a given chord.For each discharge, the time traces are averaged over 0.22s, providing D n→2 as a function of the chord number.The final average emission profile is obtained by averaging the profiles from different shots.In this work, we consider three deuterium Balmer lines, D 5→2 , D 6→2 , and D 7→2 .
The main sources of uncertainty of the line brightnesses are ∆e rep , which is estimated comparing different shots, and ∆e dia , the uncertainty due to inherent characteristics of the di-agnostics, which is estimated as 15% of the measured intensity.

III. SOLPS-ITER SIMULATIONS
SOLPS-ITER (Scrape-Off Layer Plasma Simulation-ITER) is composed of the transport code B2.5 that solves the Braginskii multi-fluid equations, and the kinetic neutral Monte Carlo solver EIRENE 23 .EIRENE is coupled self consistently with B2.5 to calculate the sources and sinks due to plasma-neutral interactions.The simulations in this work consider a multicomponent plasma, including carbon impurities and kinetic neutrals and their main reactions in the plasma 24 .Drifts are not included in these simulations, since convergence in low density, high temperature plasmas could not be achieved so far for SOLPS-ITER simulations of TCV plasmas.Unlike in the previous turbulence validation, ad-hoc diffusion coefficients for cross-field heat and particle transport are used.The simulations presented here are a good testbed for the role of the neutrals in the TCV-X21 scenario, in particular, enabling the investigation of the distribution of the ionization profile across the SOL.The absence of drifts in these simulations is a limitation, but may help disentangle the origin of the flows in the divertor, i.e., whether drift or transport-driven 25 .
Fig. 1(b) shows the computational grid used in this work.The radial particle diffusion and heat conduction is described using spatially uniform anomalous diffusion coefficients.For the simulations presented in Sec.IV A, we choose a particle diffusivity of D ⊥ = 0.2 m 2 /s and a thermal diffusivity of χ e,⊥ = 1.0 m 2 /s and χ i,⊥ = 1.0 m 2 /s, which were found in previous works to result in a good upstream match for TCV L-mode plasmas 7 .A deuterium gas puff rate of Γ D 2 = 6.8 × 10 19 /s is used for a close upstream density match, determined after a scan of the gas puff Γ D 2 = {4.5, 5.6, 6.8, 7.2, 7.6, 8.0, 8.4} × 10 19 /s.The chemical sputtering coefficient of carbon impurities on the wall is assumed as Y chem = 3.5%, and the particle recycling coefficient is set to be R = 0.99.
At the divertor targets, sheath Bohm boundary conditions are applied.At the core boundary, power transferred from the core to the edge is set to be comparable to the experimental value (160 kW).Neutrals crossing the core boundary are returned as fully ionized particles.More details about the boundary conditions used can be found in Ref. 24.
In Sec.IV B, an iterative method using the overallagreement metric and the conjugate gradient is used to adapt the input parameters Γ D 2 , D ⊥ and χ e,⊥ to achieve a better overall agreement.In Sec.IV C, D ⊥ and χ e,⊥ are scanned to find a better overall agreement.

IV. VALIDATION RESULTS
The validation of the SOLPS-ITER simulations in this work follows the same methodology as used in Ref. 13.The details of the mathematical model can be found in Ref. 17, and the basic concept of this methodology is summarized as follows: The level of agreement between simulation and experiment is evaluated using a large set of observables and is quantified by the overall agreement metric χ. χ = 0 and χ = 1, means, respectively, perfect agreement and complete disagreement.The fundamental quantity used to calculate χ is the normalised simulation-experimental distance d j , which, for each observable j is defined as: where e j,i , ∆e j,i , s j,i , and ∆s j,i denote, respectively, the experimental values, their uncertainties, the simulation values, and their uncertainties, defined on a series of discrete data points i ∈ {1, 2, ..., N j }. d j → 0 means a perfect agreement between simulation and experiment for the observable j.Due to the difficulty to provide a rigorous estimate of the simulation uncertainty, we set ∆s j,i = 0 as in Ref. 13.Other important quantities used in the validation are: the sensitivity S j of an observable j, which takes values between 0 and 1, approaching 1 for very small relative uncertainties (high precision) of the experimental and simulation observables; the hierarchy weighting H j , as given in Tab.I, is a value associated with each observable j that is smaller the higher the number of model assumptions and/or measurement combinations needed to determine the observable.Based on d j and the weighting factors S j and H j , a metric χ is then evaluated to indicate the overall simulation-experiment agreement over all (or a subset of) the observables, see Appendix A for more information.In addition, the quality Q is evaluated, which denotes the quality of the validation.This value is higher when a higher number of more directly computed, higher precision observables are included in the validation, see Appendix A.
A. Standard approach: matching upstream profiles The validation results for the standard approach of matching upstream profiles are given in Tab.II, in the columns labeled Standard (Forw) and Standard (Rev).Since the present simulations do not include drifts, the results are compared to experimental data in both field directions.Some selected profiles showing simulation-experiment comparisons are given in Fig. 2, where the radial coordinate r − r sep denotes the distance from the separatrix, after mapping along the magnetic flux surface to the outboard midplane.
The comparison at the outboard midplane and the divertor entrance are shown for the electron density n e and the electron temperature T e in Fig. 2(a) to 2(d).An appreciable match is found as expected from the values in Tab.II, and since the simulations are tuned for a good upstream match.As in the simulations, the J sat is estimated from n e and T e , its good outboard midplane agreement in Tab.II is also expected.The good agreement observed for the upstream V pl in Tab.II is a consequence mostly of large experimental error bars, as can be seen from the lower value of S j compared to that of n e and T e .On the other hand, the floating potential V f l and Mach number M ∥ at the outboard midplane show considerable discrepancy between simulation and experiment as indicated by the d j values in Tab.II.The value of M ∥ is approximately zero, different from the experiment and simulations in Ref. 13  (not shown).This is attributed to the absence of drifts, which are the basic mechanism of the Pfirsch-Schlüter flows at the outboard midplane 25,26 .
In the divertor volume (Fig. 2(e) to 2(h)), the simulation roughly reproduces the radial shape of the experimental n e and T e profiles.However, the simulated n e profile peaks at Fig. 3 shows the comparison of the Balmer line intensity profiles measured by DSS, where the DSS chords, see Fig. 1 a), are labeled with increasing number from the bottom up to the X-point.The measured intensity in the forward field direction is systematically larger than the reversed field direction.For all three Balmer lines, the simulation successfully reproduces the profile shape, with the two peaks of the intensity, corresponding to the LFS target and the X-point/high-fieldside (HFS) target.In the simulation, the location of the second peak (with the larger chord number) matches well with the experiment, while the first peak is shifted.Generally, the simulation underestimates the intensity, especially in the region between the two peaks, i.e., along the leg.At the first peak, the underestimation is small for D 5→2 and for D 6→2 , while large for D 7→2 .At the second peak, the value in the simulation is much smaller than that in the forward field measurement, while closer to that measured in the reversed field experiments.
Fig. 4 presents the experimental and simulation results for the neutral pressure in the two locations in the divertor.Both p div and p tmp in the forward field direction are higher than in the reversed field direction.The simulation result for p div in the standard approach is within the range of uncertainty of the forward field case, while the simulated p tmp matches within experimental uncertainty of the reversed field case.Compared to the full-field TCV SOLPS simulations studied in Ref. 7, where the simulated p div systematically exceeded the measured value by a factor ∼ 400%, here in the reversed field, p div is overestimated only by ∼ 50% at most.In summary, the validation in this case gives an overall agreement metric χ = 0.770 for the forward field direction, and χ = 0.807 for the reversed field direction.Good agreement is found for the n e , T e , J sat , and V pl profiles at the outboard midplane and the n e and T e profile in the divertor entrance, for both field directions.V f l in the midplane and p div in the forward field direction, and p tmp in the reversed field direction also show good quantitative agreement.Poorer matches are found in M ∥ at the outboard midplane, M ∥ and V f l in the divertor volume, q ∥ at the low field side target, and V pl at the high field side target, in both field directions, and in V pl and V f l at the low field side target in the reversed field direction.
To gain insight on role of the neutrals in the TCV-X21 scenario, we plot in Fig. 5 the simulated total ionization sources (more precisely, what is shown is the source of electrons due to ionization with the generation of both D + and D 2 + ).We can observe that in the SOL, most of the ionization happens along the separatrix, especially at the two targets.According to the SOLPS simulation, the ionization in the SOL (regions 2 ⃝, 3 ⃝ and 4 ⃝ in Fig. 1(b)) accounts for ∼ 65.0% of the total ionization.The HFS divertor region ( 3 ⃝ in Fig. 1(b)) accounts for ∼ 14.9% and the LFS divertor region ( 2 ⃝ in Fig. 1(b)) accounts for ∼ 31.2%.Another difference between these SOLPS-ITER simulations and the turbulence simulations in 13 is the inclusion of impurity species.Here the SOLPS simulation includes carbon impurities.Their radiation is relatively weak in this scenario, ∼ 15% of the total input power.

B. Conjugate gradient method
We explore here a systematic method to determine SOLPS-ITER input parameter values that lead to an overall improvement of the agreement metric χ.This is done by minimizing the multi-variable function χ (we recall that χ = 0 indicates perfect experiment-simulation agreement: considering its dependence on the gas puff rate Γ D 2 , the particle diffusivity D ⊥ , and the electron thermal diffusivity χ e,⊥ .As a proof of principle, only these three input parameters are used in this test, but one could, in principle, use all input parameters of the simulation subjected to tuning.For this task, we apply the conjugate gradient method, an algorithm used to solve unconstrained optimization problems 28,29 .The main advantage of this method, compared to the gradient descent method, is the avoidance of oscillating behaviors when calculating the gradient directions in the iterative minimization 29 .The algorithm to determine the input parameters can be briefly described as follows: Step 1: The index r indicates the iteration step number, with r ∈ {0, 1, 2, ...}.In the first iteration, r = 0, we compute the gradient g 0 = ∇χ(x 0 ) at the starting point x 0 = (Γ D 2 ,0 , D ⊥,0 , χ e,⊥,0 ).This is done using finite differences between χ(x 0 ) and three neighbouring simulations.Then, the minimization direction is set as s 0 = −g 0 .
Step 2: For the current r, perform simulations along the direction s r using the parameters determined as x ′ r+1 = x r + λ s r , with the values of λ chosen appropriately.
Step 3: Evaluate χ for all simulations in Step 2 and determine the local minimum, x r+1 .TABLE II.Quantitative validation result for SOLPS-ITER simulations.For each observable, the normalized distance d j and the sensitivity weighting S j are listed.The composite agreement metric χ and the overall quality Q are given, for each diagnostic and for the overall validation."Standard", "Gradient", and "D ⊥ χ e,⊥ scan" refer to the simulation with the standard upstream matching approach (Sec.IV A), the best agreement point in the first 1-D search of the gradient method (Sec.IV B), and the best agreement point in the grid scan of D ⊥ and χ e (Sec.IV C), respectively.In the standard approach, the simulation is validated against the two field directions in the TCV-X21 dataset 13  Step 4: If an asymptotic convergence is not achieved, the scan continues using the new gradient g r+1 = ∇χ(x r+1 ) , and the new scan direction (the conjugate gradient direction) is set by Step 5: Once x r+1 and s r+1 are determined in Step 3 and Step 4, restart from Step 2 with r = r + 1 for the next iteration.
In this way, x r is expected to converge to the minimum possible of x, i.e., maximized agreement.Fig. 6 shows the results of performing two iterations of the conjugate gradient method.The agreement metric is computed with respect to the reversed field case.In this demonstration, the parameters were first normalized to their values at the starting point, such that all quantities have a similar order of magnitude, and then the conjugate gradient method was applied.Finally, the parameters were denormalized back to real values.In each step, in order to find an approximate minimum along s r , we performed simulations for three positive values of λ and treated the one with minimum χ as the minimum point.The gradient needed to determine s was calculated numerically, using the metric differences between the point of interest and three neighboring points.
The starting point was chosen from the simulation in Sec.IV A, with Γ D 2 = 6.8 × 10 19 m/s, D ⊥ = 0.2 m 2 /s, and χ e,⊥ = 1.0 m 2 /s.In the first iteration (Fig. 6, red line), the metric decreases first, then increases, revealing a minimum.We select this lowest point as the starting point of the second iteration (Fig. 6, magenta line), where we observe an alternating increase and decrease, without exhibiting a clear asymptotic convergence or improvement of the agreement metric.Overall, the best case (lowest χ) is obtained in the first step, where Γ D 2 = 5.7 × 10 19 m/s, D ⊥ = 0.22 m 2 /s, and χ e,⊥ = 1.06 m 2 /s.The procedure is stopped after these two iterations, as no improvement was achieved in the second step and due to the involved simulation costs (12 simulations were needed in total for these two minimization steps).
In Tab.II, the d j and χ values for the simulations with the smallest overall metric χ found with the gradient method are presented in column "Gradient(Rev)".In general, the two step demonstration gives an overall improvement of the agreement by 5.5%.The majority (21 out of 32) of the observables have been improved as indicated by the decrease of the corresponding d j .Among them, 11 observables for the gradient method show a decrease of the d j value larger than 10%.Several observables have been significantly improved (from disagreement to agreement), for example, the n e and T e profiles measured at the LFS target by the LPs, T e of RDPA, and p div .

C. D ⊥ and χ e,⊥ scans
As an alternative of the conjugate gradient method, we also conduct here a scan of D ⊥ and χ e,⊥ , on the grids spanned by D ⊥ = {0.15,0.2, 0.25} m 2 /s and χ e,⊥ = {1.0,1.5, 2.0} m 2 /s, which are plotted as blue triangles in Fig. 6.For these nine simulations, the deuterium gas puff is fixed to be Γ D 2 = 6.8 × 10 19 /s.In Fig. 6, we find that the simulation with the largest D ⊥ and χ e,⊥ value, D ⊥ = 0.25 m 2 /s and χ e,⊥ = 2.0 m 2 /s gives the best agreement.
In Tab.II, the d j and χ values for the simulation with the smallest overall metric χ in the D ⊥ and χ e,⊥ scan are given in column "D ⊥ χ e,⊥ scan (Rev)".The overall increase of agreement given by χ is ∼ 10.4%.23 out of 32 observables have been improved as indicated by the decrease of their d j .Among them, 14 observables show a decrease of their d j value by more than 10%.Several observables have been significantly improved (from disagreement to agreement), for example, the T e profile from LPs at the LFS and HFS targets, and T e from RDPA.

V. DISCUSSION
In Sec.IV A, a SOLPS-ITER simulation without drifts and with constant particle and energy diffusivity D ⊥ and χ e/i,⊥ , was validated against the TCV-X21 reference case.The standard TCV L-mode values for D ⊥ , χ e,⊥ and χ i,⊥ were used 7 and the gas puff rate was manually tuned to have a good upstream density match.The simulation shows good agreement with experimental data at the outboard midplane and the divertor entrance, with the exception of the parallel Mach number.The agreement is especially good for the electron density n e and temperature T e .When approaching the divertor targets, a less good agreement is found, but the simulation captures the shape and order of magnitude of the n e and T e profiles.The quantitative match with the floating potential V f l , plasma potential V pl , parallel Mach number M ∥ , and parallel heat flux q ∥ , are not satisfactory.The large experimental-simulation distance for q ∥ can be mainly attributed to the small experimental uncertainties and a shift of the profile peak positions (Not shown).The mismatch of V f l and V pl with experimental data is mainly attributed to the omission of drifts in the simulations.
In Sec.IV B and Sec.IV C, using two different methods, we obtained improved quantitative agreement compared with the standard approach in Sec.IV A. The best agreement case in the conjugate gradient method of Sec.IV B features a decrease of gas puff, and slight increases in the two transport coefficients.From Fig. 7(a), we can observe that the significant improvement in the d j of the outer target n e is mainly due to the decrease of the peak value, as a result of a reduced gas puff.The improved match of p div in the divertor region for the reversed field case is also related to the decrease of the gas puff, resulting in a decrease of neutral pressure (Fig. 4).The best agreement case in the D ⊥ and χ e,⊥ scan features a significant increase of the transport coefficients, which leads to an increase of the fall-off length of n e and T e at both targets and improves the match with the experiment profile.Similarly, T e in the divertor volume displays a shallower radial decrease, which agrees better with the RDPA measurements.This may indicate that the experimental case has a higher perpendicular transport compared to what is assumed in the simulation with the standard approach.FIG. 6. Validation metric of simulations using the gradient method and parameter grid scan.The χ of each simulation is plotted as a function of the gas puff Γ D 2 , the particle diffusivity D ⊥ , and the electron thermal diffusivity χ e,⊥ .The gas puff scan of section IV A is shown by the green circles, while the D ⊥ and χ e,⊥ scan is plotted as blue triangles.The simulations of the first and second iteration of the conjugate gradient method are plotted in red and magenta squares, respectively.The arrows denote the directions of searching.The conjugate gradient method does not give the overall best agreement (lowest metric χ) although it does improve the result marginally.This can be attributed to several possible factors.First, the conjugate gradient method does not guarantee the increase of the agreement level in every iteration step.As shown in an example in Fig. 8, which aims to reproduce a similar behavior as observed in Sec.IV B, two iteration steps are not enough to reach the minimum point and non-monotonic behaviors can be found in the second iteration.To go towards the minimum point, more iteration steps may be needed.Second, the finite difference method used for the gradient calculation might introduce non-negligible errors.In the calculation of the gradient in iteration 2, we assumed the gradient along the first direction of the 1-D search to be zero, which is also an approximation.Third, the step length of each 1-D search is limited by the numerical cost of each simulation.Therefore, the 1-D minimization can only be estimated approximately, with an error of the order of the step length.Possible solutions to these problems could be using other methods for the numerical differentiation, for example the central difference method, to get a better estimate of the gradient; or trying other minimization methods independent of the evaluation of the gradients, for example, the multi-dimensional simplex method 29 .Compared to the improvement it brought, the conjugate gradient method tested here is found to be a computationally too expensive approach.
The ionization profile given by the SOLPS-ITER simulation (Fig. 5) is clearly different to what was assumed in the turbulence simulations in Ref. 13, where the ionization sources were assumed to be localized in the outer region of the confined plasma.This motivates further studies to add more realistic neutral particle profiles, or self-consistent inclusion of neutrals, to the turbulence simulations.
To explore the effect of neutrals in the flows in the divertor region, we plotted the parallel Mach number profile along a magnetic flux surface at r − r sep = 0.5 cm in Fig. 9.The GBS turbulence simulation without neutrals from Ref. 13  (positive Mach number) and reveal a significant flow all along the divertor leg, being somewhat weaker in SOLPS-ITER.Instead, the Mach number measured with RDPA is much small, close to zero.Both GBS and SOLPS simulations show an increasing Mach number as approaching the divertor target, while the RDPA measurements feature a flat profile.Comparing the GBS simulations in reversed and forward toroidal field directions, we find that the effect of drifts is small compared with the difference between simulation and RDPA measurement.Ref. 13 suggested the difference between the GBS simulations and RDPA measurements to be due to the ionization source along the divertor leg and primarily to potentially be located just in front of the target.In this study using the SOLPS simulation with neutrals, the Mach number is lower than the GBS, but the flows in the simulations are still significant along the entire divertor leg and considerably higher than in the experiment.This raises questions in the Mach number measurements with RDPA in these conditions or the model used for its interpretation.Further investigation is needed in order to disentangle the differences between the flow velocities presented here.

VI. CONCLUSION
SOLPS-ITER transport fluid simulations without drifts and with uniform particle and energy cross-field transport coefficients were qualitatively compared and quantitatively vali- dated against the TCV-X21 experimental reference case from Ref. 13.Three Balmer lines measured across the outer divertor leg by DSS and two divertor neutral pressure measurements from the Baratron gauges were added to the publicly available TCV-X21 dataset.In the standard approach, where SOLPS-ITER input parameters are tuned to match upstream quantities, qualitative comparisons of profiles from upstream to the divertor were carried out.As expected, in this standard approach, agreement between simulation and experiment is found in the outer midplane and at the divertor entrance, with the exception of parallel Mach flows.Reduced agreement is found in the divertor volume and at the divertor targets, for quantities such as density, temperature, and plasma potential profiles.Despite TCV-X21 being a (near) sheathlimited plasma designed to minimize the effect of neutrals in the SOL/divertors, the simulation still finds a ∼ 65.0% of the ionization happening in the SOL, with ∼ 14.9% in the HFS divertor region and ∼ 31.2% in the LFS divertor region.The simulation also shows ∼ 15% of the input power to be radiated by carbon impurities.Using the quantitative validation metric χ, in a proof-of-principle test, we use the conjugate gradient method and D ⊥ and χ e,⊥ scans to improve the agreement level, resulting in a 5.5% and 10.4% improvement respectively compared to the results achieved using the standard approach.This suggests that a reduced gas puff and increased particle and energy transport coefficients compared to what is used when exclusively trying to match upstream profiles lead to a better match with the experimental case, mainly via decreasing the peak target density and broadening the density and temperature profiles.While the performance achieved here with the conjugate gradient method is rather modest, this method may be improved by using finer steps and bet-ter numerical differentiation methods.Other algorithms, such as multidimensional simplex method may, however, be better a solution for the problems presented here in the iterative method to determine the input parameters resulting in an optimal match with the experiment.
The SOLPS simulations with neutrals show a significant portion of neutral ionization to occur in the SOL, a major difference compared with the assumption used in the first turbulence code validation in the TCV-X21 validation case in Ref. 13.The parallel flows in the divertor observed in SOLPS and GBS turbulence simulations from Ref. 13 are similar in shape, with the GBS divertor flows systematically larger in comparison to the SOLPS flows.This suggests some flow reduction in the divertor by the neutrals.The parallel Mach numbers from SOLPS-ITER are, however, still substantially larger than those measured with RDPA, raising questions on the latter that will be further explored in the future.The results in this work provide useful information for future turbulence simulations of TCV with neutrals, while suggesting that the contribution of the neutrals to the flow velocity and, therefore, to the parallel heat flow towards the targets, should be further investigated.

FIG. 1 .
FIG. 1.(a) Poloidal view of the magnetic reconstruction (dark blue lines) of the TCV-X21 reference case from LIUQE reconstruction 18 , also shown are: the Langmuir probes (LP, blue circles), the Reciprocating Divertor Probe Array (RDPA, black L-shaped structure) and its covered area (transparent black rectangle), the Thomson Scattering system (TS, red square array), the Fast Horizontal Reciprocating Probe (FHRP, purple rectangle) and its covered area (transparent purple), the Far InfraRed interferometer (FIR, vertical cyan line), the area of sight of the Vertical Infrared Camera (IR, green transparent patched area), the position of the top gas fueling valve (magenta rectangle), the ports of the baratron pressure gauges (dark green boxes on the wall), and the DSS viewing chords (radial orange lines in the divertor).(b) Computational grid used for the SOLPS ITER plasma model.Regions marked by different colors are 1 ⃝: the core; 2 ⃝: the low-field-side (LFS) divertor region; 3 ⃝: the high-field-side (HFS) divertor region; and 4 ⃝: other region of the SOL.

FIG. 2 .
FIG. 2. Comparison of the profiles in the standard approach of matching upstream profiles.The 1-D experimental profiles in forward field (blue lines) and reversed field (purple lines) and in the SOLPS-ITER simulation (orange lines) for electron density (left column) and electron temperature (right column).The measurements are from the FHRP (see Fig. 1(a)) at the outboard midplane ((a) and (b)), from TS at the divertor entrance ((c) and (d)), and from LPs at the low-field-side divertor target ((i) and (j)), also shown are (e) the 2-D profile of electron density, and f) temperature of the simulation and the corresponding reversed field experiments (g) and (h) measured by RDPA in the divertor volume.(k) and (l) show the 2D Mach number measured in the experiment (reversed field) and from the SOLPS simulations.

FIG. 5 .
FIG. 5. Simulated total ionization source.In (a) the values at the two targets go up to ∼ 10 23 m −3 s −1 , which is saturated in the color code, and are thus shown in the two sub figures (b) and (c).

FIG. 7 .
FIG. 7. Comparison of profiles from the standard simulation in Sec.IV A (starting point), the best case from gradient method, and the best case from the D ⊥ and χ e,⊥ scan.We show the 1-D experimental profiles of reversed field (purple lines) and the SOLPS-ITER simulation profiles for (a) electron density and (b) electron temperature measured by the LPs at the low-field-side divertor target, and (c) electron density and (d) electron temperature measured by the LPs at the high-field-side divertor target.We also show (e) the 2-D profile of temperature measured by RDPA, and (f)-(h) the corresponding simulation results.

FIG. 8 .
FIG. 8.A simple example that could give a similar behavior for the conjugate gradient method as in our study.For illustration, we consider a 2-D minimization problem with the filled contour plot.The first and second 1-D search are denoted by the red dash line and magenta dash line.The red arrow represents the direction of the first gradient (opposite to the first 1-D search direction), while the magenta arrow represents the second gradient.

FIG. 9 .
FIG.9.The Mach number along the separatrix to the LFS target.These 1-D plots consider the parallel Mach number along a flux surface (at r − r sep = 0.5 cm) in the SOL near the LFS separatrix.The RDPA measurements are plotted by solid lines with shaded errorbars.The SOLPS simulations from Sec. IV A with three gas puff values are given together with GBS simulations in Ref.13.
VII. ACKNOWLEDGEMENTThis work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 -EUROfusion).The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI).Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI.Neither the European Union nor the European Commission nor SERI can be held responsible for them.This work was supported in part by the Swiss National Science Foundation.

TABLE I .
Observables and comparison hierarchies for validation.Observables constituting the TCV-X21 dataset but not used in the present validation are shown in gray.The hierarchy weighting, H j , is a constant value depending on the number of assumptions and/or models used to obtain the observable, H j = [h sim + h exp − 1] −1 , see Sec.IV and Appendix A. The data from the Divertor Spectroscopy System (DSS) and the Baratron gauges (BAR) are newly added into the TCV-X21 dataset. .