The fluid transport code [trans-electric field (Er) module] under the BOUT++ framework has been used to simulate divertor heat flux width and boundary Er with all drifts and the sheath potential in the scrape-off layer. The calculated steady state radial Er in the pedestal region has been compared with that of experimental measurements from the Alcator C-Mod tokamak. The magnitude and shape of Er are similar to those of the experimental data. In order to understand the relative role of cross-field drifts vs turbulent transport in setting the heat flux width, four C-Mod enhanced D_{α} H-mode discharges with a lower single null divertor configuration should be simulated. BOUT++ transport simulations with cross-field drifts included yield similar heat flux width *λ*_{q} to that of experimental measurements (within a factor of 2) from both the probe and the surface thermocouple diagnostics and show a similar trend with plasma current to that of the Eich experimental scaling. The simulations show that both drifts and turbulent transport compete to determine the heat flux width. The magnetic drifts play a dominant role in setting the divertor heat-flux width, while the E × B drift decreases the heat flux width by 10%–25%, leading to improved agreement with the experiment relative to Goldston’s model. A turbulence diffusivity scan (*χ*_{⊥}) identifies two distinct regimes: a drift dominant regime when *χ*_{⊥} is small and a turbulence dominant regime when *χ*_{⊥} is large. The Goldston heuristic drift model yields a lower limit of the width *λ*_{q}.

## I. INTRODUCTION

The high confinement mode (H-mode) features a steep pressure gradient at the plasma edge, called the pedestal region. The H-mode is the baseline operation regime for the high-performance requirements of future steady state tokamaks such as ITER.^{1,2} Diffusive and convective transport of plasma from the strong gradient region of the edge pedestal into the scrape-off-layer (SOL) can impact power exhaust, which will be a critical issue for future magnetic confinement fusion devices. From power balance in the steady state, the power across the separatrix must be equal to the input power, plus fusion alpha particle power, and minus core radiated power. As the alpha power in burning plasmas is large, this can lead to large power fluxes on the divertor or other material surfaces, leading to the damage of the plasma-facing components (PFCs).^{2} The width of the divertor heat flux footprint, or so called scrape-off layer (SOL) width *λ*_{q}, and the peak of the heat-flux density on divertor targets are critical parameters for characterizing the particle and power exhaust.^{3} The size and scaling of the heat-flux width have been the subject of recent experimental investigations, and an inverse scaling dependence of the heat flux width on plasma current I_{p} is observed from a dataset combining multiple machines,^{4–8} including NSTX, ASDEX Upgrade, Alcator C-Mod, JET, and DIII-D H-mode discharges. A database from these experiments gave rise to a scaling law by Eich, in which *λ*_{q} is nearly inversely proportional to the poloidal magnetic field at the outer midplane with weak dependencies on any other key SOL parameters.^{9} A similar dependence is also found in the EAST tokamak.^{10} From this scaling, a very narrow SOL inter-ELM (Edge Localized Mode) Heat flux width, *λ*_{q} ∼ 1 mm, is projected for the ITER baseline burning plasma scenario (Ip = 15 MA). The divertor target power loads with such a narrow width would in reality be impossible to handle, which raises a serious concern.^{2,11} Even though the divertor detachment is the most promising means for mitigating target heat fluxes, the narrow *λ*_{q} will make the plasma harder to detach.^{12} To date, no consensus has been reached for the magnitude of the heat flux width for ITER and future reactors. Additional recent results on C-Mod have shown consistency with the Eich scaling up to ITER-relevant B_{p} of 1.3 T.^{13} However, gyrokinetic simulations and turbulence simulations at the ITER scale indicate *λ*_{q} ∼ 5–6 mm.^{14–16} No matter which prediction is right, there is no doubt that the heat flux width definitely needs to be broadened by any applicable techniques.

For a given input power, the peak of the heat flux is inversely proportional to a critically important parameter—the scrape-off-layer heat flux width. It is generally thought that *λ*_{q} is determined by the competition between perpendicular transport governed by drifts and turbulence and parallel transport in the SOL. In order to explain the physics behind the inverse scaling of *λ*_{q} with I_{p}, Goldston proposed a heuristic drift-based (HD) model, which is consistent with the international multimachine scaling. The hypothesis in the HD model is that magnetic drifts balance the nearly sonic parallel flow in setting the ion particle flux width and cross-field transport at the separatrix is dominated by classical drifts.^{17} In the HD model, the power across the separatrix is empirical and the E × B drift is not considered. Recently, electrostatic kinetic simulations with Particle-In-Cell (PIC) code XGC reproduced the trends in the Eich scaling and HD model for existing tokamaks.^{14,18,19} In addition, the turbulence code (6-field two-fluid module)^{20–23} under the BOUT++ framework has been developed to study the edge physics in tokamaks, including ELMs and the divertor heat fluxes during the ELMy H-mode. The role of turbulence in setting the electron heat flux width has been studied for the inter-ELM phase of DIII-D,^{24} EAST ELMing H-mode discharges,^{25} and C-Mod enhanced D_{α} (EDA) H-mode discharges.^{26,27} The turbulence simulations of C-Mod discharges further demonstrate that reducing the radial electric field Er, which is calculated from the transport module under the BOUT++ framework without the cross-field drifts, can effectively enhance the radial turbulence transport, thus increasing the heat flux width.^{27}

