Numerical Study of Divergence Cleaning and Coronal Heating/Acceleration Methods in the 3D COIN-TVD MHD Model

In the solar coronal numerical simulation, the coronal heating/acceleration and the magnetic divergence cleaning techniques are very important. The coronal–interplanetary total variation diminishing (COIN-TVD) magnetohydrodynamic (MHD) model is developed in recent years that can effectively realize the coronal–interplanetary three-dimensional (3D) solar wind simulation. In this study, we focus on the 3D coronal solar wind simulation by using the COIN-TVD MHD model. In order to simulate the heating and acceleration of solar wind in the coronal region, the volume heating term in the model is improved efficiently. Then, the influence of the different methods to reduce the ∇ ⋅ B constraint error on the coronal solar wind structure is discussed. Here, we choose Carrington Rotation (CR) 2199 as a study case and try to make a comparison of the simulation results among the different magnetic divergence cleaning methods, including the diffusive method, the Powell method, and the composite diffusive/Powell method, by using the 3D COIN-TVD MHD model. Our simulation results show that with the different magnetic divergence cleaning methods, the ∇ ⋅ B error can be reduced in different levels during the solar wind simulation. Among the three divergence cleaning methods we used, the composite diffusive/Powell method can maintain the divergence cleaning constraint better to a certain extent, and the relative magnetic field divergence error can be controlled in the order of 10−9. Although these numerical simulations are performed for the background solar corona, these methods are also suitable for the simulation of CME initiation and propagation.


