Effect of scrape-off-layer current on reconstructed tokamak equilibrium

Understanding of tokamak plasmas has greatly benefited from the separation of time scales inherent in strongly magnetized plasmas. All orderings of strongly magnetized plasma equations show that the magnetic pressure and line-bending contributions from $J×B$ and the force-density from $∇p$ are the largest terms in the center-of-mass momentum equation. Assuming symmetry and nested flux surfaces allows one to derive the Grad-Shafranov equation¹ that describes the lowest-order steady-state fields. This solution is conventionally described as the plasma equilibrium. Perturbations that drive the plasma away from equilibrium launch Alfvén waves that quickly act to restore the plasma to equilibrium.² Transport occurs on slower time scales, while symmetry-breaking macroscopic instabilities occur on intermediate time scales between the transport and the stiff Alfvén time scale.

This paradigm well describes the experimental phenomenology, and solutions to the Grad-Shafranov equation are routinely calculated hundreds of times per plasma discharge to provide a measure of the location of magnetic-flux surfaces. As such, these solutions are indispensable to controlling the plasma and interpreting diagnostic data. Thus, Grad-Shafranov theory may be considered as one of the most successful applications of tokamak-plasma theory. Measurement of both the macro- and micro-scopic perturbations confirms the small-fluctuation assumption embedded in the theory: tearing modes, for example, have perturbed magnetic energies four orders of magnitude below the equilibrium stored magnetic energy,³ and ion-temperature-gradient-driven turbulence produces density fluctuations less than 1% of the equilibrium density.⁴ The small-fluctuation hierarchy coupled with slow temporal evolution permits the common approach of extended-MHD modeling about a Grad-Shafranov equilibrium (e.g., Ref. 5), as opposed to modeling the symmetric fields with the full extended-MHD equations (e.g., Ref. 6). The effects encompassed by extended-MHD are variable: extended-MHD refers to models beyond resistive MHD that include some combination of anisotropic thermal conduction and stresses,⁷ two-fluid evolution,⁸ finite-Larmor-radius closures,^9,10 and/or advanced drift-kinetic-equation closures.¹¹ 

One challenge to simulating small fluctuations about an equilibrium state is that errors in the equilibrium can be on the order of the perturbation magnitude. This issue is recognized by Grimm et al.¹² with the first “mapping code,” a code that employs numerical maps to transfer fields from the spatial discretization of an equilibrium code onto the spatial discretization of another code, where the second code typically assesses stability or the evolution of nonlinear perturbations. Despite considerable effort to increase the accuracy of mapping codes, mapping errors are difficult, if not impossible, to completely eliminate. In practice, extended-MHD codes, such as NIMROD^7,8 and M3D-C1,^13,14 recompute the Grad-Shafranov equilibrium with their native spatial discretizations to circumvent these errors. Burke et al.¹⁵ discussed an example of the impact of mapping errors. In that work, the NIMROD extended-MHD code is benchmarked with the linear-MHD codes GATO¹⁶ and ELITE¹⁷ on peeling-ballooning modes (PBMs) in tokamak equilibria. GATO, like other global linear-MHD codes, uses mapped equilibria while ELITE uses Miller equilibria¹⁸ that effectively act to provide the same accuracy as re-solving the Grad-Shafranov equation. In NIMROD calculations with mapped equilibria, better agreement with GATO is obtained, whereas when the NIMEQ code¹⁹ is used to recompute the Grad-Shafranov equation for the NIMROD initial condition, better agreement with ELITE is obtained.

While this example demonstrates the importance of high-accuracy Grad-Shafranov solutions that are recomputed within a code's native spatial discretization, a further example is provided with Fig. 1. The figure shows a comparison of the toroidal current density from both a mapped equilibrium along with a case where the Grad-Shafranov solution is recomputed with NIMEQ. This high-resolution computation is performed with a 72 × 512 finite-element grid with bi-quartic elements. The bumps in the mapped solution roughly correspond to the resolution of the 129 × 129 grid used during the reconstruction. Computations of edge-mode growth rates with the mapped equilibria give markedly different results relative to using the solution recomputed with NIMEQ.

