Studying the process of divertor detachment and the associated complex interplay of plasma dynamics and atomic physics processes is of utmost importance for future fusion reactors. While simplified analytical models exist to interpret the general features of detachment, they are limited in their predictive power, and complex twodimensional (2D) or even threedimensional (3D) codes are generally required to provide a selfconsistent picture of the divertor. As an intermediate step, onedimensional (1D) models of the scrapeoff layer (SOL) can be particularly insightful as the dynamics are greatly simplified, while still selfconsistently including various source and sink terms at play, as well as additional important effects such as flows. These codes can be used to shed light on the physics at play, to perform fast parameter scans, or to interpret experiments. In this paper, we introduce the SPLEND1D (Simulator of PLasma ENabling Detachment in 1D) code: a fast and versatile 1D SOL model. We present in detail the model that is implemented in SPLEND1D. We then employ the code to explore various elements of detachment physics for parameters typical of the Tokamak à Configuration Variable, including the atomic physics and other processes behind power and momentum losses, and explore the various hypotheses and free parameters of the model.
I. INTRODUCTION
Plasma power and particle exhaust is a crucial issue for future fusion reactors. If unmitigated, the target heat fluxes in ITER and DEMO are expected to greatly exceed the $10\u2009MW/m2$ limit that is considered necessary for tolerable steadystate conditions. It will thus be necessary to operate in an, at least partially, detached regime.^{1–3} In such a regime, the total plasma pressure develops strong parallel gradients along the field lines in the scrapeoff layer (SOL), driven by volumetric and perpendicular momentum losses,^{4,5} and a significant fraction of the plasma power is dissipated by volumetric power sinks. This results in a reduction of the target heat and particle flux densities, as well as target temperature, which is important to reduce target erosion.^{6,7} The detachment process typically sets in at a target temperature below 5 eV,^{5,8,9} where plasmaneutral interactions are further enhanced, further reducing target heat flux, temperature, and ion flux.
While simplified analytical models exist to interpret the general features of detachment, they are limited in their predictive power, and complex 2D or even 3D codes are generally required to provide a selfconsistent picture of the divertor. However, the results of these codes can be difficult to interpret. Furthermore, the numerical cost of running such simulations generally prohibits large parameter scans that could be useful to shed light on the role of the various terms entering the equations. Therefore, as an intermediate step, onedimensional (1D) models of the scrapeoff layer (SOL) can be particularly insightful as the dynamics are greatly simplified, while still selfconsistently including various source and sink terms at play, as well as additional important effects such as flows. Typically based on a projection of the Braginskii's equations^{10} along a field line, these models simplify the perpendicular dynamics either by ignoring it, or by simplifying it with some simple ansatz. Such models have been typically employed in the context of a specific study.^{11–13} More recently, efforts have focused on implementing dedicated 1D codes that can be shared, reused and coupled with other codes or integrated into workflows for integrated modeling. Examples of such codes are SOLF1D,^{14} SD1D,^{15} and DIV1D.^{16} In this, paper, we introduce the Simulator of PLasma ENabling Detachment in 1D (SPLEND1D) code, a fast and versatile 1D model used to explore the complex interplay of atomic physics and plasma dynamics underlying the detachment process. The purpose of this paper is twofold. First, to present in great detail the model implemented in the SPLEND1D code and demonstrate the speed and robustness of the code. Second, to apply the code for the study of a base case scenario for parameters typical for the Tokamak à Configuration Variable (TCV).^{17} The neutral particle source is increased in this scenario to reach the onset of detachment, similarly to the experimental onset of detachment through increased plasma fueling. The different processes are subsequently investigated, thus providing a first illustration of the SPLEND1D capabilities and possible applications. This paper is organized as follows. In Sec. II, the derivation of the models for the charged species (plasma) and the neutral species is presented, as well as the numerical methods, boundary conditions, and implemented source terms. Section III presents a reference base case that is used to demonstrate the capabilities of the SPLEND1D code in terms of ease of use and result interpretation. We also present some measurements of the code accuracy and runtime. In Sec. IV, we simulate detachment of the base case through upstream density ramps, and highlight the mechanisms at play in the model that enable the onset of momentum and energy losses. Section V investigates the role of some free parameters of the model, such as the impurity concentration, the neutral confinement time, the heatflux limiters, and the choice of boundary condition for the parallel velocity. Finally, in Sec. VI, we present some advanced studies that have been enabled by SPLEND1D, investigating the choice of the neutral model, the effect of separating ion and electron energy equations, and the results of timedependent simulations. These advanced studies demonstrate SPLEND1D's aptness in interpreting current TCV experiments, for example to study the role of connection length on divertor detachment.^{18} Comparisons to other 1D codes in the community are discussed throughout the text.
II. MODEL
In this first section, we describe the model implemented in SPLEND1D. The equations solved by SPLEND1D are based on the Braginskii equations,^{10} which are typically used to describe the evolution of the plasma in the SOL in 1D, 2D or 3D codes, using a fluid approximation. In particular, we present in detail the assumptions used to derive the SPLEND1D code.
A. Geometry
B. Plasma model
1. Equations
2. Heat fluxes
3. Viscosity tensor
C. Fluid model for the neutrals
This section presents the model retained to describe the dynamics of the neutrals. As the neutrals are not affected by the magnetic field, their dynamics are intrinsically 3D. Furthermore, in typical divertor conditions, the mean free path of the neutrals can be large compared to the system size, leading to a high Knudsen number that in principle requires a kinetic description of the neutral dynamics. This is the approach of EIRENE, one of the main workhorse for neutral dynamics simulations in the fusion community,^{26} where a Monte Carlo method is used to simulate the behavior of the neutrals. However, such methods are computationally expensive and subject to statistical noise. In recent years, to alleviate these costs, there has been significant work devoted to the development of advanced fluid neutral models, or hybrid neutral models, mixing a kinetic and a fluid description, Refs. 27–29 and references therein. As 1D models aspire for simplicity, in this work, we use fluid neutral models. SPLEND1D implements two different neutral models, the choice of which to use is made by the user: a diffusive neutral model, in which the neutrals diffuse along the θdirection; and an advective one, where neutrals move in the $\theta \u2212\phi $ plane (that is, a flux surface), with a velocity that is not necessarily parallel to the magnetic field, similar to the model presented in Refs. 27 and 30. In the following, we further assume that the vectors describing the neutral dynamics (velocities, heat fluxes) do not have a component in the radial direction ( $\Psi $). Therefore, while neutrals are not confined to a particular flux tube, they remain confined to a particular fluxsurface. This assumption, despite being clearly a simplification, provides the highest fidelity level available for the neutral dynamics description, while maintaining a 1D fluid picture and not requiring further information about the geometry of the problem or the plasma quantities in the divertor. In its current version, SPLEND1D models only a single population of neutral atoms, and does not include contributions from molecules. It is, however, well known from experimental and simulation works that molecules are very strong contributors to the divertor dynamics,^{31–33} in particular through elastic collisions with the main ion species. Future works on the SPLEND1D neutral model will therefore encompass the addition of neutralion collisions, as well as the addition of molecular effects.
1. Continuity equation for the neutrals
2. Diffusive neutral model
D. Atomic source terms
Now, we introduce the source terms resulting from atomic interactions between the neutral and plasma populations. We restrict ourselves here to the atomic processes of ionization, charge exchange, excitation, and recombination since we model a singleion population and consider only atomic neutrals. We remark here that ionneutral elastic collisions are not considered, as we assume that the neutral and ion populations interact only through these atomic processes (ionization, recombination, chargeexchange). Inclusion of ionneutral elastic collisions is left for future work.
1. Rates
SPLEND1D's rate coefficients are either obtained from a bilinear interpolation of the openADAS^{34} database, or evaluated using the fit provided in the AMJUEL manual,^{35} as decided by the user. Table I presents the sources used for the ionization, excitation, recombination and charge exchange rates as functions of n_{e} and T_{e}. The chargeexchange reaction is discussed further below. In the case of AMJUEL data, the excitation rate Q_{exc} includes both the power radiated by excitation and the potential energy cost E_{ion} in the case of ionization ( $Eion=13.6\u2009eV$ in the case of hydrogen).^{35} In the case of openADAS data, the latter is added through an additional sink term proportional to E_{ion} and the ionization rate.
Reaction .  Rate coefficient .  AMJUEL^{35} .  OpenADAS^{34} . 

