A Quantitative Model for Heat Pulse Propagation across Large Helical Device Plasmas

Copyright and reuse: The Warwick Research Archive Portal (WRAP) makes this work of researchers of the University of Warwick available open access under the following conditions. This article is made available under the Creative Commons Attribution-3.0 Unported (CC BY 3.0) license and may be reused according to the conditions of the license. For more details see A note on versions: The version presented in WRAP is the published version, or, version of record, and may be cited as it appears here. It is known that rapid edge cooling of magnetically confined plasmas can trigger heat pulses that propagate rapidly inward. These can result in large excursion, either positive or negative, in the electron temperature at the core. A set of particularly detailed measurements was obtained in Large Helical Device (LHD) plasmas [S. Inagaki et al., Plasma Phys. Controlled Fusion 52, 075002 (2010)], which are considered here. By applying a travelling wave transformation, we extend the model of Dendy et al. successfully describes the local time-evolution of heat pulses in these plasmas, to include also spatial dependence. The new extended model comprises two coupled nonlinear first order differential equations for the (x, t) evolution of the deviation from steady state of two independent variables: the excess electron temperature gradient and the excess heat flux, both of which are measured in the LHD experiments. The mathematical structure of the model equations implies a formula for the pulse velocity, defined in terms of plasma quantities, which aligns with empirical expectations and is within a factor of two of the measured values. We thus model spatio-temporal pulse evolution, from first principles, in a way which yields as output the spatiotemporal evolution of the electron temperature, which is also measured in detail in the experiments. We compare the model results against LHD datasets using appropriate initial and boundary conditions. Sensitivity of this nonlinear model with respect to plasma parameters, initial conditions, and boundary conditions is also investigated. We conclude that this model is able to match experimental data for the spatio-temporal evolution of the temperature profiles of these pulses, and their propagation velocities, across a broad radial range from r=a ' 0:5 to the plasma core. The model further implies that the heat pulse may be related mathematically to soliton solutions of the Korteweg-de Vries-Burgers equation. V C 2015 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution 3.0 …


I. INTRODUCTION
The transport of energy across magnetically confined fusion plasmas, and the storage of energy within them, reflects a wide range of turbulent and nonlinear phenomenology; see, for example, Refs.1-32.There is extensive experimental evidence for transport phenomena that are non-diffusive and may be non-local.5][26][27] A broad range of techniques for data analysis, see, for example, Refs.12, 14, and 28, have been used to identify various forms of perturbation of heat and particle fluxes from their steady states, see the reviews. 29,30Measurements of the spatio-temporal propagation of strongly nonlinear localised heat pulses provide a particularly interesting, and potentially fruitful, challenge to theoretical understanding and models.Here, we focus on cold pulse experiments, see, for example, Refs. 18, 19, and 24.These cold pulses are typically initiated by injection into the edge plasma of ice pellets or supersonic molecule beams, by gas puffing, and by laser ablation.The zero-dimensional model of Ref. 25, which incorporates only time dependence, is successful in quantitatively capturing the local time evolution of drT e and dq e at a specific radius during two well diagnosed cold heat pulse experiments 24 using pellet injection in LHD; we refer, in particular, to Figs. 3 and 5 of Ref. 25.This motivates the following physical conjecture, which we test in the present paper.The zero-dimensional model is known to work at the best diagnosed spatial location, capturing the time evolution of the pulse there, where t ¼ 0 is defined to be the local arrival time of the initial impulsive perturbation.Therefore, we conjecture that the zero-dimensional model ought to work at each location across the radial domain of the plasma within which the same physical processes determine the behaviour of the pulse.From this, we infer that the model ought to apply in a frame that is co-moving radially with the heat pulse across this region.This final step in the conjecture provides a simple path to construct a spatio-temporally dependent model from the tested time-dependent-only model of Ref. 25.We apply a travelling wave transformation by replacing t in Eqs.(1-3) of Ref. 25 by n ¼ x þ v 0 t; here, v 0 is to be considered as a proxy of the pulse propagation velocity in the radial direction x, and the sign convention adopted for t assists consideration of inward propagation.Importantly, as we shall show, the mathematical structure of the resultant model equations, combined with the choice of physical model parameters that carries over from Ref. 25, yields a formula for v 0 that aligns with prior empirical expectations.