The work of Burke et al. considered the simple case of an artificial peeling-ballooning equilibrium with exclusively closed and nested-flux surfaces without diverted magnetic topology. A diagram of the topology associated with cases considered in this work, a diverted lower-single-null shot within the DIII-D tokamak, is shown in Fig. 2 along with labels of the terminology used for the topological regions and contours. Operation with a divertor implies the presence of a separatrix field line that separates the open- and closed-flux regions. The surface immediately inside the separatrix is the last-closed-flux surface (LCFS). The separatrix and LCFS are indistinguishable from a practical standpoint, aside from the separatrix field lines that extend to the divertor strike point. The open field lines between the LCFS and the first wall define the halo region,^20,21 and in the modeling of Burke et al., this region was approximated as a closed-field-line region. Traditionally in linear-MHD modeling, this halo region is treated as a vacuum, i.e., it has neither current nor plasma density. In the extended-MHD modeling, the term halo region is used instead to denote that it is a cold-plasma region, capable of containing current, surrounding the hot, closed-field-line region.

The goal of this work is to relax the vacuum constraint of the halo region while remaining consistent with the measurements that are used to produce the reconstruction. For model fidelity to experiments, it is important to carefully consider the vacuum assumption in the halo region relative to the dynamics of study. Edge modeling is likely sensitive to current and flows present in the scrape-off-layer (SOL) region,²² a subregion of the halo that directs the energy and particle exhaust from the hot plasma onto the divertor. Most published simulations to date are initialized from equilibrium that do not have current or flow in the halo region. This is largely because this current is typically not included in reconstructed equilibria from the highly successful, workhorse EFIT code.^23,24 Without SOL current, reconstructions must have one of two undesirable properties: (1) either there is an artificial constraint on the current which must smoothly vanish at the separatrix or (2) there is a current discontinuity at the separatrix. The former constraint leads to the incorrect edge profiles, whereas the latter impacts convergence—particularly for high-order spatial discretizations such as those employed by NIMROD and M3D-C1. Importantly, cases with either constraint fail to include the measured profiles outside the LCFS.

In this paper, we discuss our method for adding SOL current to equilibrium reconstructions generated by EFIT. One notable case that uses this method is the DIII-D QH-mode modeling in Ref. 5. The difficulties with the current and flow discontinuities and a desire for more accurate modeling of the QH-mode instabilities motivated the development of our methods that include SOL profile gradients. This paper begins by reviewing the details of the EFIT reconstructions, and discussing important metrics for determining the quality of the reconstructions in Sec. II. The methods for recomputing equilibria with separatrices along with the details on how we add SOL current to these new equilibria are described in Sec. III. Example cases that use these methods are introduced in Sec. IV. We then compare the new and original equilibria using synthetic diagnostics that are similar to those used to constrain the original reconstructed equilibrium in Sec. V. Finally, we examine the minor modification to the linear growth rates of peeling-ballooning modes (PBMs) by the addition of SOL current in Sec. VI before making concluding remarks in Sec. VII.

Determination of the experimental configuration of tokamak plasmas has become essential for understanding and optimizing stability and confinement in fusion research devices. Reconstruction of the experimental equilibrium from a combined set of magnetic, temperature, and density measurements is computed by minimizing the error between modeled and observed signals. This technique, now routine, was pioneered with the EFIT code which originally used only magnetic diagnostics external to the first wall as a constraint.^23,24 Later, the motional Stark effect diagnostic^25,26 greatly improved the accuracy of the reconstructions by providing internal measurements of the magnetic and electric fields. More recently, measurements of the density and temperature profiles via Thompson scattering^27,28 and charge exchange recombination (CER)^29,30 spectroscopy provide further constraints on the pressure profiles.³¹ 

To quantify a reconstruction's accuracy, results from EFIT Grad-Shafranov solves are applied to a χ-squared test against the experimental measurements. In other words, the goal is to find a Grad-Shafranov solution that minimizes

where M_i is the measured signal, C_i is the computed signal, and σ_i is the measurement uncertainty. In this paper, we perform this χ-squared calculation to quantify the change that arises from the numerical errors associated with mapping, our recomputation of the Grad-Shafranov solution, and the addition of the SOL pressure profiles and associated current during the recomputation.

