THE INFLUENCE OF NUMERICAL RESOLUTION ON CORONAL DENSITY IN HYDRODYNAMIC MODELS OF IMPULSIVE HEATING

Figure 1 shows the temporal evolution of the spatially averaged density (n), electron temperature (T_e), total pressure (P), and the corresponding T_e versus n phase-space plot for the numerical runs comprising groups 2 (upper four panels) and 3 (lower four panels), as described in Table 1. These correspond to roughly an active region nanoflare and a small flare, respectively. The spatial average was calculated over the uppermost 50% of the loop, which is comprised solely of coronal plasma. Each curve corresponds to a different maximum refinement level for the adaptive grid. Taking group 2 as an example, the lowest solid line in the upper-left panel for Figure 1 shows the average density for a refinement level of zero (no refinement, a uniform grid of 4 × 10⁷ cm grid cells), the dotted line shows the average density for a refinement level of RL = 2, and so forth. The curves for refinement levels of 10 (long dashes) and 12 (upper solid line) are essentially identical, which indicates that a maximum spatial resolution of order 100 m is required for the density solution to converge for this group of runs.

Zoom In Zoom Out Reset image size

It is striking that the average electron temperatures are only weakly dependent on the resolution and converge at RL = 6: even coarser grids perform adequately. The peak temperature is similar for all resolutions, though underresolved loops cool somewhat more rapidly than the resolved ones. Since the peak temperature is reached when heat input and coronal energy loss by thermal conduction are in balance (Bradshaw & Cargill 2006), one would expect this to require adequate spatial resolution of the temperature gradients. However, we are only considering the average coronal temperatures here and even our coarsest grid is sufficient to resolve coronal scales (of order 1000 km), which is all that is needed for thermal conduction to efficiently remove the excess energy deposited in the corona. This shows that inadequate spatial resolution does not "bottle up" the heat flux and thereby slow coronal cooling.

On the other hand, the density shows a major difference between low and high resolutions. For these two cases with small τ_H and L, the density reaches a peak a few tens of seconds after the temperature and then oscillates as the plasma sloshes to and fro within the loop. For both cases in Figure 1, there is a difference of 80% between the initial peaks for RL = 0 and RL = 12, and a larger difference in the density decay phase, with the poorly resolved cases also showing little in the way of oscillations. Group 2 has fairly moderate heating with a peak electron temperature at just above 6 MK, and yet at 500 s there is a factor of greater than three difference in the density between the coarsest and most refined solutions. Even with 25 km grid cells (RL = 4) the density is about 40% lower than the converged value. Group 3, for which an order of magnitude more energy was released over the same timescale (60 s), shows greater differences for different values of RL than group 2, though the runs for RL = [10, 12] again demonstrate successful convergence. Indeed, this group of runs demonstrates that as one approaches flare territory the demands made of spatial resolution become more extreme and the need to satisfy them to derive accurate solutions becomes more critical. It is evident, for example, that grid cells greater than 1 km in width are not sufficient. The final two columns in Table 1 show the resolution required for the density to exceed 90% and 75% of its fully resolved value. These numbers should be treated as approximate, and are evaluated around the second and third density peaks.

The magnitude of the density oscillations is also diminished for poorer resolution, indicative of faster damping. These should damp conductively and we find that while the damping time of the resolved cases is broadly consistent with the calculated conductive damping time, the fast damping of the unresolved cases is not. The conductive timescale scales as n/T^5/2, but the lower density in the unresolved cases is compensated by a somewhat smaller temperature so that the conductive timescale changes little. So it appears that the lack of resolution also enhances the damping, a result of importance for coronal seismologists who may use 1D HD models to interpret their results (for a further discussion of wave damping in the context of coronal seismology, see Morton et al. 2010).

Figure 2 and Table 1 show the set of results for the runs comprising group 9, another very impulsive heating case corresponding to a strong nanoflare in a long (180 Mm) loop. The pattern followed by the density and temperature evolution, and seen in the phase-space plot, is broadly the same as in Figure 1. However, even though this group is subject to the most extreme volumetric heating rate for the 180 Mm runs, the density difference (e.g., at 1000 s) between the coarsest and most refined grids is only about a factor of two. The reason for this is that longer loops require shallower temperature gradients and consequently the coarser grids are more equal to the task of providing adequate resolution. For example, RL = 6 does a reasonable job in this case, though one should bear in mind that this still corresponds to a grid cell width of just over 5 km for a loop of 180 Mm total length, reaching a peak electron temperature in the region of 12 MK. (We note from the phase-space plot that the RL = [0, 2] curves do not recover the initial equilibrium, but find another nearby in temperature but of about half the density. The curves for RL > 2 can all be seen to recover the initial equilibrium. The reason for this is that the initial conditions fed to HYDRAD are calculated on an extremely high-resolution grid with cell widths as small as 1 cm; these are then interpolated onto the much coarser grid used for the time-dependent calculation, which is then adaptively refined as the calculation proceeds. The grids of RL = [0, 2] are not sufficient to resolve even the initial conditions in the case of group 9 and so these runs ultimately return to a different equilibrium state. The grids of RL > 2 do not have this problem and converge on the initial equilibrium.)