Ionization  $\u27e8\sigma v\u27e9ion$  Reaction 2.1.5, section 4  SCD 
Excitation  Q_{exc}  Reaction 2.1.5, section 10.2  PLT 
Recombination  $\u27e8\sigma v\u27e9rec$  2.1.8, section 4  ACD 
Recombination cooling rate  Q_{rec}  2.1.8, section 10.4  PRB 
Charge exchange  $\u27e8\sigma v\u27e9CX$  2.1.9, section 3.1.8  CCD 
Reaction .  Rate coefficient .  AMJUEL^{35} .  OpenADAS^{34} . 

Ionization  $\u27e8\sigma v\u27e9ion$  Reaction 2.1.5, section 4  SCD 
Excitation  Q_{exc}  Reaction 2.1.5, section 10.2  PLT 
Recombination  $\u27e8\sigma v\u27e9rec$  2.1.8, section 4  ACD 
Recombination cooling rate  Q_{rec}  2.1.8, section 10.4  PRB 
Charge exchange  $\u27e8\sigma v\u27e9CX$  2.1.9, section 3.1.8  CCD 
Since these rates are estimated from either tabulated data (openADAS) or from fitted expressions, they are only defined over a certain range of validity. If, during a SPLEND1D simulation, the temperature or density were to go outside of these ranges, the rates are computed using the minimum (or maximum) values for which they are defined.
2. Sources and sinks
E. Impurity radiation
We note here that such impurity model, limited to constant impurity fraction, does not account for the impact of impurity transport and the subsequent occurrence of noncoronal effects. These could be implemented through an ad hoc (free) parameter, τ, a residence time that would allow to form the noncoronal parameter $ne\tau $ as done in Ref. 40, together with tabulated values of $Lz(Te,ne\tau )$ produced with the collisional radiative model (CRM),^{39} or with the inclusion within SPLEND1D itself of a collisional radiative model for impurities, evolving in time the concentration of impurities' charge states, as done for instance in Ref. 41.
F. Boundary conditions
To solve the equations describing the evolution of the plasma parameters [Eqs. (2)–(5) or (10)–(12)] and the evolution of the fluid neutral parameters [Eqs. (28), (C3), (C4), and (C10)], a set of boundary conditions is required. We distinguish two types of boundaries, depending on whether the considered boundary is a symmetry plane or a “target.” Two kinds of configurations can be simulated with SPLEND1D (Fig. 2):

