Transition region adaptive conduction (TRAC) in multidimensional magnetohydrodynamic simulations

In solar physics, a severe numerical challenge for modern simulations is properly representing a transition region between the million-degree hot corona and a much cooler plasma of about 10000 K (e.g., the upper chromosphere or a prominence). In previous 1D hydrodynamic simulations, the transition region adaptive conduction (TRAC) method has been proven to capture aspects better that are related to mass evaporation and energy exchange. We aim to extend this method to fully multidimensional magnetohydrodynamic (MHD) settings, as required for any realistic application in the solar atmosphere. Because modern MHD simulation tools efficiently exploit parallel supercomputers and can handle automated grid refinement, we design strategies for any-dimensional block grid-adaptive MHD simulations. We propose two different strategies and demonstrate their working with our open-source MPI-AMRVAC code. We benchmark both strategies on 2D prominence formation based on the evaporation--condensation scenario, where chromospheric plasma is evaporated through the transition region and then is collected and ultimately condenses in the corona. A field-line-based TRACL method and a block-based TRACB method are introduced and compared in block grid-adaptive 2D MHD simulations. Both methods yield similar results and are shown to satisfactorily correct the underestimated chromospheric evaporation, which comes from a poor spatial resolution in the transition region. Because fully resolving the transition region in multidimensional MHD settings is virtually impossible, TRACB or TRACL methods will be needed in any 2D or 3D simulations involving transition region physics.