Zoom In Zoom Out Reset image size

Figures 3 and 4 show the moderate and extreme volumetric heating cases for loops of total length 60 (groups 5 and 6) and 180 Mm (groups 11 and 12), respectively, with the energy released over an extended period of 600 s, more akin to flaring timescales. Order of magnitude differences between the density curves for RL = [0, 12] can be seen for group 6. Even at refinement levels of RL = 6 (≈6 km grid cells) there are differences of 30% (group 5) and 20% (group 6) between the peak densities attained and the fully converged peak densities. The peak density also shifts to earlier times by a couple of hundred seconds with changing RL, which is particularly noticeable for group 5. The density oscillations seen in the earlier resolved results have also disappeared for groups 5 and 6 due to the short conductive damping timescale.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The differences between the density curves for groups 11 and 12 are less severe, indicating once again that the shallower gradients associated with longer loops, for a given peak temperature, can be satisfactorily handled by coarser grids on the order of at least km resolution (cell widths of 6 km in the case of group 11 and about 2 km for group 12). The phase-space plots corresponding to groups 5 and 6, in particular, provide rather dramatic examples of the need to provide adequate spatial resolution in order to accurately follow the complete coronal heating and cooling cycle. The density oscillations have reappeared due to the increase in the conductive damping time in the longer loops and there is a little evidence for the oscillations persisting in the unresolved long loops.

In summary, analysis of groups 1–12 shows the following: (1) shorter loops require greater spatial resolution than longer loops for a given peak temperature, (2) the coronal density is far more strongly resolution dependent than the temperature, (3) grid cell widths of less than 1 km are required for fully converged solutions in the case of 60 Mm loops and peak temperatures greater than 6 MK, (4) cell widths of no more than 5 km are required for reasonably accurate solutions (densities within a few percent of converged solutions) in the case of 180 Mm loops and peak temperatures exceeding 6 MK, and (5) obtaining a physically realistic evolution through the complete coronal heating and cooling cycle demands adequate spatial resolution at all times.

We can thus conclude that (1) the mechanism responsible for driving the ablation of material into the corona during heating is strongly inhibited by a lack of spatial resolution, and we will examine the physics and the numerics of this process in greater detail in Section 3.3. The phase-space plot tracks the solutions all the way back to the recovery of the initial conditions and reveals significant differences in the evolutionary paths taken by the individual runs. (2) The underresolved cases have rapid damping of the density oscillations. An analysis of the density and temperature in the RL = 0 and RL = 12 cases shows that there is little change in the conductive cooling times, with the smaller values of T and n in the unresolved case offsetting each other. We discuss the reason for this rapid damping in the underresolved cases in Section 3.3.

The first question that needs to be addressed in understanding these results is whether the low coronal density for poorly resolved cases is caused by a lack of energy conservation due to numerical underresolution. Figures 1–4 show the temporal evolution of the average total pressure (electron + ion) calculated in the same way as the average density and electron temperature. The total pressure can be considered a proxy for the total thermal energy in the loop, and the amount of energy ultimately transported into the corona as the system evolves, since the kinetic energy component is in general relatively small. It is clear that the total pressure, and hence the total energy and energy transported, is strongly dependent on the spatial resolution. The total energy released is the same within each group of runs, so an increasing fraction of this energy does not appear in the corona as the spatial resolution is decreased. To take an extreme example, the lower four panels in Figure 3, corresponding to group 6, show differences in the total pressure of a factor of almost six between RL = 0 and RL > 8. There are also significant differences between intermediate refinement levels.

To assess the ability of the code to conserve energy, we check whether the following condition is satisfied:

$$\begin{equation} \int _0^{2L} \frac{\partial E}{\partial t} - [-E_R + E_H - \rho v g_{||}] ds = 0, \end{equation} \tag{ 1 }$$

where E is the sum of the thermal and kinetic energy, E_R and E_H are the volumetric radiation and heating rates, respectively, and ρvg_|| is the rate of change of gravitational potential energy. The other terms involved in the energy equation that determine ∂E/∂t are flux divergencies (e.g., enthalpy and thermal conduction) which spatially integrate to zero over the entire loop due to there being no mass or conductive flux through the lower boundaries. Thus, they serve only to redistribute energy within the system. Therefore, if all of the energy sources, sinks, and fluxes are correctly treated and accounted for then the condition described by Equation (1) will be satisfied.