The BOUT++ turbulence code mainly focuses on fast turbulence dynamics and a scale of up to a few thousand Alfvén times. Due to slow ion parallel dynamics, even though the turbulence simulations are able to reach a steady state solution at the midplane for both ions and electrons and electron divertor heat flux is able to reach a well saturated state, the ion heat flux on the divertor targets has not been easily characterized for the BOUT++ turbulence simulations. For quick scoping studies, a new transport module, including all drifts, which mainly focuses on 2D plasma evolution to the steady state on a longer transport time scale, is developed under the BOUT++ framework^{28,29} to get a steady state solution of the divertor heat flux for both electrons and ions. In this work, we will use the new transport code and the same 4 discharges from C-Mod as in Refs. 26 and 27 (1) to study the relative role of cross-field E × B and magnetic drifts vs turbulent transport in setting the heat flux width and (2) to demonstrate the transition from drift to turbulence dominant regime. In addition, we will make the comparison of electric field between simulations and experimental measurements with or without magnetic drifts. These studies feature the following significant new results:

Compared the simulated electric field Er with that of experimental measurements, as shown in Fig. 5.

Compared the simulated divertor heat flux width with that of measurements as shown in Figs. 8(a) and 8(b). Most importantly for the first time, the comparisons are made between simulations and experimental data from both probe and the surface thermocouple (STC) measurements.

Compared transport simulation results with turbulence simulation results, HD model, C-Mod experimental data, and experimental scaling law, as shown in Fig. 9.

Performed sensitivity studies of divertor heat flux vs separatrix temperature, as shown in Fig. 10.

Confirmed that the magnetic drift has a significant influence on the heat-flux width, while the E × B drift decreases the heat flux width by 10%–25%, which gives improved agreement with the experiment relative to Goldston’s model, as shown in Fig. 11.

Identified two distinct regimes via a turbulence diffusivity scan (

*χ*_{⊥}): a drift dominant regime when*χ*_{⊥}is small and a turbulence dominant regime when*χ*_{⊥}is large. Magnetic drifts give a lower limit of the width*λ*_{q}, similar to Goldston, as shown in Fig. 12.

These studies suggest that both drifts and turbulent transport compete to determine the divertor heat flux width in C-Mod EDA H-mode discharges.

The rest of this paper is organized as follows. Section II introduces the transport physics model with all drifts and the setup of the heat flux width simulations. The steady state Er has been compared with that of experimental measurements of a C-Mod discharge in Sec. III. The simulation results of the heat flux for C-Mod EDA H-mode discharges are described in Sec. IV. Importantly, the heat flux width calculated from the transport code (trans-Er module) has also been compared with the turbulence code (6-field module) results under the BOUT++ framework, experimental measurements, and the Goldston HD model. The competition between drifts and turbulence in setting the divertor heat flux is shown in Sec. IV C. Section V summarizes the results.

## II. PHYSICAL MODEL AND SETUP OF SIMULATIONS

In this paper, the transport model (trans-Er) with all drifts under the BOUT++ framework has been used to simulate the plasma transport and power deposition in C-Mod discharges.^{30} Starting from the Braginskii equations,^{20,31,32} the evolving variables are the vorticity (*ϖ*), electrostatic potential (*ϕ*), ion density (*N*_{i}), parallel ion velocity ($V\u2225$_{i}), and electron and ion temperature (*T*_{e} and *T*_{i}, respectively). The set of transport equations has been given in Ref. 30. For completeness, we write the equations here again,

Here, the E × B velocity is $VE\xd7B=1B0b0\xd7\u2207\varphi $ and the ion diamagnetic velocity is $Vdiai=1ZieNiB0b0\xd7\u2207Pi$, where *P*_{i} = *kN*_{i}*T*_{i} is the ion pressure and *B*_{0} is the magnetic field. Symbols *κ*$\u2009\u2009\u2225$_{i} and *κ*$\u2009\u2009\u2225$_{e} are the flux limiting parallel thermal conduction coefficients for ions and electrons, respectively. The expressions are the same as Ref. 27. The flux limited parallel thermal conduction is dominant by competition between the Spitzer-Harm and free-streaming regime. The effective parallel thermal conduction coefficients are set as^{27}

The parallel heat conduction coefficients $\kappa \u2225iSH$ and $\kappa \u2225eSH$ are set in the form of Braginskii expressions, $\kappa \u2225iSH=3.9NiVth,i2/\nu i$ and $\kappa \u2225eSH=3.2NeVth,e2/\nu e$. Here, *V*_{th,i} and *V*_{th,e} are ion thermal velocity and electron thermal velocity, respectively. Symbols *ν*_{i} and *ν*_{e} are the collision rate for ions and electrons, respectively. The free streaming heat conduction for each species is $\kappa \u2225ifs=NiVth,iqR0\alpha j$ and $\kappa \u2225efs=NeVth,eqR0\alpha j$, where the kinetic effects are introduced through *α*_{j} and *α*_{j} = 1.0 is used in the simulation. Here, *q* is the safety factor and *R*_{0} is the major radius.