The fluxtube is assumed to be symmetric around s = 0. Symmetric boundary conditions are applied on the left boundary, while “target” boundary conditions are applied to the right boundary [Fig. 2(a)].

The flux tube connects two targets together, and no symmetry is assumed. In this case, “target” boundary conditions are applied to the left and right boundaries [Fig. 2(b)].
The second case can be used, for instance, to study the power sharing between the two sides of the flux tube, which depends on the ratio of the connection lengths.^{42}
G. Numerical implementation
The code is mostly intended to be used as an initial value problem (IVP) solver. Starting from an arbitrary initial solution, the code evolves the system of equations in time, and converges toward a steadystate, if such a solution exists. Convergence is determined by the user based on the time evolution of various macroscopic quantities, as well as by the norm of the equations' residuals. The code can also directly search for a steadystate by setting the timederivatives to 0, through a nonlinear Newton solver. This is typically run only after a temporal evolution of the equations, to provide the solver a reasonable first guess for the solution.
Several timestepping schemes are implemented. SPLEND1D can either rely on the timestepping algorithms implemented in the TS environment of the Portable, Extensible Toolkit for Scientific Computation (PETSc) library,^{47,48} using by default the fully implicit CrankNicolson method. In schemes that requires the computation of the Jacobian of the system, such as the CrankNicolson method, this operation is performed by PETSc using finitedifferences with coloring. The Jacobian is typically recomputed every 10–40 iterations, a number that is set at runtime by the user. The time step is free to evolve between userprescribed minimum and maximum values, with PETSc taking care of the time step adaptivity. Alternatively to the use of PETSc, and not demonstrated in this paper, SPLEND1D is equipped with a “linearized” IMEX (IMplicitEXplicit) scheme, where the advection and source terms are treated explicitly, and the diffusion and viscosity operators are handled implicitly, using the nonlinear transport coefficients of the previous time step. The implicit part can then be rewritten as a succession of tridiagonal matrix inversions, performed using the Thomas algorithm, which scales linearly with the number of cells. While this scheme comes with a CourantFriedrichsLewy (CFL) condition that restricts the time step, its linear scaling with the resolution and number of equations makes it an interesting solver for cases with a large number of cells (provided the cells are not so small that the CFL condition becomes very strict). Due to the relatively small size of the typical problem solved by SPLEND1D (around ∼500–1000 cells for a maximum of 8 equations, that is, ∼4000–8000 degrees of freedom), the code is currently not parallelized, although extension to an Message Passing Interface (MPI)OpenMP parallelized code would be relatively straightforward, leveraging the capabilities of PETSc for the MPI parallelization. In Sec. III, after presenting a base case used to demonstrate various features of the code, we will briefly present the performance of SPLEND1D in terms of convergence and speed.
III. PRESENTATION OF BASE CASE
A. Parameters
This section presents the base case investigated in this paper to illustrate in more detail the capabilities of the code, to highlight the role of various modeling parameters and to demonstrate the code's capabilities to unravel the physics behind plasma detachment. The simplified model given by Eqs. (10)–(12) is used, assuming T_{i} = T_{e} as opposed to the twofluid model. The diffusive neutral model is employed [Eqs. (33) and (32)], and we impose T_{n} = T_{i}. Symmetric boundary conditions are applied, and the velocity is imposed greater or equal than the sound speed at the right boundary. Table II summarizes the values of the various parameters in these simulations. R, γ_{i}, γ_{e}, and τ_{n} are based on reasonable values, whereas $L$, c_{c}, and α are chosen based on typical TCV^{17} values. 500 grid cells are simulated, with the grid accumulation chosen so that the width of the cell in front of the target is approximately 0.8 mm in the parallel direction. These choices will be further discussed in Sec. III B.
Parameter .  Definition .  Value . 

$ln\u2009\Lambda $  Coulomb logarithm  Equation (A4) 
$\tau \xaf$  $Ti=\tau \xafTe$  1 
γ_{i}  Ion sheath transmission coefficient  3.5 
γ_{e}  Electron sheath transmission coefficient  5.5 
α_{i}  Ion heat flux limiter  0.6 
α_{e}  Electron heat flux limiter  0.6 
c_{c}  Carbon concentration  2% 
R  Recycling rate  99% 
τ_{n}  Neutral confinement time  $0.05\u2009ms$ 
B  Magnetic field  Constant 
α  Fieldline angle  $4\xb0$ 
$L$  Parallel connection length  $25\u2009m$ 
ν_{num}  Numerical viscosity  0 
Parameter .  Definition .  Value . 