Introduction
Numerical tools for solving magnetohydrodynamic (MHD) equations play an important role in solar and astrophysics today. In the past decades, MHD simulations achieved remarkable progress through developments in high-performance computing and in the design of efficient and accurate shock-capturing solvers. All types of phenomena such as sunspots (Rempel et al. 2009;Rempel & Schlichenmaier 2011), solar prominences (Xia et al. 2012;Xia & Keppens 2016), coronal rain (Fang et al. 2013(Fang et al. , 2015Xia et al. 2017), and solar flares (Cheung et al. 2019;Ruan et al. 2020) can be simulated successfully.
In these examples, multidimensional MHD simulations account for anisotropic thermal conduction and radiative cooling in the upper tenuous solar atmosphere, and incorporating these effects requires specific algorithmic approaches. The numerical strategy chosen to handle these effects also translates into constraints that limit affordable resolutions. At the same time, incorporating these effects is expected to be more realistic or reliable than assuming simpler zero-beta, isothermal, or purely ideal MHD conditions.
To handle more realistic dynamics in the solar atmosphere from the upper chromosphere to the corona, the governing equations of ideal MHD must at least be extended by additional source and sink terms for field-aligned thermal conduction and for (optically thin) radiative losses. In this approach, they appear only in the (total or internal) energy balance. It is diffi-cult to include these terms in simulations because these source terms, especially the thermal conduction term, alter the purely hyperbolic character of the overall set of partial differential equations. The ideal MHD set only needs to account for the fastest wave speed, that is, the fast magnetosonic wave, in setting a time-step limit directly based on the chosen grid spacing. The thermal conduction introduces a parabolic diffusion-like source term, while optically thin radiative losses in essence introduce the light speed (or an infinite speed) in the system. They have to be solved in an alternative way in order to prevent a too stringent limit on the allowed time step. Many efforts have been introduced to accelerate and stabilize their calculation while ensuring a correct result. Because it is an explicit approach, the super-time-stepping (STS) scheme is favored in handling thermal conduction, especially for large-scale simulations of clusters. This STS method allows the thermal conduction term to advance forward by several super time-steps in a single MHD time step. The restrictive Courant-Friedrichs-Lewy (CFL) condition can be relaxed in these super time-steps. There are two variants of the STS method, the Runge-Kutta-Chebyshev STS (RKC-STS) method (van Der Houwen & Sommeijer 1980;Verwer et al. 1990;Alexiades et al. 1996), and the recently proposed Runge-Kutta-Legendre STS (RKL-STS) method (Meyer et al. 2012(Meyer et al. , 2014. The latter is considered to be more attractive and stable for simulating anisotropic terms on clusters (Vaidya et al. 2017). Currently, STS methods have been implemented into many MHD codes, such as Athena++ (Stone et al. 2020 Bifrost (Nóbrega-Siverio et al. 2020), and our open-source 1 MPI-AMRVAC code (Xia et al. 2018). For the optically thin radiative loss treatment, the exact integration scheme (Townsend 2009) proved to be both efficient and robust compared with traditional fully explicit or (semi-)implicit schemes .
However, even with the modern advanced algorithmic methods for dealing with radiative losses and conduction, we still have to be very careful when these source terms need to be treated. This is because they are still difficult to include, especially when a steep temperature gradient is involved. In solar physics, the transition region (TR), a thin layer between the cool chromosphere and the hot corona, might be the most concerned. Similarly, whenever a coronal condensation forms by means of thermal instability (e.g., Fang et al. 2015;Xia & Keppens 2016;Claes et al. 2020), the dynamically created boundary between the prominence or coronal rain blob and the corona is yet another example of a TR.
An approximate way to estimate the thermal conduction length scale for a static coronal loop can follow Eq. (1) in Johnston et al. (2017a). With typical solar transition region parameters, the resulting length scale would be ∼ 30 km or even shorter. As pointed out by Bradshaw & Cargill (2013), when the grid cell size is larger than this length scale, for a typical energetic event in the solar atmosphere, the downward enthalpy flux cannot be handled correctly. Part of the energy might then skip the TR and be transported directly down to the chromosphere where plasma conditions are much cooler and denser. The consequence is that part of this energy reaches too low in the atmosphere and is radiated away by the strong optically thin radiation in the chromosphere. This in turn implies that the evaporation from the upper chromosphere or TR cannot be triggered effectively, which will finally affect the density and temperature in the solar corona. Bradshaw & Cargill (2013) performed a series of onedimensional (1D) hydrodynamic simulations with different grid cell sizes to demonstrate this problem. Their result showed that the coronal density of the most refined solution (reaching a grid size of ∼ 100 m) is several times higher than that of the coarsest solution (with grid cell sizes of ∼ 400 km), indicating that in order to fully resolve the thermal conduction in the TR, a length scale of about 1 km or even 100 m is required.
Considering the current achievable computational resources in the context of modern solar physics, a grid cell of about 100 m is only possible for long-term 1D simulations, assuming the length of the loop is typically 10 Mm, even with the help of parallel computing and adaptive mesh refinement (AMR). Depending on the size of the simulation domain, a grid cell of about 1 km is very difficult or almost impossible for long-term simulations that follow dynamics in large coronal volumes for twodimensional (2D) or three-dimensional (3D) simulations. An approximate method is therefore required to compensate or correct for the underestimated evaporation process.
Several approaches have been proposed. Johnston et al. (2017a,b) treated the unresolved TR as a discontinuity. As mentioned, the evaporation is normally underestimated. This means that the velocity at the top of the unresolved TR is incorrect. To correct this velocity, a jump condition is derived from the conservation of energy throughout this discontinuity. This jump condition is solved approximately to obtain a corrected velocity at the top of the unresolved TR. In this way, the underestimated enthalpy flux is fixed. In their tests, the results of this artificial fixing strategy agree well with their benchmark test, which was performed with a highly refined grid. However, the accuracy of this method is limited and requires a certain range of the grid cell size, typically 100-200 km. Moreover, to solve the nonlinear jump condition, a Newton-Raphson or similar iterator is necessary, which might affect the stability and efficiency of this method.
In contrast, another method proposed by Linker et al. (2001) that was later developed and implemented by Lionello et al. (2009) and Mikić et al. (2013), appears to be more robust. This heuristic approach introduces a fixed cutoff temperature (T c ) and then modifies the temperature dependence of thermal conduction and radiative cooling terms where the temperature is lower than T c . The TR would then be artificially broadened, which has the advantage that it can be numerically resolved. This method can also help to compensate for the missing enthalpy flux, either for a static loop (Lionello et al. 2009) or for a dynamically evolving loop (Mikić et al. 2013). In addition, this method could be easy implemented not only for a 1D system, but also for multidimensional systems. However, in some cases, especially for some strong impulsive heating problems, this method is likely to fail. This problem comes from the use of a fixed T c , which cannot always capture the correct position of the TR in such a quickly evolving system.
More recently, Johnston & Bradshaw (2019) combined ideas from the above two methods. This new method is called the transition region adaptive conduction (TRAC) method. In the TRAC method, a cutoff temperature T c is defined as well, but his time, T c is no longer a fixed value. It is calculated every time step, using the criterion defined in their previous works (e.g., Johnston et al. (2017a,b)). Then, the authors adopted the same heuristic modification for thermal conduction and radiative cooling terms as in Lionello et al. (2009);Mikić et al. (2013). This TRAC method is expected to be both efficient and accurate, eliminating the shortcomings of the previous two methods. Johnston et al. (2020) gave a detailed explanation and demonstration of the physical idea behind this method. Based on this, an additional modification to the heating term was proposed. At the same time, a small improvement was made to remove some abnormal jumps in the cutoff temperature between different time steps.
One topic intimately related to the evaporation process that the TRAC approach wishes to improve is the chromospheric evaporation-coronal condensation scenario (Antiochos et al. 2000). The evaporation will trigger thermal nonequilibrium (TNE) in the magnetic loop, which can explain the omnipresent coronal rain or prominence formation (Klimchuk 2019;Froment et al. 2020;Antolin 2020). Succesful condensations also involve thermal instability (TI), a well-known linear instability (Parker 1953;Field 1965) that affects the otherwise marginal entropy mode (Claes & Keppens 2019;Claes et al. 2020). In the evaporation-condensation scenario, a physically correct evaporation is crucial to either TNE cycles and/or to encounter conditions ripe for a condensation into solar prominence or rain due to in situ linear TI. If the TR cannot be resolved correctly, it is possible that numerical simulations cannot produce a correct TNE cycle or prominence formation, as shown by .
The TRAC method has so far been implemented and used only in 1D hydrodynamic simulations. Here, we extend and apply the basic TRAC idea in 2D or 3D MHD simulations that already allow dynamic AMR. We show a 2D simulation of prominence formation by applying two new variants of the TRAC method, which we call TRACL and TRACB. Both of them can be extended directly to 3D simulations. The new 2D numerical model as well as a review of the 1D TRAC method is presented in Sect. 2. The result of applying our new TRACB or TRACL methods, which are directly suitable for multidimensional MHD simulations, is addressed in Sect. 3 in representative 2D AMR-MHD runs. Discussions and conclusions are made in Sect. 4.

Simulation setup
As shown in Figure 1(a), the computational domain is x len = 80 Mm in width and y len = 50 Mm in height along the y-direction. The initial force-free magnetic field is given by where B 0 is 30 G, k = π/x len , θ = π/3 and j = k cos (θ). Black lines in Fig.1(a) show the configuration of this magnetic field. We used the temperature profile obtained by the 1D RADYN code (Hong et al. 2017), which is relaxed from the well-known VAL-C model (Vernazza et al. 1981), as our initial temperature distribution. Then, given the number density 4.6 × 10 14 cm −3 at the bottom y = 0 boundary, the whole atmosphere can be established according to the hydrostatic equilibrium condition. This initial temperature and density profile is shown in Figure  1(b). The equations in our simulation are as follows: All the symbols have their usual meaning: density ρ, velocity v, magnetic field B, total pressure p tot = p + B 2 /2, and temperature T . We assumed that the whole atmosphere is fully ionized and that the ratio between the number density of helium and hydrogen, n He /n H , is 0.1. The strength of the gravitational acceleration g is a constant 274 m/s 2 , and gravity acts downward. On the right-hand side of equation (6), the second term is heat conduction, where we adopted the Spitzer-type purely field-aligned heat conductivity κ = 10 −6 T 5/2 erg cm −1 s −1 K −1 . The third term is optically thin radiative cooling, where the cooling table Λ(T ) is interpolated from the data given in Colgan et al. (2008). The last term H is the assumed artificial heating. Following our previous 2D simulation works on prominence formation (Xia et al. 2012;Zhou et al. 2020), this term is composed of two parts: H bgr and H loc . H bgr , the background heating, is used to maintain the coronal temperature, where H 0 = 10 −4 erg cm −3 s −1 and l 0 = 50 Mm (Withbroe & Noyes 1977;Withbroe 1988), while H loc is the localized heating, used to trigger strong evaporation to induce the condensations, We adopted H 1 = 2 × 10 −2 erg cm −3 s −1 , which is two orders higher than the background heating (Withbroe & Noyes 1977;Aschwanden 2001). For all other parameters, we took We first relaxed the system, during which time we only had H bgr active, for t relax = 71.5 min. To this thermodynamically adjusted equilibrium configuration, we then gradually added H loc within t ramp = 500 s. The simulations were made with our block-based numerical code MPI-AMRVAC Porth et al. 2014;Xia et al. 2018). For the base level, we used 128 grid cells in the x direction and 80 grid cells in the y direction. Allowing a fourlevel AMR leads to an effective resolution of 78 km. A three-step Runge-Kutta time discretization (Shu & Osher 1988) and HLL scheme (Harten et al. 1983) with the second-order van Leer limiter (van Leer 1974) was used to solve equations (5)-(8). The STS method was used to solve the thermal conduction term, discretized as explained in Xia et al. (2018), and the exact integration scheme was applied to solve the radiative cooling term, as mentioned in Section 1. The background field splitting method in combination with an eight-wave method (Powell et al. 1999) were used to treat the magnetic field and control any numerical magnetic monopole errors.
For the boundary conditions, all the physics quantities were fixed at the bottom boundary while the top boundary was open. Reflective boundary conditions were imposed on the two-side boundaries.

1D TRAC method and test
Before we applied the TRAC method to our 2D simulation, we briefly reviewed the original 1D version. We tested this 1D version in our code to have an overall impression of the difference in our setup with and without the TRAC method.
The 1D TRAC method was used to treat 1D hydrodynamic equations,  (14) where s is the 1D coordinate and g is the field-aligned gravity. These are the hydrodynamic variants of our full MHD system, where E = ρv 2 /2 + p/(γ − 1), in contrast to the total energy density e = E + B 2 /2 from Eq. (6). The ratio of specific heats is γ = 5/3. According to Johnston & Bradshaw (2019) and Johnston et al. (2020), we defined the temperature length scale L T as When we set L R to be the size of the local grid cell (in our gridadaptive runs, this is the smallest cell size in use in the grid hierarchy), for example, we can define a cutoff temperature T c as Here, δ is a free parameter that should be smaller than 1/2. Because we used four grid points, that is, a fourth-order center difference method, to calculate the gradient of the temperature, our δ was chosen as 1/4. To restrict this cutoff temperature to a reasonable range, an upper bound, 0.2 max (T (s)), and a lower bound, 2×10 4 K, were set, as recommended by Johnston & Bradshaw (2019). Other restrictions to prevent some sudden jumps can be found in the Appendix A of Johnston et al. (2020). For the region in which the temperature is lower than T c , the TRAC method applies the following modifications: The left-hand side symbols with prime are used to substitute the original quantities in equation (14). In order to test the 1D TRAC method, we picked out a field line at t = 0 in our 2D setup, which is then a coronal loop, starting from point x = −18 Mm at the bottom boundary. This field line is shown by the dashed line in Fig. 1(a). The total length of this loop is 47 Mm. We now discretized it with 160 grid cells at the base level, using a three-level AMR, leading to an effective resolution of 74 Mm, which is similar to the size of our 2D grid cells (78 Mm). All the other settings were adopted directly from the 2D setting described above. Assuming that the shape of this magnetic loop remains invariant, we now ran the equivalent 1D hydro-simulation for over 9 h. The result is shown in Fig. 2.
The upper panels show the result of this 1D hydro-run without the TRAC method, and the lower three panels show the results with the TRAC method. Panels (a) and (d) show the temperature evolutions for both runs. When no TRAC method is used, only a hot coronal loop with gradually increasing temperature is visible; the increase occurs very slowly. However, when the TRAC method is applied, thermal instability occurs after we add the localized heating for about half an hour, so that condensation sets in: a prominence therefore forms in this case. Panels (b) and (e) show the number density distribution along the loop at t = 179 min. The density of the prominence is nearly two orders higher than the coronal density.
In our previous works (Xia et al. 2011;Zhou et al. 2014), this type of 1D simulations was performed in a dipped flux tube configuration, so that the formed condensation was gravitationally trapped and remained very stable for days or even weeks. In the 2D interpretation of such a rigid magnetic field 1D run, the dip corresponds to an upward Lorentz force. In the situation of Fig. 2, however, the configuration is a simple arcade, which cannot provide upward support. Because the flux tube is treated as a rigid body in a 1D simulation, gravity from the heavy prominence cannot change the shape of the prescribed flux tube. The prominence therefore eventually slides down to the chromosphere. This occurs at about t = 250 min. Afterward, evaporation Panels (c) and (f) show the evolution of the averaged coronal number density for both runs. Here, the corona is simply defined as the region between s = 13 Mm to s = 34 Mm. Clearly, the TRAC method improves the evaporation and results in a higher coronal density. The averaged density is typically several times the density found without the TRAC method.
Only in a higher coronal density situation can the thermal instability (this hinges on the quadratic density dependence of the optically thin radiative loss term) as well as the condensation occur. Therefore we conclude that for this type of prominence formation simulation, the TRAC method is very important. This is all clearly very sensitive to the employed heating prescription, but when the TRAC method is used, we can meaningfully explore a relatively large parameter space, beyond what is possible without using TRAC, at much more affordable grid resolutions.

Applying the TRAC method in the 2D setup
Because thermal conduction in multidimensional MHD is aligned with the possibly relocating field lines, the key point of applying a TRAC method in a 2D MHD simulation is the requirement of tracing magnetic field lines from time step to time step. We have recently (Ruan et al. 2020) developed a magnetic field-line-tracing module for MPI-AMRVAC to solve the fast electron heating problem in the context of solar flares, where fieldaligned fast electron beams deposit their energy in the denser atmospheric regions after gaining energy in a coronal reconnection site. Based on the experience there, we now propose two different approaches for applying the TRAC idea in multidimensional MHD settings. While details are provided in the appendix, we briefly introduce the main ideasin the following subsections.

TRACL
A straightforward field-line based or TRACL method is directly using the magnetic field-tracing module from our previous work. For this simulation, we have 128 × 80 grid cells with four AMR levels, which means that we effectively have 1024 × 640 grid cells. Because the arcade configuration is not expected to undergo dramatic changes, we chose to trace 256 field lines whose starting points are uniformly distributed on the lower y = 0 boundary. By tracing the field lines at every time step, we can calculate the cutoff temperature T c for each of them. This T c then essentially changes from field line to field line, and it also changes with time. Then, by interpolating T c from the field lines back onto the AMR grid cells, we obtain the distribution of the cutoff temperature for all grid cells. Consequently, we can applied equations (18)-(20) to modify the local evaluations of these thermal terms, similar to the previous 1D situation. This TRAC method is based on directly integrating field lines, and is appropriately named TRACL.
To demonstrate its effect, we ran two cases, one without any TRAC method and another with this TRACL method. The results are shown in Fig. 3. Panels (a) and (b) are snapshots at t = 175 min for the two different cases. Similarly to what happened in the 1D case, when we did not use the TRACL method, only a hot loop develops in panel (a). This shows that without TRACL, there is still no condensation even at the end of our simulation after t = 250 min. In panel (b) the same setting with the TRACL method is shown. The formation of a solar prominence whose typical lower temperature is 2 × 10 4 K is clearly visble, corresponding to the lower boundary of our employed cooling table.
The distributions of the cutoff temperature T c at different times are shown in the lower row of Fig. 3. Panel (c) shows its distribution at t = 0. The lower and shorter loops have a lower cutoff temperature, and the higher and longer loops tend to have a higher cutoff temperature. The shapes of these initial cutoff temperature contours are clearly consistent with the magnetic field lines, which shows that our magnetic field-linetracing method is robust. After relaxation at t = 72 minutes, the pattern adjusts, as shown in Fig. 3, panel (d). The lower cutoff temperature region near the center has grown: although the temperature in this region increases, the resulting cutoff temperature decreases. In the last panel (e), after some time of evaporation, we show the cutoff temperature distribution at t = 175 minutes, when condensation has begun. The cutoff temperature of the field lines, which went through the heated loop that then collapsed into a filament, is slightly higher than in the outside loops. The field lines have (slightly) adjusted during the entire run, and especially show some minor central dipped portion that develops at this end time.

TRACB
From the above numerical experiment, we conclude that the TR-ACL method can successfully correct for underestimated evaporation. Nevertheless, integrating many magnetic field lines at every time step, the core of the TRACL method, is a time-consuming process. In a modern parallelized simulation, global communication is required every time the magnetic field line goes across processors, and this happens several times for each field line (Ruan et al. 2020). Considering that the precise value of T c is not necessarily highly accurate, we here advocate an approximate and fast way to trace the magnetic field lines in setups such as the one studied here.
Our MPI-AMRVAC code is a block-based AMR code, which means that the computational domain is first discretized into blocks, and these blocks all have an equal amount of grid cells. Each block is composed of several grid cells at the same AMR level, and for this 2D simulation, we chose the block size to be 16 × 16 (MPI-AMRVAC allows changing this setting at runtime). Our effective resolution of 1024 × 640 grid cells means that we handle 64×40 blocks at most. We may therefore calculate the local cutoff temperature T c,local using equation (17) for each block. Then, the maximum operator in equation (17) is first done locally in each block, and only this T c,local , as well as the unit vector of the magnetic field at the center of this block, denoted as b local , is then broadcasted to all the processors. By sharing this 64 × 40 global array consisting of magnetic field unit vectors and cutoff temperatures, the field line can be traced independently in each processor using this information. Again, the cutoff temperature along the (coarsely represented) field lines is then interpolated back onto the AMR grid cells, so that we have the distribution of T c throughout the whole computational domain. Because this method is based on blocks, we call it the TRACB method. Additional global communication is required only once for each time step, so tat the computational efficiency is expected to be higher than for the TRACL method, at the cost of a coarser evaluation of the cutoff temperature distribution.
We ran the same simulation as in the previous subsection, but applied the TRACB method instead. This result is shown in Fig. 4: as in Fig. 3, panel (a) is the result without the TRAC method, for comparison. With TRACB, panel (b) shows that the condensation indeed occurs, similar to the TRACL case. The atmosphere is similar to that in Fig. 3(b), and the condensation occurs at the same position, as expected. However, because this method only employs a fairly coarse cutoff temperature distribution, we find that this time, the prominence is not as stable, and especially the left-right symmetry is easily broken. The newly formed prominence therefore has already started to slide down to one side, similar to what we obtained in the 1D simulation, for instance, after t = 250 min in Fig. 2(d).
Again, in the lower row of Fig. 4, we show the distribution of T c at t = 0, 72, and 275 minutes. For the initial snapshot in panel (c), the distribution is similar to Fig. 3, but clearly not as smoothly represented as with the TRACL method. The cutoff distribution is only aware of magnetic field lines as traced by the block-center magnetic field. For this initial state, only using the block center magnetic field can also trace the magnetic field lines and give a reliable result. However, after relaxation, as shown in panel (d), the result differs slightly from the TRACL method. This is mainly because equation (17) relies on a constant δ parameter. Such a fixed value throughout the whole domain, in combination with the hierarchical AMR grid structure, will sometimes result in a hierarchical T c structure, as shown in Fig. 3(d), which might be considered a drawback of the TRACB method. In the TRACB method, the final distribution of T c relies more heavily on interpolation. Then, for the final panel at the time t = 175 minutes, the distribution is again similar with Fig. 3(e), especially when we recognize that the atmosphere itself is already slightly different. In the last two panels (d) and (e), some regions near the upper corner are not appropriately covered by the field lines we selected because the AMR structure uses a coarser grid there. In most applications, the cutoff temperature at the upper corners does not affect the result, we simply filled them with some arbitrary values.

Comparing TRACL and TRACB
To be more precise in the comparison of the two new TRAC methods, we cut a slice along the horizontal line y = 2 Mm and quantified the time evolution of T c . The result is shown in Fig. 5. The left panel shows the TRACL and the right panel the TRACB method. The short loops near the center always have the lowest temperature, which means that along these corresponding field lines, it is not necessary to modify the TRAC from the criterion expressed in equation (17). However, on either side of this central region, the typical T c is about 2 × 10 5 -3 × 10 5 K for a steady loop. After the localized heating is introduced, T c near the heating source increases quickly and maintains a high temperature value of over 4 × 10 5 K. Both methods show a similar pattern in the temporal evolution of the cutoff temperature at this height, although there are noticable differences, mostly confirming that TRACB yields a more coarser representation of the actual evolution.

From 1D hydro-to multidimensional MHD
Although the TR is a very thin layer in the solar atmosphere, it plays a crucial role in many phenomena occurring in the solar atmosphere. Traditionally, the TR is considered to be the region in which temperature is about 10 5 K, between the chromosphere and the corona. Many of our own previous numerical simulations that related to TR evaporation, such as prominence formation (Xia et al. 2012;Xia & Keppens 2016) or coronal rain formation (Fang et al. 2013(Fang et al. , 2015Xia et al. 2017), did not fully resolve the TR, so that the evaporation was underestimated. The combination of modern shockcapturing schemes, employing limited linear or higher order reconstructions in our finite volume code, did not preclude a numerically stable representation of TRs, but this at the expense of an altered energy balance from the true TR regime. Then, with insufficient evaporation, the density in the solar corona would be severely affected, as shown by Bradshaw & Cargill (2013) in their 1D hydro-simulation. This is mainly because the enthalpy flux in the transition region is underestimated in an unresolved TR. For 1D problems, brute force high resolution that is for example efficiently achieved with AMR could be used to solve this problem. However, for multidimensional problems, a similarly high resolution is beyond computational reach, especially because the requirement might be even higher than estimated from 1D simulations. The 1D case shown in Sec. 2 had an effective resolution of 74 km, similar to the 78 km, that is, the effective resolution of the 2D case in Sec. 3, which used four AMR levels. In the 1D case, we find that when we add just one more grid level, that is, resolve to 37 km, the condensation starts to occur in a non-TRAC simulation, at nearly the same time as the 1D TRAC case. However, for the 2D case, even when we increased to seven AMR levels and run the simulation for 3.5 h, we still did not see any condensation when we did not use TRACB or TRACL.
There are two distinct possibilities: (a) the resolution required to correctly represent the TR is more strict in 2D simulations, or (b) the TRACL/TRACB methods in 2D give unphysical results. To judge this aspect and argue in favor of option (a), Fig. 6 shows the time-evolution averaged coronal density for different cases. For the same number of four AMR levels, the run without the TRAC method solid line is much lower than the two runs with our new TRACB and TRACL methods dotted and dashed lines. This is expected from the design of the TRAC method, and is meant to capture the effect of a higher resolution that is needed to obtain the proper evaporation. However, the dash-dotted line represents the average coronal density evolution for a high-resolution run (employing seven AMR levels) without any TRAC method. Before the localized heating is added (before 75 minutes) and when the localized heating starts to work, this reference result is very close to the two dashed lines. However, after about t = 100 minutes, it deviates from the two dashed lines, well before the condensations that occur in the last two cases (at about t = 120 minutes). Based on this result, we argue that multidimensional simulations require an even higher resolution when the intricacies of anisotropic thermal conduction along changing and evolving magnetic field lines are taken into account.

Efficiency considerations
The 1D TRAC method is a fast and accurate method, and in our MPI-AMRVAC implementation, the computational time for a 1D TRAC simulation is almost the same as without the TRAC method. Considering that the modified temperature prescription in the thermal conduction term allows us to use a larger time step in each STS substep, the TRAC method sometimes is even faster. However, applying it in a multidimensional simulation requires additional time spent on tracing the magnetic field lines. Ruan et al. (2020) developed and used a generic magnetic fieldline-tracing module for our open-source code. In this application, we did not (yet) optimize the field-line-tracing module, and hence it was very time-consuming and performed poorly when the number of processors was increased (because the many field lines were processed serially). For the results presented in this paper, we optimized this module by parallelizing it and allowing a reliable result using the TRACL method, obtained within a reasonable time span.
Considering the block-based structure of our code, we implemented and presented an alternative block-based TRACB method, keeping the essential functionality of estimating the T c distribution over the domain. This TRACB method is mostly performed locally in each processor, with only one additional global communication every time step. This is expected to improve the overall efficiency.
We have also experimented with yet another way to gain efficiency by noting that the TRAC method only modifies the region in which the temperature is lower than the cutoff temperature, which is typically the region below the TR. This means that integrating the magnetic field lines throughout the whole computational domain is slightly in vain as long as no coronal condensations have formed. In practice, we can therefore only integrate field lines up to a certain height, for example, 8 Mm or 10 Mm, and mask out the higher coronal part. This will save memory space as well as computational resources, and although we did not show this here, the results are quite similar to the unmasked version that we used throughout her (in both TRACL and TRACB). This may be a crucial ingredient for actual 3D simulations using TRACL or TRACB variants.
InJohnston & Bradshaw (2019) andJohnston et al. (2020, the upper value of T c was set to T ub = 0.2 max(T (s)). In our applications, we also gave this upper boundary to T c when each magnetic field line was traces. When we choose to integrate only part of the field lines, as suggested in the previous paragraph, we cannot obtain the maximum temperature for each field line. However, this upper bound is not so frequently invoked in the current setup. It is therefore not necessary to distinguish T ub from field line to field line, and the overall maximum temperature can instead be used over the full grid as a good substitute. For instance, for flare-related simulations that reach truly high temperatures, setting this global T ub might also be a good choice. Table 1 quantifies the computational time used for the simulation described in Section 2, as run by the different methods mentioned above. When field lines are only integrated in a region below 8 Mm, we refer to our result as TRACL mask and TRACB mask . For a fair comparison, we ran the simulation always from t = 0 to t = 71.5 min, that is, before we added the localized heating, so that the number of blocks on different AMR levels is the same or at least very similar. The simulations were run on  From these timing results, we conclude that for the TRACL method, masking the coronal part may well be necessary for large-scale 2D and ultimately 3D problems. If a full-domain integration is necessary for the problem at hand, the TRACB method could be used, while the most efficient implementation is our TRACB mask approach.

Broadening of the TR
As discussed in our introduction, before the TRAC method, several other methods aimed to solve the unresolved TR problem. The method by Linker et al. (2001) and Lionello et al. (2009), using a globally fixed cutoff temperature T c , could be considered as the original idea behind the TRAC method. An actual unresolved TR is then artificially broadened when the thermal terms are considered, and hence becomes a resolved TR. This is achieved by keeping the product of κΛ in equations (18)-(19) unchanged. However, a globally fixed cutoff temperature cannot always capture the correct position of the TR because, as proposed by Bradshaw & Cargill (2010a,b), from the physical point of view, the interface between the TR and the corona should be the point where the enthalpy flux changes from a sink into a source. However, this criterion is not so easily implemented into a numerical MHD code. Johnston et al. (2017a,b) therefore defined the top of the TR by an artificial correction formula, which was then used in the TRAC method; this is our equation (17) here.
The effect of the TRAC method on the height of the TR is shown in Fig. 7. The left column shows the run without the TRAC method, the middle column shows the TRACL method, and the right column shows the result with the TRACB method. The upper row shows the time evolution of the height of the top of the TR, h tra . Here, the top of the TR is simply defined as the place at which the temperature gradient along the vertical direction (y-direction) is maximum. Directly using equation (17) is not appropriate because for the short loops at the center, it is possible that no points on the loop fulfill this threshold. The time evolution of the horizontally (x-direction) averaged transition region height h tra is plotted in the lower row, as well as the averaged height of the bottom of the TR denoted with h chr . This h chr is defined as the first point from the bottom at which the temperature exceeds 2 × 10 4 K. When no TRAC is applied, the TR does not have any obvious response to the strong localized heating. After we introduced the localized heating, h tra near the source region decreased slightly while the center part increased slightly, but the overall change is minor without TRAC. In contrast, when the TRAC method was used, the TR broadened, especially after we introduced the localized heating. The width of the TR was four to six times broader. The evolution of the TR-ACL method and the TRACB method is basically very similar, with some deviations near the end of the simulation.
In this example and in other simulations that follow the evaporation-condensation scenario, we are interested in the time and location at which condensations occur in the solar corona. With that in mind, having the correct densitytemperature throughout the corona is more important than truly resolving a 'correct' TR. A broadened TR is acceptable in this type of simulations. However, when condensations do occur (directly related to TI onset), another "TR" develops between the prominence or coronal rain blob and the solar corona. This is the prominence-corona transition region (PCTR). When equation (17) is used, this PCTR also broadens artificially. The PCTR broadening is not as strong as the broadening in the TR because a much weaker heating acts at these heights. When we instead require that the PCTR not be modified at all, as mentioned in Sec. 4.2, we can add the mask and use TRACB mask or TRACL mask approaches.

Conclusion
Chromospheric evaporation is an important source of mass and energy for the solar corona. However, the TR between the chromosphere and corona has a very steep temperature gradient, and this implies a real challenge for numerical simulations, especially when we consider the anisotropic nature of thermal conduction. Earlier 1D hydrodynamic settings have shown that using a too low grid resolution, that is, one that cannot resolve the TR gradients fully, might finally underestimate the evaporation, resulting in too low coronal densities. The TRAC method (Johnston & Bradshaw 2019;Johnston et al. 2020) is an ad hoc way to correct for these problems, and this method has been investigated in 1D simulations so far.
In this work, we introduced two novel means to apply the TRAC method in multidimensional simulations that employ block-grid-adaptive techniques to efficiently resolve chromosphere to coronal dynamics. Both involve the use of a dynamic field line tracing from the chromosphere to the corona, along which an essentially 1D approach is then employed. Because modern simulations use parallel computers where the computational domain is split into blocks that are spread over the CPUs, such field line tracing must handle the complexities of field line relocations and evolutions that involve multiple processors. In the most modern applications, this parallelism is combined with AMR, introducing the dynamic creation and destruction of hierarchical grid blocks. Our TRACL and TRACB methods are fully compatible with AMR, but can also easily be used in pure domain-decomposed MHD simulations that use a simple structured grid throughout. Both methods were tested here and were implemented in the open-source MPI-AMRVAC code. They are available and integrated in our open-source software.
Although the example presented here is 2D, it is trivial to extend the method from 2D to 3D, and it is also already implemented in the dimension-independent MPI-AMRVAC code. In 3D, the TRACB method shows a higher efficiency than the TR-ACL method. The magnetic configuration in this example 2D scenario is a relatively simple arcade structure. No reconnection occurs during the times that are simulated. However, the magnetic field-line-tracing method is robust and can be applied to more complex configurations such as reconnection events, as shown in Ruan et al. (2020), who followed a flare event. One notable future improvement might be the choice of the initial points used to integrate field lines. Currently, the initial points are dis-tributed near the bottom boundary, which is applicable in most cases. This might not be the best choice for some configurations, especially near side boundaries in a complex realistic field.