A Fast and Accurate Method to Capture the Solar Corona/Transition Region Enthalpy Exchange

TL;DR

This paper introduces TRAC, a new adaptive thermal conduction method that accurately models the solar corona and transition region interactions, significantly reducing errors and computation time in simulations.

Contribution

The paper presents TRAC, a novel adaptive thermal conduction approach that resolves the transition region more efficiently and accurately in solar corona simulations.

Findings

TRAC reduces coronal density errors to less than 3% with coarse resolutions.

The method is three orders of magnitude faster than fully resolved models.

TRAC maintains high accuracy comparable to fully resolved simulations.

Abstract

The brightness of the emission from coronal loops in the solar atmosphere is strongly dependent on the temperature and density of the confined plasma. After a release of energy, these loops undergo a heating and upflow phase, followed by a cooling and downflow cycle. Throughout, there are significant variations in the properties of the coronal plasma. In particular, the increased coronal temperature leads to an excess downward heat flux that the transition region is unable to radiate. This generates an enthalpy flux from the transition region to the corona, increasing the coronal density. The enthalpy exchange is highly sensitive to the transition region resolution in numerical simulations. With a numerically under-resolved transition region, major errors occur in simulating the coronal density evolution and, thus, the predicted loop emission. This Letter presents a new method that addresses the difficulty of obtaining the correct interaction between the corona and corona/chromosphere interface. In the transition region, an adaptive thermal conduction approach is used that broadens any unresolved parts of the atmosphere. We show that this approach, referred to as TRAC, successfully removes the influence of numerical resolution on the coronal density response to heating while maintaining high levels of agreement with fully resolved models. When employed with coarse spatial resolutions, typically achieved in multi-dimensional MHD codes, the peak density errors are less than 3% and the computation time is three orders of magnitude faster than fully resolved field-aligned models. The advantages of using TRAC in field-aligned hydrodynamic and multi-dimensional magnetohydrodynamic simulations are discussed.

Equations

$$L_T(T(s))=\frac{T}{|dT/ds|},$$

$$L_R(s)=\Delta s,$$

$$\displaystyle T_{c}=\textrm{max}(T(s))\,\Big{|}\,\frac{L_{R}(s)}{L_{T}(s)}>\delta=\frac{1}{2},$$

$$T_{\text{max}}=0.2T_{\text{peak}},$$

$$T_{\text{min}}=T_{\text{chrom}}.$$

$$T_{\text{min}}\le T_c\le T_{\text{max}}.$$

$${\kappa }_{\parallel }(T)$$

$$\Lambda (T)$$

Full text

1. Introduction

One of the main difficulties encountered when using field-aligned hydrodynamic models to study the physics of magnetically closed coronal loops is the need to implement a grid that fully resolves the steep and dynamic transition region (TR) (see e.g. Bradshaw & Cargill, 2013, hereafter BC13). In a static coronal loop, the energy balance in the TR is between a downward heat flux from the corona and optically thin radiative losses to space. Defining the temperature length scale as $L_{T}=T/|dT/ds|$, where $s$ is a coordinate along the magnetic field, such static loops can have $L_{T}$ of order 1 km in the TR. During impulsive heating, the heat flux from the corona can be enhanced by orders of magnitude, so that $L_{T}$ becomes very small, of order 100 m and less in some cases.

Properly resolving these small length scales is essential in order to obtain the correct coronal density in response to heating (BC13), otherwise the downward heat flux ‘jumps’ over an under-resolved TR to the chromosphere, where the incoming energy is then strongly radiated. BC13 showed that major errors in simulating the coronal density evolution were likely with a lack of spatial resolution.

However, since the thermal conduction time scale across a grid cell scales as $(\Delta s)^{2}$ and $T^{-5/2}$ (where $\Delta s$ is the grid cell width), cells in the cooler lower TR must become dramatically smaller to resolve the temperature gradient that supports the incoming heat flux. In addition, numerical stability requires that the minimum time scale across the entire grid is temporally resolved. Thus, locally satisfying $\Delta s<L_{T}$ in the TR implies long computation times for fully resolved field-aligned simulations. The problem is more severe in three-dimensional (3D) magnetohydrodynamic (MHD) models where computational resources place significant constraints on the achievable resolutions. Therefore, other approaches need to be found.