We have calculated the deviation from zero of Equation (1) as a percentage of the total energy stored in the transition region and corona (T > 0.02 MK) for the most extreme cases, which are groups [3, 6, 9, 12]. We chose to omit the chromospheric part of the loop from the calculation in order that the large energy content of the chromosphere could not mask any significant lack of conservation in the overlying atmosphere. Between t = 0 s and the time of peak density for each group of runs (as seen in Figures 1–4), we found the maximum percentage deviations from perfect conservation to be: group 3, 0.36% (RL = 0), 0.20% (RL = 12); group 6, 0.69% (RL = 0), 0.09% (RL = 12); group 9, 0.05% (RL = 0), 0.12% (RL = 12); group 12, 0.31% (RL = 0), 0.26% (RL = 12). In addition, the temperature versus density phase-space diagrams in Figures 1–4 show that following the conclusion of the impulsive heating, when only the background heating used to find the initial hydrostatic conditions remains, the code does a good job of recovering the initial conditions by returning the solution to the starting point in the phase space. Thereafter, the atmosphere is held in hydrostatic equilibrium with only residual background flows having a Mach number M ≪ 1.

We find that energy is conserved to better than 1%, even for RL = 0, and so the large discrepancies in the total coronal thermal energy indicated by Figures 1–4 are not due to a lack of spatial resolution directly (because the code is evidently doing a good job of accounting for all energy sources, sinks, and fluxes), but the lack of spatial resolution may be driving the solutions into different physical regimes (e.g., strong versus weak ablation), as we now discuss. Whether energy equations implemented in 3D codes that do not use the conservative form do so well is unclear.

We now turn to a consideration of the different physical regimes that may dominate the numerical solutions for the range of RL studied here. Since the major differences between these runs are a consequence of the transition region resolution then it makes sense to focus our investigation on this part of the domain. Furthermore, we demonstrated in our earlier work (Bradshaw & Cargill 2010a) that it is the demands of the transition region which dominate the energetics and the dynamics of the loop during all phases of its evolution. Figure 5 shows the average heat flux (left-hand column) and the average enthalpy flux (right-hand column) through the transition region for groups 3, 6, 9, and 12. They were calculated by spatially integrating the heat flux from the location of the transition region/corona interface down to the footpoint (the first grid cell in which the temperature reaches or falls below 20,000 K) and then dividing by the transition region width. The location of the transition region/corona interface can be found in a reasonably straightforward way during this period by identifying the grid cell in which the divergence of the heat flux transitions from positive (corona) to negative (transition region), indicating that energy is being removed from the heat flux to power the transition region radiation. These average quantities were calculated for the right-hand leg of the loop and so positive fluxes are downflowing (toward the footpoint) and negative fluxes are upflowing (toward the apex). Their temporal evolution is shown for the entire heating phase and part of the ablative phase (coronal cooling by thermal conduction). At later times, when the ablated material from opposite legs of the loop collides at the apex and triggers acoustic waves (the density oscillations already noted, with accompanying localized temperature perturbations) then this interface becomes far more difficult to locate.

Zoom In Zoom Out Reset image size

An examination of the left-hand column of Figure 5 shows that the ability of the transition region to support a substantial heat flux is almost independent of the spatial resolution, particularly in the case of the groups of runs for the longer loops (groups 9 and 12) where the gradients are shallower. Recall that these results are for the most extreme heating scenarios chosen for our study. Even in the case of short loops (groups 3 and 6) the heat fluxes for the RL = [6, 12] runs are in very close agreement, given the orders of magnitude difference in the spatial resolution (a factor of 64). Once again, the discrepancies between the coronal density and the total energy for different values of RL are not due to any significant heat flux limiting caused by inadequate spatial resolution in the transition region.

The right-hand column of Figure 5 reveals a clue to what is actually behind the discrepancies. The common assumption is that the heat flux drives chromospheric ablation, where the excess energy carried by the heat flux that is not radiated by the transition region heats the chromosphere (a poor radiator), which then expands upwards and carries with it an enthalpy flux in the opposite direction to the heat flux. However, the return enthalpy flux is severely diminished at low transition region spatial resolution, which is indicative of an alternative scenario. If it were purely chromospheric ablation then, since most of the heat flux is accounted for even at small RL, the resulting ablation and enthalpy flux should reflect the agreement between the heat fluxes for different values of RL. Evidently, this is not the case and a significant amount of energy is somehow lost before it is converted into enthalpy at low spatial resolution.

