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

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 [1], a key element of the ITER design [2, 3]. In a diverted 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 [3] and SPARC [4], and even more so of a DEMO [5] 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 [6, 7]. In this approach, the largest possible number of experimental observables are quantitatively 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 [8], SOLEDGE2D [9], UEDGE [10] which solve the mean field steady state fluid equations without turbulence, but with inclusion of neutral dynamics, often in a kinetic description. The effect of fluctuation induced transport is represented in these codes by effective diffusion coefficients, set by the users. Alternatively, the modelling is often conducted with edge turbulence codes, which self-consistently evolve turbulence dynamics of the plasma but have limited ability to include neutral dynamics and are usually computationally more expensive. An example for the latter approach is the recent multi-code validation involving GBS [11], GRILLIX [12], and TOKAM3X [13], which presented the first full tokamak size edge turbulence simulations [14] of a diverted single-null plasma in the Tokamak à Configuration Variable (TCV) [15]. The scenario simulated by the turbulence boundary codes in this study, referred to as TCV-X21 reference case, is a L-mode 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 [14].

In this work, we validate the Scrape-Off Layer Plasma Simulation-ITER (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 available [14] TCV-X21 dataset (table 1) with five new observables for the neutral dynamics: Balmer line intensities measured by the divertor spectroscopy system (DSS) [16] and neutral pressure measurements from the Baratron gauges [17]. We conduct a similar quantitative validation as in [14], using the methodology of [18] (briefly reviewed in section 4). 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.

Table 1. 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_\text{sim}+h_\text{exp}-1]^{-1}$, see section 4 and appendix. The data from the divertor spectroscopy system (DSS) and the baratron gauges (BAR) are newly added into the TCV-X21 dataset.

Diagnostic	Observable	Hierarchy
		hExp	hSim	Hj
Wall Langmuir Probes (LP) at the low-field-side (LFS) and high-field-side targets	ne , Te , Vpl	2	1	1 / 2
	Jsat, $\sigma (J_{\text{sat}})$	1	2	1 / 2
	$\mathrm{skew}(J_{\text{sat}})$, $\mathrm{kurt}(J_{\text{sat}})$	1	2	1 / 2
	$J_{\parallel}$, $\sigma (J_{\parallel})$	1	1	1
	Vfl , $\sigma (V_{fl})$	1	2	1 / 2
Infrared camera (IR) for LFS target	$q_\parallel$	2	2	1 / 3
Reciprocating divertor probe array (RDPA) for divertor volume	ne , Te , Vpl	2	1	1 / 2
	$M_{\parallel}$	2	2	1 / 3
	Jsat, $\color{gray}\sigma (J_{\text{sat}})$	1	2	1 / 2
	$\mathrm{skew}(J_{\text{sat}})$, $\mathrm{kurt}(J_{\text{sat}})$	1	2	1 / 2
	Vfl , $\sigma (V_{fl})$	1	2	1 / 2
Thomson scattering (TS) for divertor entrance	ne , Te	2	1	1 / 2
Fast horizontally-reciprocating probe (FHRP) for outboard midplane	ne , Te , Vpl	2	1	1 / 2
	$M_{\parallel}$	2	2	1 / 3
	Jsat, $\sigma (J_{\text{sat}})$	1	2	1 / 2
	$\mathrm{skew}(J_{\text{sat}})$, $\mathrm{kurt}(J_{\text{sat}})$	1	2	1 / 2
	Vfl , $\sigma (V_{fl})$	1	2	1 / 2
Divertor spectroscopy system (DSS)	$D_{5\rightarrow2}$, $D_{6\rightarrow2}$, $D_{7\rightarrow2}$	2	3	1 / 4
Baratron gauges (BAR)	$p_\textrm{tmp}$	2	2	1 / 3
	$p_\textrm{div}$	1	2	1 / 2

The paper is organized as follows: the TCV-X21 experiment reference case and the associated dataset are introduced and extended in section 2. The SOLPS-ITER simulations are described in section 3. Then in section 4, 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 section 5, 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 section 6.

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. We also provide a static repository for the version used in this paper at https://doi.org/10.5281/zenodo.1084117. This scenario is a lower single null L-mode Ohmic plasma in the TCV tokamak [15] with a toroidal field of $B_{\phi, \textrm{axis}}\simeq 0.95\mathrm{T}$, a plasma current of $I_\textrm{p}\simeq 165\,\mathrm{kA}$ and an electron line-average density $\langle n_e \rangle \sim 2.5\times 10^{19} \mathrm{m^{-3}}$ measured by the far infrared camera (FIR, figure 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 $\nabla B$ drift ($\mathbf{B} \times \nabla B$) points downwards, from the core towards the X-point and 'reversed' (Rev) when it points upwards [14].