In the ion density and ion temperature equations, E × B drift and ion diamagnetic drift are considered. The E × B and electron diamagnetic drifts are included in the electron equation. The ion polarization velocity is a small correction to the E × B and diamagnetic velocities for transport problems, which is neglected in Eqs. (1)–(4) but is kept in Eq. (5) of the vorticity equation. This transport model with all drifts and radial electric field calculation are described explicitly, as in Ref. 30. It is worth noting that in the following simulations, the effective radial velocity from particle diffusion is included in the velocity, temperature, and vorticity equations as part of convection, while the source terms associated with neutral particles are not included, such as ionization, recombination, charge exchange, and radiation. The plasma fluid equations are discretized in a 2D orthogonal flux-surface mesh using a finite-difference algorithm.

Using a method similar to that in Ref. 30, for a steady state solution without volume sources, drifts, and flow, the ion continuity equation (1) becomes

Given the core boundary particle flux and experimentally measured density profile inside the separatrix, the particle diffusivity D_{⊥} can be obtained by $D\u22a5\u2202Ni\u2202x=\u2212\Gamma |x=0$. Here, the core boundary flux is determined by nominal boundary particle diffusivity *D*_{⊥} = 1 m^{2}/s and the experimentally measured density profile at the core boundary.

With the effective radial velocity from particle diffusion in the temperature equations, for a steady state solution without volume sources, drifts, and flow, the electron and ion temperature equations [Eqs. (3) and (4)] become

The heat transport coefficients *χ*_{⊥i} and *χ*_{⊥e} can be calculated by

if a boundary flux $\Gamma i,e|x=0=Ne\chi \u22a5e\u2202Te\u2202x|x=0$ is given. In our simulations, at the core boundary x = 0, the gradients are calculated from experimental profiles and a nominal thermal diffusivity *χ*_{⊥i,e} = 1 m^{2}/s is given.

The sheath boundary conditions (SBCs) are applied on the divertor target surface which we assume that the ion parallel velocity matches the sound speed: $V\u2225i=Cs=kTi+TeMi,$ where *M*_{i} represents ion mass. The same SBCs are used as Ref. 30, which the total parallel current through the sheath is described as

The sheath potential can be derived as

In our simulations, the cross-field drifts (E × B and magnetic drifts) are included in the plasma transport equations, which are important in calculating the radial electric field (Er) and the parallel heat flux on the divertor targets.^{30} To compare the Er calculation from our transport code with that of experimental results, a C-Mod EDA H-mode discharge (No. 1080321020) with lower single null divertor geometry is simulated. Figure 1 shows the simulation domain, which includes the plasma edge, the SOL region, and the private flux region (PFR). In the radial direction, the normalized poloidal flux is from ψ = 0.9 to ψ = 1.05, where ψ = 1.0 is the location of the magnetic separatrix. The grid generated from the equilibrium file (EFIT g-file) and the spatial resolution of the mesh is 260 grid points in the radial direction and 64 grid points in the poloidal direction. Section III shows the results of this comparison.

Next, in order to explore the physics behind the Eich scaling of the heat flux width *λ*_{q} vs I_{p}, such as drifts vs turbulent transport, a set of four C-Mod EDA H-mode discharges with lower single null divertor geometry are simulated. Table I shows the detailed key plasma parameters for 4 different currents with shot No. 1100223012, No. 1100212023, No. 1100303017, and No. 1160729008. The simulation domain is the same as for the discharge in Fig. 1. The initial plasma profiles inside the separatrix of the four shots used in these simulations are fitted onto a radial coordinate of the normalized poloidal magnetic flux from the experimental data with a modified tanh function. As an example, the squares in Fig. 2 are the experimental data points and the red curves are the fitting profiles of electron density and temperature for the Ip = 1.0 MA case. In the simulations, outside the separatrix ψ > 1.0, the initial plasma profiles are assumed to decrease linearly from separatrix to outer boundary. The transport coefficients are calculated based on the initial experimental profiles. The flux-limited parallel thermal transport is implemented in the simulations. The initial profiles for four different currents in the simulations are shown in Fig. 3. The pedestal pressure profiles show a clear trend for these four shots as shown in Fig. 3(a). With the plasma current I_{p} increasing, the pedestal height increases or the scale length decreases, as seen in previous empirical results.^{33–35}

C-Mod shot . | Time (ms) . | B_{t} (T)
. | B_{pol,OM} (T)
. | T_{sep} (eV)
. | I_{p} (MA)
. |
---|---|---|---|---|---|

1100223012 | 1149 | 5.4 | 0.67 | 50.23 | 0.83 |

1100212023 | 1236 | 5.4 | 0.81 | 60.81 | 0.95 |

1100303017 | 1033 | 5.4 | 0.89 | 55.10 | 1.04 |

1160729008 | 0970 | 7.8 | 1.11 | 108.00 | 1.42 |

C-Mod shot . | Time (ms) . | B_{t} (T)
. | B_{pol,OM} (T)
. | T_{sep} (eV)
. | I_{p} (MA)
. |
---|---|---|---|---|---|

1100223012 | 1149 | 5.4 | 0.67 | 50.23 | 0.83 |

1100212023 | 1236 | 5.4 | 0.81 | 60.81 | 0.95 |