INTRODUCTION
The 3D COIN-TVD MHD model which was proposed in [1][2][3] and was improved in [4][5][6][7][8] in recent years can effectively realize the coronal-interplanetary 3D solar wind simulation. This model uses the TVD Lax-Friedrichs (TVD-LF) scheme uniformly in the corona region and the interplanetary space region, and a combination of Open Multi-Processing (OpenMP) based on shared memory and Message Passing Interface (MPI) based on distributed memory has been successfully used to study the solar wind background from the corona to the interplanetary space.
The solar energy is stored in the solar nucleus, and the generated radiant energy spreads from the inside to the outside. The solar temperature should theoretically decrease with the increase of the heliocentric distance. However, the temperature of the upper atmosphere corona is much higher than that of the lower atmosphere (photosphere). The reason for the abnormal warming of the atmosphere has not yet been investigated. Therefore, coronal heating/acceleration is a central issue in the solar coronal simulation and has been discussed by many researchers (e.g., [9][10][11][12][13][14][15][16]). Parker proposed a basic theory for the problem of heating an expanding solar corona [17][18][19]. Later, various methods for solar wind acceleration and coronal heating have been developed. For example, the Alfvén wave heating method (AHM) can accelerate solar wind through the exchange of momentum and energy between large-scale Alfvén wave turbulence and solar wind plasma [10]. The turbulent heating method (THM) assumes that the turbulent free energy is transformed into the energy accelerated by the solar wind when the turbulent free energy changes with the heliocentric distance [10]; By adding momentum and energy source terms to the MHD equations [16], the volume heating method (VHM) has been widely used in solar wind simulation (e.g., [15,20,21]. In the MHD simulation, the divergence of the magnetic field should be strictly controlled to zero. The nonzero divergence of the magnetic field can lead to the ∇ · B error during the calculation. When this occurs, numerical instability may develop and the simulation can break down. Therefore, scientists have proposed many methods to control the divergence of the magnetic field, such as the generalized Lagrange multiplier (GLM) method [22][23][24], the CT method [21,[25][26][27], the projection method [28], the vector potential method [29,30], the Powell method [31,32], the diffusion method [7] and the globally solenoidality-preserving (GSP) method [33].
In this study, we adopt the COIN-TVD model to simulate the coronal solar wind. Similar to [20,21], we use the volume heating sources to model the solar wind heating/acceleration process in the simulation.
In Governing Equations of Coronal Interplanetary-Total Variation Diminishing Model, we introduce the equations of the COIN-TVD MHD model. Mesh Grid System and Numerical Scheme describes mesh grid system and boundary conditions. Volume Heating Method and Magnetic Field Divergence Cleaning Methods presents the VHM method and three magnetic field divergence processing methods. Numerical Results shows the results of numerical simulation and comparisons of three methods for processing magnetic field divergence. In Conclusions and Discussions, we make the conclusion and discussion.

GOVERNING EQUATIONS OF CORONAL-INTERPLANETARY TOTAL VARIATION DIMINISHING MODEL
The ideal MHD equations are used to simulate the coronal solar wind. Under the Corotating coordinate system, equations can be written as: where ρ is the mass density, v is the plasma velocity, B is the magnetic field, P is pressure, μ 0 is the magnetic permeability of free space, I is the unit tensor, G is the gravitational constant, M s is the solar mass, f −ω × (ω × r + 2ω × v) is the additional fictitious force densities, in which ω is the angular velocity of the rotation, and c is the polytrophic index, which is set to be 1.05 in this study.

Mesh Grid System
In the spherical coordinate, the range of the calculation area is expressed as 1Rs ≤ r ≤ 22.5Rs, − π 2 ≤ θ ≤ π 2 , and 0 ≤ ∅ ≤ 2π, where r is the radial distance from the solar center to the solar surface, θ is latitude, and ∅ is longitude. To avoid the singularity, the computation domain is divided into six identical component meshes to envelop a spherical surface with partial overlap on their boundaries [34]. The following grid partitions are employed; the grid mesh is built in the form of 224(r) × 180(θ) × 360(∅). The radial direction uses a proportional grid, the radial step length increases from 0.0161R S at the inner boundary of 1R S to 0.3636R S at the outer boundary near 22.5R S , and the total number of grids at r-direction is 224. In the latitudinal and longitudinal directions, the grid resolution is Δθ Δ∅ 1 + .

Numerical Scheme
In the COIN-TVD model, all of the physical quantities are computed from the TVD-LF numerical scheme in a facecentered grid structure (e.g., [7,8]). And this scheme is performed in the six-component mesh grid system.
The Carrington Rotation (CR) 2199 is chosen for background establishment. The initial magnetic field B 0 is given by using the potential field source surface (PFSS) model [35,36], the spherical harmonics coefficients were used to obtain the initial PFSS solution is 6. And other initial parameters, such as plasma density ρ 0 , temperature T 0 , and velocity v, are calculated by Parker's solar wind flow solution [17]. The temperature and the number density on the solar surface are set to be 1.5 × 10 6 K and 1.67 × 10 8 cm −3 , respectively. The boundary condition of the magnetic field at the inner surface also remains fixed all through the simulation.
The parameters at the outer boundary are set according to the projected characteristic boundary conditions e.g., [32,37,38].

VOLUME HEATING METHOD AND MAGNETIC FIELD DIVERGENCE CLEANING METHODS
In this section, we introduce the numerical schemes of the volume heating method and three methods to constrain ∇ · B in MHD simulation.

Volume Heating Method
Due to the limitations of observation and theory, there is no mature theoretical model to describe the mechanism of coronal heating and solar wind acceleration. Here, we use the volume heating method to solve the issue of coronal heating and solar wind acceleration. We add the source terms of momentum S M and energy Q E to the MHD Eq. 1, Eq. 2, Eq. 3, and Eq. 4 as follows: According to the work in [7,[39][40][41], we set energy and momentum source terms as follows: Here, γ 1.05, which is the polytrophic index. In the calculation region, the polytropic index γ need not be set very large. γ 1.05 can heat the corona and accelerate solar wind.
Here, r is the heliocentric distance, Q 1 and L Q1 are the intensity and attenuation length of heating, and S 1 and L M are the intensity and decay length of the momentum addition. The parameters L Q 1 and L M are set to be 1, Q 1 Q 0 C a , and S 1 S 0 C a . To test the influence of the parameter of energy and momentum source terms, we set two groups of different parameters for comparison. In model A, we set: and find that the coronal heating and solar wind acceleration were not obvious. In model B, we adjust the parameters based on [33], which are . Here, f s R s R ss 2 BR s B Rss is the expansion factor, where Rs is the solar radius, R SS 2.5R S , and B R S and B R SS are magnetic field strength at the solar surface and at R SS , respectively. Inspired by the Wang-Sheeley-Arge (WSA) model [42,43], the solar wind speed is related to the magnetic field expansion factor f s and the minimum angular distance θ b . As f s increases, the speed decreases, the high-speed stream originating from the center of the open field region always has large θ b , and the lowspeed stream from the coronal hole boundary has a relatively small θ b .
Following [20,44,45], the source term Q E also contains a heat conduction term, the expression of the heat conduction term is · B, ξ is the collisional thermal conductivity parallel to the magnetic field as given in [46] and the proton and electron temperatures are equal to T. If we add the heat conduction term in the Q E , the partial differential in the formula decreases the calculation accuracy. And after the research in [45], many works (e.g., Feng, 2012, [21]; 2017 [33]) verify that without adding heat conduction item, the coronal solar wind can also be accelerated and heated.

Powell Method
The Powell method to maintain the magnetic divergence cleaning constraint is given as follows.
Two divergence source terms, −(∇ · B)B and −(∇ · B)v, are added separately on the right side of Eq. (6) and Eq. (7) to get the following MHD equations: In this way, the divergence of the magnetic field can be propagated to the boundary to reduce the numerical error of ∇ · B in the computational region [31]. From Eq. 10 with the source term, the quantity ∇ · B/ρ satisfies the advection equation, which is, This means that the ∇ · B must be transported by the plasma motions when Powell correction is applied, since the initial and boundary conditions satisfy ∇ · B 0, and the ∇ · B will be near zero for all later times throughout the simulation.

Diffusive Method
The diffusive method is proposed to reduce the error of the magnetic divergence, in which an artificial diffusivity is added at each time step as zB zt η∇(∇ · B). Under the condition of Δt ≤ μ (Δx) 2 η , where μ ∈ (0, 2), the error of the magnetic divergence is diffused away at the maximal rate allowed by iterating: Here, (Δx) 2 Frontiers in Physics | www.frontiersin.org July 2021 | Volume 9 | Article 705744 For satisfying the condition max ⎛ ⎝ Bnds |Bn|ds ⎞ ⎠ ≤ 10 −2 , less than ten iterations are needed at each time step. This method does not violate shock capturing property, at least to the second-order accuracy in smooth regions [8,47].

Composite Diffusive/Powell Method
We combined the Powell method and the diffusive method together in the MHD calculation in the composite diffusive/ Powell method for the first time, and this method can further control the error of the magnetic field divergence. The composite diffusive/Powell method adds two divergence source terms, −(∇ · B)B and −(∇ · B)v, to Eq. 6 and Eq. 7 to get Eq. 10. The divergence of the magnetic field can be propagated to the boundary, and the quantity ∇ · B/ρ satisfies the advection equation (Eq. 11). When solving the equation, the error of magnetic divergence is diffused away at the maximal rate allowed by iterating Eq. 12.

NUMERICAL RESULTS
In this section, we show the numerical results of the solar coronal simulation from 1R S to 22.5R S for CR2199, which are obtained by executing the methods introduced in Volume Heating Method and Magnetic Field Divergence Cleaning Methods.
It takes about 100 h in physical time to obtain the steady state in our simulation. Figures 1,2 present the distribution of the magnetic field lines, the radial velocity, the number density and the temperature on the meridional plane at Φ 180°-0°from model A and model B, respectively. From these figures, it can be seen that the high latitude areas always have fast speed, high temperature, and low density. On the contrary, the radial speed is slower, the temperature is lower, and the number density is higher at lower latitudes around the heliospheric current sheet (HCS), and this is the characteristic feature of the solar wind in the corona [47]. Model B is successful in simulating the acceleration and heating of the solar wind in the corona, as shown in Figure 2. Compared with Figure 1, we can find that both the radial speed and temperature in Figure 2 are higher than those in Figure 1 obviously. This result indicates that the VHM can accelerate and heat the coronal solar wind, and the parameters S 0 , Q 0 , and C a in VHM can affect the coronal heating and solar wind acceleration process significantly.
Then, we present the simulation results of the coronal solar wind with three magnetic divergence cleaning methods. Figures  3-5, respectively, show the variation in the radial speed, the number density, and the temperature along heliocentric distance from 1 to 22.5 Rs with different latitudes of θ −80°and θ −10°a t the same longitude of Φ 0°, where θ −80°locates at the open field region and θ −10°locates at the HCS region. Comparing the three figures, we can find that the radial speed in the open field region is larger than that in the HCS region, the temperature is higher in the open field, and the number density is smaller in the high latitude region.
The composite diffusive/Powell method which combines the diffusive method and the Powell method is our new try to handle the ∇ · B constraint. From Figures 3-5, we can also see that the curve from the composite diffusive/Powell method is always in the middle, so it can generate a stable solar wind structure like the other two methods.
To quantitatively see how ∇ · B evolves, we define the relative divergence error [48] as follows: Here, Δh is the characteristic length of the mesh element.
To investigate how the three magnetic divergence cleaning methods control the ∇ · B error quantitaitvely, we make a numerical comparison for the Error(B) among the three methods. Figures 6,7 show the distributions of the Error(B) on the different meridional planes of Φ 180°-0°and Φ 270°-90°, respectively, for the steady-state solar wind. The three panels in Figures 6, 7 present the results from the composite diffusive/ Powell method, the diffusive method and the Powell method, from left to right, respectively. It is obvious that the Error(B) deduced from the composite diffusive/Powell method is lower than that from the other two methods, on both meridional planes. This indicates that the composite diffusive/Powell method is the most effective method among the three methods in dealing with the magnetic field divergence.
Here, we use the following metric for measuring divergence, which was also adopted by other research studies (e.g., [33,47]): where M is the total number of grid points in the computational domain. We know that there are other metrics that can be used to measure the divergence. As pointed in [49], the metric defined by Eq. 14 may rely on the spatial resolution. However, in this simulation, we make the comparison among the three cases with the same mesh system and the same metric definition; therefore, the influence of the spatial resolution on the comparison of the metric by Eq. 14 can be ignored. Figure 8 shows the evolution of the Error(B) ave with time deduced from the three methods. It can be recognized that the value of the Error(B) ave from the composite diffusive/Powell method is around 10 -8.7 -10 -8.5 , from the diffusive method is around 10 -8.6 -10 -8.2 , and from the Powell method is around 10 -8.6 -10 -7.1 . The composite diffusive/Powell method has the smallest Error(B) ave , and this method is a new try to maintain the magnetic divergencefree constraint. From Figure 8, we can also find that the Error(B) ave from the composite diffusive/Powell method and diffusive method is smaller than that from the Powell method obviously. Moreover, the Error(B) ave from the composite diffusive/Powell method keeps on decreasing after 60 h and is significantly smaller than that from the diffusive method near 100 h. Overall, we can find that all the divergence cleaning methods can keep the related errors under control, though the divergence errors of the Powell method are larger than those of the other methods, the divergence errors shown in Figures 6-8 are indeed small, and the largest worst number is 10 -7 ,

CONCLUSIONS AND DISCUSSIONS
In this study, by using the 3D COIN-TVD MHD model, we simulate the solar wind in the coronal region, in which the divergence cleaning and coronal heating/acceleration methods are included. The volume heating method is an effective way for coronal heating, in which the parameters can be adjusted according to the WSA model in the simulation of the coronal solar wind. In the COIN-TVD MHD model, increasing the parameters S 0 and Q 0 of the energy and momentum source terms can make the solar wind accelerate more obviously. For the divergence cleaning methods, here we choose the diffusive method, the Powell method and the composite diffusive/Powell method. We compared the numerical characteristics of the combination of each method for handling the divergence of the magnetic field and the COIN-TVD MHD model in the solar coronal simulation. The numerical results show that all of them can produce large-scale structured solar wind and reduce the divergence of the magnetic field more or less. The difference between the three divergence cleaning methods is summarized as follows:  1) The Powell method is relatively simple to apply. It only needs to add two items to the source term of the MHD equations. In this study, the Powell method can reduce the error of the relative magnetic field divergence, but it is less effective than the other two methods in dealing with magnetic divergence. 2) The diffusive method also has a good effect on reducing magnetic field divergence error in this study. It reduces the error of divergence by adding a source term in the induction equation and the ∇ · B error is diffused away by iterating B k+1 B k + μ(Δx) 2 ∇∇ · B k . If it is coupled with different numerical schemes, the effects of controlling divergence error are different. In this study, the diffusive method is not as good as the composite diffusive/Powell method in controlling the divergence of the magnetic field, but better than the Powell method.
3) The composite diffusive/Powell method is a preliminary new try in this study, and it combines the Powell method and the diffusive method during the simulation. It has been proven   Frontiers in Physics | www.frontiersin.org July 2021 | Volume 9 | Article 705744 that this composite method is the most efficient way to reduce the relative divergence errors among the three methods we used. Moreover, it also ensures the conservation of the MHD equations during the simulation.
In addition to the methods we mentioned, there are many other methods to simulate the coronal heating and the solar wind acceleration process and to control the divergence of the magnetic field. For example, both the Alfvén wave heating method and the turbulent heating method are effective for coronal heating. The Powell method can also company with other methods to control the magnetic divergence, which may be implemented in the future. Moreover, although these simulations are performed for the background solar corona, these methods can also be used for the simulation of CME initiation and propagation in the interplanetary space.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/supplementary files, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
FS provided the thesis research topic, FS and YL provided the code for the three-dimensional solar wind numerical simulation. CL modified the code, ran the program code, and drew pictures based on the data. FS, CL, and MZ participated in the analysis of the results and the writing of the manuscript. XL modified the manuscript.