Zoom In Zoom Out Reset image size

Table 1 lists the diagnostics, observables and their respective hierarchies (used in the validation metric) of the TCV-X21 dataset. Figure 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 [14]. 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.

2.1. Baratron pressure gauges

2.2. Divertor spectroscopy system

SOLPS-ITER is composed of the transport code B2.5 that solves the Braginskii multi-fluid equations, and the kinetic neutral Monte Carlo solver EIRENE [25]. Information on atomic and molecular processes considered and reaction rate models used in EIRENE can be found in [26]. 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 multi-component plasma, including carbon impurities and kinetic neutrals and their main reactions in the plasma [27]. 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.

Figure 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, a standard approach adopted for SOLPS-ITER modelling of TCV. Without a justified model, varying the transport coefficients spatially might introduce too many additional degrees of freedom and make the simulation results complicated to interpret. In terms of the optimization, this may introduce some overfitting. For the simulations presented in section 4.1, we choose a particle diffusivity of $D_{\perp} = 0.2 \ \mathrm{m^2\,s^{-1}}$ and a thermal diffusivity of $\chi_{e,\perp} = 1.0 \ \mathrm{m^2\,s^{-1}}$ and $\chi_{i,\perp} = 1.0\ \mathrm{m^2\,s^{-1}}$, which were found in previous works to result in a good upstream match for TCV L-mode plasmas [8]. A deuterium gas puff rate of $\Gamma_\mathrm{D_2} = 6.8\times 10^{19}\,\mathrm{s^{-1}}$ is used for a close upstream density match, determined after a scan of the gas puff $\Gamma_\mathrm{D_2} = \{4.5, 5.6, 6.8, 7.2, 7.6, 8.0, 8.4\}\times 10^{19}\,\mathrm{s^{-1}}$. The chemical sputtering coefficient of carbon impurities on the wall is assumed as $Y_{\text{chem}} = 3.5\%$, and the particle recycling coefficient is set to be R = 0.99. An alternative choice would be to choose R as a variable input. In [27], it was shown that, within a reasonable range, vary $1-R$ has a very similar effect on the plasma as vary $\Gamma_\mathrm{D_2}$.

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 \ \mathrm{kW}$) and is set to be equally distributed between electrons and ions. Neutrals crossing the core boundary are returned as fully ionized particles. More details about the boundary conditions used can be found in [27].

In section 4.2, an iterative method using the overall-agreement metric and the conjugate gradient is used to adapt the input parameters $\Gamma_\mathrm{D_2}$, $D_{\perp}$ and $\chi_{e,\perp}$ to achieve a better overall agreement. In section 4.3, $D_\perp$ and $\chi_\mathrm{e, \perp}$ are scanned to find a better overall agreement.

The validation of the SOLPS-ITER simulations in this work follows the same methodology as used in [14]. The details of the mathematical model can be found in [18], 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:

$$\begin{equation} d_j = \left[ \frac{1}{N_j}\sum_{i = 1}^{N_j} \frac{\left(s_{j,i} - e_{j,i}\right)^2} {\Delta s_{j,i}^2 + \Delta e_{j,i}^2} \right]^{1/2} \end{equation} \tag{ 1 }$$

where $e_{j,i}$, $\Delta e_{j,i}$, $s_{j,i}$, and $\Delta 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\in\{1, 2, \ldots , N_j \}$. $d_j\rightarrow0$ 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 $\Delta s_{j,i} = 0$ as in [14]. 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 table 1, 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 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.

4.1. Standard approach: matching upstream profiles

The validation results for the standard approach of matching upstream profiles are given in table 2, 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 figure 2, where the radial coordinate $r-r_{\text{sep}}$ denotes the distance from the separatrix, after mapping along the magnetic flux surface to the outboard midplane.

Zoom In Zoom Out Reset image size

Table 2. 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_\perp$ $\chi_\mathrm{e, \perp}$ scan' refer to the simulation with the standard upstream matching approach (section 4.1), the best agreement point in the first 1D search of the gradient method (section 4.2), and the best agreement point in the grid scan of $D_\perp$ and χ_e (section 4.3), respectively. In the standard approach, the simulation is validated against the two field directions in the TCV-X21 dataset [14].

		Standard (Forw)		Standard (Rev)		Gradient (Rev)		$D_\perp$ $\chi_\mathrm{e, \perp}$ scan(Rev)