$ln\u2009\Lambda $  Coulomb logarithm  Equation (A4) 
$\tau \xaf$  $Ti=\tau \xafTe$  1 
γ_{i}  Ion sheath transmission coefficient  3.5 
γ_{e}  Electron sheath transmission coefficient  5.5 
α_{i}  Ion heat flux limiter  0.6 
α_{e}  Electron heat flux limiter  0.6 
c_{c}  Carbon concentration  2% 
R  Recycling rate  99% 
τ_{n}  Neutral confinement time  $0.05\u2009ms$ 
B  Magnetic field  Constant 
α  Fieldline angle  $4\xb0$ 
$L$  Parallel connection length  $25\u2009m$ 
ν_{num}  Numerical viscosity  0 
This results in an upstream density of $1.5\xd71019$ m^{−3} and upstream temperature of $\u224840$ eV. The full temperature and density profiles are shown in Fig. 3 in green, alongside the velocity and Mach number profiles. A target temperature of 11.6 eV suggests that the base case represents an attached plasma regime. In this scenario, increasing the upstream plasma density via an increase in the particle source H_{n} leads to a rollover in the target ion flux $\Gamma ,t$, Fig. 4, indicating the onset of detachment. Example profiles for detached and strongly detached cases are shown in orange and red in Fig. 3.
We note that in Fig. 4, each point is the steadystate result of simulations with different values of H_{n}. They were obtained by simulating $1\u2009s$ of plasma dynamics, and then applying the steadystate solver discussed in Sec. II G. This will be the case in all simulation results presented in this paper, except the timedependent simulation presented in Sec. VI B.
B. SPLEND1D performance in the base case
We next explore the performance, in term of computing time, of SPLEND1D, across a density ramp with the same input parameters as that of Fig. 4. For this, we evaluate the simulation time, toward steadystate, for different values of $H\xafn$, applying three different strategies:

Each simulation is run independently, starting from flat, unphysical profiles ( $n=2.0\xd71019\u2009m\u22123,\u2009Te=Ti=20\u2009eV,u=0\u2009ms\u22121$).

Each simulation is run independently, starting from “physical” profiles obtained from a steadystate simulation with $H\xafn=1.58\xd71021\u2009m\u22123s\u22121$.