The true physics of ablation is as follows. At low spatial resolution, where the transition region is effectively a discontinuity, the heat flux jumps from the corona to the chromosphere where the incoming energy is strongly radiated (the chromosphere is an inefficient radiator relative to its density, but it is dense), leaving little left over to heat the plasma and drive ablation. At high spatial resolution, the heat flux passes through the transition region in a series of steps and at each step some energy is radiated, some goes into local heating and expansion/ablation, and the remainder continues on. More and more energy is radiated with increasing depth in the transition region since the volumetric loss rate scales as n². Ultimately, depending on the magnitude of the initial heat flux, a proportion of the original heat flux will make it to the chromosphere, but a significant contribution to the upward material and enthalpy flux will have been made by transition region plasma. This contribution is neglected in low-resolution runs, with the resultant discrepancies between the coronal density and the total energy. It may also partially explain why observational evidence of the upflows that must somehow fill the corona has been so difficult to provide; there is no single, dominant source of plasma and instead the corona is filled in a more piece-meal fashion.

Figure 6 confirms this conclusion. The upper-left panel shows the electron temperature profile in the footpoint region for the RL = [0, 6, 12] runs of group 3, approximately halfway through the ablation phase. Group 3 has been selected for detailed analysis because it imposes very severe constraints on the spatial resolution (see Table 1 and Figure 1). The temperature gradient is extremely poorly resolved in the RL = 0 case to the extent that there are no data points for temperatures of O(10⁵) K. The upper-right panel shows that for RL = 0 a large amount of energy is deposited by the incoming heat flux in the grid cell at s = 6.2 × 10⁸ cm, corresponding to T = 10^4.5 K, and then it steeply falls away. The lower-right panel shows strong radiation from the same grid cell and the lower-left panel shows that at this critical grid cell there is no energy left over to drive ablation. For RL = [6, 12], it is clear that there is a net heating throughout the region of the steep temperature gradient and so energy is available to ablate the increasingly dense plasma found there, driving it upward and leading to enhanced coronal densities compared with RL = 0. The transition region/corona interface is properly described as the location where the heat flux transitions from a sink (the corona) to a source (the transition region) as defined by Bradshaw & Cargill (2010a). The upper-right panel of Figure 6 shows that this extends beyond s = 10⁹ cm for the RL = [6, 12] runs (actually to 1.5 × 10⁹ cm) and to about 9 × 10⁸ cm for the RL = 0 run. The steep temperature gradient may then be thought of as the lower transition region and, since the energy deposited by the heat flux in the upper transition region is quite similar for RL = [0, 6, 12], then the dominant source of coronal plasma is the lower transition region. It is also interesting to note that underresolving the transition region evidently leads to an underestimate of its spatial scale, despite the similarities in the time evolution of the average temperature between the runs at this time. Other groups with severe resolution constraints show the same result. Thus, the mechanics of ablation is far more dependent on resolving the lower transition region and treating its radiative properties correctly than being overly concerned with the energetics of the chromosphere. This is the underlying reason why relatively simple treatments of the chromosphere, particularly in flare modeling, can yield similar results to those that are more sophisticated; the chromosphere does not make the dominant contribution. Of course, the more sophisticated treatments are necessary for detailed spectral modeling of chromospheric emission.

Zoom In Zoom Out Reset image size

The nature of the interaction between the corona and the chromosphere, mediated by the transition region, also explains why the density oscillations, caused by ablated material sloshing back and forth in the corona, are rapidly damped in the underresolved cases. When there is inadequate spatial resolution the wave energy falls over the transition region cliff to be radiated away by the chromosphere, whereas in the well-resolved cases the effective refractive index of the transition region is greater, leading to more of the wave energy being reflected from the steep density gradient and back into the corona.

We have demonstrated in this paper that one can run 1D HD models using highly refined adaptive grids (tens of meters) for several hours real (solar) time. The situation with multi-dimensional MHD models is more complicated for several reasons:

• 1.  
  A limit on resolution is imposed by the scale of the simulation being run. A contemporary massively parallel code such as Bifrost (Gudiksen et al. 2011) can run something of the order 512³ cells (depending on available resources) and the grid can be stretched in one dimension, permitting a fixed concentration of points where any transition region is expected. (Adaptive grids are more complex to implement in 3D MHD since there are many places that one might require high resolution, such as current sheets, in addition to the transition region.) As examples of contemporary 3D codes, Gudiksen & Nordlund (2005a) have cells of 150 km in the vertical direction, Bingert & Peter (2011) 230 km, which was decreased subsequently by Zacharias et al. (2011) to 75 km. These correspond to RL = 2–4. All these papers discuss relatively cool loops, with temperatures of order 1 MK, and so are closest to our group 1 (weakest heating for short loops) where we find RL = 4 is needed (Table 1). This suggests that the resolution is adequate for these particular cases, though even marginally higher temperatures will result in inadequate resolution. The Bifrost code solves for thermal conduction on a sub-grid with the potential for good resolution equivalent to RL = 12; however, at the present time a resolution of up to 25 km is used (B. V. Gudiksen 2013, private communication): we discuss other aspects of their implementation technique below.