Typically, EFIT reconstructions do not include SOL current and treat the region between the separatrix and wall, i.e., the halo region, as current free. This implies that if the plasma discharge has a finite current on the separatrix, which is often the case for H-mode discharges that contain large bootstrap and Pfirsch-Schlüter currents, there will be a discontinuity at the LCFS. Equivalently within the context of the Grad-Shafranov equation, the EFIT profiles are constrained such that the pressure and toroidal magnetic field are constant in the halo region. With finite gradients at the LCFS, this leads to discontinuities in the first derivatives of these profiles. These inconsistencies lead to subtle but significant issues when evolving tokamak-edge unstable cases initialized from EFIT reconstructions. For example, in codes that retain the magneto-sonic wave physics both the SOL current and profile gradients must be included such that the equilibrium is consistent with force balance. Furthermore during nonlinear computations, any inconsistencies to the Grad-Shafranov equation, including discretization errors from discontinuous fields, can launch spurious magneto-sonic waves.

It must be noted that EFIT has the ability to include force-free, poloidal current in the SOL through finite gradients in the toroidal magnetic flux.^32,33 However, this capability is rarely exercised and there is no published work on the inclusion of finite pressure gradients in the SOL.

Initializing a NIMROD computation from a reconstruction is a two-step process that involves mapping from the reconstructed equilibrium, and then recomputation of the equilibrium. The mapping code FLUXGRID creates both a partially flux-aligned finite-element mesh and maps the reconstructed-equilibrium fields onto this mesh. This mapped solution is refined through the solution of the Grad-Shafranov equation with the NIMEQ code.¹⁹ The solution can then be iteratively passed between FLUXGRID and NIMEQ for further grid refinement and Grad-Shafranov solves, as needed. The Grad-Shafranov solve is formulated as a boundary-value problem where the boundary condition is specified by the value of ψ from the mapped reconstruction. This boundary condition constrains both the fields from the external coils and those generated by the internal plasma that are also determined by the pressure and toroidal-flux profiles. With this method, the recomputed fields closely resemble the reconstructed versions; however, mapping errors are largely eliminated and the Grad-Shafranov equation is satisfied up to an input tolerance.

Relative to the methods described in Ref. 19, we employ extensions that identify the open- and closed-flux regions of the domain in order to apply the appropriate fitted form of the pressure and toroidal flux profiles. Within the closed-flux regions, these profiles are specified by EFIT reconstructions as a function of normalized flux (⁠$ψn=(ψ−ψo)/(ψx−ψo)$ which is zero at the O-point where $ψ=ψo$ and unity on the separatrix where $ψ=ψx$⁠). With the exception of private-flux regions, regions with $ψn<1$ contain closed field lines and regions with $ψn>1$ contain open-field lines. Applying these profiles to the domain requires identification of the closed-flux region. In practice, we find the bounding LCFS contour and use a simple algorithm that counts the number of crossings of a line that extends from a finite-element node to the boundary to determine if each node is enclosed by the contour. The separatrix contour and location of the extrema of ψ in the core are allowed to change and must be recomputed after each iteration of the Grad-Shafranov solve. As the X-point on the separatrix contour associated with diverted magnetic topology consists of a stagnation point for field-line tracing, we instead elect to find a contour vanishingly close to the separatrix which numerically approximates the LCFS. For our purposes, we only need to bound the finite-element nodes that are within the separatrix. We employ two methods that use field-line tracing of the poloidal field to find this approximate LCFS contour. The first method uses bisection of the domain to determine where the field lines transition from closed to open within a specified tolerance. For the parallel implementation, it becomes an N-section method where each core is assigned a seed point between the known closed- and open-field-line locations. The second method is to use the OCULUS code to find the saddle point in ψ that is associated with the X-point, and a field-line seeded with a vanishingly small offset towards the O-point is used as the LCFS in order to avoid stagnation of the field-line tracing near the X-point. The latter method has the advantage of being somewhat more robust for high-resolution cases, but the disadvantage of not being amenable to a straight-forward parallel implementation.

As discussed earlier, EFIT solutions have zero profile gradients outside the LCFS. A cubic-spline fit of this data is used to evaluate these fields within the LCFS. Equilibrium generated using these spline fits with constant profile values outside the LCFS closely match the solution given by EFIT while largely eliminating mapping errors. One goal of this work is to contrast this recomputed solution with a solution that contains SOL-profile gradients.