Diagnostics	Observables	dj	Sj	dj	Sj	dj	Sj	dj	Sj
Fast horizontally reciprocating probe (FHRP) for outboard midplane	ne	0.800	0.835	1.275	0.889	0.809	0.880	0.683	0.884
	Te	0.359	0.733	0.907	0.791	1.151	0.799	0.567	0.792
	Jsat	1.523	0.905	1.845	0.894	1.780	0.890	1.559	0.891
	Vpl	0.797	0.695	1.386	0.734	1.248	0.743	1.179	0.745
	Vfl	1.749	0.703	4.217	0.835	3.851	0.846	4.091	0.830
	$M_{\parallel}$	7.295	0.847	6.879	0.876	6.881	0.876	6.946	0.875
	(χ; Q)	(0.310; 2.218)		(0.501; 2.363)		(0.459; 2.371)		(0.388; 2.362)
Reciprocating divertor probe array (RDPA) for divertor volume	ne	2.559	0.871	2.335	0.888	2.640	0.881	2.467	0.888
	Te	2.041	0.899	1.929	0.901	1.642	0.906	1.196	0.909
	Jsat	4.191	0.909	3.722	0.917	3.506	0.913	4.164	0.917
	Vpl	2.588	0.869	3.748	0.884	3.367	0.890	3.087	0.892
	Vfl	10.122	0.745	23.186	0.816	21.370	0.828	22.232	0.818
	$M_{\parallel}$	11.944	0.898	13.329	0.901	12.870	0.899	12.374	0.895
	(χ; Q)	(0.979; 2.446)		(0.966; 2.503)		(0.915; 2.509)		(0.829; 2.510)
Thomson scattering (TS) for divertor entrance	ne	0.848	0.876	1.148	0.898	1.393	0.891	1.153	0.894
	Te	0.818	0.893	1.095	0.901	1.342	0.905	0.908	0.900
	(χ; Q)	(0.004; 0.884)		(0.045; 0.899)		(0.190; 0.898)		(0.031; 0.897)
Infrared camera (IR) for low-field-side target	$q_\parallel$	6.350	0.911	6.322	0.941	6.437	0.943	6.588	0.943
	(χ; Q)	(1.000; 0.304)		(1.000; 0.314)		(1.000; 0.314)		(1.000; 0.314)
Wall Langmuir probes for low-field-side target	ne	3.222	0.879	2.046	0.880	1.667	0.866	1.953	0.881
	Te	2.170	0.891	1.907	0.866	1.599	0.876	1.293	0.880
	Jsat	3.688	0.909	2.342	0.908	2.066	0.902	2.347	0.909
	Vpl	3.441	0.888	4.982	0.887	4.483	0.893	4.031	0.897
	Vfl	2.184	0.639	4.356	0.732	4.237	0.741	4.278	0.733
	$J_{\parallel}$	3.451	0.746	3.015	0.768	3.028	0.767	2.903	0.755
	(χ; Q)	(0.985; 2.849)		(0.954; 2.904)		(0.842; 2.906)		(0.841; 2.906)
Wall Langmuir probes for high-field-side target	ne	3.360	0.886	4.269	0.901	3.046	0.888	3.602	0.901
	Te	3.136	0.922	2.750	0.892	2.248	0.899	1.531	0.902
	Jsat	2.079	0.878	3.702	0.900	3.054	0.892	2.992	0.900
	Vpl	5.234	0.915	4.756	0.884	4.490	0.890	3.972	0.895
	Vfl	2.907	0.673	3.583	0.687	3.537	0.696	3.730	0.687
	$J_{\parallel}$	3.617	0.788	3.648	0.781	3.643	0.780	3.846	0.777
	(χ; Q)	(0.987; 2.925)		(0.999; 2.913)		(0.994; 2.912)		(0.904; 2.920)
Divertor spectroscopy system (DSS) for divertor volume	$D_{5\rightarrow2}$	3.460	0.869	3.014	0.870	3.240	0.862	2.830	0.876
	$D_{6\rightarrow2}$	2.218	0.753	2.627	0.825	2.774	0.818	2.481	0.835
	$D_{7\rightarrow2}$	3.467	0.818	2.896	0.826	3.053	0.821	2.641	0.835
	(χ; Q)	(0.986; 0.610)		(0.997; 0.630)		(0.998; 0.625)		(0.993; 0.637)