One such approach has been developed by Linker et al. (2001), Lionello et al. (2009) and Mikić et al. (2013) for use in global 3D MHD simulations. This technique modifies the temperature dependence of the parallel thermal conductivity $\kappa_{\parallel}(T)$ and radiative emissivity $\Lambda(T)$ below a fixed temperature value ($T_{c}$) in the TR. It has been demonstrated that such modifications artificially broaden the TR for $T<T_{c}$ without significantly changing the coronal properties of the loop. While this is a particularly attractive property the method still requires sufficiently high resolution to properly resolve the small length scales in the TR above $T_{c}$. Furthermore, the characteristic temperature of the TR can change dramatically during the heating and cooling cycle, making it unclear how suitable a fixed value of $T_{c}$ is for capturing the dynamic evolution of coronal loops that are impulsively heated.

Another approach has been developed by Johnston et al. (2017a, b), who treated the (unresolved) TR as a discontinuity across which energy is conserved through the imposition of a jump condition. No attempt is made to properly resolve the TR. Instead, the method is employed with coarse spatial resolutions that are supported by dynamically locating the top of the unresolved transition region (UTR). This approach gives good agreement with fully resolved HYDRAD simulations (Bradshaw & Mason, 2003; Bradshaw & Cargill, 2006, BC13) whilst significantly speeding up the computation time (e.g. Johnston et al., submitted).

This Letter proposes a new approach for modelling the TR that incorporates the strengths of both approaches. Principally, the ability to dynamically identify the UTR to prescribe an adaptive $T_{c}$ and then broaden the length scales in this unresolved region only in order to eliminate the need for highly resolved numerical grids. The outcome is an extremely powerful method that (1) removes the influence of numerical resolution on the coronal response to heating and (2) accurately predicts the results of properly resolved field-aligned models (e.g. BC13) when employed with the coarser spatial resolutions achieved by multi-dimensional MHD codes.

2. Numerical Model and Experiments

2.1. Numerical Model

We solve the field-aligned hydrodynamic equations using the HYDRAD code (Bradshaw & Mason, 2003; Bradshaw & Cargill, 2006, BC13) run in single fluid mode. HYDRAD uses adaptive regridding to ensure adequate spatial resolution, with the grid being refined such that cell-to-cell changes in the temperature and density are kept between user defined values (taken as 5% and 10% here), where possible. The largest grid cell in all of our calculations has a width of $10^{6}$ m (1,000 km) and each successive refinement splits the cell into two. Thus, a refinement level of RL leads to cell widths decreased by $1/2^{\textrm{RL}}$. The maximum value of RL is limited to 14, corresponding to a grid cell width of 61 m in the most highly resolved parts of the TR.

2.2. Transition Region Adaptive Conduction

While it is possible to employ such high resolution grids to properly resolve the TR in field-aligned codes, this approach is unlikely to be a viable way of running multi-dimensional MHD simulations. Therefore, we have developed a new method that addresses this difficulty of obtaining adequate spatial resolution in numerical simulations by modelling the transition region using adaptive conduction (TRAC) coefficients that act to broaden any unresolved parts of the atmosphere. This treatment of thermal conduction in the TR, referred to as the TRAC method, relies on a dynamic capability to track resolved and modify unresolved conductive fluxes.

2.2.1 Identification of an adaptive cutoff temperature

The implementation of TRAC is comprised of two main parts. The first part is to identify the maximum temperature of any unresolved grid cells in the numerical domain and use the calculated value to prescribe an adaptive cutoff temperature ($T_{c}$). This is done by using an algorithm that is based on the method employed by Johnston et al. (2017a, b) for locating the top of the UTR in the jump condition method as follows.

We define the temperature length scale as,

[TABLE]

and the resolution in a simulation is given by the local grid cell width,

[TABLE]

where $s$ is the spatial coordinate along the magnetic field.

Using these definitions, the cutoff temperature is defined as the maximum temperature that violates the resolution criteria of Johnston et al. (2017a, b),

[TABLE]

which corresponds to not having multiple grid cells across the temperature length scale (i.e. unresolved temperature gradients).

An upper bound for the cutoff temperature is set as 20% of the peak coronal temperature in the loop at the time when $T_{c}$ is evaluated (though the results are only weakly dependent on the maximum temperature fraction),

[TABLE]

and a lower bound set as the temperature value of the isothermal chromosphere,

[TABLE]

In this Letter the lower bound is taken as $T_{\textrm{chrom}}=2\times 10^{4}$ K. Employing these definitions, we dynamically adjust $T_{c}$ with the criteria that it should satisfy,

[TABLE]

Hence, the cutoff temperature is adaptive in identifying the UTR, with the value of $T_{c}$ used in the method changing in response to coronal heating and cooling.

2.2.2 Broadening the unresolved transition region

The second part of TRAC is to broaden the steep temperature and density gradients in the UTR. This is achieved using the approach developed by Linker et al. (2001), Lionello et al. (2009) and Mikić et al. (2013)