SOL-profile gradients and associated currents are included by defining bounding contours and normalizing the flux as an extension to the methods used to determine the open- and closed-flux regions. The SOL region is defined as the region between the LCFS contour and contours with $ψn=ψsol$ and $ψn=ψpf$ where ψ_sol (>1) defines a contour(s) at the edge of the SOL region and ψ_pf (<1) defines a contour(s) in the private-flux region(s) (see Fig. 2). The algorithm that determines these contours is able to handle an arbitrary number of intersections of this region and the computational boundary of the domain and thus diverted topology such as double null configurations is tractable. The algorithm works as follows: The initial SOL bounding region is assumed to be the computational boundary. The algorithm checks the value of ψ_n at each finite-element node location outside the LCFS but contained within this “working” SOL region. If the condition $ψpf<ψn<ψsol$ is not satisfied, the location between the wall and the O-point where $ψn=ψsol$ or $ψn=ψpf$ is identified, and a new contour is traced. In the absence of integration error, the new contour is terminated at the boundary. The section of the domain that does not encircle the O-point is removed from the “working” SOL region until only nodes with the property $ψpf<ψn<ψsol$ remain. This new SOL contour, in addition to the LCFS contour, bound the SOL region.

The appropriate pressure and toroidal magnetic flux profiles are defined within the SOL region via two methods: either by fits from the experimental data (if available) or through modified bump function fits. Two fits are performed for each field: one fit for the SOL region immediately outside the separatrix with $ψn>ψx$⁠, and one for the private-flux SOL region with $ψn<ψx$⁠. The modified bump function uses the form,

This function has vanishing derivatives of all orders at its endpoint and thus ensures that the current goes smoothly to zero at the transition contour between the SOL region and the current-free region. With this form, five free parameters are available (ψ_sol, ψ_pf, f₀, f_c, and Δ). The values of ψ_sol and ψ_pf are inputs that influence the width of the SOL and may be inferred from experimental measurements. The other three parameters may be determined either by requiring C² continuity at the LCFS or by setting the functional value at ψ_sol and enforcing C¹ continuity at the LCFS. The former constraint has the advantage of producing a C¹ smooth current profile, whereas the latter method allows specification of the density and/or temperature in the current-free regions. If the fit requires that the function first reverse the sign of its derivative, a truncated Gaussian is fit to half of the domain followed by the bump function as is used in Ref. 34.

In order to demonstrate and quantify the impact of our methods that recompute the Grad-Shafranov solution and add SOL-profile gradients, we choose two specific cases to study in detail: a low SOL-current case and a high SOL-current case. The first of these cases is an EFIT reconstruction of DIII-D shot 160414 at 3025 ms with profiles as shown in Fig. 3. Profiles are specified as a function of normalized flux. The profiles shown are as included from the EFIT reconstruction and for fits with SOL-profile gradients (NIMEQ-SOL) along with the experimental Thomson and CER measurements.⁴² This shot is from an experiment of lithium pellet injection for edge-localized-mode pacing and the specific time represents the last 20% of the inter-ELM period. The reconstruction contains a relatively cold plasma at the LCFS, and thus a small pressure gradient in the fitted profiles within the SOL region. This implies modest SOL current and modification to the resulting equilibria when the SOL profiles are included in the Grad-Shafranov solve. In particular, at the LCFS, T_e = 72 eV, T_i = 470 eV, and $ne=8.6×1018$ $m−3$⁠; and in the current-free region, T_e = 20 eV, T_i = 50 eV, and $ne=2.1×1018$ $m−3$ where the latter values only apply to equilibria with SOL-profile gradients.

