Extended MHD modeling of the steady solar corona and the solar wind

The history and present state of large-scale magnetohydrodynamic (MHD) modeling of the solar corona and the solar wind with steady or quasi-steady coronal physics is reviewed. We put the evolution of ideas leading to the recognition of the existence of an expanding solar atmosphere into historical context. The development and main features of the first generation of global corona and solar wind models are described in detail. This historical perspective is also applied to the present suite of global corona and solar wind models. We discuss the evolution of new ideas and their implementation into numerical simulation codes. We point out the scientific and computational challenges facing these models and discuss the ways various groups tried to overcome these challenges. Next, we discuss the latest, state-of-the art models and point to the expected next steps in modeling the corona and the interplanetary medium.


Introduction
It was realized thousands of years ago that the space between heavenly objects must be filled by something that is much lighter than the materials found on Earth. Aether (ancient Greek for light) was one of the primordial deities in Greek mythology. He was the personification of the upper air. Later the name "aether" was used by Greek philosophers to describe a very light fifth element. Plato mentioned that "there is the most translucent kind which is called by the name of aether." Aristotle introduced a new "first" element to the system of the classical elements (Earth, Fire, Air, Water). He noted that the four terrestrial classical elements were subject to change and naturally moved linearly. The first element however, located in the celestial regions and heavenly bodies, moved circularly and had none of the qualities the terrestrial classical elements had. It was neither hot nor cold, neither wet nor dry. With this addition the system of elements was extended to five and later commentators started referring to the new first one as the fifth and also called it aether. Medieval scholastic philosophers granted aether changes of density, in which the bodies of the planets were considered to be more dense than the medium which filled the rest of the universe. Later, scientists speculated about the existence of aether to explain light and gravity. Isaac Newton also suggested the existence of an aether (Newton 1718).
He wrote: "From a physical point of view it is most probable that solar rays are neither exclusively negative nor positive rays, but of both kinds." From his geomagnetic surveys, Birkeland (1908) realized that auroral activity was nearly uninterrupted, and he concluded that the Earth was being continually bombarded by "rays of electric corpuscles emitted by the Sun." His work was, however, pretty much ignored at the time. Some seventy years later, Dessler (1984) wrote an intriguing article arguing that personality conflicts between Sidney Chapman and Scandinavian scientists, as well as British imperial arrogance, were partly responsible for ignoring the pioneering works of Birkeland, Enskog, Størmer and Alfvén. In the early 1930s, Chapman and Ferraro (1931a, b, 1932a, b, 1933 published the first quantitative model of an infinitely conducting quasi-neutral plasma beam and its interaction with a magnetic dipole. Chapman was fundamentally a mathematician and he always tried to simplify physics phenomena to treatable mathematical problems. Ferraro was a PhD student looking for a thesis topic. Akasofu (1995) quotes Chapman on why he did not consider a continuous plasma stream from the Sun: "He [Lindeman] said it [a stream of gas] must consist of charges of opposite signs in practically equal numbers, so that it could hold together. Lindemann never tried to develop what would be the consequences on the Earth of the impact of such a stream of gas. I made an attempt at that while I was Professor at Manchester in 1919Manchester in -1924 unfortunately I started at the wrong end; I tried to find out what would be the steady state, as if the stream had been going on forever. It didn't work out; so I was still wanting to find out what would happen, and this was the subject I proposed for Ferraro." In plain English, Chapman did not consider the case of a continuous solar wind because he could not solve this problem mathematically.

The puzzle of comet tails
The idea of a continuous corpuscular radiance emanating from the Sun did not resurface until the 1950s. Interestingly, it was the observation of comet tails that lead to the reemergence of the idea. Comets exhibit two distinct types of tails. A nice example is comet Hale-Bopp, shown in Fig. 1.
By the 1950s, it had been recognized that the broad and curved dust tails (also called Type II tails) usually lag behind the Sun-comet line (opposite to the direction of cometary orbital motion). Since the gravity of the cometary nucleus is negligible, the motion of the dust particles is controlled by an interplay between solar gravity and solar radiation pressure. Assuming that dust grains have a more or less constant mass density, ρ d , the anti-sunward radiation pressure, F rad is inversely proportional to the square of the heliocentric distance, d h , and proportional to the cross section of the dust grain, πa 2 (a is the equivalent radius of a spherical dust grain), F rad ∝ a 2 /d 2 h . The sunward pointing gravitational force is proportional to the particle mass (4πρ d a 3 /3) and inversely proportional to the heliocentric distance: F grav ∝ a 3 /d 2 h . Since the two forces point in the opposite direction, have the same heliocentric dependence, but exhibit different dependence on the grain size (F rad ∝ a 2 and F grav ∝ a 3 ), the  Hale-Bopp (1997) showing two distinct tails-a broad dust tail (white) and a narrow ion tail (blue). Source: http://www.tivas.org.uk/solsys/ tas_solsys_comet.html resulting effect is a complex "dust mass spectrometer" where each particle size is moving under the influence of its own "reduced solar gravity." This effect results in the broad and curved dust tail.
For active comets the straight, narrow plasma tails (also called Type I tails) are 10 7 -10 8 km long and, within a few degrees, always point away from the Sun. Observation of various tail structures as they moved down the Type I tail determined that the acceleration in plasma tails ranged from about 30 to 300 cm/s 2 and occasionally even larger. This value was some three orders of magnitude larger than any acceleration due to solar radiation pressure. Some major process was missing. In order to account for the observed large acceleration Biermann (1951) postulated the existence of a continuous "solar corpuscular radiation" composed of electrons and ions. Assuming that the antisunward acceleration of small irregularities in Type I comet tails was due to Coulomb collisions between electrons in a radially outward plasma flow from the Sun and newly ionized cometary particles, Biermann (1951) inferred a solar wind density and velocity of n sw ≈ 1000 cm −3 and u sw ≈ 1000 km/s that represents a particle flux ∼ 500 times larger than was later observed.
The important consequence of Biermann's (1951) idea was that if solar corpuscular radiation is responsible for the antisolar acceleration of comet tails then the Sun evidently emits solar corpuscular radiation in all directions at all times. This follows from the fact that comets "fill" the heliosphere; antisunward pointing comet tails were observed over the poles of the Sun as well as at low heliographic latitudes. In addition, comets come by as frequently at sunspot minima as at sunspot maxima. Yet none fail to show an antisolar Type I tail. This means that interplanetary space must be completely filled with solar corpuscular radiation.
Bierman's idea was pretty much ignored by the solar physics community. The fact that the Sun has a million-degree corona was first discovered by Grotian (1939) and Edlén (1941) by identifying the coronal lines as transitions from low-lying metastable levels of the ground configuration of highly ionized metals (the green Fe xiv line at 530.3 nm, but also the red line Fe x at 637.4 nm). In the mid-1950s Chapman (1957) calculated the properties of a gas at such a temperature and determined it was such a superb conductor of heat that it must extend way out into space, beyond the orbit of Earth. However, Chapman-and others-considered a static corona that was gravitationally bound by the Sun. If the Earth is moving at about 30 km/s along its orbit around the Sun, the interaction between the Earth and the stationary solar corona could only result in very minor geomagnetic disturbances. Large geomagnetic storms, however, were already well observed at that time.

Solar wind
Sometime in 1957, Ludwig Biermann visited John Simpson's group at the University of Chicago. While in Chicago he had extensive discussions with one of Simpson's postdocs who was working on the problem of cosmic ray modulation in the solar system. Biermann explained his comet tail idea to the postdoc, Eugene Parker. According his own recollection, Parker started to think about the Biermann-Chapman puzzle: how to reconcile Chapman's (1957) hot, highly conducting corona with Biermann's (1951) idea of a continuously outward streaming fast solar corpuscular radiation. Parker's solution to the Biermann-Chapman puzzle was the idea of the continuous expansion of the hot solar corona, the solar wind (Parker 1958).
In his seminal paper, Parker (1958) pointed out that the static corona has a finite pressure at infinity that exceeds the pressure of the interstellar medium by a large factor. He concluded that the solar corona cannot be static [however, as it was later pointed out by Velli (1994) and Del Zanna et al. (1998), the situation is quite complicated]. Next, he considered a steady-state spherically symmetric isothermal corona that expands with a velocity v(r ), where r is the heliocentric distance. In addition, he assumed that the electron and ion temperatures were identical and approximated the plasma pressure by p = 2nkT 0 , where n(r ) is the ion number density, T 0 is the ion temperature, and k is the Boltzmann constant. Using the conservation of particle flux and the momentum equation, he obtained the following analytic relation for the plasma flow velocity: where a is the radius of the hot coronal base (Parker used a value of a = 10 6 km), v esc is the escape velocity at radius a (v 2 esc = 2G M /a, where G is the gravitational constant and M is the solar mass), while v m is the most probable velocity (v 2 m = 2kT 0 /m p , where m p is the proton mass). In Eq. (1), the integration constant was chosen to ensure that the velocity is real and positive for all r > a values. In addition, only solutions with v 2 0 2kT 0 /m and v(r → ∞) > v m solutions were considered, since the expansion velocity at the base of the hot corona, v 0 was assumed to be very small and Parker was interested in a hydrodynamic escape solution. This solution was later called the "solar wind." Figure 2 shows the class of escaping corona solutions (Parker 1958).
It is interesting to note that in Parker's (1958) solution the expansion velocity reaches the most probable particle velocity at the heliocentric distance of r c = av 2 esc /4v 2 m . Since it is assumed that the outflow velocity is small at the coronal base, this means that hydrodynamic outflow can only take place if the v 2 m > v 2 esc /4 condi-

Fig. 2
Spherically symmetric hydrodynamic expansion velocity v(r ) of an isothermal solar corona with temperature T 0 plotted as a function of r /a, where a is the radius of the base of the solar corona and has been taken to be 10 6 km. Image reproduced with permission from Parker (1958), copyright by AAS tion is met. In plain English, one needs a very hot corona for the solar wind to exist.
In his original paper Parker (1958) also considered the effect of the general outflow of solar gas upon the solar dipole magnetic field. He assumed that there are no field-free regions in the Sun, so that each volume element of gas flowing outward from the Sun will carry the embedded magnetic field lines with it. The field lines, being embedded in both the Sun and the ejected gas, will be stretched out radially as the gas moves away from the Sun. The radial configuration will be as universal as the outward-gas motion, which is responsible for it. Parker (1958) suggested that the gas flowing out from the Sun is not field-free but carries with it magnetic lines of force originating in the Sun. Hence, he predicted a radial solar magnetic field, falling off approximately as 1/r 2 in interplanetary space. The magnetic field lines will follow an Archimedean spiral described by where Ω is the angular velocity of the Sun, φ is the solar azimuth angle and φ s is the azimuth angle of the magnetic field line at the distance R s [Parker (1958) assumed a value of R s = 5R ]. Using this simple model Parker (1958) also expressed the magnetic field vector at an arbitrary point outside the r = R s sphere: where v s is the asymptotic speed of the solar wind that is reached at a heliocentric distance, where B (θ, φ s ) is the radial component of the magnetic field at the origin of the field line (on the r = R s sphere). At large heliocentric distances (where r R s ) the radial magnetic field component decreases as B r ∝ 1/r 2 , while the azimuthal component only decreases as B φ ∝ 1/r . This means that in the interplanetary space the magnetic field lines become more and more "wound up" as one moves further and further away from the Sun. This heliospheric magnetic field configuration is now referred to as the "Parker spiral." A schematic of the interplanetary magnetic field line topology is shown in Fig. 3 (Parker 1958).
Another important consequence of Parker's (1958) solution was that, due to the conservation of mass flux, it predicted that the particle density in the interplanetary medium would decrease as n ∝ 1/r 2 , much slower than the exponential decrease obtained in the case of a hydrostatic atmosphere. This result was consistent with Biermann's (1951) conclusions, but it was very different than the accepted model predictions.
Opposition to Parker's (1958) hypothesis on the solar wind was immediate and strong. He submitted the manuscript to the Astrophysical Journal in early January of 1958. Here is how Parker remembers the events (Parker 2001): …sometime in the late Spring of 1958 the referee's report on the original paper appeared, with the suggestion that the author familiarize himself with the subject before attempting to write a scientific article on solar corpuscular radiation. There was no specific objection to the arguments or calculations in the submitted paper. The author pointed out to the editor the absence of any substantive objections by the referee. The editor, Subrahmayan Chandrasekhar, sent the paper to a second referee. Sometime in the summer the second referee responded that the paper was misguided and recommended against publication, again with no specific criticism except that the author was obviously unfamiliar with the literature.
Again, the author's response to the editor was to note that there was no substantive criticism, no specific fault pointed out. Some days later Chandrasekhar appeared in the author's office with the paper in his hand and said something along the lines of, 'Now see here, Parker, do you really want to publish this paper?' I replied that I did. Whereupon he said, 'I have sent it to two eminent experts in the field and they have both said that the paper is misguided and completely off the mark.' I replied that my problem with the referees was that they were clearly displeased, but had nothing more to say. Chandrasekhar was silent for a moment and than he said, 'All right, I will publish it.' And that is how the paper without the words 'solar wind' finally appeared in the November issue of the Astrophysical Journal.
Just about the same time when Parker (1958) published his paper solving the Biermann-Chapman puzzle, Alfvén (1957) pointed out that the plasma densities inferred by Biermann (1951) were inconsistent with observed coronal densities (assuming that the plasma density decreases with the square of the heliocentric distance as the solar corpuscular radiation moves outward). Alfvén (1957) offered an alternative explanation for the formation of cometary plasma tails that is depicted in Fig. 4. The main assumption in this model was that the solar corpuscular radiation was carrying a "frozen-in" magnetic field that "hangs up" in the high density inner coma, where the solar particles strongly interact with the cometary atmosphere and consequently the solar plasma flow considerably slows down. This interaction results in a "folding" of the magnetic field around the cometary coma that creates the long plasma tail. Disturbances along the folded magnetic field lines propagate as magnetohydrodynamic waves and they can reach velocities of 100 km/s even if the solar plasma has a density of ∼ 5 cm −3 . It is interesting to note that Alfvén accepted Biermann's (1951) idea of the continuous plasma outflow from the Sun and naturally concluded

Pa rk er
Cham berla in that this plasma outflow must carry an embedded magnetic field, but he did not worry about the origin of such a plasma flow.

Solar breeze?
Parker's (1958) paper generated swift negative reaction. The community was not ready to give up the idea of a hydrostatic corona and accept a continuously escaping solar atmosphere.
The most prominent critic of Parker was Joseph Chamberlain, who received his PhD from the University of Michigan in 1952. At that time he was on the faculty of the University of Chicago, the same institution where Chandrasekhar and Parker worked. In September 1959, less than a year after Parker's (1958) paper was published, he submitted a paper which strongly criticized Parker's (1958) solar wind solution. Chamberlain (1960) pointed out that Eq. (1) had two solutions that met the condition that the velocity was small at the coronal base. In addition to Parker's (1958) solution, a second solution described a solar corona that started to expand, but the expansion gradually stopped beyond the critical point (r c = av 2 esc /4v 2 m ). A comparison of the Parker (1958) and Chamberlain (1960) solutions is shown in Fig. 5. Chamberlain's (1960) solution provided a hydrostatic coronal density distribution and the expansion velocity returned to zero at infinity. This solution "saved" most features of the prevailing hydrostatic atmosphere model and it was greeted with great relief by most of the solar community. The disagreement between Parker and Chamberlain became quite heated. In 1959 when the debate took place Parker was a junior (untenured) faculty member at the University of Chicago and, as a result of the ongoing controversy with a more senior faculty member, his tenure case became more challenging. The controversy was eventually resolved-at least in the eyes of the space physics community-by the beginning of the space age, when in-situ observations proved the existence of the solar wind. Gringauz et al. (1960) analyzed the results of the Lunik-2 spacecraft and determined an interplanetary corpuscular particle flux of about 2 × 10 8 ions/cm 2 /s: "Starting from 9.30 hours Moscow time on 13 September 1959 up to the moment of the container of the second space rocket reaching the moon the container was recorded as passing through a positive ion flux (in all probability protons) with energies exceeding 15 eV; Φ ∼ 2 × 10 8 ions/cm 2 /s." In a May 1962 presentation in Washington, Gringauz et al. (1964) used observations of the first deep space probe (Venera-1) to estimate the speed of the solar corpuscular radiation to be about 400 km/s. In addition, measurements by Bonetti et al. (1963) on Explorer 10 confirmed these initial results. However, because the Lunik, Venera and Explorer observations were short-term, the doubters had some wiggle-room and the observations were not regarded as definitive.
It was not until the end of 1962, when the first Mariner 2 results were reported (Neugebauer and Snyder 1966), that the existence of the solar wind was widely accepted. The observed solar wind had typical proton densities of 5-20 cm −3 and velocities between 300 and 700 km/s (see Fig. 6). The observed temperature was in the range of 3 × 10 4 -3 × 10 5 K. These observations confirmed the predictions of the Biermann (1951)-Alfvén (1957)-Parker (1958 theory and proved that the concept of a static, slowly evaporating corona was incorrect (Chapman 1957;Chamberlain 1960).
It is interesting to note Chamberlain's reaction to the eventual acceptance of the solar wind concept. Even some 30 years later, he tried to show that those who advocated the solar wind concept were right for the wrong reasons and they advocated inappropriate solutions. Here is a quote from Chamberlain's (1995) paper: In the early days of solar-wind theory, Parker (1958) appeared to be influenced by two seminal hypotheses: (a) Biermann's conclusion that the behavior of comet tails was governed more by the interplanetary medium than by solar-radiation pressure, and (b) Chapman's advocacy of a static solar corona that was heated to great distances by conduction. Parker showed that a static corona was untenable and then constructed a primitive hydrodynamic model, which he labeled the solar wind, that would account for Biermann's analysis of comet tails. I describe this solar-wind model as 'primitive' because its temperature distribution, instead of being derived physically, was characterized by a polytropic index, which simplified the model enormously.
At this point I was skeptical that Parker's supersonic solutions were realistic (Chamberlain 1960) and, to investigate the problem, I wrote down the three equations -now called the solar-wind equations -for a hydrodynamic solar corona that was heated from below by thermal conduction (Chamberlain 1961). I restricted my investigation to slow (subsonic) expansion, or models in which the parameter ε (described below) was 0. ln a jocular spirit, I referred to those solutions as solar-breeze models, and I suggested that they were merely the hydrodynamic counterpart to Waterston-Jeans evaporative escape.
Shortly afterward, a solar-wind model, based on the same three equations but encompassing supersonic expansion, was developed by Scarf and Noble. (The pair of papers, Noble and Scarf (1963) and Scarf and Noble (1965), are hereinafter abbreviated by 'SN-I' and 'SN-II,' respectively.) But their model, although quite acceptable as an illustrative case, is not physically accurate, even within the constraints of its own simplifying assumptions.   Neugebauer and Snyder (1966), copyright by AGU

Numerical solution
The first numerical solution of the solar wind equations was carried out by Scarf and Noble (Noble and Scarf 1963;Scarf and Noble 1965), who solved the spherically Fig. 7 Comparison between simulated and observed electron densities in the inner corona (neglecting heat conduction and viscous effects). Image reproduced with permission from Scarf and Noble (1965), copyright by AAS symmetric Navier-Stokes equations including heat conduction and viscosity effects. Instead of considering a polytropic corona, they solved the energy equation, thus making the problem more challenging. The set of differential equations had a singularity when the local solar wind speed, v(r ), reached the value of kT (r )/m p . In order to avoid numerical problems, Noble and Scarf (1963) integrated the equations from the Earth inward. This was made possible by the fact that at that time the solar wind conditions were observed. Noble and Scarf (1963) chose 1 AU values of n = 3.4 cm −3 , v = 352 km/s, and T = 2.77 × 10 5 K. A typical solution without heat conduction and viscous effects is shown in Fig. 7. We note that even in this relatively simple case the transonic solar wind solution describes the observed electron density profile within a factor of 2 or 3, a surprisingly good agreement. Scarf and Noble (1965) also considered the "subsonic solution," representing the solar breeze (Chamberlain 1960). When they took into consideration heat conduction and viscosity, Scarf and Noble (1965) were able to get excellent agreement with observations (however, in effect, they increased the number of free parameters, so the improved agreement is not really surprising). A similar solution with heat conduction but without viscosity was also obtained by Whang and Chang (1965).

Two-fluid model
The first two-fluid (separate electron and proton equations) was published by Sturrock and Hartle (1966). They realized that as the solar wind leaves the vicinity of the Sun collisional coupling between electrons and ions becomes weak and the electron and ion temperatures (T e and T i ) can deviate from each other. They developed a twotemperature model where the plasma remains quasi-neutral (n e ≈ n i ) and current-free (u e = u i ), but the two temperatures can be different. Sturrock and Hartle (1966) still considered a spherically symmetric problem and neglected magnetic field effects, but their model represented a step forward. By combining the continuity, momentum, and energy equations they derived separate "heat equations" for electrons and ions (Sturrock and Hartle 1966): where s = e, i refers to either electrons or ions, the t subscript refers to the other species, Φ = nvr 2 is the spherical particle flux, v is the plasma flow speed, κ s is the heat conductivity of the appropriate species and ν ei is the electron-ion energy transfer collision frequency. Sturrock and Hartle (1966) used the Chapman (1954) approximation for the collision frequency (ν ei = 8.5 × 10 −2 nT −3/2 e , where n is given in units of cm −3 ) and the Spitzer (1962) Figure 8 shows the two-temperature solar wind solution (Sturrock and Hartle 1966). At Earth orbit they obtained v = 270 km/s, n = 13 cm −3 , T i = 2800 K and T e = 4.6 × 10 5 K. While these values were in the right ballpark, the wind speed was too slow, the density was too high, and the two temperatures were quite a bit off (T e too high and T i too low). The authors lamented that "We are left with the problem of temperatures as a function of radial distance. Image reproduced with permission from Sturrock and Hartle (1966), copyright by APS understanding why the solar-wind parameters do not always agree with our model." In other words, it was Nature's fault that it did not agree with the predictions of the model.

Potential magnetic field
By the late 1960s, it has become clear that the magnetic field plays a major role in determining the density as well as the velocity and temperature structure of the corona. The first models used the observed line-of-sight component of the photospheric magnetic field to determine the magnetic field of the solar corona in the current-free (or potential-field) approximation (Schatten , 1969(Schatten , 1971Schatten et al. 1969;Altschuler and Newkirk 1969;Newkirk and Altschuler 1970).
The potential field model is based on the fundamental assumption that the magnetic field above the photosphere is current free (∇ × B = 0 when r > R ) and therefore the coronal magnetic field can be represented by a magnetic scalar potential, ψ: Since there are no magnetic monopoles (∇ · B = 0) we obtain The solution of Eq. (6) can be written as an infinite series of spherical harmonics: where the coefficients g m n and h m n need to be determined from solar line-of-sight observations and P m n (θ ) are the associated Legendre polynomials with Schmidt normalization: With the help of the magnetic scalar potential, the magnetic field vector can be obtained anywhere above the photosphere: It was recognized by Schatten et al. (1969) that the coronal magnetic field follows the current-free potential solution between the photosphere and a "source surface" (located at r = R s ) where the potential vanishes and the magnetic field becomes radial. This requires a network of thin current sheets for r ≥ R s (cf. Schatten 1971).
Such a potential field can be described with the help of Legendre polynomials that define the magnetic potential between two spherical shells, located at r = R and r = R s , each containing a distribution of magnetic sources (Altschuler and Newkirk 1969): where For the optimal radius of the source surface Schatten et al. (1969) found R s = 1.6 R , while Altschuler and Newkirk (1969) recommended R s = 2.5 R . Today, most potential field source surface (PFSS) models use the R s = 2.5 R value. We note that at the source surface With the help of expression (10), the magnetic field components can be written as At the source surface the magnetic field becomes radial because f n (R s ) = 0 and g m n cos mφ + h m n sin mφ P m n (θ ) .
(16) Fig. 9 Magnetic field line map obtained with the potential field source surface (PFSS) model. Image reproduced with permission from Altschuler and Newkirk (1969), copyright by Springer Outside the source surface it is assumed that the radial flow of the solar wind carries the magnetic field outward into the heliosphere. This region is not described by the PFSS model. Between the photosphere and the source surface the magnetic field vector components can be described by expressions (13), (14) and (15). The inner boundary condition at the photosphere is obtained from the observed line-of-sight magnetic field components using a least square fit (cf. Hoeksema et al. 1982). These measurements are used to determine the expansion coefficients g m n and h m n . The outer boundary condition at the source surface is that the field is normal to the source surface, consistent with the assumption that it is then carried outward by the solar wind. This condition is automatically satisfied by our selection of the f n function [Eq. (11)]. An example of the PFSS solution is shown in Fig. 9.

Expanding magnetized corona
Nearly simultaneously with the development of the source surface models the effect of the coronal expansion on the magnetic field was also explored (Pneuman 1966(Pneuman , 1968(Pneuman , 1969Pneuman and Kopp 1971). In a series of papers, Pneuman (1966Pneuman ( , 1968Pneuman ( , 1969 analytically investigated how centrifugal, pressure gradient, and magnetic forces impact the flow of an infinitely conducting fluid where the magnetic field lines and plasma flow lines must coincide in the corotating frame. An early example of a numerical model of the corona with solar wind is by Pneuman and Kopp (1971). As in the Parker solution (Parker 1963), this model possesses the high temperature (T > 1 MK) and comparatively high density (ρ ≈ 10 −16 g/cm 3 ) plasma whose pressure cannot be held in equilibrium by solar gravity or the pressure of the interstellar medium. Consequently, the coronal plasma expands rapidly outward achieving supersonic speeds within a few solar radii, and in doing so forms the solar wind.
Early numerical models prescribed volumetric coronal heating in ways that strived to accounted for the effects thermal conduction and radiative losses, as well as sat- Fig. 10 Analytic solution of a helmet configuration. Image reproduced with permission from Pneuman (1968), copyright by Springer isfying known constraints of coronal heating. However, these works found that the observed fast solar wind (speed 700-800 km/s) cannot be produced by thermal pressure without temperatures greatly exceeding coronal values, particularly in coronal holes where the fast wind originates.
Using simplified conservation laws, Pneuman (1966Pneuman ( , 1968Pneuman ( , 1969 obtained helmetlike closed field line regions buffeted by converging open field lines (Fig. 10). Later, Pneuman and Kopp (1971) numerically solved the conservation equations for an axially symmetric steadily expanding corona with an embedded magnetic dipole, assuming North-South symmetry. In this case the solution is not only axially symmetric, but it is also symmetric with respect to the equatorial plane.
The Pneuman and Kopp (1971) solution shows several interesting features. An important aspect of the solution is the development of a neutral point at the top of the helmet where the magnetic field vanishes. A current sheet forms between the regions of opposite magnetic polarity above the neutral point. The outflowing plasma reaches the Alfvén velocity just above the neutral point (at the top of the helmet), so the top of the helmet is, in effect, the transition from sub-Alfvénic to super-Alfvénic flow. Figure 11 shows a comparison of the field lines of the numerical solution and a potential field model with the same normal component at the reference level (Pneuman and Kopp 1971). To make the comparison meaningful, the field lines for the two cases are chosen so as to be coincident at the surface. As expected, the field is everywhere stretched outward by the gas with this distention becoming large near the neutral point. In the closed region well below the cusp the difference between the two configurations is small, mainly because the pressure at the reference level is taken to be independent of latitude. The outward distention of the field in this region, as a result, is produced solely by currents generated by the expansion along open field lines. For the general case of variable surface pressure, however, a pressure and gravitational force balance cannot be satisfied normal to the field and significant j × B forces are expected to maintain the equilibrium.
In spite of its obvious limitations, the Pneuman and Kopp (1971) simulation established the usefulness of MHD simulations of the solar corona. It was the first successful attempt to apply the conservation laws of magnetohydrodynamics to explain largescale features of the solar corona. Fig. 11 A comparison of the magnetic field (solid curves) with a potential field (dashed curves) having the same normal component at the reference level. The field lines for the two configurations are chosen so as to be coincident at the surface. Image reproduced with permission from Pneuman and Kopp (1971), copyright by Springer

Steady-state models of the solar wind
We adopt the "ambient-transient" paradigm to modeling the dynamic and highly structured solar wind. By "ambient" we mean a quiet-sun-driven full 3-D structure for the interplanetary magnetic field and a 3-D distribution of the solar wind parameters. They are both close to being steady-state in the frame of reference co-rotating with the Sun, except for the highly intermittent solar wind region. This ambient solution determines the bimodal structure of the solar wind. It affects magnetic connectivity between the active regions at the Sun and the corresponding regions in the heliosphere. In turn such connectivity affects energetic particle acceleration and transport.
As we discussed in Sect. 3, the first generation of magnetohydrodynamic models of the interplanetary medium were developed in the second half of the 1970s and were used for about two decades (Steinolfson et al. 1975(Steinolfson et al. , 1978(Steinolfson et al. , 1982Pizzo 1978Pizzo , 1980Pizzo , 1982Pizzo , 1989Pizzo , 1991Pizzo , 1994aSteinolfson 1988Steinolfson , 1990Steinolfson , 1992Steinolfson , 1994Steinolfson and Hundhausen 1988;Pizzo et al. 1993;Pizzo and Gosling 1994). These models were designed to describe only large-scale bulk-average features of the plasma observed through the solar cycle. At solar minimum, these coronal structures are the following: 1. open magnetic field lines forming coronal holes; 2. closed magnetic field lines forming a streamer belt at low latitudes; 3. the bimodal nature of the solar wind is reproduced with fast wind originating from coronal holes over the poles and slow wind at low latitudes.
A thin current sheet forms at the tip of the streamer belt and separates opposite directed magnetic flux originating from the two poles. At solar minimum, the fast wind lies at 30 degrees heliographic latitude and has an average velocity of 750 km/s at distances greater than 15 solar radii, at which distance the wind has attained the majority of its terminal velocity. The slow wind, by contrast, is confined close to the global heliospheric current sheet, which lies near the equator at solar minimum. This component of the wind is highly variable, with speeds that lie between 300 and 450 km/s. The slow solar wind has been suggested by Wang and Sheeley (1990) to originate from highly expanding plasma traveling down magnetic flux tubes that originate near coronal hole boundaries. It has also been suggested that opening of closed flux tubes by interchange reconnection with open flux may release plasma to form the slow solar wind (Fisk et al. 1998;Lionello et al. 2005;Rappazzo et al. 2012). More recent theories have related the slow solar wind to complex magnetic topology flux tubes near the heliospheric current sheet, which are characterized by the squashing factor Antiochos et al. 2012). At solar maximum, the current sheet is highly inclined with smaller coronal holes forming at all latitudes, while the fast wind is largely absent. Steinolfson et al. (1975) numerically solved the MHD equations for a spherically symmetric (2D) solar corona, neglecting the polar components of velocity and magnetic field. They solved for ρ, v r , v φ , p, B r , and B φ as a function of time and radial distance. The model was applied to a forward-reverse shock pair propagating in this simplified solar wind and the solution was compared with the results of similarity theory. In the end Steinolfson et al. (1975) concluded that MHD effects might be important in the dynamical behavior of the solar wind near Earth orbit. In a follow-up paper, Steinolfson et al. (1978) considered a different situation, when the flow and magnetic field are in the meridional plane. They solved time evolution equations for ρ, v r , v θ , p, B r , and B θ as a function of polar angle and radial distance. They considered solar transients in two magnetic configurations: with the magnetic field radial (open), and with the magnetic field parallel to the solar surface (closed). The solar event is simulated by a pressure pulse at the base of an initially hydrostatic atmosphere. The pressure pulse ejects material into the low corona and produces a disturbance that propagates radially and laterally through the corona. The disturbance is preceded by waves (which may strengthen into shocks) that travel in the meridional plane with the shape of an expanding loop. The portion of the disturbance between the ejected material and the preceding waves consists entirely of coronal material whose properties have been altered by the waves. This simulation, in spite of its many limitations, was the first successful attempt to simulate coronal transients.

2D and quasi-3D models
A few years later, Steinolfson et al. (1982) revisited the steady global solar corona that was investigated a decade earlier by Pneuman and Kopp (1971). This more sophisticated simulation allowed for non-constant temperature and did not require specific conditions to be met at the cusp (Pneuman and Kopp 1971). They applied their previous model (Steinolfson et al. 1978) to study the axially symmetric global corona. Their solution essentially confirmed the configuration obtained by Pneuman and Kopp (1971) with some notable differences. Since the coronal plasma is subsonic and sub-Alfvénic near the Sun the flow is considerably more complex than predicted by Pneuman and Kopp (1971). The main feature, however, is essentially the same: Fig. 12 Steady-state coronal structure for a plasma β value of 0.5. Solid lines represent magnetic field lines, while arrows show plasma flow velocity vectors. The horizontal axis is the solar equator. Image adapted from Steinolfson et al. (1982) the coronal outflow results in closed field lines near the equator and open field lines (coronal holes) at high latitudes (see Fig. 12).
In parallel with Steinolfson's efforts, Pizzo took a very different approach to simulating the interplanetary medium (Pizzo 1978(Pizzo , 1980(Pizzo , 1982(Pizzo , 1989(Pizzo , 1991(Pizzo , 1994aPizzo et al. 1993;Pizzo and Gosling 1994). As a first step, Pizzo (1978) developed a 3D hydrodynamic model of steady corotating streams in the solar wind, assuming a supersonic, inviscid and polytropic flow beyond approximately 35 R s . This approach takes advantage of the fact that in the inertial frame the temporal and azimuthal gradients are related by ∂ t = −Ω ∂ φ , where φ is the azimuth angle and Ω is the equatorial angular velocity of the Sun.
The steady-state conservation equations were solved using an explicit Eulerian approach on a rectangular (θ, φ) grid covering the 45 • ≤ θ ≤ 135 • and −60 • ≤ φ ≤ 60 • spherical surface. Periodic boundary conditions were imposed at the azimuthal edges of the mesh, while the latitudinal boundaries are free surfaces, with the meridional derivatives approximated by one-sided differences. The density, velocity vector, and scalar pressure were defined at the inner boundary (35 R ). Knowing the solution at a heliocentric distance r , the solution was advanced to r +Δr using the conservation equations. Using a value of Δr = 30 km, Pizzo (1978) obtained a steadily corotating stream structure between 35 R and 1 AU.
Specifically, Pizzo (1978Pizzo ( , 1980 solved the governing equations that describe the dynamical evolution of 3-D corotating solar wind structures. The model is limited to those structures that are steady or nearly steady in the frame rotating with the Sun and utilizes the single-fluid, polytropic, nonlinear, 3D hydrodynamic equations to approximate the dynamics that occur in interplanetary space, where the flow is supersonic and the governing equations are hyperbolic. In the inertial frame, the equations are the following (Pizzo 1978(Pizzo , 1980: where ρ is the mass density, u is the center of mass velocity, p is the total isotropic (scalar) gas pressure, G is the gravitational constant, M is the solar mass, γ is the polytropic index and Ω is the equatorial angular rotation rate of the Sun. The independent variables are the spherical polar coordinates (r , θ, φ). Conduction, wave dissipation, differential rotation, the magnetic field, and shock heating are all neglected. Equations (17)- (19) can be rearranged to obtain a set of differential equations describing the radial evolution of the primitive variables and solved by marching in radial distance from the inner boundary outward. This method yields a 3D solution that is steadystate in the corotating frame and describes the stream structure in the interplanetary medium.
In a subsequent paper, Pizzo (1982) extended his marching method to include the interplanetary magnetic field. He considered the implication of a high-speed stream emanating from the Sun and expanding to 1 AU. The central portion of the stream is a circular plateau some 30 • in diameter where the radial velocity is a uniform 600 km/s. The speed falls off smoothly in all directions, bottoming out at the background speed of 300 km/s at a distance of 7.5 • from the periphery of the plateau. The results at 1 AU are shown in Fig. 13. The panels from top to bottom show the bulk speed, the flow angle in the inertial frame, the flow angle in the rotating frame, the number density, the single-fluid temperature, the magnetic field intensity and the total pressure as a function of azimuth. At the front of the stream (west, or left), forward (F) and reverse (R) waves have formed about an interface (I). The two waves, which propagate in opposite directions from the interface, originate in the general compression at the stream front near the Sun.
The Pizzo (1982) model was later applied to qualitatively explain Ulysses observations at larger helio-latitudes (Pizzo and Gosling 1994). Ulysses discovered that the forward-reverse shock pair structure commonly bordering corotating interaction regions beyond 1 AU near the ecliptic plane undergoes a profound change near the maximum heliographic latitude of the heliospheric current sheet. At this latitude the forward shock is observed to weaken abruptly, appearing as a broad forward wave, while the reverse shock weakens much more slowly (Gosling et al. 1993). These observations were well reproduced by the Pizzo (1982) model (Pizzo and Gosling 1994).

Connecting the corona and the heliosphere
Probably the most challenging region to model is the transition from the cold (∼ 5000 K) photosphere to the million K solar corona. This narrow, but critical region has very complex physics and challenging numerics: the plasma temperature increases by a Fig. 13 Equatorial solution for an initially steep-sided circular stream at I AU. F and R mark the forward and reverse wave fronts, which demarcate the dynamic compression ridge. The two shocks propagate in opposite directions, but both are convected outward in the bulk flow. The interface, I, is a shear layer separating fast and slow flow regimes that initially had very different fluid properties. Image reproduced with permission from Pizzo (1982), copyright by AGU factor of several hundred over ∼ 500 km (< 10 −3 R ). The first attempts to include this physics in 3D simulations are just beginning (see Sect. 5 in this paper and Sokolov et al. 2016).
In the early 1990s it was recognized that the solar wind speed at 1 AU negatively correlates with a magnetic flux tube expansion factor near the Sun (Wang and Sheeley 1990, 1992. This expansion factor describes the ratio of a given flux tube's cross sectional area at some heliocentric distance and the cross sectional area of the same flux tube at the solar surface. Wang and Sheeley (1990, 1992 found that the solar wind speed at a heliocentric distance can be expressed as where u min and u max represent the slow and fast solar wind speeds, f exp is the expansion factor at heliocentric location r, while α is an empirical factor (near unity). Later Arge and Pizzo (2000) and Arge et al. (2004) generalized the Wang and Sheeley (1990, 1992 formula and introduced several additional parameters. A comprehensive description of the Wang-Sheeley-Arge (WSA) model and its parameter values can be found in Riley et al. (2015).
The WSA formula provides an efficient and simple way to circumvent the complex physics of the low corona and transition region. It is used by a number of heliosphere models to provide inner boundary conditions beyond the Alfvén surface (where the solar wind speed exceeds the local Alfvén speed). These models typically place their inner boundaries in the 20-30 R range (e.g., Wold et al. 2018).

3D MHD heliosphere models
Pizzo's early work was followed with the development of the ENLIL heliospheric model by Odstrčil and Pizzo (1999), who applied it to study several phenomena (Odstrčil et al. 1996;Odstrčil and Karlicky 1997) including coronal mass ejection (CME) propagation (Odstrčil and Pizzo 1999;Odstrčil et al. 2004bOdstrčil et al. , 2005Siscoe and Odstrčil 2008). ENLIL has been adapted to accept inner boundary solar wind conditions from a variety of sources, including the WSA model (Wang and Sheeley 1990;Arge et al. 2004) and coronal MHD models (Odstrčil et al. 2002(Odstrčil et al. , 2004a. In the case of Hayashi (2012), the inner boundary conditions (outside the critical point) are derived from interplanetary scintillation (IPS) observations (Jackson et al. 1998). More recently, other groups also developed 3D inner heliosphere models with super-Alfvénic inner boundary conditions using the WSA approach. These include the LFM-Helio model (Merkin et al. 2011(Merkin et al. , 2016, the SUSANOO-SW code (Shiota et al. 2014;Shiota and Kataoka 2016), the MS-FLUKSS suite (Kim et al. 2016) and EUHFORIA developed by the Leuwen group (Pomoell and Poedts 2018).
WSA-like empirical inner boundary conditions were also used in the first generation of outer heliosphere models describing the interaction between the solar wind and the interstellar medium. Linde et al. (1998) published the first 3D MHD model describing the interaction of the magnetized solar wind with the magnetized interstellar medium. This simulation also took into account the presence of the neutral component of the interstellar medium and the resulting mass loading process inside the heliosphere. Figure 14 presents a 3D view of the global heliosphere.
Following the early work of Linde et al. (1998), several groups developed sophisticated 3D models of the outer heliosphere and its interaction with the magnetized interstellar medium. While details of these models go beyond the scope of this review we note the progress made by Opher et al. (2003Opher et al. ( , 2006Opher et al. ( , 2007Opher et al. ( , 2016 and the Huntsville group (Pogorelov et al. 2013(Pogorelov et al. , 2015Zirnstein et al. 2017).

3D MHD coronal models
The next level of sophistication of MHD coronal models came with reduced (Usmanov 1993;Mikić et al. 1999) and variable adiabatic index models (e.g., Wu et al. 1999;Roussev et al. 2003;Cohen et al. 2007). These models approximate coronal heating by greatly lowering the adiabatic index from 5/3 to typically 1.05 so that the plasma remains nearly isothermal as it expands and maintains the pressure necessary to drive The color code shows the log of plasma density. Yellow lines are the plasma velocity streamlines and the white lines follow the magnetic field lines. Black arrows indicate the direction of the magnetic field along a cross-tail cut placed 225 AU downstream from the Sun. Image reproduced with permission from Linde et al. (1998), copyright by AGU a fast wind. Such models have been successful in largely reproducing the observed solar wind speed at 1 AU over an entire solar cycle.
A number of coronal heating models have been published over the years using two fundamentally different approaches: (i) use a general heating function with an ad-hoc heating rate that is chosen to fit observations; or (ii) include a semi-empirical coronal heating function that is based on the physics of Alfvén waves. Examples for the first approach include papers by Groth et al. (1999aGroth et al. ( , 2000a, Lionello et al. (2001Lionello et al. ( , 2009, Riley et al. (2006), Feng et al. (2007Feng et al. ( , 2010, Nakamizo et al. (2009), Titov et al. (2008, and Downs et al. (2010). An important limitation of this approach is that models utilizing an ad-hoc approach depend on some free parameters that need to be determined for various solar conditions. While the ad-hoc heating function approach is well-suited for typical conditions, it can't properly account for unusual conditions, such as those that can occur during extreme solar events.
A major benefit of the WSA model is its simplicity and therefore it can be seamlessly integrated into global-scale space weather simulations (e.g., Odstrčil 2003;Cohen et al. 2007). In the Cohen et al. (2007) study the WSA formulae were used as the boundary condition for a large-scale 3-D MHD simulation with varied polytropic gas index distribution (see Roussev et al. 2003). These models can successfully reproduce observed solar wind parameters at 1 AU.
An example of a solar model driven by empirical heating is shown in Fig. 15, which depicts solar minimum conditions on a meridional plane close to the Sun. This result from Manchester et al. (2004) is based on a model by Groth et al. (2000a). Here the color image indicates the velocity magnitude, |u|, of the plasma while the magnetic field is represented by solid white lines. Inspection reveals a bimodal outflow pattern with slow wind leaving the Sun below 400 km/s near the equator and high-speed wind above 700 km/s found above 30 • latitude. In this model, we find that the source of the slow solar wind is plasma originating from the coronal hole boundaries that overexpands and fails to accelerate to high speed as it fills the volume of space radially above the streamer belt. This model is consistent with the empirical model of the solar wind proposed by Wang and Sheeley (1994) that explains solar wind speeds as being inversely related to the expansion of contained magnetic flux tubes. The magnetic field remains closed at low latitude close to the Sun, forming a streamer belt. At high latitude, the magnetic field is carried out with the solar wind to achieve an open configuration. Closer to the equator, closed loops are drawn out and, at a distance (r > 3R ), collapse into a field reversal layer. The resulting field configuration has a neutral line and a current sheet originating at the tip of the streamer belt. Several validation and comparison studies have been published using this approach (Owens et al. 2008;Vásquez et al. 2008;MacNeice 2009;Norquist and Meeks 2010;Gressl et al. 2014;Jian et al. 2015;Reiss et al. 2016). However, these models do not capture the physics of Alfvén wave turbulence or even neglect it altogether. Even though some of these models were designed to account for the Alfvén wave physics (e.g., Cohen et al. 2007), they do not capture many aspects of the interaction of the turbulence with the background plasma flow, which include both energy and momentum transfer from the turbulence to the solar wind plasma. Because Alfvénic turbulence effects are likely to be of great significance in the near-Sun domain, these simpler models should be used with caution for simulating the solar atmosphere.

Thermodynamic corona models
While the polytropic solar corona models were quite successful in describing the quiet low latitude corona they had major problems in accounting for the two-state (slow and fast) solar wind and large dynamical processes (such as CMEs) that can interrupt the quasi-steady state (in the corotating frame) situation. In particular, the shock thermodynamics got seriously distorted by the reduced adiabatic index and various spurious phenomena (like runaway plasma temperatures) appeared in the solutions.
In order to overcome this problem the full energy equation-with all accompanying computational challenges-needs to be considered. Mikić et al. (1999) proposed that the energy equation be solved with a realistic adiabatic index (γ = 5/3) and an empirical heating function be introduced to account for the effects of heat conduction, radiative cooling, and various heating processes. In particular, Mikić et al. (1999) introduced the following form of the heating function: where H ch is the coronal heating source, D is the Alfvén wave dissipation term, H d = η J 2 + ν∇v : ∇v represents heating due to viscous and resistive dissipation, and Q(T ) is the radiative loss function. In the collisional regime (below ∼ 10R ), the heat flux is q = b(b · ∇)T , where b is the unit vector along the magnetic field vector, and κ = 9 × 10 7 T 5/2 is the Spitzer value of the parallel thermal conductivity. The polytropic index γ is 5/3. In the collisionless regime (beyond ∼ 10R ), the heat flux is modeled by q = αn e kT v, where α is an empirical parameter. The coronal heating source is a parameterized function given in the form of where the empirical functions H 0 (θ ) and λ(θ ) express the latitudinal variation of the volumetric heating and scale length, respectively. While this "thermodynamic" approach sidesteps the underlying physics of coronal heating and solar wind acceleration it provides an adequate mathematical framework to describe the coronal processes in a way that is consistent with solar dynamical processes.
The thermodynamic coronal model has been successfully applied to simulate the solar corona for the first whole Sun month (Lionello et al. 2009). The global magnetic configuration obtained with the model is shown in Fig. 16.

Model inputs
The magnetic field is an essential component of the corona. As a matter of fact, the surface magnetic field distribution is the primary direct quantitative observable for models, physical variables associated with other observations need to be inferred. Early models used simple dipolar magnetic fields to simulate solar minimum conditions (Groth et al. 2000b;Steinolfson et al. 1982). In this case, the dipole is chosen such that the maximum field strength at the poles approximately 10 gauss. However, simulating realistic solar conditions requires the use of the observed line-of-sight global magnetic field. To provide these data, full disk magnetograms are taken as the Sun completes a full rotation, and then combined into a full-surface synoptic map. Early examples of these maps include those provided by the Stanford Wilcox and Mount Wilson observatories as well as the Global Oscillation Network Group (GONG). These data sources are further augmented by offerings from the National Solar Observatory's  (Carrington Rotation 1913 August 22 to September 18). The simulation was carried out with the thermodynamic model described in Sect. 4.5. The magnetic flux distribution was projected on the solar surface with selected magnetic field lines from the MHD solution. Image reproduced with permission from Lionello et al. (2009), copyright by AAS Synoptic Long-term Investigation of the Sun (SOLIS), and NASA's Solar Dynamics Observatory/Helioseismic and Magnetic Imager (HMI). Figure 17 provides an example input of a full-surface synoptic map based on GONG data. The use of synoptic maps in global MHD models was pioneered by Wu et al. (1999), Mikić et al. (1999) and Linker et al. (1999) and applied by many others (e.g., Roussev et al. 2003;Hayashi 2005;Feng et al. 2007;Cohen et al. 2007;Nakamizo et al. 2009;. While the line-of-sight photospheric magnetic field can be accurately measured and extrapolated to coronal heights with reasonable accuracy, the base densities and temperatures remain much more difficult to accertain, especially on the global scale required to specify boundary conditions for numerical models. In the case of van der Holst et al. (2010), the base temperature and mass density were derived empirically from the differential emission measure tomography (DEMT) technique by Vásquez et al. (2008). While promising, this approach requires an involved, time-consuming calculation not suitable for operational use. In recent years, a physics-based approach has been used to self-consistently calculate the thermodynamic properties at the base of the corona by applying field-aligned electron heat conduction and radiative processes (typically derived from CHIANTI; Dere et al. 1997). With this approach, it is possible to self-consistently reproduce the transition region. With that region, and the appropriate coronal density and temperatures, realistic heating functions can be applied (Lionello et al. 2011;Sokolov et al. 2013;van der Holst et al. 2014).

Model validation
3D coronal models applying synoptic magnetograms, in conjunction with thermodynamic processes such as field-aligned heat conduction, are capable of both predicting and interpreting the detailed magnetic and plasma structure of the corona. Successful comparisons began with the predicted appearance of the corona in Thomson-scattered white light, compared to coronagraph and eclipse images . Figure 18 provides a very impressive prediction of white-light images of a solar eclipses (Mikić et al. 2007). Models capture the low-density coronal holes and high-density helmet streamers, including plasma sheets extending into the low corona. Where models can reproduce coronal mass density and electron temperature, they can predict thermal emission in the extreme ultraviolet and X-ray spectrum, and line-of-sight integration yields synthetic images that show reasonably good agreement with observations, as seen in Fig. 19. When coronal models extend into the heliosphere, in particular beyond 1 AU, they offer predictions of the solar wind plasma parameters, including charge state composition that can be compared directly to in-situ observations. Early examples include those by Wu et al. (1999). Steady-state models can be validated by replicating solar wind observations for the synodic rotation period measured at the Earth, namely 27.27 days. Examples are shown in Cohen et al. (2008), Feng et al. (2012b, Meng et al. (2015) and Török et al. (2018).

Numerical mesh techniques
The computational domain for coronal simulations typically extends to r > 20R . Flows will be superfast at the outer boundary so that simple outflow boundary conditions will be well prescribed. Most simulations use spherical grids with fixed angular resolution, but highly stretched in the radial direction to provide high resolution to the lower corona. In the case of Groth et al. (2000b) and Manchester et al. (2004), an adaptive Cartesian grid was used with 4 × 4 × 4 blocks, with grid cells ranging in size from 1/32R to 2R . Grids are spatially positioned to highly resolve the corona as well as the heliospheric current sheet. More recent models have used spheric adaptive grids (e.g. Sokolov et al. 2013;van der Holst et al. 2014), while others have employed a cube- Fig. 18 Comparison between the MHD prediction (with magnetic field lines and polarization brightness shown in the first two columns) and eclipse observations (shown in the third column). Image reproduced with permission from Mikić et al. (2007), copyright by ASP sphere to maintain the advantage of nearly constant angular resolution while avoiding the singularity at the poles. For example, Feng et al. (2007Feng et al. ( , 2010Feng et al. ( , 2012a developed the Solar-InterPlanetary Conservation Element/Solution Element (SIP-CESE) MHD model that employs a six-component Yin-Yang grid system similar to a cubed sphere grid, with spherical shell-shaped domains extending from the low corona to 1 AU.

The role of Alfvén wave turbulence
The concept of Alfvén waves was introduced more than 70 years ago by Hannes Alfvén (1942). The importance of the role they play in the heliosphere was not immediately recognized due to the lack of relevant observations. Results from Mariner 2 allowed an observational study of a wave-related phenomena in solar wind (Coleman 1966(Coleman , 1967. This pioneering study culminated in a paper by Coleman (1968), a work that concluded that Alfvén wave turbulence has the potential to drive solar wind in a way that is consistent with observations at 1 AU.
Interest in the role Alfvén waves play in the heating and acceleration of the solar wind goes back to the early years of in-situ space exploration. Examples include papers by Belcher et al. (1969), Belcher and Davis Jr. (1971), and by Alazraki and Couturier (1971). A consistent and comprehensive theoretical description of Alfvén wave turbulence and its effect on the averaged plasma motion has been developed in a series of works, particularly by Dewar (1970) and by Jacques (1977Jacques ( , 1978 (see also references therein). More recent efforts to simulate solar wind acceleration utilize the approach developed by Usmanov et al. (2000). Currently, it is commonly accepted, that the gradient of the Alfvén wave pressure is the key driver for solar wind acceleration,  Usmanov et al. (2000), copyright by AGU at least in fast flows. It is important to emphasize, that while incorporating the Alfvén wave-driven acceleration is usually accomplished by including the wave pressure gradient in the governing equations (Jacques 1977), there is still no generally accepted approach to describe the coronal heating via the Alfvén wave turbulence cascade. Figure 20 shows the results of the first axisymmetric (2-D) simulation of the solar corona and solar wind using self-consistent Alfvén turbulence (Usmanov et al. 2000). After 64 h of relaxation time, a closed field region develops near the equator; the flow velocity is high in the open polar field region but decreases toward the equator above the region of the closed magnetic field.
Damping of Alfvén wave turbulence as a source of coronal heating has also been extensively studied from the early days of in situ solar wind observations (e.g., Barnes 1966Barnes , 1968. Later, it was demonstrated that reflection from sharp pressure gradients in the solar wind (Heinemann and Olbert 1980;Leroy 1980) is a critical component of Alfvén wave turbulence damping (Matthaeus et al. 1999;Dmitruk et al. 2002;Verdini and Velli 2007). For this reason, many numerical models explore the generation of reflected counter-propagating waves as the underlying cause of the turbulence energy cascade (e.g., Cranmer and Van Ballegooijen 2010), which transports the energy of turbulence from the large-scale motions across the inertial range of the turbulence spatial scale to short-wavelength perturbations. The latter can be efficiently damped due to wave-particle interaction. In this way, the turbulence energy is converted to random (thermal) energy.
Many recent efforts aim to develop models that include Alfvén waves as a primary driving agent for both heating and accelerating of the solar wind. Examples are papers by Hu et al. (2003), Suzuki and Inutsuka (2005), Verdini et al. (2010), Matsumoto and Suzuki (2012), and Lionello et al. (2014a, b).

Alfvén wave turbulence driven solar atmosphere model
The ad hoc elements can be eliminated from the solar corona model by assuming that the coronal plasma is heated by the dissipation of Alfvén wave turbulence . The dissipation itself is caused by the nonlinear interaction between oppositely propagating waves (e.g., Hollweg 1986).
Within coronal holes, there are no closed magnetic field lines, hence, there are no oppositely propagating waves. Instead, a weak reflection of the outward propagating waves locally generates sunward propagating waves as quantified by van der Holst et al. (2014). The small power in these locally generated (and almost immediately dissipated) inward propagating waves leads to a reduced turbulence dissipation rate in coronal holes, naturally resulting in the bimodal solar wind structure. Another consequence is that coronal holes look like cold black spots in the EUV and X-ray images, while closed field regions are hot and bright. Active regions, where the wave reflection is particularly strong, are the brightest in this model (see Sokolov et al. 2013;Oran et al. 2013;van der Holst et al. 2014).
The continuity, induction, and momentum equations used in the model are the following: where ρ is the mass density, u is the bulk velocity (u = |u| is assumed to be the same for the ions and electrons), B is the magnetic field, G is the gravitational constant, M is the solar mass, r is the position vector relative to the center of the Sun, μ 0 is the magnetic permeability of vacuum. As has been shown by Jacques (1977), the Alfvén waves exert an isotropic pressure ( p A in the momentum equation). The relation between the wave pressure and wave energy density is p A = (w + + w − )/2. Here, w ± are the energy densities for the turbulent waves propagating along the magnetic field vector (w + ) or in the opposite direction (w − ). The isotropic ion and electron pressures, p i and p e , are governed by the appropriate energy equations: ∂ ∂t where T e,i are the electron and ion temperatures, N e,i are the electron and ion number densities, and k B is the Boltzmann constant. Other newly introduced terms are explained below. The ideal equation of state, p e,i = N e,i k B T e,i , is used for both species. The polytropic index is γ = 5/3. The optically thin radiative energy loss rate in the lower corona is given by where Λ(T e ) is the radiative cooling curve taken from the CHIANTI v7.1 database , and references therein). The energy exchange rate between ions and electrons due to Coulomb collisions is defined in terms of the collision frequency where m e and e are the electron mass and charge, m p is the proton mass, ε 0 is the vacuum permittivity, b = B/B, and Λ C is the Coulomb logarithm. Finally, the electron heat flux q e is expressed in the collisional formulation of Spitzer and Härm (1953):

Transport and dissipation of Alfvén wave turbulence
Describing the dynamics of Alfvén wave turbulence and its interaction with the background plasma requires special consideration. The evolution of the Alfvén wave amplitude (velocity, δu, and magnetic field, δB) is usually treated in terms of the Elsässer (1950) variables, z ± = δu ∓ δB/ √ μ 0 ρ. The Wentzel-Kramers-Brillouin (WKB) approximation (Wentzel 1926;Kramers 1926;Brillouin 1926) is used to derive the equations that govern transport of Alfvén waves, which may be reformulated in terms of the wave energy densities, w ± = ρz 2 ± /4. Dissipation of Alfvén waves, Γ ± w ± , is the physical process that drives the solar wind and heats the coronal plasma.
Alfvén wave dissipation occurs when two counter-propagating waves interact. Alfvén wave reflection from steep density gradients is the physical process that results in local wave reflection, thus maintaining a source of both types of waves. In order to describe this wave reflection we need to go beyond the WKB approximation that assumes that the wavelength is much smaller than spatial scales of the background variations.
The equation describing the propagation, dissipation, and reflection of Alfvén turbulence has been derived in van der Holst et al. (2014): where V A = B/ √ μ 0 ρ is the Alfvén velocity, while the dissipation rate (Γ ± ) and the reflection coefficient (R) are given by and Here L ⊥ is the transverse correlation length of Alfvén waves in the plane perpendicular to the magnetic field line and I max = 2 is the maximum degree of turbulence "imbalance." If √ w ± /w ∓ < I max , then Alfvén wave reflection is neglected and R = 0. With the help of the dissipation rate of Alfvén turbulence one can express the ion and electron heating rates: where f p ≈ 0.6 is the fraction of Alfvén wave energy dissipated to the ions. Finally, to close the system of equations, we use the following boundary condition for the Poynting flux of Alfvén waves, Π A : The transverse correlation length is assumed to scale with the magnetic field magnitude (e.g., Hollweg 1986):

Chromosphere boundary conditions
A simulation model based on the Alfvén wave turbulence may be extrapolated down to the top of the chromosphere. In order to save computational resources for this physicsbased model (which would require a gigantic amount of resources), the temperature and density at the top of the chromosphere are specified as: One needs to use an innovative approach to handle the sharp density gradients that have a spatial scale length of which would greatly complicate the Alfvén wave turbulence model and introduce unmanageable wave reflection. To avoid this problem, we apply WKB Alfvén wave turbulence effects and let the Alfvén waves freely propagate through the plasma at T ≤ T ch . To both balance the radiative cooling and ensure the hydrostatic equilibrium, we apply an exponential heating function, Q h = A exp(−x/L), to maintain the analytical solution of the momentum and heat transfer equations, as follows: Here g = 274 m/s 2 is the gravity acceleration near the solar surface, the direction of this acceleration being antiparallel to the x-axis, and m i is the proton mass. The two constants in the solution, N ch and T ch , which are the boundary values for the density and temperature, respectively, are unambiguously related to the amplitude, A, of the heating function: and to the scale-length (see Eq. 38). Notice that there is a very simple relationship for the exponential scale-length for the heating function, which is half of the barometric scale-length of density variation: 2L = L g = k B (T e + T i )/(m i g). The solution satisfies the equation for the heat conduction as long as the heat transfer in the isothermic solution is absent and heating at each point exactly balances the radiation cooling. The hydrostatic equilibrium is also maintained, as long as The suggested solution does a good job describing the chromosphere. The short scale-length of the heating function, (see Eq. 38), which is equal to ≈ 0.6 Mm for T ch = 2 × 10 4 , may presumably mimic absorption of (magneto)acoustic turbulent waves, rapidly damping due to the wave-breaking effects. Physically, including this chromosphere heating function would imply that the temperature in the chromosphere is elevated compared to the photospheric temperatures due to some mechanism acting in the chromosphere itself. By no means can this energy be transported from the solar corona as long as the electron heat conduction rate at chromospheric temperatures is very low.

Transition region
One of the first successful models describing the Transition Region (TR), based on an analytical solution, was published by Lionello et al. (2001) (see also Lionello et al. 2009;Downs et al. 2010). This model treats the TR as a thin continuous layer and it does not agree well with observations. This discrepancy suggests that the TR could be more accurately described as a carpet of 1D-like "threads" (maybe spicules, see Cranmer et al. 2013). Collectively these threads give the impression of a "thin" layer described by a solely height-dependent 1D solution. To derive this solution one uses 1D governing equations to close the MHD model with the boundary condition at "low boundary", which is at the same time the top boundary for the TR model. By solving the said 1D equations, one can merge the MHD model to the chromosphere, which is the bottom of the TR.
The heat transfer equation for a steady state hydrogen plasma in a uniform magnetic field reads: Here Q h = Γ (w − + w + ) is the coronal heating function, assumed to be constant at spatial scales typical for the TR. Note that the coordinate is taken along the magnetic field line, not along the radial direction. By multiplying Eq. (42) by κ 0 T 5/2 e (∂ T e /∂s), and by integrating from the interface to the chromosphere at temperature T e , one can obtain: Here the product, N e T e , is assumed to be constant. Therefore, it is separated from the integrand. For a given T ch , the only parameter in the solution is (N e T e ). It can be expressed at any point in terms of the local value of the heating flux and the radiation loss integral: The assumption of constant (N e T e ) is fulfilled only if the effect of gravity is negligible. Quantitatively, this condition is not trivial, as long as both the barometric scale and especially the heat conduction scale are functions of temperature. The barometric scale may be approximated as L g (T e ) ≈ T e × (60 m/K). The heat conduction scale, L h , can be estimated by noticing that within a large part of the transition region the radiation losses dominate over the heating function, so they are balanced by heat conduction: κ 0 T 5/2 e × (T e/L 2 h ) ∼ Q r . Thus, the condition for neglecting gravity is: In Fig. 21 we plot temperature dependencies L h (T e ) and L g (T e ) for (N e T e ) = 10 20 K/m 3 . We see that near the chromosphere boundary the approximation given by Eq. (45) works very well as long as the temperature changes with height are very abrupt. The increase in temperature to 10 5 K occurs in less than 0.1 Mm. This estimate agrees with the temperature profile seen in the chromospheric lines (see, e.g., Figs. 2 and 8 in Avrett and Loeser 2008). However, as the temperature increases with height, the effect of gravity on the temperature and density profiles becomes more significant. It becomes comparable to the heat conduction effect at T e ≈ 4.5 × 10 5 K, which can be accepted as the coronal base temperature, so that the transition region corresponds to the temperature range from T ch ≈ 2 × 10 4 K to 4.5 × 10 5 K, with a typical width of ∼ 10 Mm ≈ R S /70. The TR solution merges to the chromosphere solution with no jump in pressure. The merging point in the chromosphere, therefore, is at the density of (N e T e )/T ch ∼ 10 16 cm −3 . The short heat conduction scale at the chromosphere temperature (see Fig. 21) ensures that the heat flux from the solar corona across the transition region does not penetrate to higher densities.
In Sect. 7 we revisit the transition region analytical model. However, in Sect. 6 we use another way to model the transition region, by artificially increasing the heat conduction in the lower temperature range (see Abbett 2007). Consider the transformation of the temperature functions shown in Eqs. (42)- (43): with a common factor, f ≥ 1. The equations do not change in this transformation and the only effect on the solution is that the temperature profile in the transition region becomes a factor of f wider. By applying the factor, f = (T m /T e ) 5/2 at T ch ≤ T e ≤ T m , the heat conduction scale in this range is almost constant and is close to ≈ 1 Mm for a choice of T m ≈ 2.5 × 10 5 K (see Fig. 21). It should be emphasized, however, that the choice of temperature range for this transformation is highly confined by the condition given in Eq. (45). If a higher value of T m is chosen, the heat conduction scale at the chromospheric temperature exceeds the barometric scale in the chromosphere, resulting in a unphysical penetration of the coronal heat into the deeper chromosphere. The global model of the solar corona, with this unphysical energy sink, suffers from reduced values for the coronal temperature and produces a visible distortion in the EUV and X-ray synthetic images. Thus, in formulating the transition region model we modify the heat conduction, the radiation loss rate, and the wave dissipation rate, and the maximum temperature for this modification does not exceed T m ≈ 2.5 × 10 5 K.

Alfvén wave-driven MHD coronal models
With the help of the mathematical formalism described in Sect. 5, Alfvén wavedriven self-consistent models of the solar atmosphere and corona can be developed. The potential role of low-frequency Alfvén waves to provide heat and momentum to accelerate the solar wind has been already recognized in the first decade of in-situ exploration of the interplanetary medium. Alfvén waves have long been measured in situ in the solar wind (Belcher and Davis Jr. 1971), and have more recently been remotely observed in the solar corona (De Pontieu et al. 2007;Cranmer et al. 2009), where their energy is sufficient to heat and accelerate the solar wind. The theoretical exploration of Alfvén waves was first suggested in early work by Hollweg (1978Hollweg ( , 1981 and Hollweg et al. (1982). Based on this early work, theories were developed that describe the evolution and transport of Alfvénic turbulence, e.g., Zank et al. (1996Zank et al. ( , 2017, Matthaeus et al. (1999) and Zank (2014). To self-consistently describe the heating and acceleration of the solar wind with Alfvénic turbulence, several extended magnetohydrodynamic (MHD) models have been developed. One-dimensional models include those by e.g., Tu and Marsch (1997), Laitinen et al. (2003), Vainio et al. (2003), Suzuki (2006), Cranmer and Van Ballegooijen (2010) and Adhikari et al. (2016), while multi-dimensional models include e.g., Usmanov et al. (2000), Suzuki and Inutsuka (2005), Cranmer et al. (2009 and Lionello et al. (2014a).
These models have many common features. First, they employ low-frequency Alfvén waves, which are assumed to dissipate below the ion cyclotron frequency. Wave amplitudes are typically prescribed at the inner boundaries to match observed wave motions in the low corona (De Pontieu et al. 2007). Wave energy propagates at the Alfvén speed along the magnetic field and drives the corona in two ways: (i) wave pressure gradient provides a volumetric force that accelerates the solar wind, while (ii) wave dissipation heats the plasma. In a variety of implementations, these Alfvén wave-driven models have been shown to self-consistently reproduce the fast/slow solar wind speed distribution (e.g., Usmanov et al. 2000;van der Holst et al. , 2014. An alternative method of coronal heating was developed by Suzuki (2002Suzuki ( , 2004. Here, the focus is on shock wave heating instead of turbulent dissipation. The assumption is that slow and fast magneto-acoustic waves are generated by small scale reconnection events. These wave steepen into shocks while propagating along the field lines into the corona to heat and accelerate the plasma. Suzuki (2004) demonstrated that the dissipation of shock trains can satisfactory reproduce the fast and slow wind speeds, except for the observed high temperatures in the slow wind, where other heating mechanisms might be needed. To the best of our knowledge, there are no 3D simulations that self-consistently include shocks and turbulence to assess which mechanism dominates in the heating and acceleration.
The recently developed Alfvén Wave Solar Model (AWSoM) van der Holst et al. 2014;Meng et al. 2015) extends the description with a threedimensional solar corona/solar wind model that self-consistently incorporates lowfrequency Alfvén wave turbulence. The model employs a phenomenological treatment of wave dissipation, with a prescribed correlation length inversely proportional to the magnetic field strength. In this case, the wave spectrum is not resolved, so only the total forward and backward propagating wave energy densities and the partitioning of dissipated wave energy between electrons and protons is fixed. Turbulence parameters are the wave energy densities, the correlation length, and the reflection rate. The wave reflection model used is essentially the same formulation as introduced by Matthaeus et al. (1999). In this case, the energy of the dominant wave is transferred to a counterpropagating minor wave, with the reflection coefficient controlled by the gradient of the Alfvén speed.
With AWSoM, Alfvén waves are represented as two discrete populations propagating parallel and antiparallel to the magnetic field, which are imposed at the inner boundary with a Poynting flux of the outbound Alfvén waves assumed to be proportional to the magnetic field strength. The wave spectrum is not resolved, so the waves are presented as frequency-integrated wave energies that propagate parallel to the magnetic field at the local Alfvén speed. The waves possess a pressure that does work and drives the expansion of the plasma. In this model, outward propagating waves experience partial reflection on field-aligned Alfvén speed gradients and the vorticity of the background. The partial reflection leads to nonlinear interaction between oppositely propagating Alfvén waves and results in an energy cascade from the large outer scale through the inertial range to the smaller perpendicular gyroradius scales, where the dissipation takes place. The apportioning of the dissipated wave energy to the isotropic electron temperature and the parallel and perpendicular proton temperatures depends on criteria for the particular kinetic instabilities that are involved (Meng et al. 2015). To apportion heating to the various ion species, we use the multispecies generalization of the stochastic heating, as described by Chandran et al. (2013).
In the AWSoM model, the partitioning strategy is based on the dissipation of kinetic Alfvén waves (KAWs) with the stochastic heating mechanism for the perpendicular proton temperature (Chandran et al. 2011). In this mechanism, the electric field fluctuations due to perpendicular turbulent cascade can disturb the proton gyro motion enough to give rise to perpendicular stochastic heating, assuming that the velocity perturbation at the proton gyroradius scale is large enough. The firehose, mirror, and ion-cyclotron instabilities, due to the developing proton temperature anisotropy, are accounted for. When the plasma is unstable because of these instabilities, the parallel and perpendicular temperatures are relaxed back toward marginal stable temperatures, with the relaxation time inversely proportional to the growth rate of these instabilities. See the work of Meng et al. (2012Meng et al. ( , 2015 for a detailed description. In this global model, excess of energy in the lower corona is transported back to the upper chromosphere via electron heat conduction where it is lost via radiative cooling. The AWSoM model is representative of the state of the art of extended MHD Alfvén wave-driven coronal models, presenting many significant advances in modeling capability. First, turbulent dissipation rates are based directly on counter propagating wave amplitudes, which are greatly enhanced by wave reflection at Alfvén speed gradients. Second, the model captures temperature anisotropies caused by preferential perpendicular heating in the fast solar wind. Third, the effects of kinetic instabilities: fire hose, mirror mode, and cyclotron instabilities limit temperature anisotropies with thresholds that are dependent on the proton temperature ratio and plasma β. Finally the three-dimensional model includes the entire structure of the corona including active regions and slow and fast streams. This is the first time such kinetic physics has been incorporated into a global numerical model of a CME propagating through the solar corona, which allows us to address both particle heating, Alfvén wave damping, and their nonlinear coupled interaction as shown in Manchester and Van Der Holst (2017).

Multi-temperature coronal models
Magnetohydrodynamic (MHD) theory is the simplest self-consistent model describing the macroscopic structure of the corona comprising the global distribution and temperature of the coronal plasma and magnetic field. Such MHD models ignore the extreme complexity of a coronal environment that is made up of many plasma species affected by a wide range of wave-particle and particle-particle interactions, where heating occurs by the dissipation of waves and time varying electric currents. Particle populations are far from equilibrium and exhibit vastly different temperatures and distribution functions with extended high energy tails, the full complexity of which can only be described by kinetic models with non-Maxwellian velocity distribution functions (Landi and Pantellini 2003). For electrons, there are two nearly isotropic populations: the thermal core and the suprathermal halo, and a field aligned strahl component (Rosenbauer et al. 1977) that travels away from the Sun. Ions are more often characterized by a population that is anisotropic with a temperature perpendicular to the magnetic field higher than that parallel to the field. Hydrogen is fully ionized, and all other atomic species are highly ionized. Protons, being almost 2000 times more massive than electrons, thermodynamically decouple at a distance of roughly 1.5 solar radii where Coulomb collisions become infrequent.
To begin to address a range of physical processes as well as reduce the number of free parameters and ad hoc assumptions, a new generation of extended MHD global coronal models were developed. First and foremost, thermodynamic processes were added, beginning with heat conduction and radiative losses, which allowed models to accurately capture the temperature structure of the lower atmosphere. The use of stretched radial grids allow these models to resolve the transition region so that they may extend down to the upper chromosphere (Lionello et al. 2009;Downs et al. 2010;Sokolov et al. 2013;van der Holst et al. 2014). The radiative looses for these models are almost universally based on CHIANTI tables (Dere et al. 1997), which specify the optically thin losses from the corona, which is dominated by line emission from heavy ions impacted by thermal electrons.
With the thermal processes captured, extended MHD simulations can successfully reproduce images of the low corona provided by extreme ultraviolet imaging telescopes, including SOHO/EIT, STEREO/EUVI, and SDO/AIA. Observations provided by Downs et al. (2010) and van der Holst et al. (2014) indicate the qualitative match to coronal temperature and density available with the new models. Figure 19 provides an example of coronal ultraviolet images from Jin et al. (2017). Here, the simulated active regions for 7 March 2011 (CR2107) naturally produce the enhanced emissions observed by SDO/AIA 211 Å, STEREO A/EUVI 171 Å, and STEREO B/EUVI 195 Å. The distinct feature in the present model is the enhanced wave reflection in the presence of strong magnetic fields, such as in close proximity to active regions that can increase the dissipation and thereby intensify the observable EUV emission.
Even more complex, nonequilibrium thermodynamics can be captured with multiple-temperature single-fluid coronal models. Two-temperature coronal models describe protons and electrons with a single fluid velocity but with individual energy equations and temperatures (e.g., Sturrock and Hartle 1966;). The impetus for the feature stems from two facts: First, protons are almost 2000 times more massive than electrons, so that at one million degrees, their respective sound speeds are 120 and 5000 km/s. Second, Coulomb collisions are so infrequent that within a fraction of a solar radius above the surface the ions and electrons thermally decouple from each other. Consequently, heat conduction in the corona is completely dominated by electrons, which is particularly conspicuous in CME-driven shocks.
The speed of fast CMEs occurs in the range where protons are shocked but not the electrons, and beyond a distance of two solar radii (R ) collisions become so infrequent that protons and electrons thermally decouple on the timescale of the shock passage. As a result, protons can be shock heated to high temperature, while in the same location electrons cool from adiabatic expansion and heat conduction. These temperature structures in CMEs were first modeled with one-dimensional two-temperature simulations by Kosovichev and Stepanova (1991) and Stepanova and Kosovichev (2000), and later in three-dimensional simulations by Manchester et al. (2012) and Jin et al. (2013).
The three-temperature thermodynamics model captures the electron temperature and resolves proton temperature into components parallel and perpendicular to the magnetic field. Such models can capture the temperature anisotropy produced by a nearly collisionless plasma heated by wave-particle interaction. A leading example is the three-temperature version of AWSoM described in van der Holst et al. (2014) and Meng et al. (2015) and shown in Fig. 22. Here, particle temperatures and heating rates are shown on the left and right, respectively, and are determined by the incorporation of a description of turbulent dissipation developed by Chandran et al. (2013). This theory describes the turbulent cascade and dissipation of kinetic Alfvén waves, providing the thermal energy partitioning between protons and electrons. As seen in Fig. 22, for regions of low plasma beta, such as coronal holes, most energy goes to perpendicular proton heating, while in high beta regions, such as the current sheet, parallel heating dominates. Electron heating dominates at intermediate beta levels found at the margins of the current sheet. In AWSoM, temperature anisotropies are limited by kinetic instabilities, which are invoked when temperature ratios surpass the instability thresholds of fire hose, mirror mode, and cyclotron kinetic instabilities. This three-temperature model has also been applied to study the thermodynamics and the interaction Alfvén turbulence with CMEs and CME-driven shocks  7 Threaded field line model

Magnetic "treads"
In the transition region the plasma temperature increases some two orders of magnitude over ∼ 10 2 km, resulting in a temperature gradient of ∼ 10 4 K/km. To resolve this gradient 3-D numerical simulations require sub-kilometer grid spacing, making these simulations computationally very expensive.
An alternative approach is to reformulate the mathematical problem in the region between the chromosphere and the corona in a way that decreases the computational cost. Instead of solving a computationally expensive 3-D problem on a very fine grid, one can reformulate it in terms of a multitude of much simpler 1-D problems along threads that allows us to map the boundary conditions from the the solar surface to the corona. This approach is called the Threaded-Field-Line Model (TFLM) (Sokolov et al. 2016).
The physics behind the reformulated problem is the assumption that between the solar surface and the top of the transition region (R < r < R b ) the magnetic field can be described with a scalar potential. A thread represents a field line of this potential field. One can introduce a 1-D problem that describes a magnetic flux tube around a given thread (Sokolov et al. 2016).
The magnetic field is divergenceless, therefore the magnetic flux remains constant along each thread: where s is the distance along the field line, and A(s) is the cross-section area of the flux tube. Other conservation laws are also greatly simplified due to the fact that in a low-beta plasma, the flow velocity is aligned with the magnetic field. Assuming steady-state, the basic conservation laws can be written as 1-D equations. Continuity equation: Conservation of momentum: here p = p i + p e , 2T = T i + T e , R TR is the radius of the bottom of the transition region (TR), and b r is the radial component of b. In this expression terms proportional to u 2 are neglected, j × B is omitted due to the fact that electric currents vanish in a potential field (j ∝ ∇ × B = 0), and the pressure of Alfvén wave turbulence is assumed to be much smaller than the thermal pressure, p A p. Conservation of energy: where the term ∂ T /∂t is retained because it is assumed that the electron heat conduction is a relatively slow process. In addition to the plasma equations, the Alfvén wave dynamics can also be reformulated. In Eq. (31), we introduce a new variable, a 2 ± : With the help of this substitution, the Alfvén wave transport equation becomes ∂a 2 ± ∂t + ∇ · ua 2 ± ± (V A · ∇) a 2 ± = ∓Ra − a + − 2 These equations can be additionally simplified since in the lower corona u V A (i.e., waves are assumed to travel fast and quickly converge to equilibrium), therefore we can neglect the ∂a 2 ± /∂t terms: Additionally, we introduce a new variable: Now the wave equations become Equation (55) describe boundary value problems and one needs to specify boundary conditions somewhere along the thread. Let ξ − denote the location of the lower boundary at the outgoing end of the thread (where the field direction points away from the Sun), and ξ + denote the lower boundary at the downward end of the thread. The boundary conditions now must specify the values of a ± at the location where the Alfvén turbulence enters into the thread: a + (ξ = ξ − ) = a + 0 and a − (ξ = ξ + ) = a − 0 . The values of a ± 0 are empirically specified.

From the transition region to the threaded field line corona
Finally, one must specify the plasma and turbulence conditions at the interface between the threaded field line region and the corona at the radial distance of r = R b . These conditions depend on the direction of the magnetic field at the interface.
If the magnetic field points outward, b r (r = R b ) > 0: If the magnetic field points inward, b r (r = R b ) < 0: ; (a + ) T F = (a + ) cor ; (a − ) cor = (a − ) T F . (57) In addition, the temperature and density need to be matched at the interface between the threaded field line and the corona. In order to achieve this it is assumed that at the interface the temperature gradient is mainly in the radial direction: The boundary condition for the density is controlled by the sign of b · u: In the last step, one needs to consider the energy balance in the transition region where two physical processes balance each other: heat conduction and radiative cooling. Assuming steady-state conditions, this energy balance can be expressed using Eq. (42) for Q h = 0, ∂ ∂s κ 0 T 5/2 ∂ T ∂s = N e N i Λ(T ), where for the field-aligned heat conduction coefficient the usual κ = κ 0 T 5/2 expression is used. The length of a given magnetic field line between the photosphere and the top of the transition region is obtained by integrating the thread: If the temperature at the top of the transition region, T TR , is known, one can obtain the heat flux and pressure from the following equations: where T ch ≈ (1 ÷ 2) × 10 4 K and Here U heat (T ) = 2 k 2 B T T ch κ 0 (T ) 1/2 Λ(T )dT .
The Λ(T ) and U heat (T ) functions can be easily tabulated using the CHIANTI database (Dere et al. 1997;Landi et al. 2013).
wind from tens to hundreds of R (cf. Pizzo 1991;Odstrčil and Pizzo 1999;Groth et al. 1999a;Linde et al. 1998;Pogorelov et al. 2013;Opher et al. 2003). However, modeling the transition from the dense, cold photosphere to the hot super-Alfvénic corona is still challenging. The first generation of coronal models used simplified energetics (cf. Usmanov 1993;Mikić et al. 1999;Cohen et al. 2007) with considerable success. "Thermodynamic" models use a realistic adiabatic index but empirical heating functions to heat and accelerate the solar wind (cf. Groth et al. 1999b;Lionello et al. 2001). The advantage of the thermodynamic approach is their ability to describe shock related phenomena. The latest generation of coronal models use Alfvén waves to heat and accelerate the solar wind (cf. Usmanov et al. 2000;van der Holst et al. 2014). While this approach has the promise to explain the origin of both the fast and slow solar wind states, it is still at a relatively early state of development and much work remains to be done. We urge the reader to stay tuned.
Though model validation goes beyond the scope of this paper, one can envision investigations to determine the veracity of the models described in this review. As it was pointed out, one can compare emission properties in the EUV and soft X-ray with observed signatures in the low corona. In the extended corona, many of Alfvénic effects may be more pronounced, such as non-thermal velocities and signatures of wave dissipation such as temperature anisotropies. The Parker Solar Probe is ideally designed to resolve these questions with high cadence observations of electromagnetic waves and particle distribution functions. One can also validate models by comparing their predictions of solar wind parameters with those from the WSA model. In addition, efforts are now underway to develop an integrated model describing the acceleration and transport of solar energetic particles (SEPs) directly coupled with an Alfwén wind turbulence based solar wind model (see Borovikov et al. 2018). Such a coupled model addresses the effects of wave turbulence on SEP transport. In particular, particle diffusion rates, and arrival times, are affected by the turbulence level in the solar corona. Therefore, SEP observations may provide additional validation capability for the coronal model.