• 2.  
  Other solutions have been proposed in order to sidestep the need for high resolution, such as an artificial reduction of κ₀. This would appear to permit significantly greater time steps, but by mitigating the cooling role of conduction, higher coronal temperatures are maintained, which will tend to push the time step back to shorter values. Gudiksen et al. (2011) also propose lowering the thermal conduction coefficient, but comment that "the results must be carefully analyzed." Linker et al. (2001) proposed an alternative formalism for κ(T):

  $$\begin{equation} \kappa (T) = \kappa _0 \left[ s T^{5/2} + (1 - s) T^\alpha T_{{\rm mod}}^{5/2 - \alpha } \right], \end{equation} \tag{ 2 }$$

  where

  $$\begin{equation} s(T) = \frac{1}{2} \left[ 1 + \tanh \left(\frac{T - T_{{\rm mod}}}{\Delta T_{{\rm mod}}} \right) \right]. \end{equation} \tag{ 3 }$$

  Here, T_mod is a typical transition region temperature (0.2 MK say), and above (below) this temperature the heat flux approaches the Spitzer value (is decreased). For typical active region temperatures, solution of the static equilibrium loop equations leads to a decrease of order two in the coronal density, as also noted by Linker et al. (2001). Subsequently, Lionello et al. (2009) argues that if the radiative losses are also diminished such that κ(T)R_L(T) is unchanged, one can recover proper coronal densities: indeed for static loops, we find the difference is <1%. However, it is unclear how successfully this formalism translates into dynamic evolution. The increase in coronal density after heating arises precisely because the transition region is not in equilibrium: the downward heat flux exceeding the radiative losses, leading to an upward enthalpy flux (see Klimchuk et al. 2008; Cargill et al. 2012). In the decay of impulsive events, the coronal properties are governed by a downward enthalpy flux responsible for the transition region radiation: conduction plays a minimal role (Bradshaw 2008; Bradshaw & Cargill 2010a, 2010b).In addition, the physical definition of the transition region is such that the transition region/corona interface is located where the supporting flux transitions between a sink and a source (Bradshaw & Cargill 2010a, 2010b). The transition region is not simply anywhere that T < 1 MK. The characteristic temperature of the transition region can change dramatically during the heating and cooling cycle; from several MK in hot AR cores and during flares, to below 1 MK during the later stages of cooling. It is therefore unclear how one should choose the quantity T_mod at each stage of the numerical calculation. Studies are clearly required to validate this approximation: its benefits are considerable in terms of the run time of a 3D code and even a 1D code (e.g., for flares); its disadvantages are that for any transition region spectral line properties are unlikely to be relevant.
• 3.  
  Another issue that may arise with 3D MHD codes is the reluctance to use the conservative form of the energy equation. This arises because in a conservative MHD code, one calculates the plasma pressure by subtracting the magnetic field energy density (and kinetic energy density) from the total energy density: for a low β plasma, typical of the corona, the pressure is then the difference between two large numbers, leading to a possibility of negative pressure. There are in fact forms of the energy equation that get around this (see Cargill et al. 2000), but these run into difficulties in modeling shocks. In underresolved non-conservative simulations, it is unclear whether the lack of conservation of energy is a problem. Some benchmarking is required.
• 4.  
  Some 3D codes use operator splitting, or time-splitting, in the interests of computational efficiency (e.g., not solving the MHD equations on a conductive timescale). Here, the various parts of the energy equation are solved separately. For example, Tam et al. (2013) solve the conduction and radiation parts of the equation separately, and then feed the new temperature back into the MHD equations. Gudiksen et al. (2011) solve the conduction equation implicitly on a fine grid, and then generate a smoothed solution to feed into the MHD equations. It is unclear to us whether such procedures guarantee, for example, the correct interaction between a downward conductive flux, and the response through an upflow. For the benign problems discussed in Bingert & Peter (2011) and Gudiksen et al. (2011), there may not be an issue. But a benchmarking of such procedures seems to be needed.To close this part of the discussion, we recommend that those using 3D MHD codes to study the interaction between corona, transition region, and chromosphere carry out a set of "resolution tests" as we have done here with more extreme coronal temperatures than those published. This can be accomplished by reducing the 3D code to a 1D one, while keeping the overall 3D grid structure, with the remaining coordinate aligned along a straight magnetic field. Such a comparison will provide confidence in how the 3D codes address this important aspect of modeling the solar atmosphere.

Our numerical experiments have employed a fully explicit 1D HD code, which forces us to resolve the principal physical timescales of the system at each stage of its evolution. Other workers have resorted to implicit solvers for their 3D MHD codes (e.g., Botha et al. 2011; Gudiksen et al. 2011; Tam et al. 2013). The advantage of, for example, an implicit Crank–Nicolson scheme is that it is guaranteed to be numerically stable. Some sort of iteration is undertaken in multi-dimensions (e.g., a successive overrelaxation (SOR) scheme; Botha et al. 2011) with a criterion for convergence being a small discrepancy in the temperature (either averaged or pointwise) between iterations. This can slow down an implicit scheme if the number of iterations is large, or, worse, if convergence never occurs.