The second case studied is an EFIT reconstruction of DIII-D shot 145098 at 1800 ms as shown in Fig. 4. Again, the profiles shown are as included from the EFIT reconstruction and for fits with SOL-profile gradients (NIMEQ-SOL) along with the experimental Thomson and CER measurements. This shot is from a DIII-D QH-mode experiment with ITER-like shaping during a period with edge harmonic oscillations (low-$nϕ$ perturbations). This shot contains a relatively large pressure gradient from the fitted profiles in the SOL region. Thus, relative to the reconstruction from shot 160414, we expect greater modifications to the equilibria resulting from the inclusion of the SOL-profile gradients in the Grad-Shafranov solve. In particular, at the LCFS T_e = 230 eV, T_i = 1020 eV, and $ne=4.8×1018$ $m−3$⁠; and in the current-free region, T_e = 50 eV, T_i = 50 eV, and $ne=2.7×1018$ $m−3$ where the latter values only apply to equilibria with SOL-profile gradients.

As is clear from Figs. 3 and 4, the NIMEQ-SOL ion-temperature profiles do not match the CER measured data. This is a consequence of a two-temperature model which over-constrains the profiles. Specifically, with a two-temperature model assuming quasi-neutrality (⁠$ne=Zini$⁠) only the electron and ion fluids contribute to the pressure

Thus, only three of the four profiles shown in the figures (n_e, p, T_e, and T_i) can be matched exactly. Initializing computations with the measured pressure profile is essential as both the Grad-Shafranov solution and ideal stability are critically dependent on this profile. Secondary to this, the electron temperature profile determines the resistivity profile and the density profile sets Alfvén speed (among other collisionality parameters).

Within the context of an extended-MHD simulation, the ion-temperature profile typically becomes significant only when a two-fluid model is evolved and thus is often allowed to vary with respect to experimental measurements. In experiment, the pressure has contributions from both impurities and non-Maxwellian, or “hot,” ion particles from neutral-beam injection

Comparing Eqs. (3) and (4), the ion temperature used to initialize our simulations includes the contributions from the other species

As expected from this relation and as shown in the plots, the NIMEQ-SOL ion temperature (equivalent to $Ti2T$ in Eq. (5)) over estimates the measured ion temperature. For these cases, Z_i is chosen such that the edge ion-temperature profile approximately matches the measured data in Figs. 3 and 4.⁴³ Thus, computations in this work use the NIMEQ-SOL fits for $ψn<1$ and only the computations with SOL-profile gradients use the NIMEQ-SOL fits where $ψn>1$⁠. Modeling that includes separate species for hot particles and/or impurities is required to eliminate the discrepancy with the measurements in the ion-temperature profile, and is thus planned for future study.

To further examine the reduction of mapping errors through the recomputation of the Grad-Shafranov solution, Figs. 5 and 6 show the resulting current-density distributions from the mapped and recomputed (GS) cases without SOL-profile gradients, along with the two cases with SOL-profile gradients that are discussed later in this section, based on reconstructions from shots 160414 and 145098, respectively. The DIII-D first wall, LCFS, and SOL contours are superimposed into the plotted computational domains. Shot 145098 uses a reversed plasma current relative to 160414, which uses the standard current orientation for DIII-D. These cases use a 72 × 64 high-order finite-element mesh with bi-quartic elements. For both reconstructions, the mapped current is clearly distorted and contains numerical oscillations. This is particularly evident for the edge current shown in the zoomed insets of the figures and outside the LCFS as shown in the inset figures that magnify the details near the X-point. The numerical bumps in the current profiles inside the LCFS are a result of the mapping and are associated with the resolution of the EFIT grid. They are not eliminated but rather only resolved with enhanced NIMEQ resolution (see Fig. 1). While it is possible to partially circumvent this issue with the closed-flux mapped current by using a high-resolution EFIT (see, for example, Refs. 35 and 36), in practice most EFIT reconstructions are generated at relatively low resolution relative to what is required for extended-MHD computations. The spatial requirements to solve for Grad-Shafranov equilibria are less stringent than those to solve for 3D MHD perturbations where, for example, the edge computations presented in Sec. VI use a finite element grid with 72 × 512 with high-order bi-quintic elements. The current outside the LCFS in the mapped case results from the representation of the discontinuous current profile on NIMEQ's C⁰ finite-element spatial discretization. The mapped plots in Figs. 5 and 6 are generated with finite-element calculations that compute the poloidal magnetic field and current from the mapped ψ and $RBΦ$ fields. Alternatively, splines may be used to map the poloidal magnetic field and current density. While this method produces smooth, but not consistent, fields, the derivatives of these fields with NIMEQ's C⁰ finite-element representation are used in extended-MHD calculations. Thus, mapped magnetic fields and current densities with a spline spatial representation in effect hide, but do not eliminate, the mapping errors.

