Dynamics of Ion Temperature Gradient Turbulence and Transport with a Static Magnetic Island
Olivier Izacard, Christopher Holland, Spencer D. James, Dylan P. Brennan
DDynamics of Ion Temperature Gradient Turbulence and Transport with a StaticMagnetic Island
O. Izacard, ∗ C. Holland, S. D. James, and D. P. Brennan Lawrence Livermore National Laboratory, 7000 East Avenue, Livermore, California 94550, USA University of California - San Diego, La Jolla, California 92093-0417, USA University of Tulsa, Tulsa, Oklahoma 74104, USA Princeton University, Princeton, New Jersey 08544, USA (Dated: September 19, 2018)Understanding the interaction mechanisms between large-scale magnetohydrodynamic instabili-ties and small-scale drift-wave microturbulence is essential for predicting and optimizing the per-formance of magnetic confinement based fusion energy experiments. We report progress on under-standing these interactions using both analytic theory and numerical simulations performed withthe BOUT++ [B. Dudson et al., Comput. Phys. Comm. , 1467 (2009)] framework. Thiswork focuses upon the dynamics of the ion temperature gradient instability in the presence of abackground static magnetic island, using a weakly electromagnetic two-dimensional five-field fluidmodel. It is found that the island width must exceed a threshold size (comparable to the turbulentcorrelation length in the no-island limit) to significantly impact the turbulence dynamics, with theprimary impact being an increase in turbulent fluctuation and heat flux amplitudes. The turbulentradial ion energy flux is shown to localize near the X-point, but does so asymmetrically in thepoloidal dimension. An effective turbulent resistivity which acts upon the island outer layer is alsocalculated, and shown to always be significantly (10x - 100x) greater than the collisional resistivityused in the simulations.
I. INTRODUCTION
One of the most crucial factors in determining the sizeof magnetic-confinement based fusion energy (MFE) re-actors is the level of plasma confinement achieved [1]. InMFE-relevant plasmas, the dominant physical processeswhich determine this confinement are a variety of insta-bilities driven by the inherent temperature, density, andcurrent gradients of the confined plasma. Two of themost important such instabilities in tokamaks are the iontemperature gradient (ITG) mode [2] and tearing modes[3]. As its name suggests, the ITG instability is driven bythe radial gradient (i.e. across nested magnetic flux sur-faces) of the equilibrium ion temperature profile, whichmust exceed a critical value (the so-called critical gradi-ent) for onset of the instability. Once the critical gra-dient is exceeded, the ITG instability leads to a broadspectrum of ion gyroradius ( ρ i ) scale fluctuations, whosecollective nonlinear behavior drives rapid cross-field par-ticle, energy, and momentum fluxes that greatly exceedcollisional transport levels [4], and often determine theoverall confinement level achieved. Since these fluctua-tions generally have correlation lengths on the order of1 − ρ i , which is much smaller than the minor radius a in MFE-relevant devices, the saturated fluctuations andtheir dynamics are often referred to as microturbulence.Magnetic islands can appear through magnetohydro-dynamic (MHD) instabilities [5, 6], such as the tear-ing mode (TM) driven by equilibrium gradients, or beexternally imposed with resonant magnetic perturba- ∗ Electronic address: [email protected] tions (RMP). These islands can reduce the confinementachieved, similar to the ITG instability, or even termi-nate confinement by leading to a disruption. It is theo-retically predicted and experimentally observed that thepresence of a magnetic island flattens the electron and iontemperature and density. This flattening occurs both in-side the island through rapid parallel equilibration, boththrough parallel sound waves [7] and rapid parallel heatconduction [8], and in time averaged profiles by enhanc-ing the effective cross-field transport. The classic pictureis that the radial heat flux increases at the X-point ofthe magnetic island due to the fast parallel dynamics inthe highly anisotropic 3D island structure [8]. Indeed,because of the direct connection of the magnetic fluxsurface at the separatrix of the island, all density andtemperature structures on one side of the island can berapidly transported to the other side of the island onshorter timescales than the timescale of the backgroundcross-field transport.Given the importance of both classes of instabilities indetermining the level of confinement achieved in magnet-ically confined plasmas, it is important to develop mod-els of their dynamics which can be used to accuratelyinterpret existing observations and confidently design fu-ture experiments and reactors. While great progress hasbeen made in recent years in understanding the nonlin-ear dynamics and scalings of each class of instabilities inisolation (often through use of massively parallel compu-tation), these instabilities frequently coexist in currenthigh-performance tokamak discharges, and are expectedto continue to do so in many future reactor scenarios.We must therefore also understand how these instabili-ties couple linearly and nonlinearly if we are to be ableto predict the overall confinement and performance of a r X i v : . [ phy s i c s . p l a s m - ph ] S e p such discharges. For example, it is now well-known thatmicroturbulence (and ITG microturbulence in particu-lar) generally saturates via the nonlinear generation ofaxisymmetric zonal flows [9] which in turn shear apartthe ITG eddies, forming a self-regulating system of tur-bulence and zonal flows. How this process is altered bythe presence of magnetic islands (which generally havetoroidal mode numbers n = 1 − m = 2 −
5) and their associated flow fields re-mains an open question. Alternatively, the saturation ofthe neoclassical tearing mode depends sensitively on thelevel of profile flattening that occurs in and near the is-land [10, 11]. As this flattening depends sensitively onthe ratio of parallel and perpendicular (microturbulence-dominated) transport, it is clear that the saturated islandwidth will depend self-consistently on the level of micro-turbulence present (which also depends on the amountof profile flattening that occurs).These two instabilities (ITG and TM) are well studiedin the literature, with many different research teams find-ing a wide array of complex nonlinear couplings betweenmicroturbulence and tearing modes. For instance, McDe-vitt and Diamond [12] demonstrated that drift-wave mi-croturbulence could nonlinearly excite an island via theReynolds stress in a fashion analogous to the formation ofzonal flows, while Sen et al [13] used a similar approach todemonstrate that electron temperature gradient (ETG)modes could nonlinearly damp the island via an effectiveturbulent resistivity. More recently, a wide range of re-searchers have utilized both fluid [14–26] and gyrokinetic[27–37] approaches to investigate couplings of various mi-croturbulence instabilities to either statically imposed is-lands or dynamically evolving tearing modes. Althoughthere is no simple and clear consensus yet, it is clear thatthe presence of the island can lead to a poloidal localiza-tion of linear microturbulence eigenmodes, and that themicroturbulence can both damp and destabilize islandsdepending on the type of microturbulence instability, andthe width and rotation speed of the island. The impactsof the turbulence on the island evolution are generallyrepresented by inclusion of turbulent polarization currentterms in generalized Rutherford equations [18, 38].Building upon these results, we present in this paperthe findings of a study examining the responses of weaklyelectromagnetic ITG turbulence to statically imposed is-lands of varying width, performed using a simple fluidmodel in slab geometry. We find that above a criticalisland width, the magnitude of both the microturbulencefluctuations and associated radial ion heat flux Q i in-crease approximately linearly with island width, bothnear and far from the island resonant surface, for val-ues of the the ion temperature gradient both close to,and well above the linear critical gradient. The criticalisland width is found to be approximately equal to theradial correlation length of the turbulence in the absenceof an island. Perhaps not surprisingly, for islands smallerthan the critical size, the microturbulence is found to beeffectively insensitive to the presence of the island. We also observe the turbulent Q i to localize near the islandX-point as expected, but somewhat surprisingly, this lo-calization is not seen to be symmetric in the poloidaldirection. To begin assessing how the turbulence may beexpected to back-react upon the island, we present calcu-lations of an effective turbulent resistivity which operateson the “outer solution” of the island structure. We showthat for the parameters considered, while this resistivityis always much larger than the collisional value, it max-imizes for small island size and decreases rapidly withincreasing island width.The paper is organized as follows. Section II details thefluid model used to study the effects of a static magneticisland on the ITG microturbulence, and includes detailsof the numerical implementation in BOUT++ (II B) andverification of the implementation using analytic growthrate calculations (II C) and a global energy balance anal-ysis suitable for nonlinear simulations (II D). Sec.III dis-cusses the results of numerical simulations quantifyingthe impact of statically imposed islands of varying widthon ITG turbulence, and Sec.IV discusses initial work ex-ploring the expected feedback of the turbulence on theisland. Conclusions and future research directions arediscussed in Sec.V II. DETAILS OF THE FIVE-FIELD FLUIDMODEL
This section describes the 2D electromagnetic five-fieldmodel used in this study, and its implementation in theBOUT++ framework [39]. Two verification methods arepresented to test the implementation. The first is a com-parison of calculated linear growth rates against analyticderivations, detailed in Sec. II C. The second is a globalenergy balance analysis approach, detailed in Sec. II D.
A. Model overview
In order to investigate the self-consistent couplings ofITG turbulence and tearing modes in a simple, tractableform, we employ a simple five-field two-fluid model [21,40, 41] in a Cartesian slab geometry [21, 41–44] simi-lar to that used in previous studies, considered togetheronly in Ishizawa et al. [21, 41]. The simulation radialand poloidal domains span the range − a ≤ x ≤ a and − b/ ≤ y < b/ x = ± a and periodicboundary conditions at y = ± b/
2. Using the standarddefinitions of (cid:126)E = − (cid:126) ∇ φ and (cid:126)B = B ( x )ˆ z + (cid:126) ∇ ψ × ˆ z , themodel describes the evolution of fluctuations in the elec-tron density ˜ n = δn/ ¯ n , vorticity ˜Ω = ∇ ⊥ ˜ φ , ion parallelvelocity ˜ v i,z = δv i,z /c s , ion temperature ˜ T i = δT i / ¯ T , andmagnetic flux ˜ ψ = (1 /β )( c s /c ) eδψ/ ¯ T , where ˜ φ = eδφ/ ¯ T is the normalized electrostatic potential fluctuation, andthe normalized axial current is given by j z = −∇ ⊥ ψ .Here, ¯ n and ¯ T are constant normalizing factors takento be equal to the values of the equilibrium density andion temperature profiles at x = 0; c s = (cid:112) ¯ T /M i and ρ s = c s / Ω ci , with Ω ci = eB (0) /M i c . In this model,we can chose a form of the equilibrium magnetic fluxthat allows us to specify both an equilibrium poloidalmagnetic field B y, ( x ) = − ∂ x ψ and a static magneticisland B x,island ( y ) = ∂ y ψ island . These relations comefrom B = ∇ ψ ( x, y ) × ˆ z and the current is given by j z = ˆ z · ∇ × B = −∇ ψ . To further simplify the problem,in the rest of this paper we restrict our attention to thecase for which the equilibrium density n is constant inradius and the equilibrium electrostatic potential φ andparallel ion flow v i,z are equal to zero. The equilibriumion temperature profile is taken to have a simple lineardependence T i ( x ) = ¯ T (1 − x/L T i ), where L T i is constantacross the radial domain. The equilibrium magnetic fluxis taken to have the form ψ eq ( x, y ) = ψ ( x )+ ψ island ( y ) = − x / (2 L S ) + ψ i cos(2 πy/b ), with ψ i = k W island / (4 a )and the associated current is j z ( x, y ) = −∇ ⊥ ψ eq where ∇ ⊥ = ∂ xx + ∂ yy . The constituent equations of the modelare the electron density, the perpendicular momentumassociated with the quasi-neutrality, the parallel momen-tum, the ion temperature conservation, and Ohm’s lawas follows: d t ˜ n = v (cid:63)B ∂ y (cid:16) ˜ φ − ˜ n (cid:17) + ∇ (cid:107) (cid:0) ˜ j z − ˜ v i,z (cid:1) + ˜ ∇ (cid:107) j z − ν s (cid:104) ˜ n (cid:105) y + D ⊥ ∇ ⊥ ˜ n (1) d t ˜Ω = ∇ (cid:107) ˜ j z + ˜ ∇ (cid:107) j z (2) − v (cid:63)B ∂ y { (1 + τ )˜ n + τ ˜ T i } + τ v (cid:63)T i ∂ y ˜Ω+ τ ∇ ⊥ · (cid:104) ∇ ⊥ ˜ φ, ˜ n + ˜ T i (cid:105) + D ⊥ ∇ ⊥ ˜Ω (3) d t ˜ v i,z = −∇ (cid:107) (cid:16) (1 + τ )˜ n + τ ˜ T i (cid:17) − ν s (cid:104) ˜ v i,z (cid:105) y + D ⊥ ∇ ⊥ ˜ v i,z (4) d t ˜ T i = − v (cid:63)T i ∂ y ˜ φ − (Γ − ∇ (cid:107) ˜ v i,z + λ v (cid:63)B (Γ − ∂ y (cid:16) ˜ φ/τ + ˜ n + ˜ T i (cid:17) − ν s (cid:68) ˜ T i (cid:69) y + D ⊥ ∇ ⊥ ˜ T i (5) β∂ t ˜ ψ = −∇ (cid:107) (cid:16) ˜ φ − ˜ n (cid:17) − η ˜ j z (6)where v (cid:63)B = a/L B and v (cid:63)T i = a/L T i are the driven im-posed gradients of the magnetic field and the ion tem-perature. In these equations x and y are normalized to ρ s and t to a/c s . Using the standard Poisson bracketnotation { f, g } = ∂f /∂x ∂g/∂y − ∂f /∂y ∂g/∂x , we writethe total time derivative as df /dt = ∂f /∂t + (cid:110) ˜ φ, f (cid:111) toinclude the E × B convection term, and the total paral-lel derivative ∇ (cid:107) f = − β (cid:110) ψ + ψ island + ˜ ψ, f (cid:111) includesderivatives along the total magnetic field (the equilib-rium field, any imposed island, and the magnetic fieldfluctuations). We also note that the use of the notation˜ ∇ (cid:107) f = − β (cid:110) ˜ ψ, f (cid:111) is associated only with the magneticfluctuations part. The coefficient Γ = 5 / τ is the ratio between the ion andelectron reference temperatures. The coefficient λ , equalto 0 or 1, is used to investigate the impact of differentclosure models in the ˜ T i equation. We have included λ in order to investigate the effects of the absence in somemodels found in the literature of the divergence of theion E × B and diamagnetic velocity terms. This term,describing finite Larmor-radius (FLR) effects, is retained(or not) in the fluid closure model used for the ion tem-perature equation. Moreover, this model contains otherFLR effects (i.e., τ ∇ ⊥ · (cid:104) ∇ ⊥ ˜ φ, ˜ n + ˜ T i (cid:105) ) [45–47] whichare not always considered in the literature but includedhere in order to take into account the ion temperature ef-fects which come from the convection of the ion pressureby the E × B drift. This term appears in the vortic-ity equation due to the divergence of the ion polariza-tion current. When the inhomogeneity of the magneticfield is kept, the models including this FLR term suchas in Refs. [40, 45, 46, 48–51] are not Hamiltonian [52]in the ideal limit ( D = η = ν s = 0). However, theideal version of our model is very close to a Hamilto-nian structure because the terms required to recover theHamiltonian property involve, for example, higher orderspatial derivatives of the magnetic field which are can-celed by the assumption of a linear background magneticgradient, or higher order spatial derivatives of the vor-ticity which are assume negligible. Moreover, the energyconservation of our simulations is systematically verifiedas well as the error between the numerical and theoret-ical computation of the energy balance. We found thatthe inclusion of the FLR term (with λ = 1) changes thegrowth rate less than 10% in the linear regime, but thisFLR term in the ion temperature equation and the onein the vorticity equation are kept in the following be-cause they contribute to the numerical stability of oursimulations for high k y modes.This set of equations becomes an electrostatic two-dimensional four-field fluid model with gradient driventurbulence when β = 0. Moreover, when β is relativelysmall (i.e., ∼ − ) the differences between the five-fieldand the four-field electrostatic model are not significant.The five-field model is preferred for consistency with fu-ture work on dynamic islands. The additional dissipativeviscosity ν s is included to prevent quasilinear relaxationand maintain input background profiles such as the back-ground ion temperature gradient. Without this forcing,the ion temperature gradient cannot be maintained in thepresence of a magnetic island and the ITG which drivesthe turbulence becomes stable. This choice is needed tomaintain on average the same level of turbulence due tothe ITG. We show below that this assumption is stillcompatible with the flattening effect of the island on thetemperature and density profiles. The numerical dissipa-tion terms proportional to D ⊥ ∇ ⊥ utilize small values of D ⊥ to damp grid-scale fluctuations without significantlyimpacting the large-scale dynamics of interest. Note thatrelative to previous tearing mode studies which have uti-lized anomalous transport coefficients, these simulationsself-consistently determine those coefficients via the tur-bulent fluxes e.g. Q i = (cid:68) ˜ T i ˜ v E × B (cid:69) . B. Numerical implementation
The BOUT++ [39] framework is used for the numer-ical computation of the five-field fluid model given byEqs. (1)-(6). Even though BOUT++ can deal with ahelical axisymmetric magnetic field, as a first step ourwork focuses on a slab geometry choosing the box ratioof the simulations in order to fix a specific regime of theinstabilities. The use of BOUT++ is an asset since thedynamical equations are written in a small script withobject oriented functions. Also, all multi-dimensionalscans and post-treatment analyses have been simplifiedby the use of the OMFIT [53] framework and has mo-tivated the creation and development of the BOUT++interface module in the OMFIT framework as explainedin Ref. [54]. As a summary, OMFIT is a software (withGraphical User Interface) for integrated studies whichcan quickly connect many codes and experimental datain a workflow and which can perform data analysis. Forinstance, OMFIT was used to handle remote execution,data transfer, and analysis of the simulations presentedhere, which were performed mainly on the TSCC com-puter at the San Diego Supercomputer Center.
C. Linear dispersion relation
In order to verify the accuracy of our numerical sim-ulations against the target fluid model, the first step isto compare the computed linear growth rates against an-alytic calculations in the limit of no background staticisland and no static magnetic shear (i.e., ψ eq = 0).The linearization of the fluid equations by ˜ f ( x, y, z, t ) =˜ f exp( ik x x + ik y y − iωt ) yields − ω ˜ n = v (cid:63)B k y (cid:16) ˜ φ − ˜ n (cid:17) , (7) ωk ⊥ ˜ φ = − v (cid:63)B k y ((1 + τ )˜ n + τ ˜ T i ) − τ v (cid:63)T i k y k ⊥ ˜ φ, (8) − ω ˜ v i,z = 0 , (9) − ω ˜ T i = − v (cid:63)T i k y ˜ φ + v (cid:63)B k y λ (Γ − (cid:16) ˜ φ/τ + ˜ n + ˜ T i (cid:17) , (10) − ωβ ˜ ψ = iηk ⊥ ˜ ψ. (11)After some analytic computation the dispersion relationis reduced to D ( ω ) = A ω + A ω + A ω + A ω + A ω + A , (12) with A = βk ⊥ , (13) A = k ⊥ ( iηk ⊥ + βC BT i ) , (14) A = − β Ω B ((1 + τ + Γ )Ω B − τ Ω T i ) + k ⊥ iηC BT i − βk ⊥ Ω B (Γ Ω B − τ (Γ − T i ) , (15) A = − iηk ⊥ ((1 + τ + Γ )Ω B − τ Ω T i Ω B − k ⊥ C BT i ) − βτ Ω T i Ω B (Γ k ⊥ + 1) , (16) A = k ⊥ iη (Γ k ⊥ + 1) τ Ω T i Ω B , (17) A = 0 , (18)where Ω B = v (cid:63)B k y , Ω T i = v (cid:63)T i k y , Γ = λ (Γ −
1) and C BT i = (Γ − B + τ Ω T i .This linear dispersion relation is obviously driven bythe imposed gradients of the ion temperature v (cid:63)T i andmagnetic field v (cid:63)B . If the latter are dropped, the disper-sion relation becomes D ( ω ) = k ⊥ ω ( βω + iηk ⊥ ) and theroots are found to be stable (imaginary part of ω , i.e.the growth rate, is 0 and the propagation frequency isequal to 0 or − k ⊥ η/β ). A scan of the ion temperatureand magnetic field gradients has shown that the lineargrowth rate is unstable beyond a threshold dependent onboth gradients. We typically perform our simulations inthe unstable regime with v (cid:63)B = 0 . v (cid:63)T i = 0 .
3, highlyrelevant to experimental values.A comparison of the computed linear growth ratesagainst the analytic linear dispersion relation (Eq. (12))is shown in Fig. 1, demonstrating excellent agreement inboth growth rate γ and real frequency ω for this basecase; similar agreement is obtained at other parameters.Error estimates on the simulation results are calculatedvia the standard deviation of the instantaneous frequencyand growth rate of the time-averaging window used (gen-erally on the order of 100’s of a/c s ). For wave num-bers with sufficiently small or zero growth rate, a well-converged value is often not found, leading to the largeerror bars shown at the highest values of k y plotted andat k y = 0. −1.0 −0.5 0.0 0.5 1.0 k y γ (a) Growth rate γ −1.0 −0.5 0.0 0.5 1.0 k y −0.2−0.10.00.10.2 ω (b) Propagation ω FIG. 1: Verification of the linear dispersion relation compu-tation between the theory (blue curves) given by the Eq. (12)and the BOUT++ simulation (dashed-black curves) of theEqs. (1)-(6). The error bars (x-markers) of the simulationare plotted. The error bar becomes artificially big when thegrowth rate is stable or very small in the unstable regime.
D. Energy evolution and energy balance
The next step of the verification process is to quan-tify the global energy balance of the simulations, whichprovides an error measure when the analytic calcula-tion of the linear dispersion relation cannot be used, dueto inclusion of a non-uniform background magnetic fieldand/or nonlinear effects. The total energy of the modelgiven by Eqs. (1)-(6) is E = 12 (cid:90) d x (cid:20) (1 + τ ) (cid:16) ˜ n + |∇ ⊥ ˜ φ | + β |∇ ⊥ ˜ ψ | (cid:17) +˜ v i,z + τ Γ − T i (cid:21) , (19)where the first term corresponds to the internal poten-tial energy of the ions and electrons, the second term tothe perpendicular kinetic energy, the third term to themagnetic fluctuation energy, the fourth term to the par-allel kinetic energy and the final term to the ion thermalenergy. After some analytic computations the time evo-lution of the total energy is d t E = E P + E M + E Q + E J − D η − D D − D ν , (20)by specifying all components of the energy balance wherethe individual terms are given by E P = λτ v (cid:63)B P, (21) E M = τ (1 + τ ) v (cid:63)B M, (22) E Q = (cid:20) { τ (1 + τ ) − λ } v (cid:63)B + τ Γ − v (cid:63)T i (cid:21) Q, (23) E J = (1 + τ ) (cid:68)(cid:16) ˜ n − ˜ φ (cid:17) ˜ ∇ (cid:107) j z ( x, y ) (cid:69) x,y (24) D η = (1 + τ ) η (cid:10) ˜ j z (cid:11) x,y , (25) D D = D ⊥ (cid:28) (1 + τ ) (cid:16) |∇ ⊥ ˜ n | + |∇ ⊥ ˜ φ | (cid:17) + |∇ ⊥ ˜ v i,z | + τ Γ − |∇ ⊥ ˜ T i | (cid:29) x,y (26) D ν = ν s (cid:18) (1 + τ ) (cid:104) ˜ n (cid:105) x,y + (cid:104) ˜ v i,z (cid:105) x,y + τ Γ − (cid:68) ˜ T i (cid:69) x,y (cid:19) , (27)with the positive driven terms related to: a pressuregradient flux P = (cid:68) ˜ T i ∂ y ˜ n (cid:69) x,y , the radial particle flux M = (cid:104) ˜ n ˜ v x (cid:105) x,y , the radial ion heat flux Q = (cid:68) ˜ T i ˜ v x (cid:69) x,y due to the E × B drift [55] where ˜ v x = − ∂ y ˜ φ , and thenegative dissipative terms related to: the resistivity η ,the numerical diffusion D ⊥ and the axisymmetric sinkterms proportional to ν s which are used for the conser-vation of the background radial profiles. The coupling ofthe fluctuations to background axial and island currentgradients E J is added in order to verify its contribution in the energy balance. However, our choice of magneticshear B y ( x ) does not contribute to E J since its secondderivative vanishes, and so for the cases considered here, E J contains only the contribution of the imposed islandcurrent gradient. In Fig. 2 and 3 below it is shown to benegligible for all island widths relative to the dominantsource and sink terms.Each of these terms can be calculated in the simula-tions and the energy conservation properties of a givensimulation can thereby be quantified in detail. The terms E f are in general positive drives ( E J can be positive ornegative), and the terms D f are dissipative contributions.As mentioned previously, we note that for the equationsfor E P and E Q , the presence of the FLR term λ = 1stabilizes the turbulence. In fact, the magnetic gradientpart of the energy balance related to the heat flux E Q is divided by 2 (when τ = 1) and Fig. 3 shows a nega-tive contribution of E P due to the presence of the FLRterm. By calculating each of these terms from the simu-lation output, a numerical error measure (cid:15) can be readilydefined as (cid:15) = d t E num − d t EE . (28)For the nonlinear simulations we are using 125 radialgrid points and 128 poloidal grid points with poloidal pe-riodicity and radial Dirichlet boundary conditions. Thesize of the box is (2 a, b ) = (65 ρ s , ρ s ) so the discretiza-tion is ( dx, dy ) ∼ (0 . ρ s , . ρ s ). There are no signifi-cant effects of increasing the number of grid points by afactor 2 in each direction. With these parameters, the mi-croturbulence is well described and the saturation stateis reached after a few hundreds of time steps ( a/c s ) asshown below.Fig. 2 shows the time evolution of the accuracy of theenergy balance between the analytic prediction given byEq. (20) and the numerical computation. All componentsof the analytic prediction of the energy balance are shownin different colors. We can clearly see that the dominantcomponent of the energy balance is the one related to theradial ion heat flux, as expected in our ITG driven model.The numerical computation of the energy balance givenby the time derivative of the numerical total energy isalmost exactly identical to the analytic prediction. Therelative error between the analytic and numerical com-putation of the energy balance (cid:15) given by Eq. (28) isnegligible with respect to the dominant term E Q . We donot show here but, as expected, the part of the energy re-lated the ion temperature fluctuations is dominant withrespect to the other terms in Eq. (19). This explains thedominance of E Q with respect to the other components E P , E M and E J . The values of the dissipative terms( D = 0 . ν s = 0 . η = 10 − , and β = 10 − ) act on thesaturation of the numerical simulations. Indeed in thisnonlinear simulation, with the chosen parameters, thesaturation is obtained after approximately 400 a/c s timeunits since the energy balance fluctuates around zero fora static island of width W island = 15 ρ s .In Fig. 3, time-averages (over the steady-phases) ofnormalized energies and energy balances d t E with vary-ing island width and the ion temperature gradient v (cid:63)T i are shown, demonstrating the maintenance of good con- t ( a/c s ) × d t E E P E M E Q E J − D η − D D − D ν d t E num † FIG. 2: Time evolution of the energy balance in arbitraryunits with a static background magnetic island of width W island = 15 ρ s . Each term of Eq.(20) is represented in coloras well as their sum d t E (plain black curve) and the numericalvalue d t E num (dashed-black curve). The relative difference (cid:15) defined by Eq. (28) between the two later is negligible (boldblack curve). W island ( ρ s ) d t E = S − D† =( d t E num − d t E ) / E SDE J E Q E Q for v Ti =0 . E Q for v Ti =0 . FIG. 3: Scaling in arbitrary units of time-averaged terms inthe energy balance analysis as a function of v (cid:63)T i = a/L Ti overthe time range t ∈ [500 , a/c s in the presence of islandsof widths W island = { , , , , , } , and magnetic islandsize at fixed v (cid:63)T i = { . , . , . } . The sum of the positivesource (respectively negative dissipative) terms of the energybalance given by Eq. (20) are noted S (respectively D ). Theerror bars correspond to the standard deviation. servation for the simulation results discussed in the fol-lowing sections. All curves correspond to an ion temper-ature gradient of v (cid:63)T i = 0 . E Q for v (cid:63)T i = 0 . .
25. The numerical erroris shown by the difference d t E = S − D between thesource S = E P + E M + E Q + E J and dissipative terms D = D η + D D + D ν which is near 0. In particular, thevalue of the standard deviation of the theoretical predic-tion of the energy balance (std( d t E )) averaged over theisland width is ≈ × − which is much smaller thanthe standard deviation of the dominant term (std( E Q ))averaged over the island width ≈ − . Moreover, thetheoretical prediction of the energy balance d t E is not0 but is around 10 − , which is much smaller than thedominant term E Q ≈ .
03. Saturated energies increasewith island width as discussed in Sec. III, however thedependence of the normalized energy and energy balanceon island width is weak, and the numerical error is low.
III. QUANTIFYING THE IMPACT OF STATICISLAND ON ITG TURBULENCE
In our simulations we force the system to maintainon average the same level of turbulence due to the ITGby keeping an average background ion temperaturegradient. The electrons are assumed to be correlatedwith the ions and are not the focus of this work. Thetime evolution of the magnetic island is known to beat least one order of magnitude slower than the timeevolution of the microturbulence driven by the ITG.An electrostatic model ( β = 0) allowing the exact forcebalance of the fluctuations in the Ohm’s law betweenthe parallel classical resistive current and the parallelgradient of the drifts (i.e., the electric and diamagneticdrifts) is not used in this paper. Instead, weaklyelectromagnetic (small β = 10 − ) cases are considered,close to the electrostatic cases but, preparing for thefull electromagnetic self-consistency simulations andallowing comparisons with other published electrostaticresults. The advantage of full electromagnetic sim-ulations is that the effects of the turbulence on theslow dynamical evolution of the magnetic island willbe self-consistently included without using the artifi-cially large numerical perpendicular transport coefficient.In Fig. 4 the saturated energy E sat averaged on thetime range t ∈ [500; 3000] a/c s is plotted against thesize of the static magnetic island as well as its stan-dard deviation. We observe that beyond a thresholdof the island size (around 5 ρ s ) the saturated energy in-creases significantly. Further analysis has shown that itis due to an increase of the ion temperature fluctuationsand not to the island magnetic energy. Indeed, as men-tioned previously, the dominant part of the energy givenby Eq. (19) is the fluctuations of the ion temperature1 / (cid:82) τ ˜ T i / (Γ − T i fluctuate around W island ( ρ s ) E sat E sat for v Ti =0 . E sat for v Ti =0 . FIG. 4: Energy saturation in arbitrary units for different is-land sizes averaged over the time range t ∈ [500; 3000] a/c s for an ion temperature gradient v (cid:63)Ti = 0 . . . . . W island = 15 ρ s .The amplitude of the fluctuations are approximately 2times bigger (i.e., ∼ ± .
2) for W island = 15 ρ s with re-spect to the one with no island. Fig. 5 shows time slicesat times t = { , , , , } a/c s of the iontemperature T i + ˜ T i without ( W island = 0 ρ s ) and with( W island = 15 ρ s ) a static magnetic island. At this partic-ular time slices with a static magnetic island, Fig. 5(b),5(d), 5(f), 5(h) and 5(j) show the propagation of an iontemperature finger-like structure across the island asso-ciated with a clockwise flow inside the island. This prop-agation is due to the enhancement of the perpendicularheat flux transport close to the X-point and the clockwiseflow rotation inside the island. At other time slices, wecan observe opposite characteristics with a radial sym-metry with respect to x = 0 for lower ion temperaturestructures associated with counter-clockwise flow rota-tion inside the island. Fig 6 shows the time average ofthe ion temperature over t ∈ [500 , a/c s . We ob-serve the expected flattening of the ion temperature dueto the effect of the island. However, this flattening isnot poloidally invariant. The flattening seems to be ex-tended toward the top of the island. This is due to thedirection of the poloidal zonal flow which goes along de-creasing poloidal coordinates at the middle of the box(i.e., | x | <
10) and along increasing poloidal coordinatesat the edge of the box (i.e., | x | > < T i + ˜ T i > y,t of the total ion temperature as well as the profiles acrossthe X-point < T i + ˜ T i ( y = b/ > t and the O-point < T i + ˜ T i ( y = 0) > t . As a result, the flattening of theion temperature is dominant on the profile across the O-point ( y = 0). As explained with the previous figure, thisis due to the fact that all positive fluctuations of the ion
20 0 20 x ( ρ s ) y ( ρ s ) (a) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (b) W island = 15 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) (c) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (d) W island = 15 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) (e) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (f) W island = 15 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) (g) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (h) W island = 15 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) (i) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (j) W island = 15 ρ s FIG. 5: Snapshots of the total ion temperature T i = T i ( x ) +˜ T i ( x, y ) at t = { , , , , } a/c s respectivelyfor (a)-(b), (c)-(d), (e)-(f), (g)-(h) and (i)-(j), obtained fromtwo simulations with W island = 0 ρ s for (a), (c), (e), (g), (i)and with W island = 15 ρ s for (b), (d), (f), (h), (j). The samecolorbar limits (i.e., [0 .
55; 1 .
45] and [0 .
1; 1 . W island = 0 ρ s and W island = 15 ρ s ).
20 0 20 x ( ρ s ) y ( ρ s ) (a) W island = 0
20 0 20 x ( ρ s ) y ( ρ s ) (b) W island = 15 ρ s FIG. 6: Time average on t ∈ [500; 3000] a/c s of the total iontemperature T i = T i + ˜ T i for (a) W island = 0 ρ s and (b) W island = 15 ρ s .
20 0 20 x ( ρ s )
0; 10; 20 } ρ s respectively for the figures (a),(b) and (c). With no static island (a), the flux is concen-trated around the resonant surface at x = 0. This timeaveraged radial ion heat flux is almost homogeneous inthe poloidal direction. By increasing the size of the staticmagnetic island, we observe a compression of the time av-eraged ion heat flux toward the X-point but, more pre-cisely toward the top of the island. Indeed, the poloidalprofile is more peaked near y ∼ ρ s and the maximumvalue of this poloidal profile is approximately 3 timeshigher for W island = 20 ρ s than for W island = 10 ρ s andapproximately 10 times higher than for the case withno island. We observe that the shape of the radial ionheat flux is consistent with the previous observation ofthe flattening of the ion temperature across the O-point.We also remark that the radial ion heat flux average ispositive because the radial transport is such that ˜ v x isnegative (respectively positive) for negative (respectivelypositive) ion temperature fluctuations ˜ T i . To obtain moredetailed understanding of the global effects of the islandon the turbulence and transport, we need to compare theradial ion heat flux by averaging Fig 8, over the poloidaldimension y .We use a parameterization of the poloidal average (cid:104)· · · (cid:105) y of the time-averaged ion heat fluxes of Fig. 8. The pa-rameterization is done using a fitting of the profiles withthe following Gaussian f ( x ) = f + f exp (cid:32) − x W f (cid:33) . (29)The variation of the parameter f describes the effects ofthe island outside the island (i.e., | x | > W island ), f + f the effects around the position x = 0 inside the island,and W f the characteristic width of the Gaussian pro-file. Fig. 9 represents the variations of these parametersused by the Gaussian fitting of the radial ion heat flux Q = ˜ v x ˜ T i corresponding to the transport property andthe squared electric potential ˜ φ corresponding to theturbulence property.We observe the fact that the width of the Gaussian pro-file fits as well as the maximum value at x = 0 increasefor the transport and the turbulence beyond a thresholdisland size. The threshold is around 5 ρ s for the trans-port and between 5 ρ s and 10 ρ s for the turbulence. These
20 0 20 x ( ρ s ) y ( ρ s ) × (a) Radial heat flux for W island = 0 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) × (b) Radial heat flux for W island = 10 ρ s
20 0 20 x ( ρ s ) y ( ρ s ) × (c) Radial heat flux for W island = 20 ρ s FIG. 8: Comparison of the time average of the radial heatflux Q = ˜ v x ˜ T i over the time range t ∈ [500 , a/c s between the following sizes of the static magnetic island W island ∈ {
0; 10; 20 } ρ s . The white plain (resp. dashed)curves are the island separatrix (resp. 15 iso-contour mag-netic flux surfaces). The radial ion heat flux is peaked insidethe island between the X-point and the O-point. W island ( ρ s ) W Q (a) Width W Q of the fittingGaussian profile W island ( ρ s ) × Q Q + Q (b) Max Q + Q and min Q of the fitting Gaussian profile W island ( ρ s ) W φ (c) Width W φ of the fittingGaussian profile W island ( ρ s ) φ φ + φ (d) Max φ + φ and min φ ofthe fitting Gaussian profile FIG. 9: Comparison of the Gaussian fitting parametersagainst the size of the static magnetic island W island ∈{
0; 1; 2; 5; 10; 15; 20 } ρ s for (a)-(c) the radial ion heat flux andfor (b)-(d) the square of the electrostatic potential. The timerange is t ∈ [500 , a/c s . Beyond a threshold around W island ∼ ρ s , the static island globally enhances the trans-port and turbulence even outside the island. thresholds are about the same order as the characteristiclength of the turbulence 5 ρ s to 10 ρ s of our simulations.Moreover, the parameters Q and φ correspond to theeffect of the island respectively on the transport and theturbulence outside the island. Even if it is a small vari-ation, we observe an increase of the transport and theturbulence far away from the island. This observationmeans that there is an enhancement of the turbulenceand the transport independent of a simple displacementby the presence of the island.Finally, these results given by a simple poloidal averagecan be verified by using another average which conservesthe structure of the static magnetic island. To meet thisgoal, a last investigation is to try to recover the pre-vious results by using the flux surface average (cid:104)· · · (cid:105) ψ eq instead of the poloidal average (cid:104)· · · (cid:105) y . Fig. 10 representsthe flux surface averages of the time-averaged radial ionheat flux for different sizes of the magnetic island. Thefirst figure corresponds to profiles of the radial ion heatflux inside the island versus the flux surface coordinate( ψ − ψ ) / ( ψ sep − ψ ) for different sizes of the magneticisland and the second figure to profiles outside the islandversus ( ψ − ψ sep ) / ( ψ a − ψ sep ). The values ( ψ , ψ sep , ψ a )of the magnetic flux correspond respectively to the O-point, the separatrix, and the last closed flux surface inour simulation box. Beyond the last closed flux surface,all integrals along each opened flux surface linearly de-crease to 0, which is their minimum values at the corner0of the simulation box of the flux surface correspondingto ψ eq ( − a, b/ ( ψ − ψ ) / ( ψ sep − ψ ) W island =5 ρ s W island =10 ρ s W island =15 ρ s W island =20 ρ s (a) Flux surface average of the heat flux inside the island. ( ψ − ψ sep ) / ( ψ a − ψ sep ) W island =5 ρ s W island =10 ρ s W island =15 ρ s W island =20 ρ s (b) Flux surface average of the heat flux outside the island. FIG. 10: Comparison of the time and flux surface averagesof the radial ion heat flux for the following sizes of the staticmagnetic island W island ∈ {
5; 10; 15; 20 } ρ s . The effect of thepoloidal asymmetry inside the island and the threshold of W island ∼ ρ s are consistent with these flux surface averages. the island (Fig. 10(a)) an increase of the shift of the ra-dial ion heat flux toward the separatrix instead of beinghomogeneous everywhere inside the island. However, theasymmetry between the top and the bottom of the islandcannot be recovered in comparison to what have beenshown above with the 2D slices. With the observation ofthe radial ion heat flux outside the island (Fig. 10(b)), thebackground transport increases beyond a threshold sizeof the magnetic island between 5 ρ s and 10 ρ s , related to the characteristic length of the transport with no island.This enhancement of the transport outside the island isdue to the filaments which appear at different time slicesas described previously. As an example with a staticmagnetic island W island = 20 ρ s , we typically observe fil-aments of the positive fluctuations of the ion temperaturestarting around the position ( x, y ) = ( − , x, y ) = (20 ,
0) and finally diffusing outsidethe island in the region of lower ion temperature.
20 0 20 x ( ρ s ) y ( ρ s ) (a) (cid:104) φ (cid:105) t for W island = 0 ρ s .
20 0 20 x ( ρ s ) y ( ρ s ) (b) (cid:104) φ (cid:105) t for W island = 15 ρ s . FIG. 11: Comparison of the time averaged electric potential (cid:104) φ (cid:105) t on t ∈ [500; 3000] a/c s for a static island of width (a) W island = 0 ρ s and (b) W island = 15 ρ s .
20 0 20 x ( ρ s ) y ( ρ s ) (a) Mode k y = 0 for W island = 0 ρ s .
20 0 20 x ( ρ s ) y ( ρ s ) (b) Mode k y = 0 for W island = 15 ρ s .
20 0 20 x ( ρ s ) y ( ρ s ) (c) Mode k y = ± π/b for W island = 0 ρ s .
20 0 20 x ( ρ s ) y ( ρ s ) (d) Mode k y = ± π/b for W island = 15 ρ s . FIG. 12: Reconstruction of the 2D electric potential from oneselected poloidal mode: (a)-(b) from k y = 0 and (c)-(d) from k y = ± π/b with (a)-(c) no island and (b)-(d) W island = 15 ρ s . It is important to note that in addition to the turbu-lent heat flux, there is also transport due to static E × B flows which arise due to the presence of the island, thatwill also drive transport. Fig. 11 shows the time-averagedelectrostatic potential for the cases of W island = 0 ρ s and15 ρ s , respectively. In the no-island case, the mean po-tential corresponds to a k y = 0 zonal flow with radial1wavenumber approximately equal to π/a . However, inthe finite island case (Fig. 11(b)), the mean potentialhas both significant k y = 0 and k y = 2 π/b components,with comparable magnitude to the no-island case. Thesecomponents are separately visualized in the (x,y) plane inFig. 12. Note in particular that the finite k y equilibriumflow in the finite island case is not symmetric about theO-point at y = 0, and thus contributes to the asymmetryin radial heat flux shown in Fig. 8.As a conclusion of our investigation of the effects of astatic island on turbulence and transport, we observedthe flattening of the ion temperature profile across theO-point ( y = 0) due to the island as well as an increaseof its gradient across the X-point. The increase of theion temperature gradient around the X-point enhancesthe ITG instability and drives more turbulence in thisregion. Due to the radial transport, positive or negativestructures of the ion temperature enter inside the islandand flow along the magnetic flux surfaces. Inside the is-land, the time average of the radial ion heat flux is thenpeaked on the top of the island due to the fact that theflow naturally goes along decreasing poloidal coordinatesaround the radial position x = 0 (as observed with noisland). This privileged direction of the flow is due toour choice of the sign of the equilibrium poloidal mag-netic flux (i.e., the sign of the magnetic shear). Finally,because of the presence of the island, the turbulence andtransport are significantly increased when the size of thestatic magnetic island is bigger than the characteristicsize of the turbulence and transport.To make progress toward understanding the self-consistent interaction, some results of the feedback ofthe modified turbulence on the dynamical evolution ofthe island is investigated below. IV. FIRST RESULTS TOWARD THEFEEDBACK ON THE ISLAND
One of our interests is to understand the self-consistentinteraction between fast dynamics of microturbulenceand the slow dynamics of a magnetic island. In the pre-vious section we focused on the effect of a static magneticisland on transport and ITG microturbulence. Fromthese results, we report here a quantitative study of theeffect of microturbulence on the slow dynamics of themagnetic island.
A. Formulation of a turbulent resistivity
We begin by adopting a mean-field type approach toexploit the separation of time and spatial scales betweenthe turbulence and the island. We can then rewrite themagnetic flux evolution equation as β∂ t ψ I = − B y, ∂ y φ I − η cl j z,I + β (cid:104) ˜ ψ, ˜ φ − ˜ n (cid:105) , (30) where B y, ∂ y = ∇ (cid:107) = ∇ (cid:107) − ˜ ∇ (cid:107) , the second term is theusual classical resistive current and the last term corre-sponds to the average of the Ohm’s law stress. This lastterm can be written as a coherent part − η turb j z,I anda non-coherent part − η NC j NC with respect to the islandcurrent where the effective turbulent resistivity describesthe self-consistent feedback of the turbulence on the dy-namics of the island η turb = − (cid:68)(cid:104) ˜ ψ, ˜ φ − ˜ n (cid:105) j z,I (cid:69) x,y (cid:68) j z,I (cid:69) x,y . (31)Moreover, in addition to the dynamical evolution of themagnetic flux of the island ψ I , a dynamical evolutionof the electric potential φ I associated with the magneticisland has also been developed. Here, we do not focus onthis part because in our previous simulations we imposedonly a poloidal magnetic flux and we omit the presenceof an associated electric potential. Indeed, the associatedelectric potential φ I of the island is a very slow dynamicalquantity in comparison to the fast evolution of ˜ φ and φ I has almost no effect on the dynamical evolution of themagnetic island poloidal flux. B. Evaluation of turbulent resistivity in simulation
For the dynamical equation of the island, simulationscan be challenging since we need to resolve the fasttime scales of the microturbulence as well as the slowtime scales of the evolution of the magnetic island andthe tearing mode. In this work, we have exploited thetimescale separation between the island and turbulencedynamics to approximate the island as fixed, consistentwith previous studies by other groups. However, we canstill begin to probe how we might expect the turbulenceto act on the island in a fully dynamic model by calculat-ing the time-averaged forcings of the turbulence on theimposed island structure.Fig. 13 represents the Ohm’s Law stresses (cid:104) ˜ ψ, ˜ φ − ˜ n (cid:105) averaged over t ∈ [500; 3000] a/c s for different sizes of thestatic magnetic island W island ∈ {
1; 2; 5; 10; 15; 20 } ρ s .There is an obvious link between the structures of theseOhm’s Law stresses and those of the radial ion heat fluxesshown in Fig. 8. The effective turbulent resistivity η turb given by the Eq. (31) is then computed from the extrac-tion of the coherent part of the Ohm’s law stress with re-spect to the static current associated with the magneticisland. The ratio between the effective turbulent resis-tivity and the classic resistivity η cl is shown in Fig. 14.The first result is the expected fact that the effectiveturbulent resistivity, which quantifies the effects of thefeedback of the turbulence on the dynamical evolution ofthe magnetic island, decreases by 1 order of magnitudeas the island size increases. The second result is the am-plitude scale between the effective turbulent resistivityand the classic resistivity. Indeed, for all island sizes, the2 y ( ρ s ) × × y ( ρ s ) × ×
20 0 20 x ( ρ s ) y ( ρ s ) ×
20 0 20 x ( ρ s ) × FIG. 13: Ohm’s law stresses (cid:104) ˜ ψ, ˜ φ − ˜ n (cid:105) averaged over t ∈ [2250; 3000] a/c s for different sizes of the static magnetic is-land W island ∈ {
1; 2; 5; 10; 15; 20 } ρ s (from left to right and topto bottom). The black plain (resp. dashed) curves are the is-land separatrix (resp. 15 iso-contour magnetic flux surfaces).The structures of the Ohm’s law stresses are consistent withthe radial ion heat fluxes. W island ( ρ s ) η turb /η cl FIG. 14: Ratio of the effective turbulent resistivity η turb givenby Eq.(31) over the classic resistivity η cl versus the size of theisland W island . effective turbulent resistivity is at least 1 order of mag-nitude higher than the classic resistivity. However, evenwith a difference of 2 orders of magnitude for small is-lands, the total effect of the self-consistent feedback ofthe turbulence on the island is given by the multiplica-tion of the effective turbulent resistivity η turb with the current generated by the magnetic island j z,I . The lateris approximately 10 times smaller for small island sizesthan for larger island sizes. This investigation clarifiessome expected mechanisms toward the self-consistent in-teraction between the ITG microturbulence and the dy-namical evolution of a magnetic island. Perspectives arediscussed below. V. DISCUSSIONS AND FUTURE WORK
In this work, we have studied the dynamics and asso-ciated thermal transport due to ITG microturbulence inthe presence of static magnetic islands of varying width.We find that the magnetic island is responsible for en-hancement of the turbulent fluctuation amplitudes andradial heat flux, which increases rapidly above a thresh-old island width comparable to the correlation length ofthe turbulence in the no-island limit. Broadly consistentwith theoretical expectations, we observe a flattening ofion temperature inside the island, and localization of theradial heat flux near (but not symmetric about) the is-land X-point. It should be noted that this flatteningoccurs due to turbulent transport processes and in-planeflow along flux surfaces, as the simulations presented heredo include explicit conductive parallel heat fluxes of theform q i, (cid:107) = − k (cid:107) ∇ (cid:107) T i that are often used to representboth parallel collisional dynamics and Landau-dampingeffects. Thus, the combination of self-consistently calcu-lated radial and poloidal turbulent heat fluxes in the pres-ence of the imposed island (along with the self-consistentelectrostatic flow which also forms) are sufficient to leadto temperature flattening and equilibration on the totalmagnetic flux surfaces. A novel observation has beenmade that a asymmetric shift in location of the peak ra-dial heat flux relative to the X-point is due to inclusionof self-consistent turbulence effects (such as the inherentphase velocity of the turbulence) which are not includedin theoretical studies which simply approximate the tur-bulent fluxes in terms of anomalous diffusion coefficients.In addition to calculating the response of the turbulenceto the presence of the island, we also investigated theback-reaction of the turbulence on the island, using amean-field approach to calculate an effective turbulentresistivity which would act on the island structure. Forthe parameters considered, our simulations predict thatwhile this turbulent resistivity is always at least an or-der of magnitude larger than the collisional value used,it is most effective at small island widths and decreasesquickly with increasing island size.There are a wide variety of future directions to be pur-sued in future work, building upon these results. Fore-most are inclusion of parallel heat flux due to both col-lisions and Landau damping, as well as fully 3D effects,to make closer connection with experiments. Equally im-portant in our opinion is to transition from statically im-posed islands to self-consistently evolving islands drivenby tearing unstable current profiles. An important nu-3merical challenge to be solved for these simulations in3D is the implementation of preconditoners and implicitadvance schemes suitable for long-time integration of theslowly evolving island in the presence of fast ITG turbu-lence and even faster damped Alfv´en waves in 3D geom-etry. Greater understanding the self-consistent transportenhancements and profile flattening effects, as well as tur-bulent forcings of the island are needed to develop a fullypredictive model of coupled tearing mode and turbulencedynamics. We can then begin to include more realisticgeometric effects in the equilibrium structure, to prop-erly describe the coupled turbulence and reconnection processes which drive observed magnetic islands. Acknowledgements
This work was supported by U.S. DoE Contract No.DE-SC0007783 and DE-SC0010520. O.I. would like tothank Orso Meneghini for his help with OMFIT. Com-puting time was provided on the Triton Shared Com-puting Cluster (TSCC) at the San Diego SupercomputerCenter (SDSC). [1] E.J. Doyle, et al.,
Nucl. Fusion , S18 (2007). dx.doi.org/10.1088/0029-5515/47/6/S02 . Nucl. Fusion , 099801 (2008). dx.doi.ord/10.1088/0029-5515/48/9/099801 [2] W. Horton, Rev. Mod. Phys. , 735 (1999). dx.doi.org/10.1103/RevModPhys.71.735 [3] H.P. Furth, Phys. Fluids , 48 (1963). dx.doi.org/10.1063/1.1724507 [4] P. Helander and D.J. Sigmar, Collisional Trans-port in Magnetized Plasmas , Cambridge UniversityPress, Cambridge, UK (2002). adsabs.harvard.edu/abs/2002ctmp.book.....H [5] J.P. Freidberg,
Ideal magnetohydrodynamics , PlenumPress, New York, NY (1987).[6] Biskamp,
Phys. Fluids B , 3893 (1993). dx.doi.org/10.1063/1.860612 [7] B.D. Scott, A.B. Hassam and J.F. Drake, Phys. Fluids , 275 (1985). dx.doi.org/10.1063/1.865197 [8] R. Fitzpatrick, Phys. Plasmas , 825 (1995). dx.doi.org/10.1063/1.871434 [9] P.H. Diamond, S-I Itoh, K Itoh and T.S. Hahm, PlasmaPhys. Control. Fusion , R35 (2005). dx.doi.org/10.1088/0741-3335/47/5/R01 [10] Z. Chang, J.D. Callen, E.D. Fredrickson, R.V. Budny,C.C. Hegna, K.M. McGuire, M.C. Zarnstorff, and TFTRgroup, Phys. Rev. Lett. , 4663 (1995). dx.doi.org/http://dx.doi.org/10.1103/PhysRevLett.74.4663 [11] R. Carrera, R.D. Hazeltine and M. Kotschenreuther, Phys. Fluids , 899 (1986). dx.doi.org/10.1063/1.865682 [12] C.J. McDevitt and P.H. Diamond, Phys. Plasmas ,032302 (2006). dx.doi.org/10.1063/1.2177585 [13] A. Sen, R. Singh, D. Chandra, P. Kaw and D. Raju, Nucl. Fusion , 115012 (2009). dx.doi.org/10.1088/0029-5515/49/11/115012 [14] A. Ishizawa and N. Nakajima, Phys. Plasmas , 040702(2007). dx.doi.org/10.1063/1.2716669 [15] A. Ishizawa and N. Nakajima, Nucl. Fusion , 1540(2007). dx.doi.org/10.1088/0029-5515/47/11/016 [16] F. Militello, F.L. Waelbroeck, R. Fitzpatrick, andW. Horton, Phys. Plasmas , 050701 (2008). dx.doi.org/10.1063/1.2917915 [17] F.L. Waelbroeck, F. Militello, R. Fitzpatrick, andW. Horton, Plasma Phys. Cont. Fusion , 015015(2009). dx.doi.org/10.1088/0741-3335/51/1/015015 [18] F.L. Waelbroeck, Nucl. Fusion , 104025 (2009). dx. doi.org/10.1088/0029-5515/49/10/104025 [19] M. Muraglia, O. Agullo, S. Benkadda, X. Garbet,P. Beyer and A. Sen, Phys. Rev. Lett , 145001 (2009). dx.doi.org/10.1103/PhysRevLett.103.145001 [20] M. Muraglia, O. Agullo, S. Benkadda, M. Yagi,X. Garbet and A. Sen,
Phys. Rev. Lett , 095003(2011). dx.doi.org/http://dx.doi.org/10.1103/PhysRevLett.107.095003 [21] A. Ishizawa and N. Nakajima,
Phys. Plasmas , 072308(2010). dx.doi.org/10.1063/1.3463435 [22] A. Ishizawa and P.H. Diamond, Phys. Plasmas ,074503 (2010). dx.doi.org/10.1063/1.3460346 [23] A. Ishizawa, S. Maeyama, T.H. Watanabe, H. Sugamaand N. Nakajima, Nucl. Fusion dx.doi.org/10.1088/0029-5515/53/5/053007 [24] O. Agullo, M. Muraglia, A. Poy, S. Benkadda, M. Yagi,X. Garbet, and A. Sen, Phys. Plasmas , 092303(2014). dx.doi.org/10.1063/1.4894699 [25] Z.Q. Hu, Z.X. Wang, L. Wei, J.Q. Li and Y. Kishimoto, Nucl. Fusion , 123018 (2014). dx.doi.org/10.1088/0029-5515/54/12/123018 [26] P. Hill, F. Hariri, and M. Ottaviani, Phys. Plasmas ,042308 (2015). dx.doi.org10.1063/1.4919031 [27] H.R. Wilson and J.W. Connor, Plasma Phys. Con-trol. Fusion , 115007 (2009). dx.doi.org/10.1088/0741-3335/51/11/115007 [28] A.G. Peeters, Y. Camenen, F.J. Casson,W.A. Hornsby,A.P. Snodin, D. Strintzi andG. Szepes, Comput. Phys. Comm. , 2650 (2009). dx.doi.org/10.1016/j.cpc.2009.07.001 [29] E. Poli, A. Bottino and A.G. Peeters,
Nucl. Fusion , 075010 (2009). dx.doi.org/10.1088/0029-5515/49/7/075010 [30] E. Poli, A. Bottino, W.A. Hornsby, A.G. Peeters,T. Ribeiro, B.D. Scott and M. Siccinio, Plasma Phys.Control. Fusion , 124021 (2010). dx.doi.org/10.1088/0741-3335/52/12/124021 [31] W.A. Hornsby, M. Siccinio, A.G. Peeters, E. Poli,A.P. Snodin, F.J. Casson, Y. Camenen and G. Szepesi, Plasma Phys. Control. Fusion , 054008 (2011). dx.doi.org/10.1088/0741-3335/53/5/054008 [32] M. Siccinio, E. Poli, F.J. Casson, W.A. Hornsby andA.G. Peeters, Phys. Plasmas , 122506 (2011). dx.doi.org/10.1063/1.3671964 [33] R.E. Waltz and F.L.Waelbroeck, Phys. Plasmas ,032508 (2012). dx.doi.org/10.1063/1.3692222 [34] Y.-H. Liu, W. Daughton, H. Karimabadi, H. Li, andS.P. Gary, Phys. Plasmas , 022113 (2014). dx.doi.org/10.1063/1.4865579 [35] D. Zarzoso, F.J. Casson, W.A. Hornsby, E. Poli, andA.G. Peeters, Phys. Plasmas , 022127 (2015). dx.doi.org/10.1063/1.4908549 [36] W.A. Hornsby, P. Migliano, R. Buchholz, L. Kroenert,A. Weikl, A.G. Peeters, D. Zarzoso, E. Poli and F.J. Cas-son, Phys. Plasmas , 022118 (2015). dx.doi.org/10.1063/1.4907900 [37] W.A. Hornsby, P. Migliano, R. Buchholz, D. Zarzoso,F.J. Casson, E. Poli and A.G. Peeters, Plasma Phys.Control. Fusion , 054018 (2015). dx.doi.org/10.1088/0741-3335/57/5/054018 [38] P.H. Rutherford, Phys. Fluids , 1903 (1973). dx.doi.org/10.1063/1.1694232 [39] B. Dudson, M.V. Umansky, X.Q. Xu, P.B. Snyder andH.R. Wilson, Comput. Phys. Comm. , 1467 (2009). dx.doi.org/10.1016/j.cpc.2009.03.008 [40] R.D. Hazeltine, M. Kotschenreuther and P.J. Morrison,