Such implicit schemes are useful for a slowly evolving system or one seeking equilibria, where long wavelength behavior dominates, and one is not concerned with temperature fluctuations over one or two grid points. However, for general applications, one must always bear in mind the various physical timescales that are of interest. The most important one here is the time it takes a conduction front to cross a grid cell, where the cell is narrow enough to resolve any temperature gradient. For a dynamic situation, one needs a "Courant condition" for a conduction front: Δt(implicit) < min(Δs/V_c), where V_c is the velocity of the conduction front. For an analytic solution such as the Zel'dovich thermal wave, V_c is known. However, in a simulation, one needs to calculate both Δs and V_c in real time. The former is accomplished through an adaptive grid. The latter is more challenging. An alternative way to impose this condition would be to require that the change in temperature in any cell be limited to some reasonably small fraction (e.g., 20%), and adjust the time step accordingly. Such a secondary iteration (in addition to the one needed to solve the diffusion equation itself) would need careful implementation to avoid crippling computational costs.

Besides the obvious implications for calculating the response of the lower layers of the atmosphere to heating and the subsequent energy transport, these results are collectively significant for spectral modeling and emission measure calculations. The differences between the electron temperature curves, though small, may lead to very different predictions for ionization states and spectral line intensities because the temperature-dependent contribution functions tend to be very sharply peaked. Consequently, small temperature differences matter. The much greater differences between the density curves are even more problematic. The emission measure and spectral line intensity scale as n², so differences are nonlinearly amplified, and time-dependent calculations of the ionization state are likely to be strongly affected due to the density dependence of the collisional ionization/recombination timescales. It is not difficult to identify the potential for forward modeling studies, in which numerical predictions are compared with real observational data, to yield misinterpretations and erroneous conclusions.

The evolution of the HYDRAD code used in this paper has been described in various publications (Bradshaw & Mason 2003a, 2003b; Bradshaw & Cargill 2006; Bradshaw & Klimchuk 2011). Given its now widespread use, it is useful to describe all the developments of HYDRAD and the current capabilities of the code in one location. HYDRAD solves the 1D, multi-fluid HD equations for a full loop of arbitrary geometry. The mass and momentum equations reduce to the following conservative form:

$$\begin{equation} \frac{\partial \rho }{\partial t} + \frac{\partial }{\partial s}\left(\rho v\right) = 0, \end{equation} \tag{ A1 }$$

$$\begin{eqnarray} &&\frac{\partial }{\partial t}\left(\rho v\right) + \frac{\partial }{\partial s} (\rho v^2) = -\frac{\partial }{\partial s} (P_e + P_i) + \frac{\partial }{\partial s}\left(\frac{4}{3} mu_i \frac{\partial v}{\partial s}\right) - \rho g_{||}, \nonumber\\ \end{eqnarray} \tag{ A2 }$$

where current and charge neutrality have been imposed, t and s are time and the field-aligned spatial coordinate, ρ is the ion mass density (ρ = m_in: m_e ≪ m_i), v is the bulk velocity of the flow, and g_|| is the field-aligned gravitational acceleration. Here, m_i is the average ion mass, and accounts for the fact that while hydrogen is by far the most abundant element the presence of heavy metals has a non-negligible influence on the mass density and the bulk flow. Ion viscosity is included using the classical Spitzer viscosity coefficient. Electron viscosity is negligible. There are separate energy equations for the electrons (subscript "e") and ions (subscript "i"):

$$\begin{eqnarray} && \frac{\partial E_e}{\partial t} + \frac{\partial }{\partial s} [(E_e + P_e)v] = -\frac{\partial F_{ce}}{\partial s} + v\frac{\partial P_e}{\partial s} \nonumber\\ &&\quad + \frac{k_B n}{\gamma -1}\nu _{ie}(T_i-T_e) - R + H, \end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray} \frac{\partial E_i}{\partial t} + \frac{\partial }{\partial s} [(E_i + P_i)v] &=& -\frac{\partial F_{ci}}{\partial s} - v\frac{\partial P_e}{\partial s} + \frac{k_B n}{\gamma -1}\nu _{ie} (T_e-T_i)\nonumber\\ && + \frac{\partial }{\partial s}\left(\frac{4}{3} mu_i v \frac{\partial v}{\partial s}\right) + \rho v g_{||}, \end{eqnarray} \tag{ A4 }$$

where

$$\begin{equation} E_e = \frac{P_e}{\gamma -1},\quad E_i = \frac{P_i}{\gamma -1} + \frac{1}{2}\rho v^2 \end{equation} \tag{ A5 }$$

and closure is imposed by

$$\begin{equation} P_e = k_B n T_e,\quad P_i = k_B n T_i, \end{equation} \tag{ A6 }$$