For the cases with SOL current, we use $ψn,sol=1.1$ and $ψn,pf=0.96$ in the fits to the electron density and temperature profiles from Thomson scattering measurements and extrapolate the ion temperature to an assumed 50 eV. The resulting profiles are shown as the dashed purple line in Figs. 3 and 4. The half width of the electron pressure profiles are roughly $3.3 mm$ and $2.7 mm$ at the outboard mid-plane and $7.3 cm$ and $4.7 cm$ at the divertor plate for the cases from shots 160414 and 145098, respectively. This results in SOL widths that are roughly consistent with the measured half width of the heat-flux during the later half of the inter-ELM period of DIII-D ELMy H-mode discharges in Ref. 37.

Currently, we do not have diagnostic information about appropriate profiles for the private-flux region. Two methods employing the bump function extrapolations are compared with profiles shown in Fig. 7. In the first method, C² continuity is enforced at $ψn=1$ resulting in a high-pressure, high-density private-flux region. With the second method, C¹ continuity is enforced at $ψn=1$ and the values at ψ_pf are set to the same as those at ψ_sol. In the figures and tables, these cases are referred to as GS-SOL and GS-SOLpf, respectively. The profiles generated with the second method (GS-SOLpf) roughly match the measured heat-flux profile.^34,37

In addition to the mapped and recomputed cases without SOL current, Figs. 5 and 6 plot cases with SOL currents. With low pressure in the private-flux region the toroidal current density reverses locally near the divertor strike points after the Grad-Shafranov solve as shown in the inset figure for the GS-SOLpf case. Currents near and outside the limiter are likely an artifact of our method. In this cold-plasma region, additional terms in the momentum equation become large (e.g., interactions with neutrals and the first wall), which are outside the scope of the Grad-Shafranov equation. As such, the approximation that quantities are functions of the flux surfaces breaks down. However, at present our focus is on the effects of the SOL-profile gradients near the LCFS where there is the potential for interaction with perturbations that originate from inside the separatrix. With SOL-profile gradients, there is a modest modification (less than 1%) to the plasma current. For the cases shown here, we re-normalize the total current to match the value from the EFIT reconstruction using the method of Ref. 38. Cases with and without current re-normalization produce similar results, where the $χ2$ values discussed in the Sec. V are slightly smaller for cases with re-normalization.

The poloidal current in the SOL that flows into and out of the divertor plate are solely an effect of the profile gradient in the toroidal magnetic flux within the SOL. For these cases, a bump-function fit that enforces C² continuity at $ψn=1$ determines the toroidal-magnetic-flux profiles in the SOL. The resulting poloidal currents have a maximum value on the divertor plate of 4000 and 500 $A/m2$ for the cases from shots 160414 and 145098, respectively. This current must be less than the ion saturation current (⁠$Jmax=neecs$⁠). With a Deuteron ion species, a conservative calculation of the ion saturation current (using the vacuum temperatures and densities from each case) is over $2×104$ $A/m2$ for both cases.

Analysis between each nonlinear Grad-Shafranov-solve iteration provides a confirmation that the macroscopic quantities of the equilibrium (e.g., total current, toroidal flux, and internal energy) are invariant between the mapped equilibrium and the new solution. In the interest of showing that our methods which recompute the Grad-Shafranov equilibrium and add SOL profiles and current only minimally impact the relative agreement with measurements, we use a more advanced comparison that computes a $χ2$ value. The Python code NIMNOSTICS is used to calculate a $χ2$ from the shots and sub-cases previously discussed. NIMNOSTICS models the magnetic coils, MSE, and Thomson scattering through local evaluations of the fields where linear interpolation is used between finite-element nodes. Thus, the boundary of the domain is chosen as the approximate vacuum vessel instead of the limiter in order to encompass the magnetic coil locations for comparison.