Baratron pressure gauges (BAR) for divertor volume	pdiv	0.469	0.877	3.336	0.941	1.193	0.933	2.871	0.940
	ptmp	3.685	0.940	0.845	0.936	2.099	0.931	1.512	0.933
	(χ; Q)	(0.417; 0.752)		(0.603; 0.783)		(0.411; 0.777)		(0.742; 0.781)
Overall	(χ; Q)	(0.770; 12.988)		(0.807; 13.309)		(0.763; 13.312)		(0.723; 13.326)

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 figures 2(a)–(d). An appreciable match is found as expected from the values in table 2, and since the simulations are tuned for a good upstream match. Though most of the density and temperature profile is within the uncertainty range, we also note that near the separatrix, there is a slightly narrower simulated profile in figures 2(b) and (d), compared to the experimental data. As in the simulations, the J_sat is estimated from the product of n_e and sound speed (which depends on T_e and T_i ), its good outboard midplane agreement in table 2 is also expected. The good agreement observed for the upstream V_pl in table 2 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_fl and Mach number $M_{\parallel}$ at the outboard midplane show considerable discrepancy between simulation and experiment as indicated by the d_j values in table 2. The value of $M_{\parallel}$ is approximately zero, different from the experiment and simulations in [14] (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 [28, 29], and an additional effect attributed to the poloidal distribution of the cross-field transport in the transport codes [30].

In the divertor volume (figures 2(e)–(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 the target, while in the experiments it peaks at the X-point. For T_e , both simulation and experimental peaks are close to the X-point, with the simulation showing a stronger poloidal gradient. The experimental n_e and T_e profiles show overall similar trends in this region in both field directions, although the density profile is shifted more towards the private flux region in forward field. For the Mach number in the divertor volume (figures 2(k) and (l)), the SOLPS-ITER simulation predicts high values in the region $r-r_\mathrm{sep}\gt 0$, where RDPA measured Mach numbers close to zero.

At the LFS target (figures 2(i) and (j)), the simulation reproduces the peak magnitude of T_e , but gives a larger peak n_e value. The overestimation of n_e was also observed in the reversed field direction in previous SOLPS-ITER simulations of TCV L-mode plasma at $1.45\,\mathrm{T}$ [8]. The simulation profile of n_e and T_e in figures 2(i) and (j) are narrower compared with the experiment in both field directions. It is worth noting that in the forward field direction, the n_e experimental profile shows a double peak structure not present in the simulations. Such profile shape is usually attributed to the effect of drifts [31], which are not included in this simulation.

Figure 3 shows the comparison of the Balmer line intensity profiles measured by DSS, where the DSS chords, see figure 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/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\rightarrow 2}$ and for ${D}_{6\rightarrow 2}$, while large for ${D}_{7\rightarrow 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.

Zoom In Zoom Out Reset image size

Figure 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-ITER simulations studied in [8], where the simulated p_div systematically exceeded the measured value by a factor ${\sim}400\%$, here in the reversed field, p_div is overestimated only by ${\sim}50\%$ at most.

Zoom In Zoom Out Reset image size

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_fl in the midplane and $p_\mathrm{div}$ in the forward field direction, and $p_\mathrm{tmp}$ in the reversed field direction also show good quantitative agreement. Poorer matches are found in $M_\parallel$ at the outboard midplane, $M_\parallel$ and V_fl in the divertor volume, $q_\parallel$ 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_fl 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 figure 5 the simulated total ionization sources (more precisely, what is shown is the source of electrons due to ionization with the generation of both $\mathrm{D}^+$ and $\mathrm{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-ITER simulation, the ionization in the SOL (regions ②, ③ and ④ in figure 1(b)) accounts for ${\sim}65.0\%$ of the total ionization. The HFS divertor region (③ in figure 1(b)) accounts for ${\sim}14.9\%$ and the LFS divertor region (② in figure 1(b)) accounts for ${\sim}31.2\%$. Another difference between these SOLPS-ITER simulations and the turbulence simulations in [14] is the inclusion of impurity species. Here the SOLPS-ITER simulation includes carbon impurities. Their radiation is relatively weak in this scenario, ${\sim}15\%$ of the total input power. The effect of impurities is expected to be weak in these simulations.

Zoom In Zoom Out Reset image size

4.2. 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:

$$\begin{equation} \chi\left(\mathbf{x}\right), \quad \mathbf{x} = \left(\Gamma_\mathrm{D_2}, D_\perp, \chi_{e,\perp}\right), \end{equation} \tag{ 2 }$$

considering its dependence on the gas puff rate $\Gamma_{\mathrm{D_2}}$, the particle diffusivity $D_\perp$, and the electron thermal diffusivity $\chi_{e, \perp}$. 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 [32, 33]. 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 [33].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\in\{0,1,2,\ldots \}$. In the first iteration, r = 0, we compute the gradient $\mathbf{g}_0 = \nabla \chi(\mathbf{x_0})$ at the starting point $\mathbf{x}_0 = (\Gamma_\mathrm{D_2,0}, D_{\perp,0}, \chi_{e,\perp,0})$. This is done using finite differences between $\chi(\mathbf{x_0})$ and three neighbouring simulations. Then, the minimization direction is set as $\mathbf{s}_0 = -\mathbf{g}_0$.

Step 2: For the current r, perform simulations along the direction s_r using the parameters determined as $\mathbf{x}_{r+1}^{^{\prime}} = \mathbf{x}_r + \lambda \mathbf{s}_r$, with the values of λ chosen appropriately to cover the range of minimization and resolve typical differences of the data.

Step 3: Evaluate χ for all simulations in Step 2 and determine the local minimum, $\mathbf{x_{r+1}}$.

Step 4: If an asymptotic convergence is not achieved, the scan continues using the new gradient $\mathbf{g}_{r+1} = \nabla \chi(\mathbf{x_{r+1}})$, and the new scan direction (the conjugate gradient direction) is set by

$$\begin{equation} \mathbf{s}_{r+1} = -\mathbf{g}_{r+1} + \frac{g_{r+1}^2}{g_{r}^2}\mathbf{s}_r. \end{equation} \tag{ 3 }$$

Step 5: Once $\mathbf{x}_{r+1}$ and $\mathbf{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.

Figure 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 λ such that $\lambda \mathbf{s}_r = \{0.1,0.2,0.3\} \mathbf{s}_r/|\mathbf{s}_r|$, 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.

Zoom In Zoom Out Reset image size

The starting point was chosen from the simulation in section 4.1, with $\Gamma_\mathrm{D_2} = 6.8\times 10^{19}\,\mathrm{m\,s^{-1}}$, $D_{\perp} = 0.2 \ \mathrm{m^2\,s^{-1}}$, and $\chi_{e,\perp} = 1.0 \ \mathrm{m^2\,s^{-1}}$. In the first iteration (figure 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 (figure 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 $\Gamma_\mathrm{D_2} = 5.7\times 10^{19}\,\mathrm{m\,s^{-1}}$, $D_\perp = 0.22\ \mathrm{m^2\,s^{-1}}$, and $\chi_{e,\perp} = 1.06\ \mathrm{m^2\,s^{-1}}$. 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 table 2, 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_\mathrm{div}$.

4.3.  $D_\perp$ and $\chi_{e, \perp}$ scans

In section 4.1, a SOLPS-ITER simulation without drifts and with constant particle and energy diffusivity $D_\perp$ and $\chi_{e/i,\perp}$, was validated against the TCV-X21 reference case. The standard TCV L-mode values for $D_\perp$, $\chi_{e,\perp}$ and $\chi_{i,\perp}$ were used [8] 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_fl , plasma potential V_pl , parallel Mach number $M_\parallel$, and parallel heat flux $q_\parallel$, are not satisfactory. The large experimental-simulation distance for $q_\parallel$ in reversed field can be mainly attributed to the small experimental uncertainties and a shift of the profile peak positions as shown in figure 7. In the forward field case, the mismatch also comes from an overestimation of the magnitude. The mismatch of V_fl and V_pl with experimental data is mainly attributed to the omission of drifts in the simulations.

Zoom In Zoom Out Reset image size

In sections 4.2 and 4.3, using two different methods, we obtained improved quantitative agreement compared with the standard approach in section 4.1. The best agreement case in the conjugate gradient method of section 4.2 features a decrease of gas puff, and slight increases in the two transport coefficients. From figure 8(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_\mathrm{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 (figure 4). The best agreement case in the $D_\perp$ and $\chi_{e, \perp}$ 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.

Zoom In Zoom Out Reset image size

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 figure 9, which aims to reproduce a similar behavior as observed in section 4.2, 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 1D search to be zero, which is also an approximation. Third, the step length of each 1D search is limited by the numerical cost of each simulation. Therefore, the 1D 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, or the adjoint algorithmic differentiation sensitivity calculation [34] 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 [33]. Compared to the improvement it brought, the conjugate gradient method tested here is found to be a computationally too expensive approach.

Zoom In Zoom Out Reset image size

5.1. Discussion on the ionization profile and parallel flow

The ionization profile given by the SOLPS-ITER simulation (figure 5) is clearly different to what was assumed in the turbulence simulations in [14], 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_{\mathrm{sep}} = 0.5\,\mathrm{cm}$ in figure 10. The GBS turbulence simulation without neutrals from [14] and the SOLPS-ITER simulation results are of the same order of magnitude. The simulated parallel flows point towards the target (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-ITER 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. Oliveira, Body et al [14] 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-ITER simulation with neutrals, the Mach numbers in the divertor are lower than the results given by GBS, a turbulence code with drifts included, but the divertor flows in the present simulations are still considerably higher than in the experiment. This raises questions on 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. Possible controlled simulations, with different neutral density distributions could be interesting future directions to reveal the net contribution of the neutrals.

Zoom In Zoom Out Reset image size

SOLPS-ITER transport fluid simulations without drifts and with uniform particle and energy cross-field transport coefficients were qualitatively compared and quantitatively validated against the TCV-X21 experimental reference case from [14]. 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. For a number of observables, the difference between experiments and simulations is larger than that between forward and reversed field experimental data, which suggests that there are important ingredients other than drifts needed to aim for a closer quantitative match. An example is the use of spatially uniform transport coefficients in the present study. Despite TCV-X21 being a (near) sheath-limited plasma designed to minimize the effect of neutrals in the SOL/divertors, the simulation finds that ${\sim}65.0\%$ of the ionization occur in the SOL, with ${\sim}14.9\%$ in the HFS divertor region and ${\sim}31.2\%$ in the LFS divertor region. The simulation also shows ${\sim}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_\perp$ and $\chi_{e, \perp}$ 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 better numerical differentiation methods. Other algorithms, such as multidimensional simplex methods may, however, be a better 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-ITER simulations with neutrals show a significant portion of neutral ionization to occur in the SOL. This is a major difference compared with the assumption used in the first turbulence code validation in the TCV-X21 validation case in [14], motivating the self-consistent inclusion of neutrals in future TCV-X21 turbulence studies. The parallel flows in the divertor observed in SOLPS-ITER and in GBS turbulence simulations from [14] are similar in shape, with the GBS divertor flows systematically larger in comparison to the SOLPS-ITER 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.

This 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.

This appendix summarizes the quantities used in the validation procedure, following the same method as in [14].

The sensitivity S_j , for the observable j, is given by

$$\begin{equation} S_{j} = \exp \left(-\frac{\sum_{i} \Delta e_{j, i}+\sum_{i} \Delta s_{j, i}}{\sum_{i}\left|e_{j, i}\right|+\sum_{i}\left|s_{j, i}\right|}\right), \end{equation} \tag{ A1 }$$

where $e_{j,i}$, $\Delta e_{j,i}$, $s_{j,i}$, and $\Delta 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\in\{1, 2, \ldots , N_j \}$.

The hierarchy weighting H_j is defined as

$$\begin{equation} H_j = \left(h_\text{sim}+h_\text{exp}-1\right)^{-1}, \end{equation} \tag{ A2 }$$

where h_sim and h_exp are the simulation and experimental primacy hierarchy level for an observable, being higher the higher number of assumptions and/or measurement combinations used to obtain the observable. We note that the hierarchy weighting does not include the dimension or number of channels of an observable. The idea is that each observable should be determined in the best possible way, in terms of resolution etc. However, its dimensionality does not directly enter in the metric.

The overall agreement metric χ and the quality Q of a set of observables are obtained by:

$$\begin{equation} \chi = \frac{\sum_{j} R_j\left(d_{j}\right) H_{j} S_{j}}{\sum_{j} H_{j} S_{j}} ,\qquad Q = \sum_j H_j S_j, \end{equation} \tag{ A3 }$$

where the level-of-agreement function $R_j(d_j)$ is a nonlinear, increasing function of the normalised simulation-experimental distance d_j (equation (1)), with sharp transition when d_j is larger than certain threshold, defined as

$$\begin{equation} R_j\left(d_{j}\right) = \frac{\tanh \left[\left(d_{j}-1 / d_{j}-d_{0}\right) / \lambda\right]+1}{2}. \end{equation} \tag{ A4 }$$

In this work we set d₀ = 1, λ = 0.5, as in [18]. R_j takes values between 0 and 1. It is used to unify the distance to a level of agreement with fixed range, from perfect agreement (0) to complete disagreement within errorbars (1).