1100303017 | 1033 | 5.4 | 0.89 | 55.10 | 1.04 |

1160729008 | 0970 | 7.8 | 1.11 | 108.00 | 1.42 |

One should note that the sheath boundary conditions (SBCs) are implemented at the divertor targets in the simulations as in Ref. 30. The effects of cross-field drifts on the heat flux width are investigated by turning on or off the E × B drift and magnetic drift terms. In general, from the trans-Er module, we can derive (1) the radial electric field, (2) radial and poloidal plasma profiles in the SOL, (3) radial and parallel heat fluxes profiles, and (4) divertor heat flux magnitude and widths. The quantities that can be compared with experiments are the Er profile, the divertor heat flux magnitude, and widths.

## III. Er COMPARISON WITH EXPERIMENT

The onset of H-mode is characterized by a sharp decrease in edge fluctuations and the formation of edge transport barriers (pedestals) in both temperature and density profiles. These barriers lead to significant increases in core confinement as a result of temperature and density profile stiffness. It is commonly accepted that the E × B velocity shear is responsible for the suppression of edge turbulence, which reduces the losses of both energy and particles across the magnetic field. Classical particle motion from E × B and magnetic drifts is believed to be important for understanding tokamak edge/scrape-off-layer (SOL) transport. In this section, we perform transport calculations using the trans-Er module^{30} under the BOUT++ framework. We also do the comparison of the radial electric field with that of the experimental measurement for a C-Mod EDA H-mode discharge (No. 1080321020) with a lower single null divertor configuration. The major parameters for this EDA H-mode discharge are the following: the toroidal magnetic field *B*_{t} = 5.4 T, the separatrix temperature *T*_{e,sep} = 63.5 eV, the separatrix poloidal magnetic field at outer midplane *B*_{p,OMP} = 0.66 T, and the plasma current *I*_{p} = 0.81 MA.

Figure 4 shows the initial density and temperature profiles we used in the simulations. For all C-Mod cases in this paper, ion temperature is assumed to be the same as electron temperature (Ti = Te) due to high collisionality. For the specific discharge considered in this section (No. 1080321020), measurements from pedestal charge exchange recombination spectroscopy (CXRS) confirm the equivalence of Te and Ti within experimental errors. The profiles in Fig. 4 are fits to experimentally measured data inside the separatrix, while in the SOL, they follow a constant absolute value decay from the separatrix. The radial electric field is calculated by solving the potential *ϕ* [Eq. (7)] and vorticity *ϖ* [Eq. (5)] equations for two cases with and without cross-field drifts. In the simulation, *μ*_{i,⊥} is the anomalous perpendicular viscosity and *μ*_{i},$\u2009\u2009\u2225$ is the classical parallel viscosity. The effect of *μ*_{i,⊥}and *μ*_{i,}$\u2009\u2009\u2225$ on the linear growth rate of the peeling-ballooning model^{36} and the steady state solution of the radial electric field^{37} has been discussed in the ELM simulation and Er calculation under the BOUT++ framework. The constant value of *μ*_{i,⊥} = 5 m^{2}/s and *μ*_{i,}$\u2009\u2009\u2225$ = 1.0 × 10^{6} m^{2}/s is selected based on Refs. 35 and 36. The parallel current *J*$\u2009\u2009\u2225$ is calculated by Ohm’s law as Eq. (6). The curvature drift *∇* · (2*P∇* × ** b**) and E × B convection

*V*_{E×B}·

*∇ϖ*are included in the vorticity equation.

The calculated steady state radial electric field (Er) has been compared with that of experimental measurements in discharge 1080321020 using the pedestal CXRS diagnostic.^{38} The simulated Er is shown in Fig. 5 for two cases with and without cross-field drifts. Er is larger in the low field side (LFS) while is smaller in the high field side (HFS) both with and without cross-field drift effects included, which agrees generally with experimental results.^{39} Without magnetic drift effects, inside the separatrix, the E × B flow is balanced by the diamagnetic flow within the pedestal. The Er is mostly determined by ion diamagnetic term $E\u2207P=1ZieNi\u2207\u22a5Pi$ in the pedestal with no net flow. While in the SOL region, a positive Er is formed due to the SBCs as shown by the blue curves. The magnetic drift has a significant effect on Er through their contributions to the net flow due to the charge separation.^{30} In the simulation of discharge No. 1080321020, the radial electric field forms a deep (∼150 kV/m) negative well inside the separatrix due to the steep gradient of density and temperature profiles. The simulation results have similar magnitude and profile shape when compared to experimental measurements shown by the green squares in Fig. 5(b).

## IV. SIMULATIONS FOR DIVERTOR HEAT FLUX WIDTH

### A. The heat flux simulation results

The turbulence (6-field module) simulations under the BOUT++ framework show that in the radial direction, the particle and heat flux are turbulently transported by fluctuations across the separatrix into the SOL and that parallel transport in the SOL then causes particles and heat to flow from the outer midplane toward the divertor targets. From 6-field simulations, we can get the well saturated perpendicular particle and heat flux for both ion and electron at the outer midplane. However, due to the multiscale nature of boundary plasma turbulent transport, the turbulence simulations are able to reach well saturated electron divertor heat fluxes, but not for the ion divertor heat flux.^{27} Therefore, we developed a plasma fluid transport model with all drifts included to arrive at a steady state solution for both electrons and ions more efficiently. Figure 6 shows that both ion and electron parallel heat fluxes reach a steady state at the outer divertor target.