Below the adaptive cutoff temperature ($T_{c}$), the parallel thermal conductivity ($\kappa_{\parallel}$) is set to a constant value,

[TABLE]

and the radiative loss rate ($\Lambda$) is modified to preserve $\Lambda(T)\kappa_{\parallel}(T)$,

[TABLE]

Increasing the parallel thermal conductivity and decreasing the radiative loss rate, at temperatures below $T_{c}$, has the desired effect of broadening the temperature length scales in the UTR. This helps TRAC prevent the heat flux jumping across the unresolved region while maintaining accuracy in the properly resolved parts of the atmosphere. (The broadening effect will be explained in greater detail in forthcoming work.) Furthermore, we note that the formulation of TRAC (1) makes no assumptions about the spatial resolution in a simulation and (2) the parallel thermal conductivity reduces to the classical Spitzer-Harm value when the TR is properly resolved.

2.3. Validation Experiments

The effectiveness of the TRAC method to obtain the correct interaction between the corona and TR, in response to rapid heating events, is investigated by considering two impulsive coronal heating events, comprising short and long pulses that last for a total duration of 60 s and 600 s, respectively. The temporal profile of the heating pulses is triangular with a peak value of $Q_{\textrm{H}}=2\times 10^{-2}$ Jm*-3s-1*, while the spatial profile is uniform along the loop. We release the energy in a coronal loop of total length 100 Mm, which includes a 10 Mm chromosphere attached at the base of each TR. Thus, the total energy injected into the coronal part of the loop is $4.8\times 10^{7}$ Jm*-2* ($4.8\times 10^{8}$ Jm*-2*) in the 60 s (600 s) heating pulse simulations. These heating conditions are representative of reasonably powerful flares and present a challenge to resolving the TR.

For each simulation, the main assessment of the performance of TRAC is a comparison with the results from two alternative methods that are commonly used to treat thermal conduction. Each method is applied with spatial resolutions that cover several orders of magnitude, ranging from those required for 3D MHD codes to run in a realistic time to those that fully resolve the TR but are achievable only in field-aligned models. The first of these methods is the classical Spitzer-Harm heat flux formulation, while the second is the approach developed by Lionello et al. (2009), which artificially broadens the TR below a fixed specified temperature. For brevity and clarity we define these two conduction methods as SH and L09, respectively. We note that the L09 method uses the broadening technique outlined in Section 2.2.2 but the modifications are applied below a fixed cutoff temperature taken as $T_{c}=250,000$ K (which is a typical value used by Lionello et al. (2009) and Mikić et al. (2013)).

Using the same set of parameter values but employing the three different conduction methods (SH, L09 and TRAC), we repeated each simulation for RL = [0, 1, 3, 5, 7, 9, 11, 13, 14] refinement levels to create a group of simulations run for each conduction method.

All of the simulations start from the same initial conditions. These are calculated using the SH parallel thermal conductivity and unmodified radiative loss rate. At $t=0$ (the initial conditions) a small spatially uniform background heating term ($Q_{\textrm{bg}}$) is present. This gives a starting temperature of order 1 MK. However, we note that $Q_{\textrm{bg}}$ is switched off thereafter and the total energy in the initial conditions is negligible compared with the energy released into the loop during the heating pulses.

3. Results

3.1. Coronal Response to Heating

Figure 1 shows the temporal evolution of the coronal averaged temperature ($T$) and density ($n$) and the corresponding $T$ versus $n$ phase space plot (shown on a log-log scale) for the 60 s heating pulse simulation. The coronal averages are calculated by spatially averaging over the uppermost 50% of the loop. Each conduction method is identified with a specific colour and their results are shown in separate rows. The red curves correspond to the SH method (first row), the purple curves the L09 method (second row) and the blue curves represent the TRAC method (third row).

In the panels of each method, each curve corresponds to a simulation run with a different value of maximum refinement level. Taking the SH simulations as an example, the line styles associated with the different refinement levels are shown in the figure legend on the temperature plot, and indicate that RL increases as one moves upwards from the lowest curve (RL=0) to the highest curve (RL=14) in the density plot. These simulations are identical in all respects except for the value of RL. Note that the results for RL=11 and 13 are not shown in the SH panels. In addition, RL=11, 13 and 14 are not shown in the L09 and TRAC panels, where instead the SH simulation run with RL=14 (solid red curves) is shown. This SH simulation (RL=14) is properly resolved and used as a benchmark solution.