II. MODEL DESCRIPTION
The normalised zero-dimensional model 25 examined below is constructed in terms of the three key physical quantities that were measured 24 so as to characterise pulse propagation in LHD.These quantities are the deviation from steady state of the electron temperature gradient drT e , the excess turbulent heat flux dq e , and the deviation dT e of electron temperature from its steady-state value.The dimensionless counterparts of these variables are denoted by x 1 , x 2 , and x 3 , respectively, as defined in Ref. 25.The model of Ref. 25  shows quantitative agreement between its outputs and experimental measurements of the time evolution of these variables at fixed locations in the LHD plasma, after rapid cooling at the edge.It is successful both when core electron temperature rises on arrival of the pulse and when it drops.However, for given parameters and initial conditions, the model only simulates the time evolution of a passing heat pulse at a specific radius, for example, r=a ¼ 0:19 in Ref. 25. Spatial dependence is eliminated from the physical picture underlying the model in Ref. 25 by using the parameter 1=L c as a proxy for divergence in the heat flux energy conservation equation; here, L c is the characteristic scale-length of steady-state turbulent transport.The model equations ( 5) to (7) of Ref. 25 are as follows: (1) (2) Physically, Eqs.(1) to (3) embody on assumption of dominant strong coupling between drT e and dq e , while the background transport acts to dissipate the pulse.The physical significance of the various coefficients is described in the discussion of equations (1)-(4) of Ref. 25.We note that the mathematical structure of Eq. (3) above implies that x 3 is slave to x 2 .The change in the numerical value of x 3 at each time step is calculated directly from Eq. (3) using the timeevolving value of x 2 from the two-field system of Eqs.(1) and ( 2), in which x 3 acts as an x 2 -slaved coefficient.As noted, physically, the value of x 3 maps to the deviation in electron temperature over time, which is a primary experimentally measured time-dependent variable for the model comparisons reported below.Central to the present paper is the adoption of travelling wave transformations as a method for generating (x, t)-dependence from the t-dependent model of Eqs.(1) to (3).We assume that Here, x and t are considered to be independent variables, and we refer to v 0 as the pseudo-velocity of the pulse.This pseudo-velocity is expected to be similar, but not identical, to the real measured velocity of the pulse.It then follows that Due to the fact that y 3 is a dependent variable, we choose to simplify by neglecting all y 3 related terms in Eqs.
(1) to ( 3), and substitute Eq. ( 5) into Eqs.(1) and ( 2), yielding Our new model [Eqs.( 6) and ( 7)] comprises only two independent variables: y 3 values are deduced directly from the time series of y 2 , using Eq. ( 3) above.Let us now operate on Eq. ( 6) with d=dn and multiply Eq. ( 7) by j T0 , then eliminate dy 2 =dn by substitutions.This yields the following equation after leading order approximation and transpositions: The second coupled nonlinear ordinary differential equation is derived by applying similar procedures Equations ( 8) and ( 9) comprise our mathematical model for the propagating heat pulse.Its physical motivation is that of Ref. 25, combined with the co-moving conjecture outlined above and expressed in Eq. ( 4).The mathematical structure of Eq. ( 8) can be linked, in certain approximations, to the Korteweg-de Vries-Burgers (KdV-Burgers) equation, 31 and Eq. ( 9) can be linked to a damped wave equation in the comoving frame.In particular, the left hand side of Eq. ( 8) corresponds to that in the formulation of the KdV-Burgers equation in Eq. ( 5.3) of Ref. 31.This nonlinear differential operator, acting on y 1 , is known to support soliton pulses, and furthermore these pulses move at a speed determined by the coefficient of the term which multiplies the linear term in y 1 . 31The expression for this coefficient in Eq. ( 8) motivates the following conjecture regarding heat pulse velocity in our model: This scaling has previously been noted empirically from heat pulse experiments. 29,32It is possible that there will in future arise a convergence between the mathematical structure of the present model, empirical scalings, 29,32 and other theoretical work, 33 which encompasses Eq. ( 10).An important avenue to explore is a potential link between the present approach and the theory of turbulence spreading, 33 where the propagation velocity of a turbulence front scales as the square root of the product of a linear growth rate with heat diffusivity.