The perpendicular heat flux is determined by the cross-field physics governed by drifts and turbulent transport,

where the first terms represent the turbulent flux and the second terms are the drifts flux. *V*_{exb⊥} and *V*_{dia⊥} are E × B velocity and diamagnetic velocity, respectively. $VD\u22a5\u2207N$ is the effective radial velocity from particle diffusion. Figure 7 shows an example of the perpendicular transport coefficients for the low current case (Ip = 0.8 MA). Inside the separatrix, the perpendicular transport coefficients of *χ*_{i⊥} and *χ*_{e⊥} are prescribed in order to match the experimental profiles as described in Eqs. (11) and (12) of Sec. II and in Ref. 30, and the inferred separatrix value is used for the whole SOL region. Because a core boundary flux $\Gamma i,e|x=0=Ne\chi \u22a5e\u2202Te\u2202x|x=0$ is unknown, in our simulations, the boundary gradients are calculated from experimental profiles and a nominal thermal diffusivity *χ*_{⊥i,e} = 1 m^{2}/s is given. As shown in Table II, the electron heat transport coefficient near the separatrix is comparable to that calculated from turbulence simulations. In the simulation, the perpendicular heat transport coefficients are the same for ions and electrons (*χ*_{i⊥} = *χ*_{e⊥} = *χ*_{⊥}) due to the same ion and electron temperature profiles assumed for C-Mod discharges.

. | Transport
. | Turbulence
. | ||
---|---|---|---|---|

Cases (MA) . | χ_{⊥sep} (m^{2}/s)
. | λ_{q} (mm)
. | χ_{⊥sep} (m^{2}/s)
. | λ_{q} (mm)
. |

0.8 | 0.41 | 1.27 | 0.23 | 1.04 |

0.9 | 0.24 | 1.05 | 0.15 | 0.94 |

1.0 | 0.27 | 0.79 | 0.05 | 0.67 |

1.4 | 0.21 | 1.2 | 0.16 | 1.09 |

. | Transport
. | Turbulence
. | ||
---|---|---|---|---|

Cases (MA) . | χ_{⊥sep} (m^{2}/s)
. | λ_{q} (mm)
. | χ_{⊥sep} (m^{2}/s)
. | λ_{q} (mm)
. |

0.8 | 0.41 | 1.27 | 0.23 | 1.04 |

0.9 | 0.24 | 1.05 | 0.15 | 0.94 |

1.0 | 0.27 | 0.79 | 0.05 | 0.67 |

1.4 | 0.21 | 1.2 | 0.16 | 1.09 |

We impose a sheath boundary condition by assuming that the ion velocity matches the sound speed *C*_{s} at the divertor targets. The sheath heat flux calculations for both ions and electrons are

where *γ*_{e} ≈ 4.8 and *γ*_{i} ≈ 2.5 are electron and ion sheath heat transmission factors, respectively. The quasineutral assumption (*N*_{e} = *Z*_{i}*N*_{i}) is used in the simulations. The steady state solutions of the heat flux for both ions and electrons at the outer divertor target are shown in Table III. The peak of the total heat flux q_{t} = q_{i} + q_{e} is close to that measured on C-Mod. The heat flux increases with the current. The electron heat fluxes are dominant at the lower values of B_{p} (or I_{p}). However, the ion heat flux becomes larger with the current increase, approaching to the electron heat flux for the high current case (I_{p} = 1.4 MA).

. | . | Experiment . | BOUT++ transport (MW/m^{2})
. | ||||
---|---|---|---|---|---|---|---|

I_{p} (MA)
. | B_{t} (T)
. | q_{t} (MW/m^{2})
. | λ_{q} (mm)
. | q_{t} (MW/m^{2})
. | λ_{q} (mm)
. | q_{i} (MW/m^{2})
. | q_{e} (MW/m^{2})
. |

0.8 | 5.4 | 163 | 0.76 | 184 | 1.27 | 51 | 133 |

0.9 | 5.4 | 253 | 0.63 | 233 | 1.05 | 108 | 125 |

1.0 | 5.4 | 298 | 0.97 | 294.7 | 0.79 | 111.7 | 183 |

1.4 | 7.8 | 800 | 0.6 | 880.3 | 1.2 | 434 | 456 |

. | . | Experiment . | BOUT++ transport (MW/m^{2})
. | ||||
---|---|---|---|---|---|---|---|

I_{p} (MA)
. | B_{t} (T)
. | q_{t} (MW/m^{2})
. | λ_{q} (mm)
. | q_{t} (MW/m^{2})
. | λ_{q} (mm)
. | q_{i} (MW/m^{2})
. | q_{e} (MW/m^{2})
. |

0.8 | 5.4 | 163 | 0.76 | 184 | 1.27 | 51 | 133 |

0.9 | 5.4 | 253 | 0.63 | 233 | 1.05 | 108 | 125 |

1.0 | 5.4 | 298 | 0.97 | 294.7 | 0.79 | 111.7 | 183 |