The coronal response of the benchmark simulation (SH with RL=14), to the 60 s heating pulse, follows the standard evolution of an impulsively heated loop as described in the literature (Bradshaw & Cargill, 2006; Cargill et al., 2012a, b, 2015; Bradshaw & Klimchuk, 2015; Reale, 2016). The rise in temperature is followed by the density increasing due to the ablation process (Antiochos & Sturrock, 1978; Klimchuk et al., 2008), then, after the time of the density peak, the loop cools by radiation and drains by a downward enthalpy flux to the TR (Bradshaw & Cargill, 2010a, b). The temperature cools below the starting value of 1 MK at around 2100s, and then continues to decrease towards the chromospheric value of $2\times 10^{4}$ K while the density is evacuated from the loop.

Consistent with BC13, the coronal density evolution in the group of SH simulations is strongly dependent on the spatial resolution, requiring grid cell widths of at least 488 m (RL$\geq 11$) for convergence. These resolution requirements are slightly weaker in the L09 simulations, with convergence seen for minimum grid cell widths of 7.81 km or finer (RL$\geq 7$). However, the L09 method achieves this relatively modest relaxation of the resolution dependence at the expense of accurately modelling the loop evolution below the specified cutoff temperature ($T_{c}=250,000$ K).

It is therefore striking that the TRAC solutions are only weakly dependent on the spatial resolution, while maintaining high levels of accuracy throughout the full coronal response range. Grid cell widths of 125 km (RL=3, medium dashed curves) are sufficient to observe convergence to the properly resolved SH solution which employed TR grid cell widths of 61 m. The TRAC solution computed with RL=3 correctly captures the interaction between the corona and chromosphere during the ablation phase, up to the time of peak density where the error is less than 3% and the subsequent decay phase. Even the density oscillations (that are characteristic of the short heating pulse imposed e.g Reale (2016)) and the global cooling of the loop down to $2\times 10^{4}$ K, are both correctly captured. This demonstrates that adaptively broadening unresolved parts of the TR removes the influence of numerical resolution on the coronal density response to impulsive heating.

Figure 2 shows the results for the 600 s heating pulse simulations. Releasing the energy for an extended period of time but with the same peak heating rate has two main consequences. Firstly, the density oscillations associated with the sloshing of the coronal plasma disappear. Secondly, the differences between the density curves that represent the various SH and L09 simulations become more severe.

Considering the refinement level of RL=3 as an example, these conduction methods have peak density errors of 70% (SH) and 20% (L09) when benchmarked against the properly resolved simulation (SH with RL=14), increased from discrepancies of 55% (SH) and 15% (L09) for the shorter heating pulse. Consequently, the convergence requirement on the spatial resolution becomes stricter for the L09 method and remains as severe for the SH method. Grid cell widths of 488 m (RL$\geq 11$) and 1.95 km (RL$\geq 9$) were necessary for the solutions to converge in the SH and L09 simulations, respectively.

On the other hand, the TRAC solutions once again converge at RL=3 (125 km grid cells) with a peak density error of less than 2%, indicating that this treatment of the TR gives rise to a coronal response that can accurately follow the complete coronal heating and cooling cycle with coarse spatial resolutions, irrespective of the heating duration and resulting density regime.

In particular, TRAC eliminates the need for the very short time steps that are imposed by a highly resolved numerical grid. This enables the method to accurately capture the corona/TR enthalpy exchange significantly faster than fully resolved models. A typical speed-up is between two and three orders of magnitude. For example, the computation time of the properly resolved SH simulation (RL=14) is 3.5 days for the 600 s heating pulse, while the run time for the TRAC simulation computed with RL=3 is only 6 minutes.

3.2. Global Evolution

The main question that needs to be addressed in understanding these results is why are the TRAC simulations so successful in describing the coronal response to heating with such large grid cell widths? A comprehensive answer requires a detailed assessment of how the TRAC modifications to the parallel thermal conductivity and radiative loss rate affect the local energy balance and subsequent dynamics in the TR, coupled together with the resulting global evolution of the loop. Here we concentrate just on the latter. The former will be presented in future work.

Figure 3 shows the time evolution of the global temperature, density, pressure, and velocity profiles for the 600 s heating pulse simulations. Each solid curve represents a different snapshot from the properly resolved SH simulation (RL=14) and the dashed blue curves imposed on top are the corresponding snapshots from the TRAC simulation computed with 125 km grid cells (RL=3). (The comparison with an under-resolved SH solution (e.g. SH run with RL=3) has been previously discussed in BC13 and Johnston et al. (2017a) and will not be considered further.)