Summaries of the $χ2$ value by diagnostic and case for the reconstructions from shots 160414 and 145098 are shown in Tables I and II, respectively. In both shots, the mapped and recomputed equilibria with SOL-profile gradients (GS cases) are roughly equivalent indicating that the equilibrium resulting from the recomputation is substantially similar. For the cases with modest SOL-profile gradients and currents (shot 160414), the $χ2$ values for a given measurement either decreases and the deviations of the LCFS contour are modest or are roughly the same as the mapped and GS cases. However, cases from the shot with large SOL current and profile gradients exhibit mixed results with some measurements decreasing in $χ2$ (Thomson electron temperature) and others increasing in $χ2$ (Thomson electron density, MSE and coils) when comparing relative to the mapped and GS cases. The deviation of the LCFS contour is also as much as 2 cm on the upper side of the contour. This deviation is apparent in the inset figures of Fig. 6.

In order to examine the source of the changes in $χ2$ in detail, Figs. 8 and 9 show the local values of $χ2$ from Thomson electron density, MSE, and magnetic coils measurements for our four different cases on shots 160414 and 145098, respectively. With shot 160414, the aggregate $χ2$ values are improved or comparable for every diagnostic when comparing the cases with SOL gradients to those without (as seen in Table I). We note that $χ2$ values are particularly improved with SOL-profile gradients for the Thomson measurements, and Fig. 8 shows that this is a result of improved agreement in the SOL region. The results are mixed for shot 145098 (as seen in Table II). While the aggregate $χ2$ value for the Thompson electron temperature profile improves when SOL-profile gradients are included, the $χ2$ value for the Thomson electron density profile is degraded. Examination of Fig. 9 shows that while the $χ2$ values in the SOL region are smaller when the SOL-profile gradients are included, the values near the LCFS become large consistent with the approximately 2 cm movement of the separatrix line relative to the cases without SOL-profile gradients. Additionally, the $χ2$ values for the magnetic-coil measurements become marginally larger near the divertor region as these values are affected by the inclusion of toroidal current in this region.

As an additional test of the quality of the equilibrium, we calculate the flux-surface-averaged geometrical quantity discussed in Ref. 39,

where $b=B/|B|$⁠. Figure 10 shows the result of this calculation and the decomposition of the expression in Eq. (6) into separate contributions from the two integrals in the expression (referred to as intA and intB, respectively) for both shots examined in our studies. This integral is known to vanish analytically;⁴⁰ however, this behavior is only reproduced numerically with recomputation of the Grad-Shafranov solution. Consistent with the relatively large motion of the flux surfaces (⁠$2–3 cm$⁠) in the large SOL-current case (145098), the contributing integrals for this case differ slightly between the GS and GS-SOL cases where only the latter includes the SOL current. The contributing integrals lie on top of each other for the low SOL-current case (160414). The improved equilibria provided by recomputation of the Grad-Shafranov solution are critical in NIMROD drift-kinetic computations.¹¹ For example, an accurate account of $b·∇ln|B|$ in latter integral of Eq. (6) (intB) is essential when assessing how trapped particles affect parallel closures for NIMROD's fluid system.

In order to assess the impact, if any, on linear stability, we examine the toroidal-mode-number (⁠$nϕ$⁠) growth-rate spectrum for shot 145098 at 1800 ms with and without SOL-profile gradients. The reconstruction from shot 160414 is during a stable inter-ELM period and thus it is not considered. As seen in Fig. 11, the growth rates are only minimally impacted by the inclusion of SOL-profile gradients where the difference between the growth rates is at most 5% (at $nϕ=6$⁠). A 72 × 512 mesh with bi-quintic elements is used for these computations. The growth rates with the SOL current are larger than those without for $nϕ<30$ and 1% smaller at $nϕ=30$⁠. This is consistent with the effect of enhanced resistivity and decreased density outside the LCFS as associated with the SOL-profile gradients that leads the dynamics in this region to produce a more vacuum-like response (see Refs. 14 and 15) whereby the low-$nϕ$ modes are destabilized and the high-$nϕ$ modes are stabilized. In order to investigate this response further, we examine a case with a low electron temperature, 1 eV, at the edge of the SOL region (annotated as GS-SOLpf-lt in the figure). The stabilizing effect at low $nϕ$ and destabilization at high $nϕ$ from the enhanced vacuum-like response is more apparent for this case. Relative to the case without SOL current, the growth rates for this case are 7% larger at $nϕ=6$ and 8% smaller at $nϕ=30$⁠. Thus, the effect of the SOL-profile gradients is modest compared to the effect of drift-stabilization (see, e.g., Refs. 36 and 41)—a result that is consistent with prior PBM calculations.¹⁷ 