1.4 | 7.8 | 800 | 0.6 | 880.3 | 1.2 | 434 | 456 |

The power crossing the separatrix is the integral of the perpendicular heat flux over the separatrix surface and the power at the divertor is the integral of the parallel heat flux over the divertor targets. The last two columns in Table IV show the simulation results of the perpendicular power across the separatrix Q_{sep} and the parallel power at the divertor targets Q_{target.} As the current increases, the power across the separatrix increases, which agrees with the experiments as shown in the second column. The simulated power at the divertor targets is also consistent with the simulated power across the separatrix. The value of simulated power at the divertor targets is a little bit smaller than that across the separatrix due to the energy lost to the chamber wall (<11.5%). The total simulated divertor power is ∼2× the experimentally inferred power on the outer divertor. The inner divertor power is not measured for these discharges. Radiative losses in the divertor and SOL are neither measured nor simulated here.

. | Experiment: power across . | Experiment: power on . | BOUT++ transport: power . | BOUT++ transport: power . |
---|---|---|---|---|

w/drifts (MA) . | separatrix estimated Q_{sol} (MW)
. | divertor Q_{target} (MW)
. | across separatrix Q_{sep} (MW)
. | on divertors Q_{target} (MW)
. |

I_{p} = 0.8 | 1.38 | 0.722 | 1.35 | 1.31 |

I_{p} = 0.9 | 2.27 | 1.163 | 1.56 | 1.38 |

I_{p} = 1.0 | 2.25 | … | 2.09 | 1.94 |

I_{p} = 1.4 | 2.93 | 1.77 | 3.09 | 3.08 |

. | Experiment: power across . | Experiment: power on . | BOUT++ transport: power . | BOUT++ transport: power . |
---|---|---|---|---|

w/drifts (MA) . | separatrix estimated Q_{sol} (MW)
. | divertor Q_{target} (MW)
. | across separatrix Q_{sep} (MW)
. | on divertors Q_{target} (MW)
. |

I_{p} = 0.8 | 1.38 | 0.722 | 1.35 | 1.31 |

I_{p} = 0.9 | 2.27 | 1.163 | 1.56 | 1.38 |

I_{p} = 1.0 | 2.25 | … | 2.09 | 1.94 |

I_{p} = 1.4 | 2.93 | 1.77 | 3.09 | 3.08 |

In order to compare the BOUT++ transport simulation results of *λ*_{q} with experimental results, we use the same fitting method as that used to populate the experimental heat flux database,^{5} as follows:

and $s\xaf=(Rsep\u2212R)\u2009\u22c5\u2009fx$, where R is the major radius, *R*_{sep} is the separatrix position at the outer midplane, and *f*_{x} is the effective flux expansion. The quantity *S* is the divertor power spreading parameter, which depends on local divertor plasma parameters and geometry, *q*_{BG} is the background heat flux, and *q*_{0} is the peak heat flux. The fitting function is a convolution from an exponential decay of the outer midplane parallel heat flux profile with a Gaussian function of width *S*. This formula will be used to calculate the heat flux width *λ*_{q} from experimental data and simulated data.

Direct comparisons of simulated heat flux profiles with experimental profiles are shown in Fig. 8. Simulations of lower current discharges show a sharp drop in q_{||} in the private flux region (PFR) as shown by black squares in Fig. 8(a). This does not generally agree with experimental measurements of the divertor parallel heat flux as shown by gray squares in Fig. 8(a), which were made with infrared imaging (IR),^{40} which has an estimated spatial resolution in this coordinate system of approximately 1 mm. In order to compare the simulated heat flux profile with the measurement, a convolution with an instrument function of 1 mm width is applied to the simulated data. Figure 8(b) shows the overlay of the same IR data with synthetic data points (green squares) by applying the instrument function convolution to the simulation data. This shows a better comparison between the IR measured heat flux profile and synthetic reconstruction based on the BOUT++ simulation. Figure 8(c) compares BOUT++ transport simulation and experimental data for the high current case (I_{p} = 1.4 MA). The black squares in Fig. 8(c) are the calculated parallel heat flux from the BOUT++ simulation, and the light cyan squares and gray squares in Fig. 8(c) are the measured outer divertor parallel heat flux, mapped to the outer midplane. The experimental data come from a suite of divertor diagnostics with high spatial resolution,^{13} such as the surface thermocouples^{41} for the light cyan squares and Langmuir probes (LP)^{42} for gray squares. There is the natural scatter that is always present in experimental data. There are clear outliers in the profile and even negative values of the heat flux due to both real fluctuations, instrumental reporting error, and uncertainty in the magnetics reconstruction. Here, a radial shift of 1 mm is applied to bring the peak of the profile close to the strike point R-R_{sep} = 0. The red curve and blue curve in Fig. 8(c) are fits to the profiles using Eq. (19). The sharp falloff in q_{||} for R-R_{sep} < 0 is reproduced well by the simulation, while the simulated heat flux channel for R-R_{sep} > 0 is wider than that in the experiment.