First we focus on the evolution of the pressure. Of particular importance is the formation of the pressure gradients in the TR. Adaptively broadening the temperature and density profiles in the unresolved parts of the TR prevents the downward heat flux jumping across the unresolved region. This ensures that the incoming energy goes into increasing the gas pressure locally, rather than being lost due to artificially high radiative losses (e.g. BC13; Johnston et al., 2017a, b). The scale of the broadening is relatively small but Figure 3 shows that the resulting TRAC approximations of the SH pressure gradients are remarkably good. These pressure gradients are then responsible for driving the flows.

To that end it is clear the TRAC solution correctly captures the global velocity evolution of the properly resolved SH simulation, during both the ablation and decay phases in response to the 600 s heating pulses. (The 60 s heating pulse simulations show the same fundamental properties). This ensures that TRAC captures the mass and energy exchange that takes place between the chromosphere, TR and corona correctly, enabling the method to maintain high levels of accuracy in the coronal temperature and density evolution throughout the full heating and upflow followed by cooling and downflow cycle.

4. Discussion

The difficulty of obtaining adequate spatial resolution in numerical models of the outer solar atmosphere has been a long-standing problem. BC13 demonstrated that lack of adequate spatial resolution during impulsive heating events led to coronal densities that are erroneously small. Johnston et al. (submitted) then went on to show that one consequence of these artificially low coronal densities is to (artificially) suppress thermal non-equilibrium (TNE) in coronal loops. Thus, under-resolving the TR in numerical simulations has very significant implications for (1) the resultant loop dynamics and (2) any comparisons between model predictions and observations.

Several different approaches have been proposed in order to side-step the need for highly resolved numerical grids and the commensurate very short time steps that are required for numerical stability, yet no fully satisfactory solution is available to date. For example, Lionello et al. (2009) artificially broaden the TR below a fixed cutoff temperature ($T_{c}$), while Johnston et al. (2017a, b) dynamically locate the top of an unresolved TR and impose a jump condition that is derived from an integrated form of energy conservation. However, when employed in simulations with coarse spatial resolutions, the latter approach suffers from overestimating the coronal density response to heating while the former, in contrast, underestimates the density.

This Letter has presented a new approach, referred to as the TRAC method, that, for the first time, has successfully removed the influence of under-resolving the TR on the coronal density response to heating while retaining remarkable levels of accuracy compared with fully resolved models. The new method combines the basic ideas from the approaches developed previously by Lionello et al. (2009) and Johnston et al. (2017a, b).

We have considered only impulsive heating events where the spatial distribution of the energy deposition is uniform along the loop. However, the TRAC methodology is designed to deal with steep transition regions whenever they arise, independent of the nature of the heating. A detailed investigation of different forms of heating (e.g. footpoint heating) will be presented in future work.

TRAC accurately captures the coronal response of properly resolved field-aligned models (e.g. BC13) when employed with spatial resolutions that were up to three orders of magnitude coarser. For example, with TRAC, grid cell widths of order 100 km were sufficient to observe convergence to the fully resolved field-aligned model, which required grid cell widths of order 100 m. The peak density errors were less than 3% and the relaxation of the TR resolution culminated in computation times that were three orders of magnitude faster (e.g. 6 minute run times instead of 3.5 days).

The advantages of this new approach are multiple. For field-aligned hydrodynamic simulations of the coronal response to heating (see e.g. Reale, 2014, for a review), the short computation time means that (1) simulations of coronal heating events can be run quickly, permitting extensive surveys of large parameter spaces (e.g. as done by Froment et al. (2018) to study the occurrence of TNE in coronal loops) to be completed significantly faster at a fraction of the computational cost and (2) simulations of multiple loop strands (thousands or more) that comprise an entire active region (e.g. experiments seeking to reconcile heating models with the Hi-C observational data), can be performed with relative ease and high accuracy for the coronal emission. However, full numerical resolution is still required to deduce the details of the emission in the TR at temperatures below the adaptive cutoff temperature.

In 3D MHD codes, the method can be included without the need for high spatial resolution in the TR and a corresponding extended computation time, “freeing up” resources to resolve better the current sheets responsible for the heating. Indeed, our results suggest that high levels of accuracy can be obtained with grid cell widths of order 100 km, which is achievable in current 3D MHD simulations.

This research has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 647214). S.J.B. is grateful to the National Science Foundation for supporting this work through CAREER award AGS-1450230. C.D.J. acknowledges support from the International Space Science Institute (ISSI), Bern, Switzerland to the International Team 401 “Observed Multi-Scale Variability of Coronal Loops as a Probe of Coronal Heating”. We also thank Dr Z. Mikić for providing important clarifications on the implementation of the L09 method.