where T_i and T_e are the ion and electron temperatures, and k_B is Boltzmann's constant. The first term on the right-hand side of Equations (A3) and (A4), represents the thermally conducted flux carried by each particle species. The electron flux is written as

$$\begin{equation} F_{ce} = - \kappa _{0e} T_e^{5/2} \frac{\partial T_e}{\partial s} \end{equation} \tag{ A7 }$$

and the ion flux in an analogous form. The Spitzer coefficient for the electron thermal conduction is κ_0e = 7.8 × 10⁻⁷ and we assume that protons dominate the ion thermal conduction, giving κ_0i = 3.2 × 10⁻⁸. A flux limiter is built into the numerical code to ensure that the free-streaming limit is not exceeded and is discussed below (see also Patsourakos & Klimchuk 2005; Bradshaw & Cargill 2006). The second term on the right-hand side of Equations (A3) and (A4), P_e, represents the electric potential, the energy given to/removed from the particle species by the electric field when the separation between them becomes too large/small. The third term represents the thermal equilibration between the particle species via collisions, where ν_ie is the Coulomb collision frequency:

$$\begin{equation} \nu _{ie} = \frac{16\sqrt{\pi }}{3} \frac{e^4}{m_e m_i} \left(\frac{2k_B T_e}{m_e}\right)^{-{3}/{2}} n \left(\ln \Lambda \right), \end{equation} \tag{ A8 }$$

where ln Λ is the Coulomb logarithm. The gravitational potential energy has been included for the ions and neglected for the electrons. R is the optically thin radiative emission, which can be calculated using a set of piecewise power laws, a fit to radiative loss curves calculated by the Chianti atomic database for equilibrium ionization or a full non-equilibrium treatment of the radiating ions (thermal bremmstrahlung is also included in the calculation). R can also be calculated for optically thick conditions in the lower atmosphere using the treatment given by Carlsson & Leenaarts (2012) and in this case the presence of neutral species is accounted for (e.g., number density of neutrals, ions, and electrons, inter-species collisions, thermal conduction of neutral hydrogen). However, in the present work the chromosphere is maintained at a uniform temperature (20,000 K) by reducing the optically thin radiative losses to zero over a specified temperature interval (100 K) above the chromospheric temperature (Klimchuk et al. 1987). H is a parameterized, volumetric heating function that depends on the field-aligned coordinate and/or time.

The transport terms of the HD equations are calculated using Barton's method for monotonic transport and a staggered leapfrog algorithm to perform the time integration. Centrella & Wilson (1984) provide a detailed discussion of Barton's method, which is a simpler and more rugged version of the method devised by van Leer (1977), and which works well even on highly discontinuous HD problems. For a uniform grid, the method is second order in space away from the discontinuities and second order in time. At discontinuities, it reduces to first order in space to maintain fidelity to the nature of the solution in that region. All other ordinary and source/sink terms are spatially differenced and time integrated to second order using standard Runge–Kutta methods (Press et al. 1992).

For a uniform grid, the conductive flux F_c is calculated at the edge of cell j by adopting the following discretization:

$$\begin{equation} F_{c,j-1/2} = - \kappa _0 T_{j-1/2}^{5/2} \frac{T_j - T_{j-1}}{\Delta s_j}, \end{equation} \tag{ A9 }$$

where

$$\begin{equation} T_{j-1/2}^{5/2} = [(T_j + T_{j-1})/2]^{5/2}. \end{equation} \tag{ A10 }$$

The volumetric rate of energy change in the grid cell due to thermal conduction is then given by

$$\begin{equation} {-} \frac{\Delta F_{c,j}}{\Delta s_j} = {-}\frac{F_{c,j+1/2} - F_{c,j-1/2}}{\Delta s_j}. \end{equation} \tag{ A11 }$$

For a non-uniform grid, the conductive flux F_c is calculated at the edge of cell j in such a way as to limit the potential loss of spatial order in the derivatives that is known to arise when there are variations in the widths of adjacent grid cells (Oran & Boris 2000). The temperature T_{j − 1/2} is first found by interpolation using a four-point stencil (j − 2, j − 1, j, j + 1); however, if this third-order interpolation introduces a new local maximum or minimum then T_{j − 1/2} is calculated using a linear interpolation (j − 1, j). Using the same four-point stencil, the temperature T* is found at the position s_j − Δs_j (again, reduced to a linear interpolation if a new local max. or min. is introduced). The heat flux at the cell boundary is then given by

$$\begin{equation} F_{c,j-1/2} = - \kappa _0 T_{j-1/2}^{5/2} \frac{T_j - T^*}{\Delta s_j}. \end{equation} \tag{ A12 }$$

This procedure is effectively equivalent to interpolating the temperature onto a locally uniform grid, centered at j − 1/2, where the cell width is Δs_j. Note that this procedure is used by HYDRAD to calculate all spatial gradients (except for the transport terms) on a non-uniform grid (e.g., dP/ds, etc.).