As seen in Fig. 12, the mode structure is localized within the region just inside the LCFS. Convergence is affected by discontinuity in the current-density profile when SOL-profile gradients are not included even though the mode is localized away from this discontinuity. NIMROD does not enforce the $∇·B=0$ constraint through the spatial discretization, but rather converges to a solution with small $∇·B$ error. The $ℓ2$-norm of $∇·B$ is reduced by approximately 50% in the cases with SOL-profile gradients relative to cases without as a result of the continuous equilibrium profiles.

Relative to the effect on linear computations (Sec. VI), the inclusion of SOL-profile gradients and associated continuous current-density profiles has a greater impact on nonlinear modeling with perturbations that are advected across the separatrix such as the studies of QH-mode evolution of Ref. 5. Nonlinear computations can be affected by the inclusion of SOL-profile gradients in multiple ways: the spatial resolution required to converge on the dynamics at the LCFS is less with a continuous current profile as the dynamics are affected by discontinuities in the current and by the changes in the vacuum-like response as described in Sec. VI.

Perhaps more importantly, the methods described to extrapolate the thermodynamic profiles in the SOL can be applied to modeling with flows. Typically, the measured flows do not vanish at the LCFS and can be extrapolated to zero in the SOL region. Including this extrapolation both affects the dynamics of the perturbations as they cross the LCFS and prevents the computationally pathological case where perturbations may be advected quickly inside the LCFS and not at all outside. Again, an example that applies these methods to modeling with flows is found in Ref. 5.

The inclusion of SOL-profile gradients may have an effect outside the context of modeling with initial value codes. The last-closed flux-surface locations are shifted by up to $3 cm$ in the cases included in this study. While this shift may not greatly affect MHD stability, it could have an impact on the methods that are predicated on highly accurate reconstructions such as RF injection for current drive and/or tearing mode stabilization. These considerations may motivate the inclusion of the SOL-profile gradients within the reconstruction itself.

One limitation to the methods described here is that the flux from the plasma that penetrates through the walls remains fixed. The flux can be decomposed into plasma and external-coil contributions (⁠$ψ=ψplasma+ψext.coil$⁠). A potential extension to this work is to perform a free-boundary computation where $ψext.coil$ is fixed but ψ_plasma is allowed to vary. Additionally, $ψext.coil$ could be extended to include contributions from return currents flowing through the wall. Ultimately, as these computations become more sophisticated it may be better to include the SOL-profile gradients in the $χ2$ minimization performed during the reconstruction.

Even with this caveat, our methods represent substantial progress on the initial condition for edge modeling. In particular, we have developed a workflow whereby SOL-profile gradients and current can be included in the initial condition for NIMROD even if they are not included in the reconstruction. Using both global (e.g., total current) and local (e.g., the separatrix location) metrics as well as a $χ2$ test, we quantify the impact of our methods on the accuracy of the initial condition after inclusion of SOL-profile gradients. We find that this impact is small and that the modified initial condition closely resembles the state found by the reconstruction. While linear stability is modestly impacted by the inclusion of the SOL-profile gradients through an enhancement of the vacuum-like response, we argue that our methods are more important for nonlinear modeling of dynamics across the separatrix.

The authors would like to thank Alessandro Bortolon for generation of the EFIT reconstruction of shot 160414 at 3025 ms, Jim Myra and Dan D'Ippolito (deceased) for discussions of scrape-off-layer physics and ideas on how to extend the profiles, Carl Sovinec and Eric Howell for the development of NIMEQ and feedback on this development, and Stuart Hudson for the use of the OCULUS code and assistance on integration of the code into the NIMROD code base. This material is based on the work supported by U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences under Award Nos. DE-FC02-06ER54875,¹ DE-FC02-08ER54972¹ (Tech-X collaborators), DE-FC02-04ER54698² (General Atomics collaborators), and DE-FG02-03ER54692³ (Auburn collaborators) and used the resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.