III. COMPARISON OF MODEL RESULTS WITH LHD EXPERIMENTAL DATA
We recall that y 1 , y 2 , and y 3 are the dimensionless counterparts of drT e ; dq e , and dT e respectively.Solution of our model, embodied in the two coupled nonlinear ordinary differential equations ( 8) and (9), proceeds as follows.The numerical values for j T ; j Q and their derivatives, together with the numerical values of c L1 ; c L2 ; g; s c , and v 0 , carry over identically from Ref. 25, to which we refer for the experimental motivation for these values; see Table I for detailed information.
The new model embodied in Eqs. ( 8) and ( 9) is strongly nonlinear.In solving it, we have three kinds of degrees of freedom.First, there are the initial values of y 0 i ðn ¼ 0Þ.Second, there is the value of the pseudo-velocity v 0 .Third, it is possible to apply fixed horizontal and vertical shifts in the values of the electron temperature and time, provided these shifts are applied uniformly to all outputs of a simulation, as a way of reducing systematic errors introduced by parameter choices.We optimise model outputs, for comparison with the experimental measurements, in these three ways.These model outputs (time traces of y 1 , y 2 , and y 3 ) reflect the underlying phase space structure of the solutions of the nonlinear system of equations.This structure is robust, in the sense discussed in Ref. 25, and as demonstrated for the present system of equation in Fig. 7 below, for example.We recall also that the great majority of model parameter values, see, e.g., Table I, are determined by experimental measurements.
Figure 1 compares time traces of the evolving electron temperature at multiple radial locations, obtained from the model and from experimental data(# Te49708) for the R case.Several representative radii are marked by arrows on the right hand side.The model results are able to match experimental data from q ¼ 0:450 inward to the core, if we uniformly apply horizontal ( þ 0.01) and vertical (þ0.20) shifts, suggesting that the model applies over this broad radial range.It is also clear that electron temperature profiles from q ¼ 0:546 outward to q ¼ 0:703 are not simulated by the model.This suggests that different physics dominates heat pulse propagation in the outer region of this plasma.Figure 2 demonstrates the comparison of model results and experimental data(# Te49719) for the D case.In common with the R case shown in Fig. 1, the model results are good from q ¼ 0:450 inward to the core, with uniformly applied horizontal (þ0.04) and vertical (þ0.07) shifts, but not in the outer region of this plasma.Figs. 3 and 4 show zoomed-in electron temperature pulse plots at three specific radii, selected from Figs. 1 and 2, respectively.
Empirically, we approach v 0 by identifying a best fit straight line linking the positions of the peak or trough of the heat pulse temperature profile in the two cases, see Figs. 5  and 6.The effect of data noise on peak location is minimised by taking a five per cent running window.The width of the  filter window determines the horizontal error bars.assume, for simplicity and to match our travelling-wave model, that the inward radial propagation velocity of the heat pulse is invariant across the plasma volume of interest.
A straight line fit to the radial location of the pulse peak versus time is thus applied.The error of this fitting may come from the radial dependence of the pulse velocity, which tends to increase in the central region.We find from Figs. 5 and 6 that v 0 ¼ 15 for the temperature rise case and v 0 ¼ 30 FIG. 1.Time evolution of electron temperature at multiple radial locations, derived from LHD data (blue) and the model (red) for the core temperature rise (R) heat pulse propagation experiment in plasma 49708.Radial locations range from edge(q ¼ 0:703) to core(q ¼ 0:015), where q ¼ r=a, r is the radial co-ordinate, and a $ 0.6 m is minor radius of LHD.Model results match experimental data well from q ¼ 0:450 inwards to the plasma core, especially amplitudes and the time structure of pulse decay.The amplitudes of model time traces increase from edge to core, as in the measured electron temperature profiles.Model results do not fit experimental data outwards from q ¼ 0:546 to q ¼ 0:703, implying that different physics applies in the outer LHD plasma.
FIG. 2. Time evolution of electron temperature at multiple radial locations, derived from LHD data (blue) and the model (red) for the core temperature drop (D) heat pulse propagation experiment in plasma 49719.Radial locations range from edge(q ¼ 0:703) to core(q ¼ 0:015), where q ¼ r=a, r is the radial co-ordinate, and a $ 0.6 m is minor radius of LHD.As in Fig. 1, model results match experimental data well from q ¼ 0:450 inwards to the plasma core.Again, model results do not fit experimental data outwards from q ¼ 0:546 to q ¼ 0:703, reinforcing that different physics dominates in the outer LHD plasma.(10).We have referring also to Table I, where superscripts R and D denote the core T e rise and drop cases, respectively.We note that the empirically determined velocities v R 0 and v D 0 of heat pulse propagation need not necessarily coincide with the optimal value of the mathematical transformation parameter v 0 in Eq. ( 4).Mathematically, we may treat v 0 as a free parameter in Eqs. ( 8) and ( 9), which we have labelled the pseudovelocity.We solve Eqs. ( 8) and ( 9) repeatedly for different values of this pseudo-velocity, see, for example, Fig. 7 for the R case, and identify the best fit value.To test sensitivity with respect to n, three other cases ðn ¼ 5; n ¼ 15; n ¼ 20Þ have been examined.All three of these test options exhibit the same properties as the case of n ¼ 10.Various combinations of initial conditions and pseudo-velocities were also tested, showing that the combination specified above provides the best fit to the experimental data.For the phase plots, see, for example, the right-hand panel of Fig. 7, circulation directions are identical with those of the experimental data for both the R and D cases.
The fitting operation can be assisted by shifting the origin of co-ordinates.The robustness, with respect to variation of initial conditions, of this mathematical approach to identifying pseudo-velocity has been tested.This robustness reflects the structure of the underlying attractor in phase space, projected in ðy 1 ; y 2 Þ coordinates in these two figures; compare also Figs. 6 and 7 of Ref. 25.

IV. CONCLUSIONS
We have derived from first principles a time-dependent model in one spatial dimension, which is able to describe quantitatively the radial inward propagation of heat pulses in the core of two plasmas 24 in the LHD.In one plasma the central electron temperature rises, in the other it falls.This new model is derived from a travelling wave transformation of the zero-dimensional model of Ref. 25, which is known to capture the time-evolution of the heat pulse as it passes a fixed radial location in these two plasmas.From the experimental data, we infer that the velocity of the propagating pulse is constant in both the electron temperature rise and drop cases.The pulse velocity in the electron temperature rise case is smaller than in the drop case by a factor ' 2. This aligns with Eq. (10), which reaches back into the mathematical structure of our model, and also coincides with empirical expectations given the values of the heat conduction parameters for these two plasmas.A pseudo-velocity parameter is introduced in the travelling wave transformations, in order to model heat pulse propagation across spatial location as well as in time.From numerical tests, we discover that real pulse velocity is about two times the best estimate of the travelling-wave transformation parameter v 0 , referred to as the pseudo-velocity.Comparison between model outputs and raw experimental data suggests that our model is able to FIG. 5. Data analysis underpinning calculation of pulse velocity from the experimental data, which requires a statistically robust identification of the time of the pulse peak from the noisy data at each radius 0:450 !q !0:015.Blue lines show timeseries of electron temperature data versus time from the R case (Te49708).Red lines denote timeseries smoothed over a window whose span is 5% of the total sample points, so that approximately 50 sample points generate the moving average.Black dots mark the maximum values, at each radius, of each smoothed timeevolving electron temperature pulse.The width of the horizontal error bars is defined by span of the moving window.The black dash line is the best fit straight line joining all the peaks.From it, we infer the pulse propagation speed, which is nearly independent of radius, to be ð32:6269:89Þ ms À1 .FIG. 6.As Fig. 5, demonstrating the calculation of the pulse velocity for the electron temperature data for the D case.Despite data which are more noisy than in Fig. 5, the pulse propagation velocity calculated from the dashed line is approximately constant across all radii in this region 0:450 !q !0:015 of LHD plasma 49719.We infer a pulse velocity ð53:50620:97Þ ms À1 .
Phys.Plasmas 22, 062308 (2015)   This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: describe heat pulse propagation well, within a broad radial range of the LHD core plasma from r=a ' 0:5 to the centre.The results of the present paper provide additional support to the physical proposals, described in Sec. 2 of Ref. 25, which motivate the simple model equations reproduced as Eqs.(1) to (3) above.Central to these proposals is the conjecture that heat pulses are structures which involve plasma physics processes which are so strongly nonlinear that heat pulse evolution is primarily determined by the reactions of the perturbed heat flux and the perturbed temperature gradient on each other.In the present model, this is the dominant element of nonlinear physics, whereas conventional turbulent transport plays a relatively minor dissipative role.It is this mutually coupled interaction between perturbed heat flux and temperature gradient that governs the local plasma dynamics of the heat pulse in space, equivalent to the local up-and-down dynamics of a water wave under gravity.We have shown in this paper that this coupling model lends itself readily to a travelling wave transformation, yielding spatiotemporal pulse propagation, and that the pulse velocity that emerges mathematically provides an adequate match to empirical results and expectations.This aspect of the analysis also provides guidance on a previously unanswered question, namely, the generic character of the heat pulse: we have shown that it may be closely related to a Korteweg-de Vries-Burgers soliton.Two avenues of investigation would repay immediate attention.First, this model has been tested on, and motivated by, measurements from only two plasmas in LHD-albeit plasmas with exceptionally high quality measurements.The simplicity of the model encourages one to hope that it may be more widely applicable and clearly it should now be tested on a broader range of heat pulse experimental datasets, provided that they possess the required spatial and temporal resolution and that the other relevant plasma parameters are well diagnosed, as in LHD.Second, while gyrokinetic or other computationally intensive transport simulations have not yet (to the authors' knowledge) been applied to heat pulse experiments, the outputs of any future simulations of heat pulses could be tested directly against the analytical model presented here, which is constructed in terms of variables that are directly measurable, as distinct from the first-principles particle and field distributions of gyrokinetic theory.