The simulations are run sequentially, each starting from the steadystate obtained from the previous simulations.
In all simulations presented in this section, the same convergence parameters for the PETSc solver are used. The time step $\Delta t$ is allowed to vary between $\Delta t=5.71\xd710\u221211\u2009s$ and $\Delta t=5.71\xd710\u22124\u2009s$. It is automatically adapted by PETSc based on convergence and error estimates. The simulations are run for $1\u2009\u2009s$ of simulated plasma time, after which the steadystate solver is applied to find an exact (within nonlinear solver tolerance) steadystate solution. All simulations were performed on a typical laptop, with an Intel^{®} Core™i78565U CPU. These three simulation strategies provide identical output profiles (within the nonlinear solver tolerance), Fig. 6(a). Figure 6(b) reports the simulation runtimes. Starting from flat, unphysical profiles (strategy 1), all simulations converge within $\u2248$ 20 s, except for a few outliers that require up to 40 s. Starting from physical profiles of an attached case (strategy 2) yields even lower computation times, approximately tens, except at low densities (lower than the initial simulation), where simulations can take relatively long times or even fail to converge in less than 10 min (after which the simulations were stopped). A similar observation is done for the sequential scan of $H\xafn$ (strategy 3), restarting from the previous simulation (with a slightly different value of $H\xafn$), where the simulation time can drop to $\u2248$ 6 s, except at low density where simulations can struggle to converge in reasonable time. In conclusion, this section demonstrates the numerical performances of SPLEND1D, which can achieve convergence toward a steadystate in less than 30 s across a wide range of regimes, including the strongly detached one.
We start by considering the influence of the number of grid cells, N, keeping $hend=0.8\u2009mm$. N is varied from 100 to 1000 by increments of 50, and then from 1000 to 10 000 by increments of 500, Fig. 7. All simulations converged. We first use a simulation with 31 250 cells of 0.8 mm each as a reference for the evaluation of p_{ref}. As N is increased, we observe that both the upstream and downstream values of δ_{p} decrease quadratically with N, as one could expect from a second order scheme. ϵ_{tot} remains largely unaffected. Furthermore, using now a high resolution (125 000 cells, $0.2\u2009mm$ each) simulation as reference, we find that δ_{p} stagnates when the number of cells is higher than 2000. This is because, as will be shown in the next paragraph, the accuracy of SPLEND1D is largely dictated by the width of the cells close to the target, kept fixed in this scan. Hence, increasing further the number of cells does not lead to an increase in the accuracy. Furthermore, and interestingly, even at fairly low resolution ( $\u223c100$ cells), δ_{p} for the upstream and target pressure remain lower than 1%.
We now keep the number of grid cells constant (500 cells) while changing h_{end}. As h_{end} is reduced, ϵ_{tot} and δ_{p} (both upstream and downstream) strongly decrease, Fig. 8. This indicates that the accuracy of SPLEND1D largely depends on the width of cells near the target. In particular, ϵ_{tot} scales approximately as $hend2$.
Taken together, the results of Figs. 7 and 8 show that the accuracy of SPLEND1D is largely dictated by the grid resolution near the targets, with the total number of grid cells playing a lesser role. However, as the number of grid cells is increased, or the width of the last cell is decreased, the computational cost increases. From Figs. 6–8, the resolution chosen for the base case (500 grid cells, 0.8 mm target cell width) appears to be a good tradeoff between accuracy and numerical cost, and is retained for the rest of the simulations presented in this paper.
C. The twopoint model formatting
To evaluate the importance of physical momentum and power loss processes in the base case simulations, we first introduce briefly the twopoint model (TPM) formulation. It is a model that does not consider the spatial distribution of plasma parameters along the SOL, but considers quantities at only two points: upstream and target, where upstream can be any point along the flux tube. Furthermore, it considers a steady state situation.
IV. DETACHMENT ONSET VIA DENSITY RAMP IN THE BASE CASE
A. Observation of a target ion flux rollover and onset of a total pressure drop
As mentioned earlier, in the base case presented in Sec. III, increasing the upstream plasma density via an increase in neutral particle source H_{n} leads to a rollover in the target ion flux, Fig. 4, indicating the onset of detachment. Before the rollover, the target ion flux increases somewhat linearly with increasing n_{u}, whereas the twopoint model would predict $\Gamma ,t\u221dnu3$. Such behavior is reminiscent of the “power starvation” mechanisms, observed in previous simulation works^{15} and experiments,^{5} in which the linear evolution of the ion flux and its subsequent rollover are driven by a reduction of the power that is available for ionization. This is clearly seen in the simulations presented in this paper, where the rollover coincides with a drop of the ionization source, whereas recombination only plays a role in the reduction of the ion flux at a later stage, Fig. 9. The rollover is accompanied by a reduction of target temperature to less than $1\u2009eV$ [Fig. 10(a)]. In contrast, the upstream temperature T_{u} is much less sensitive, and only at the highest degree of detachment does it start to degrade significantly, dropping from approximately $\u224840$ eV to $\u224828$ eV. The particle flux rollover is accompanied by a total pressure drop, as a signature of detachment, as shown by evaluating f_{mom} [Eq. (57)], which increases from 0 (no pressure drop) at low density to 0.99 at the highest density achieved in the present simulations, Fig. 10(b). The plasma power loss factor, f_{pow} [Eq. (60)], exhibits a similar behavior as the momentum loss with increasing density, but reaches saturation ( $fpow\u223c1$) at a lower upstream density. This suggests that the onset of detachment is first driven by the increase in power losses, which precedes momentum losses. With the reduction of the total pressure, Fig. 10(b), a reduction of the static pressure at the target is also observed, Fig. 10(c), initially driven by the reduction of the target temperature, and, ultimately, also by a reduction of the target density.
We remark in Fig. 10(b) that, at low density, both f_{mom} and f_{pow} can be negative. For f_{mom}, as shown in Figs. 11 and 12, this is due to the action of the viscosity term, which can be relatively large in the initial, attached conditions when the gradient of velocity in front of the target is large. As for f_{pow}, this is related to the ionization process, Fig. 12, and the particular choice of neutral particle source, H_{n}, which injects in the system neutrals with a temperature T_{n} = T_{i}. Upon ionization, this thermal energy is transferred to the plasma. We note here that these observations are consequences of some of the choices made in the derivation of the model.
B. Investigate process at play
The first sign of detachment, Fig. 10, appears to be a target static pressure rollover, along with a target ion flux rollover and a target temperature decrease to below $1$ eV. With a further increase in upstream density, the target electron density rolls over, followed by a saturation in the power losses and then momentum losses. The underlying processes behind these features can be studied with SPLEND1D, which directly outputs each term contributing to the particle, momentum and energy balance equations [Eqs. (10)–(12)], such that the importance of each process can be compared. Figure 11 shows the profile of each of these terms along the flux tube for an attached, a detached, and a strongly detached case. This allows to compute the contribution of each of these terms to the momentum and power losses, Fig. 12. We must note that the choice of neutral model may affect the role of atomic processes in momentum and power losses. In particular, the simulations discussed here employ the diffusive neutral model with T_{n} = T_{i}, and so the power loss due to charge exchange may be underestimated. The role of τ_{N} is discussed in further detail in Sec. V B.
The target electron pressure rolls over as the momentum losses, which are dominated by charge exchange reactions, begin to increase significantly, resulting from a strong increase in electron density in front of the target, as well as from a drop in temperature, which favors chargeexchange reactions over ionization. This leads to the rollover of the target ion flux. The region in front of the target becomes much cooler, and a region of strong temperature gradient moves upstream, with $Te,t<1$ eV. This front movement is also seen in other atomic processes (ionization, recombination and impurity radiation, Fig. 11), and eventually in the plasma density front as strong detachment is achieved, leading to a rollover in target electron density. The impurity radiation increases as the divertor becomes cooler and denser, becoming the dominant power loss mechanism, until f_{pow} saturates. As the neutral density continues to increase in the cool divertor, charge exchange momentum losses continue to increase until f_{mom} also saturates at $fmom\u223c1$. Note that although momentum and power losses due to recombination increase with upstream density, they remain negligible compared to other atomic sources throughout the density range explored. However, recombination appears as a significant contributor of the particle balance, Fig. 11.
One of the key assumptions of the standard twopoint model (TPM) that is often used as a first model to study SOL physics^{1} is that the heat transport is mostly due to (electron) heat conduction. While convection can be enabled in the extended TPM through some ad hoc f_{cond} parameter, the TPM itself does not provide a selfconsistent way to estimate the value of f_{cond}. In SPLEND1D, we find that convective heat transport in the current simulations is small but nonnegligible, Fig. 13. The difference between the SPLEND1D results and the TPM predictions is small at low upstream density ( $nu=1.6\xd71019\u2009m\u22123$) but becomes stronger as upstream density is increased, where the convective heat flux becomes significant.
V. EXPLORING THE ROLE OF FREE PARAMETERS
As mentioned in the introduction and in the derivation of the model, 1D codes come with strong assumptions and, as illustrated by Table II, free parameters that may affect the results of the simulations. These become especially important when attempting to use the results of such a model to explain experimental observations (although, in the case of interpretative simulations, some of them can be constrained by experimental data). It can, therefore, be important to understand how sensitive the code results are to these free parameters and assumptions. In this section, we review the impact of some of the main assumptions or parameters on the observations of Sec. IV.
A. Carbon concentration
Let us first examine the role of a parameter strongly affecting the plasma energy sink: the carbon concentration c_{c}. As expected, c_{c} has a strong impact on the target parallel particle flux [Fig. 15(a)], with the upstream density required for the $\Gamma ,t$ rollover decreasing strongly with c_{c}. This is due to the strong cooling induced by increased carbon radiation, leading to lower target temperatures for a given upstream density and temperature [Fig. 15(b)], and hence a facilitated access to detachment. We note here that the observed behavior (reduction of the peak particle flux achievable before rollover) is consistent with the power starvation mechanism described in Sec. IV A. As carbon concentration, and therefore radiation losses, increases, less power is available for ionization, leading to the earlier rollover of $\Gamma ,t$.
B. Neutral confinement time
The neutral model implemented in SPLEND1D essentially features three free parameters: the neutral confinement time τ_{N} [see Eq. (28)], the choice of neutral model and neutral temperature implementation. In this section, we focus on the free parameter τ_{N}. Figure 16 plots the effect of varying τ_{N} on the inferred target ion flux and integrated neutral density in the base case density ramp. This reveals that, as expected, τ_{N} has little effect at low density, where the neutral density is low in all cases. Only at high n_{u} does the effect of τ_{N} become significant: the neutral density integrated along the flux tube increases strongly with τ_{N}, associated with a strong decrease in the target ion flux. The upstream density at which the target ion flux rolls over is not strongly affected by τ_{N}, except for the most extreme case where τ_{N} has been decreased by a factor 10.
C. Heat flux limiter
The SpitzerHarm heat flux can overestimate the physical heat flux at low plasma collisionality, when the electron mean free path, λ_{e}, is large compared to the electron temperature gradient scale length, $L\u2207T=(\u2207Te/Te)\u22121$. Figure 17 shows the ratio of these scale lengths as a function of distance along the SOL, for the base case density and temperature profiles shown in Fig. 3. At the target, the electron collisionality is high enough that the classical SpitzerHärm heat fluxes are expected to fairly accurately predict the physical heat flux, without the need for a flux limiter. However, the upstream electron mean free path exceeds the temperature gradient scale length and so the effect of enforcing a heat flux limiter should be studied.
The strength of the heat flux limiter can be controlled with the parameter $\alpha {i,e}$. Target conditions are compared here for values of $\alpha i=\alpha e=[0.06,0.3,0.6,6,6\xd71017]$, thus testing values well below and above the typical value of 0.5 used in the literature (Refs. 20 and 21). This includes $\alpha \u226b1$, approaching the situation without heat flux limiters. The results are found to be little affected by this choice. The target ion flux rollover and thus the detachment threshold is weakly affected by the choice of $\alpha {i,e}$, Fig. 18, except for the smallest value of $\alpha {i,e}=0.06$. For $\alpha {i,e}\u22650.3$, the magnitude of the target ion flux varies by approximately 10% within the range of $\alpha {i,e}$ studied, Fig. 18, and the pressure drop along the SOL to the target is negligible (not shown). The target density is largely unaffected, and the target temperature is only affected at low upstream density ( $\u223c2\xd71019m\u22123$), not shown. These results hold for a range of input power levels.
Due to the apparent insensitivity of target parameters to the heat flux limiter coefficient for the base case conditions, and for reasonable values of $\alpha {i,e}\u22650.3$, and over a large part of the base case density ramp, $\alpha {i,e}$ was set to 0.6 in all simulations presented in this paper, unless stated otherwise.
D. Bohm boundary condition
The base case, Sec. III, leads to naturally supersonic flows ( $M\u22481.7$) at the target after rollover, Fig. 19(a). When enforcing M = 1 at the target in such conditions, the code showed numerical difficulties to converge with satisfactory particle and energy balances, with ϵ_{tot} [Eq. (56)] reaching up to $\u224812%$, Fig. 19(d). This is due to the presence of a very sharp velocity gradient required to slow down the flow to M = 1. This is alleviated by the introduction of some numerical viscosity ν_{num} [Eq. (26)]. Figures 19(b) and 19(c) plot the evolution of various target quantities for three different values of ν_{num}. The effect of this artificial viscosity on $Te,t$ (not shown), n_{t} and $\Gamma t,$ is modest, but it greatly improves the code convergence, with ϵ_{tot} reduced to a maximum of $\u22481.7%$ and $\u22480.23%$ for $\nu num=1.0\xd710\u22125[a.u.]$ and $\nu num=5.0\xd710\u22125[a.u.]$, respectively, Fig. 19(d). We note here that it is not yet clear whether the occurrence of such sharp velocity gradient is a consequence of the choice of parameters for the base case, or a general observation.
For both target boundary conditions M > 1 and M = 1, $\Gamma t,$ is very similar, Fig. 19(b). Similarly, the target temperature $Te,t$ is unaffected (not shown). However, as expected, a difference arises in the target density n_{t} and parallel velocity $u,t$. The case with M = 1 features a lower $u,t$ and a higher n_{t} compared to the M > 1. Since $\Gamma t,$ remains similar across these different cases, we conclude that enforcing the strict equality $u=cs$ leads to a redistribution of $\Gamma t,$ between its velocity and density contributions. This could have some implications in simulation codes that enforce $u=cs$, by promoting higher density at the target. Since many atomic source and sink terms scale with n_{e}, or even $ne2$, this will ultimately influence the particle, momentum, and energy balance of the system.
VI. ADVANCED STUDIES ENABLED BY SPLEND1D
SPLEND1D opens up a large number of possible SOL and detachment studies, such as the investigation of the effect of the total flux expansion, of the parallel connection length, of inout power sharing, or dynamical behavior. In Secs. VI A and VI B, we present some example studies on the role of different neutral models, ion vs electron heating ratios, and heat pulses effects on the SOL plasma.
A. Independent ion and electron temperatures
SPLEND1D can either solve a single energy equation, Eq. (12), assuming a proportionality relation between the ion and electron temperatures, or two separate energy equations for the electrons and ions, Eqs. (4) and (5). The base case presented in Sec. III has so far employed the T_{e} = T_{i} assumption. In this section, we relax this constraint by enabling both temperatures to be independent. It is well known that, at low collisionality, ion and electron temperatures can be different.^{1}
To set up the simulations, we split the power source between the H_{e} and H_{i} terms [Eqs. (4) and (5)], with either H_{e} = H_{i} (50/50 split of the total input power between electrons and ions), $He=3Hi$ (75/25 split of the input power between electrons and ions), or $3He=Hi$ (25/75 split of the input power between electrons and ions). In all cases presented in this section, the total power injected in the system is kept constant, as well as all the parameters presented in Table II. Similarly to Sec. IV, we perform density ramps by scanning $Hn\xaf$. Allowing $Ti\u2260Te$ affects the rollover threshold, Fig. 20(a), and leads to typically higher T_{i} than T_{e}, Figs. 20(b) and 20(c). This is easily explained by the lower heat conduction coefficient of the ions, compared to the electrons, Eqs. (16) and (17), while the convective heatflux and the equipartition term decrease this difference in transported heat flux. If the fraction of input power carried by the electrons is increased, the differences between T_{i} and T_{e} decreases, as expected, but remain significant. As density is increased, so does the collisionality and hence the equipartition term [Eq. (6)]. The difference between T_{i} and T_{e} becomes negligible near the target (where density is high and temperature low), Fig. 20(b), and reduces at the upstream location, Fig. 20(c). This is also evident when looking at the temperature profile along the flux tubes, Figs. 20(d)–20(f), which show profiles of T_{i} and T_{e} at increasing densities. Clearly, while in the attached case [Fig. 20(d)], the two temperatures are strongly different, they get closer to each other as density increases, Figs. 20(e) and 20(f).
B. Timedependent simulations
As mentioned in Sec. III B, SPLEND1D solves the 1D Braginskii equations as a timedependent problem. Hence, in addition to finding steadystate solutions, as discussed in earlier sections, it is also possible to use SPLEND1D to explore the dynamics of the system. In this section, we briefly highlight such a possible study enabled by SPLEND1D. We perform a timedependent simulation, based on a converged, detached simulation of the basecase, with $nu=3.0\xd71019\u2009m\u22123$. We then introduce a sequence of short ( $500\u2009\mu s$) pulses during which the heatsources, H_{e} and H_{i}, are amplified by a factor 10, before being relaxed to their initial values for $10\u2009ms$, Fig. 21(a). During each pulse, we observe a strong increase in the target parallel particle flux, $\Gamma ,t$, Fig. 21(b), together with a strong increase in the target electron temperature $Te,t$ Fig. 21(c), a sign that the plasma is reattaching during these events.
The target plasma and neutral densities are also strongly affected by the heat pulses. In particular, after a short increase in the target neutral density $nn,t$ (likely due to an increased recycling caused by the increased $\Gamma ,t$), $nn,t$ drops below its steadystate value, due to the ionization of most of the neutrals present in the system, as evidenced in Fig. 21(e) by the strong decrease in the integrated neutral density in the flux tube, associated with an increase in the integrated plasma density. This leads to an increase in the plasma upstream density n_{u} and a complex dynamics of the plasma target density n_{t}, which, after an initial rise and drop, peaks again before relaxing to its steadystate value. These results highlight how SPLEND1D can be used for timedependent simulations.
VII. CONCLUSION
This paper presented the SPLEND1D code: a 1D plasma fluid model that solves the Braginskii equations projected along a flux tube, together with a simple neutral fluid model. We have highlighted the key features of SPLEND1D, including its flexibility in terms of included SOL physics (magnetic field norm and pitch angle variation along the flux tube, dependent or independent electron and ion temperatures, possibility to enable, disable, or rescale any physical term contributing to the equations, etc.), numerical accuracy, and high computational speed. We believe this makes SPLEND1D a valuable tool for studying the complex dynamics of a divertor plasma interacting with a neutral gas and in contact with a wall. The SPLEND1D code was then used to investigate the physics of detachment in a reference simulation scenario. The code outputs the individual contributions of each term constituting its particle, momentum, and energy equations, allowing for the results to be easily interpreted. The rollover of the target ion flux, and the onset of detachment, was found to be mainly due to momentum losses owing to chargeexchange reactions, and power losses owing to impurity radiation. We note, however, that the importance of atomic chargeexchange losses is due, in part, to the absence of the molecular chargeexchange loss channel, which is known to be an important player in the simulation of detachment dynamics. To address this, we therefore plan to upgrade the SPLEND1D model to account for molecular species in addition to atomic ones. As with all reduced models, SPLEND1D comes with significant assumptions and free parameters. In this paper, we have shown how they can influence simulation results. However, some important effects are not captured by the current formulation of SPLEND1D. In particular, momentum and energy losses due to radial transport are absent, as well as drifts and SOL currents. These effects, intrinsically relying on perpendicular dynamics, are difficult to capture in a 1D code, likely relying on some ansatz associated with free parameters. Such possibility is left for future work.
SPLEND1D is currently being used to interpret TCV experiments, in particular those related to alternative divertor configurations, where it is being employed to elucidate the role of parallel connection length on the onset of detachment, as well as the impact of total flux expansion on detachment threshold and SOL parallel profiles.
In the future, we plan on continuing the development of the SPLEND1D model. From a numerical perspective, it would be interesting to further improve the capabilities of SPLEND1D, for instance by adaptively refining the mesh depending on local gradients of the solution, which may prove useful in situations where the detachment “front” is moving away from the target. From a physics perspective, we plan to extend SPLEND1D to add the capability to handle multiple plasma and neutral species, in particular molecules, which are known to play a very significant role in the divertor dynamics.^{31,32} The impurity model implemented in SPLEND1D, currently limited to constant impurity fraction assuming coronal equilibrium, should also be further enhanced to properly account for the impact of impurity transport and the subsequent occurrence of noncoronal effects. This could be achieved by adding a collisional radiative model for impurities, evolving in time the concentration of impurities' charge states, either with a trace assumption, similarly to Ref. 41, or using a more advanced fluid closure such as the Zdhanov closure.
ACKNOWLEDGMENTS
This work was supported in part by the Swiss National Science Foundation. This work has been carried out within the framework of the EUROfusion Consortium, via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion) and 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.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
O. Février: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (lead); Validation (equal); Writing – original draft (equal); Writing – review & editing (equal). S. Gorno: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal). C. Theiler: Conceptualization (equal); Methodology (equal); Project administration (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (equal). M. Carpita: Validation (equal); Writing – review & editing (equal). G. DurrLegoupilNicoud: Validation (equal); Writing – review & editing (equal). M. von Allmen: Software (equal); Validation (equal); Writing – review & editing (equal).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX A: COULOMB LOGARITHM
APPENDIX B: CONTRIBUTIONS OF MOMENTUM AND POWER SOURCE TERMS TO LOSS FACTORS
APPENDIX C: ADVECTIVE NEUTRAL MODEL
APPENDIX D: PRESSUREDIFFUSION EQUATION
APPENDIX E: ON THE ROLE OF THE NEUTRAL MODEL
The choice of the neutral model, and the neutral temperature implementation, can affect the plasma dynamics, primarily through the momentum and energy volumetric source terms. The diffusive model [Eqs. (32) and (33)] can be implemented with either T_{n} = T_{i} or with a constant, imposed value of T_{n} (in this section we set $Tn=1$ eV to model cold neutrals), while the advective model [Eqs. (28), (C3), and (C4)] can be implemented with either T_{n} = T_{i} or with the evolution of T_{n} following Eq. (C10).
Figure 22 shows the main differences in the results of the base case simulations for each neutral model outlined above. In the advective model with a selfconsistent T_{n}, and for the diffusive model with cold neutrals ( $Tn=1$ eV), we see an earlier target ion current rollover, along with a lower target plasma pressure, compared to implementing T_{n} = T_{i}. The cold neutrals facilitate detachment by increasing both the power loss factor with respect to the cases with T_{n} = T_{i}, as a result of increasing the power transferred through charge exchange reactions, by increasing the energy transferred from the ions to neutrals. This, in turn, leads to an increase in the momentum losses by promoting enhanced chargeexchange reaction and reduced ionization. We remark here that the density ramps presented in this paper are performed by ramping up the value of H_{n}. For different values of H_{n}, steadystate solutions can have the same upstream density n_{u}.
The advective model, in comparison with the diffusive model, features an additional neutral transport mechanism, the advective cross“fluxtube” transport. For a fixed τ_{N}, as shown in Fig. 22, the advective model displays a much weaker rollover. The total momentum and power loss factors shown in Fig. 22 also differ between the models.