An important aspect of HYDRAD is its adaptive grid, which is flexible enough to increase the numerical resolution wherever needed and so it can resolve as many discontinuities and/or steep gradients, for example, as arise in the calculation at any time. Restriction (refinement) is carried out according to user-defined tolerances (e.g., greater than 10%) on cell-to-cell differences between some or all (again, user definable) of the quantities ρ, E_e, and E_i. When HYDRAD detects that a grid cell needs to be restricted then it is split into two. The quantities ρ, ρv, E_e, and E_i are calculated in the leftmost of the two new cells by interpolation from the original cells j − 2, j − 1, j, and j + 1. Similarly, these same quantities are calculated in the rightmost of the two new cells by interpolation from the original cells j − 1, j, j + 1, and j + 2. Cell j is then removed from the linked list that connects the grid cells and the two new cells are inserted in its place. The two new cells are also tagged with a unique identifier for each refinement level so that only the same two cells that resulted from the splitting of a single cell can be merged to form a single cell under prolongation. The user can specify the number of times that a grid cell can be split and this is the maximum refinement level. When HYDRAD detects that two grid cells (j and j + 1) with the same unique identifier satisfy the criteria to be merged into a single cell (e.g., cell-to-cell differences between quantities falls to less than 5%) then the prolongation operator calculates ρ, ρv, E_e, and E_i in the new cell by interpolation from the original cells j − 1, j, j + 1, and j + 2. Cells j and j + 1 are then deleted from the linked-list and the new cell inserted in their place. The user can select linear interpolation for the restriction and prolongation operators if they expect that sharp features such as shocks and/or steep gradients might arise in their solutions. In this case, just the inner two points of the stencils described above are used. The four-point polynomial interpolation routine also checks for new local minima or maxima and reduces to a linear interpolation if necessary. The adaptive grid algorithm employed by HYDRAD limits the variation in the widths of adjacent cells to a single refinement level as a precaution against the loss of order in calculating spatial derivatives on non-uniform grids (Oran & Boris 2000).

While Equations (A9) and (A10) are a conventional way to finite difference the conduction term, they are not unique. One option is to write Equation (A10) as

$$\begin{equation} T_{j-1/2}^{5/2} = \left[ \left(T_j^{5/2} + T_{j-1}^{5/2} \right) \left/\right. 2 \right. \end{equation} \tag{ A13 }$$

A second is to define the variable u = T^7/2 so that the conduction term is

$$\begin{equation} F_c = - \frac{2}{7} \kappa _0 \frac{\partial ^2 u}{\partial s^2} \end{equation} \tag{ A14 }$$

and finite difference this in the usual three-point manner. We have compared these with Equations (A9) and (A10) by solving the Zel'dovich problem of a propagating heat front (Zel'dovich & Raizer 1967; Reale 1995). For a resolved front, the solution agrees well with the analytical one for all three methods. For an unresolved heat front (i.e., a case where the major temperature jump at the leading edge of the front is over one cell), the three finite difference schemes give heat fluxes at the leading edge that differ by up to a factor two. While our tests show that this has little effect on the propagation of the conduction front, it is worth bearing in mind when more complex problems are solved. Also, it should be noted that the mathematically most satisfactory formalism of conduction, involving the variable u, cannot be used in implicit numerical schemes.

Heat flux saturation becomes important in many of our experiments, especially at the most energetic extremes of the heating parameter space that we explore. The flux-saturation (or free-streaming) limit describes the maximum conducted heat flux that the plasma is capable of supporting. A reasonable estimate of this limit can be obtained by assuming that it is reached when all of the particles travel in the same direction at the electron thermal speed. Bradshaw & Cargill (2006) give the saturated heat flux as

$$\begin{equation} F_{cs} = \frac{3}{2\sqrt{m_p}} n_p (k_B T_p)^{{3}/{2}}, \end{equation} \tag{ A15 }$$

where the subscript p denotes the particle species (electrons or ions). More accurate kinetic models (e.g., Fokker–Planck calculations) suggest that Equation (A15) should be multiplied by a factor 1/6. To ensure a smooth transition between heat flux regimes (saturated and unsaturated) in HYDRAD, we use the form

$$\begin{equation} F_c = \frac{F_{{\rm classical}} F_{cs}}{\sqrt{F_{{\rm classical}}^2 + F_{cs}^2}}, \end{equation} \tag{ A16 }$$

where F_classical is the heat flux calculated in the Spitzer–Härm approximation (e.g., Equation (A9) above) and F_cs is evaluated at the cell boundary, assuming that the density n remains approximately constant across the grid cell. West et al. (2008) found that heat flux saturation does not arise in loops that have reached (or are close to) hydrostatic equilibrium, but it does seem essential to include in all dynamic coronal modeling codes.