FIG. 3 .
FIG.3.Time evolution of electron temperature at three specific radii selected from Fig.1for the central temperature rise (R) case, during the heat pulse propagation experiment in LHD plasma # Te49708.Data and model output are denoted by blue and red lines, respectively.

FIG. 7 .
FIG.7.Variation and robustness of nonlinear pulse phenomenology in the model, for three different values of pseudo-velocity v 0 at three radial locations (left panel).Blue lines show experimental data for R case in plasma 49708.Solid, dash, and dot magenta lines denote simulation outputs for v 0 ¼ 15; v 0 ¼ 30, and v 0 ¼ 45, respectively.The boundary condition for R case is y 2 ð0Þ ¼ 1:5, and no horizontal or vertical shift is applied to the model outputs.(Right panel) Phase plot of model outputs from left plots, all of which lie on the same orbit.Circulation direction of this phase plot is identical with Fig.6(a) of Ref.25, where the sign of the horizontal axis is reversed.

TABLE I .
25perimentally inferred parameter values25for both T e rise and T e drop cases.

TABLE II .
Boundary conditions of both T e rise and T e drop cases.
This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: À1and (53.50 6 20.97) ms À1 , respectively, which are comparable in magnitude to measured heat pulse propagation velocities.The ratio 2 between these two velocities is broadly consistent with the value inferred by substituting experimentally measured values for transport coefficients into Eq.
2 for the central temperature drop (D) case, during the heat pulse propagation experiment in LHD plasma # Te49719.Data and model output are denoted by blue and red lines, respectively.This article is copyrighted as indicated in the article.Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 86.151.86.183On: Fri, 26 Jun 2015 15:01:39 for the temperature drop case.These dimensionless pseudovelocity values equate to the dimensional values (32.62 6 9.89) ms