The detailed comparison of simulated data with the experiment for 4 different currents is listed in Table III. Here, we assume that Eq. (19) effectively captures the instrumental resolution of the diagnostics in the S parameter and that the fitted *λ*_{q} can be compared with that fitted to BOUT++ results. Here, BOUT++ transport simulation results have an Eich heat flux width *λ*_{q} systematically larger than the Eich scaling would predict. The values are 0.8–1.8 times that of the experimental results, which exhibit significant scatter about the projection of the Eich scaling. The electron heat flux width *λ*_{q,e} is larger than that of the ions *λ*_{q,i} as shown in Fig. 9(a). The same trend is found for both ion and electron heat flux width, *λ*_{q,i} and *λ*_{q,e}, respectively. The transport simulation results for the heat flux width are consistent with the turbulence code calculation and follow the Eich scaling trend for low current cases. It should be noted that for the high current case, the calculated result does not follow the inverse I_{p} scaling but is consistent with the Goldston HD model as shown by the blue square in Fig. 9(b), possibly due to the high separatrix temperature. Further discussion will be given in the next paragraph.

To explore why the simulated heat flux width for the high current case does not continue the trend of Eich scaling, we make a comparison of the simulated divertor heat flux width to that of the Goldston HD model.^{17} The model assumes that the magnetic (∇B and curvature B) drifts balance against the near-sonic parallel flow in setting the ion heat flux width and that cross-field transport at the separatrix is dominated by classical drifts.^{17,43} The heuristic drift-based (HD) model gives a width *λ*_{q} ∼ *qρ*_{s}, where *ρ*_{s} is the ion gyroradius and *q* is the safety factor. Goldston’s formula for the heat flux width follows

For the C-Mod deuterium plasma, we approximate a = 0.22 m, R = 0.68 m, $Z\xaf=1$, and $\u0100=2$. The separatrix poloidal magnetic field at the outer midplane B_{pol,MP} (B_{p}) and the separatrix temperature T_{sep} are listed in Table I for I_{p} = 0.8–1.4 MA, respectively. The heat flux widths calculated by Goldston’s HD model are added in Fig. 9(b) as blue squares. The error bars are calculated from the uncertainty of the separatrix temperature, considering the uncertainty of the separatrix position given by the magnetic equilibrium reconstruction, *dψ* ∼ 0.01. The short vertical blue lines are the error bars, as shown in Fig. 9(b). The simulation results calculated by the transport code under the BOUT++ framework show a similar trend of the heat flux width to that of Goldston’s HD model. In the simulation, the averaged parallel flow from the outer midplane to the divertor target is calculated to be around 0.3–0.4*C*s, which is smaller compared with an estimate of 0.5*C*s in Goldston’s model. Compared to the low current case, the simulated heat flux width becomes larger at the high current. That is because the width *λ*_{q} not only depends on the poloidal magnetic field at the outer midplane B_{pol,MP} (B_{p}) but also is very sensitive to the separatrix temperature. In order to estimate the dependence of the separatrix temperature, we do a separatrix temperature scan for a high current case (1.4 MA). The results are shown in Fig. 10(a). From the transport simulations as shown in Fig. 10(b), it can be found that the heat flux width increases with the separatrix temperature, which is consistent with the Goldston formula in Eq. (16). In Fig. 9(a), the short vertical green line for B_{pol,mp} = 1.1 T is the error bar from the scanning of separatrix temperature for a high current case, as shown in Fig. 10(b). When the separatrix temperature is low, the heat flux width is close to the experimental measurement.

### B. The effects of drifts on heat flux width

Classical particle motion from E × B and magnetic drifts is important for understanding tokamak edge/scrape-off-layer (SOL) transport even in the presence of turbulent transport, and this may influence the overall behavior of the heat flux at the divertor targets.^{43} The asymmetry of the plasma density and temperature profiles changes between the inner and outer divertor regions, and the profiles are poloidally nonuniform with the drift effects.^{44–46} In this section, the impact of magnetic drift and E × B drift on the divertor heat flux width will be studied by turning on or off the options of magnetic drift and E × B drift. Based on fits to Eq. (19), the heat flux width is calculated for the three cases: without any drifts, with the magnetic drift only, and with both magnetic and E × B drifts.

The data points of the heat flux width for four C-Mod EDA H-mode discharges with different drifts are shown in Fig. 11. We find that the trend of the heat flux width with the poloidal field is similar for the cases with or without the drifts. The magnetic drift has a significant influence on the heat flux width. With the magnetic drift, the width is reduced around 2–3 times because the ∇B and curvature B drifts can cause a large parallel flow and hence reduce the SOL residence time. The values calculated are similar to Goldston’s HD model, which is also without E × B drift. Reiser and Eich conducted the studies of the drift-based scrape-off particle width in X-point geometry to confirm the HD model for the averaged density decay length.^{43} In order to know the effect of the E × B drift, the heat flux width is then calculated with both magnetic and E × B drifts as shown by the green stars in Fig. 11. The E × B drift can further increase the parallel transport and decreases the heat flux width. From the simulation results, the E × B drift decreases the heat flux width by 10%–25%, which brings the calculated heat flux width closer to the experimental values.

### C. The competition between drifts and turbulent transport for C-Mod

We expect that the heat flux width is determined by the competition between the cross-field and parallel transport.^{47} The cross-field transport is governed by drifts and turbulent transport. The effect of drifts with increasing particle turbulent diffusivity has been investigated in Ref. 48. In this section, the competition between drifts and turbulent heat transport coefficient is investigated. We perform a scan of the turbulent heat transport coefficients for the 0.8 MA C-Mod discharge. The perpendicular heat transport coefficient (*χ*_{⊥}) in the scrape-off-layer (SOL) is set constant in radius and then scanned in a series of simulations. Inside the separatrix, we use the same *χ*_{⊥} interpreted from experimental profiles. The SOL *χ*_{⊥} varies by 160× from 0.05 m^{2}/s to 8 m^{2}/s, as shown in Fig. 12(a). Figure 12(b) shows the transport simulation results of heat flux width *λ*_{q} for different *χ*_{⊥} with drifts. The scan of turbulence diffusivity *χ*_{⊥} identifies two distinct regimes. When *χ*_{⊥} < 0.5 m^{2}/s, the heat flux width *λ*_{q} is a constant, the drifts dominate cross-field transport and heat flux width is insensitive to the turbulent transport *χ*_{⊥} in this regime, which is consistent with the Goldston HD-model as shown by blue star in Fig. 12(b). The Goldston HD model yields a lower limit of the width *λ*_{q}. When *χ*_{⊥} > 0.5 m^{2}/s, heat flux width *λ*_{q} increases with *χ*_{⊥}, corresponding to turbulence dominating cross-field transport. The $\chi \u22a5c=0.5\u2009m2/s$ is the critical thermal diffusivity value for the transition from drift dominant regime to turbulent dominant regime. The green star shows the interpreted result from the experimental profile using formula (12) with *χ*_{⊥} = 0.41 m^{2}/s and *λ*_{q} = 1.27 mm. The red bullet in Fig. 12(b) is the simulation result from the turbulence code with *χ*_{⊥} = 0.23 m^{2}/s and *λ*_{q} = 1.04 mm, which is comparable to transport code results. A similar scan of thermal diffusivity has been done for the high current case (1.4 MA), and a smaller critical $\chi \u22a5c=0.3\u2009m2/s$ is found due to the decrease in drift effect. Both the thermal diffusivity calculated by turbulence code and that interpreted from the experimental data are smaller than the critical thermal diffusivity.^{15} The simulation results suggest in C-Mod EDA H-mode discharges, drifts, and turbulent transport compete to determine the characteristic heat flux width, and those drifts set a lower bound on *λ*_{q}. However, uncertainty remains, since the heat flux width calculated by all models systematically overpredicts the measured values as shown by the gray dashed-dotted horizontal line. It should be mentioned that the predictions for future machines, such as CFETR and ITER, show a lower critical thermal diffusivity $\chi \u22a5c$ and a higher width *λ*_{q}. ITER and CFETR are predicted to be in a turbulence dominant regime.^{15,16}

## V. SUMMARY

The transport code (trans-Er) with all the cross-field drifts under the BOUT++ framework is used to calculate the steady state radial electric field (Er) and to study the physics setting heat flux width. Transport coefficients are calculated from the experimental profiles inside the separatrix, then extending to the SOL. Simulation results show that Er is larger in the LFS than that in the HFS. The steady state Er has been compared favorably to experimental measurements from a C-Mod enhanced D_{α} (EDA) H-mode discharge.

In order to understand the relative role of drifts vs turbulent transport in setting the heat flux width, a set of four C-Mod EDA H-mode discharges with a lower single null divertor configuration are simulated. The simulation results yield similar *λ*_{q} to experimental measurements within a factor of 0.8–1.8 and show a trend similar to Goldston’s HD model. In simulations, the power across the separatrix increases with current for these four discharges. The power deposited on the divertor target is consistent with the power across the separatrix due to the low SOL radiation loss. The magnetic drift has a significant influence on the heat-flux width, while the E × B drift decreases the heat flux width by 10%–25%, which gives improved agreement with experiment relative to Goldston’s model. However, there are indications from modeling a separate EDA H-mode discharge that the simulation with full magnetic drifts may overestimate Er. In simulations of C-Mod discharges, the heat flux width is sensitive to the temperature near the separatrix. A turbulence diffusivity scan (*χ*_{⊥}) identifies two distinct regimes: a drift dominant regime when *χ*_{⊥} is small and a turbulence dominant regime when *χ*_{⊥} is large. Magnetic drifts give a lower limit of the width *λ*_{q}, similar to Goldston. This exercise suggests that drifts and turbulent transport compete to determine the heat flux width in C-Mod EDA H-mode discharges.

## ACKNOWLEDGMENTS

The authors would like to thank T. D. Rognlien, B. Chen, J. G. Chen, Zeyu Li, T. Y. Xia, Y. M. Wang, and T. F. Tang for useful discussions and A. Q. Kuang and D. Brunner for providing heat flux data for discharge (Grant No. 1160729008). This work was supported by the National Key R&D Program of China (Grant Nos. 2017YFE0301206 and 2017YFE0300402) and the National Natural Science Foundation of China under Grant Nos. 11675037 and 11575039. This work was also performed under the U.S. Department of Energy by Lawrence Livermore National Laboratory under Contract No. DE-AC52-07NA27344. Further support comes from US DoE (Award No. DE-SC0014264). The LLNL work was supported by the U.S. DoE, Office of Science, Office of Fusion Energy Sciences.