A novel code for numerical 3-D MHD studies of CME expansion

A recent third-order, essentially non-oscillatory central scheme to advance the equations of single-fluid magnetohydrodynamics (MHD) in time has been implemented into a new numerical code. This code operates on a 3-D Cartesian, non-staggered grid, and is able to handle shock-like gradients without producing spurious oscillations. To demonstrate the suitability of our code for the simulation of coronal mass ejections (CMEs) and similar heliospheric transients, we present selected results from test cases and perform studies of the solar wind expansion during phases of minimum solar activity. We can demonstrate convergence of the system into a stable Parker-like steady state for both hydrodynamic and MHD winds. The model is subsequently applied to expansion studies of CME-like plasma bubbles, and their evolution is monitored until a stationary state similar to the initial one is achieved. In spite of the model's (current) simplicity, we can confirm the CME's nearly self-similar evolution close to the Sun, thus highlighting the importance of detailed modelling especially at small heliospheric radii. Additionally, alternative methods to implement boundary conditions at the coronal base, as well as strategies to ensure a solenoidal magnetic field, are discussed and evaluated.


Introduction
Coronal mass ejections (CMEs) moved into the focus of several research activities during recent years. Besides a variety of observational data resulting from SOHO (Pick Correspondence to: J. Kleimann (Jens.Kleimann@uit.no) et al., 2006), SMEI (Webb et al., 2006), and, very recently, STEREO (Vourlidas et al., 2007), significant progress has also been achieved with the numerical modelling of CMEs, see, e.g., the reviews by Aschwanden et al. (2006) and Forbes et al. (2006). The motivation for the various activities is at least fourfold. First, CMEs are amongst the main mediators of the influence of the Sun on the inner heliosphere, particularly on the Earth and its environment, where they significantly co-determine the space weather conditions. There is, in view of ever advancing technology that is increasingly sensitive -if not vulnerable -to space weather effects, strong interest in an understanding of the latter. An even stronger driver of research activities is provided with the recognition that space weather phenomena offer valuable opportunities to study many aspects of plasma astrophysics in great detail (Scherer et al., 2005;Bothmer and Daglis, 2006;. Second, as a consequence of the shocks driven by CMEs, they serve as particle accelerators that do not only contribute to space weather effects, but can be used to study the actual acceleration processes (Reames, 1999;Mewaldt et al., 2005;Li et al., 2005), which are expected to occur in other astrophysical systems as well (Eichler, 2006). Third, with the recent launch of the two-spacecraft mission STEREO (Kaiser, 2005), the full three-dimensional structure of CMEs can be observed both remotely and in-situ for the first time. First results have already been reported by, e.g., Howard and Tappin (2008) and Vourlidas et al. (2007). And, fourth, the magnetohydrodynamic (MHD) modelling of CMEs provides an excellent testbed for numerical codes. Although not strongly motivated by CME physics, this is actually one of the main drivers of model development, as is manifest with numerous approaches documented in the literature. These various approaches can be ordered into three groups. There is (i) principal modelling that is either analytical and/or based on symmetry assumptions (e.g., Titov and Démoulin, 1999;Roussev et al., 2003;Schmidt and Cargill, 2003;Jacobs et al., 2005), (ii) local modelling limited to the arXiv:1204.3767v1 [astro-ph.SR] 17 Apr 2012 extended corona, i.e. a few tens of solar radii (e.g., Mikić and Linker, 1994), and (iii) global modelling covering the inner heliosphere from the solar surface out to 1 AU and beyond (e.g., Manchester et al., 2004;Odstrčil et al., 2005;Tóth et al., 2005;Riley et al., 2006). Despite these intensified efforts and activity regarding the study of CMEs, there remains both a number of unsolved problems and various modelling deficiencies. For example, the acceleration and heating processes of the plasma near the coronal base are -even nearly 50 years after the 'discovery' of the solar wind -still not known (Cranmer et al., 2007), and there is also no agreement on the processes that actually initiate CMEs . Also, their propagation and evolution in size and shape is by far not fully understood in all detail, and neither is their interaction with the background solar wind (Jacobs et al., 2007), with other CMEs (Gopalswamy et al., 2001), and with planetary magnetospheres (e.g., Groth et al., 2000;Ip and Kopp, 2002). Regarding the model formulations underlying the numerical simulation of CMEs, particularly the (non-thermal) heating of the plasma is mostly treated in a rather simplified manner via ad-hoc heating functions (e.g., Groth et al., 2000;Manchester et al., 2004), variable adiabatic indices γ = γ(r) (e.g., Lugaz et al., 2007), or phenomenological heating functions (e.g., Usmanov et al., 2000), for a discussion see Fichtner et al. (2008). With the intention to address several of the above-mentioned problems, we have applied our recently developed CWENObased MHD code (Kleimann et al., 2004), which primordially originated from that by Grauer and Marliani (2000), to the CME expansion problem. To our knowledge, this is the first published paper describing the application of a CWENO-based numerical code to an MHD problem related to space physics. In the following, we describe the model (Sect. 2) and its numerical realization (Sect. 3 to 5), present results of analyses of both the propagation of individual and the interaction of two CMEs, and suggest a possible connection of our findings to observations (Sect. 6).

Governing equations
Choosing the normalization constants summarized in Tab. 1, the set of MHD equations for mass density ρ, flow velocity u, magnetic field strength B, and gas pressure p reads (in dimensionless form): where g = −Γ/r 2 e r and (5) respectively denote gravity (with Γ := (GM )/(R c 2 s ) = 11.49 in normalized units, cf. Tab. 1), and the total energy density of a plasma with adiabatic exponent γ. Throughout this paper, · is used to denote the norm of a vector (i.e. X ≡ √ X · X for any vector X), and the symbolÎ in Eq.
Quantity Normalization (solar) mass M = 2.0 × 10 30 kg length L0 := R = 7.0 × 10 8 m temperature T T0 := 1.0 × 10 6 K number density n ≡ ρ/mp n0 := 1.0 × 10 14 m −3 plasma pressure p 2 n0 kB T0 = 2.8 × 10 −3 Pa velocity u cs := 2 kB T0/mp = 1.3 × 10 5 m s −1 time t L0/cs = 5.5 × 10 3 s energy density e mp n0 (cs) 2 = 2.8 × 10 −3 J m −3 mag. induction B cs √ µ0 mp n0 = 4.2 × 10 −5 T heating rate Q (cs) 3 /L0 = 3.0 × 10 6 W kg −1 (CME) mass Mcme mp n0 (L0) 3 = 5.7 × 10 13 kg A Parker-like solution for the solar wind is per construction isothermal, i. e. γ = 1, resulting in an adiabatic cooling for γ > 1. In reality, the decrease in temperature T ≡ p/ρ due to this adiabatic cooling of the expanding plasma is compensated by processes such as reconnective energy release and Alfvénic wave heating. A realistic inclusion of such effects, while certainly desirable, is beyond the scope of this first approach, and, thus, reserved for future refinements of our model. As an alternative, we employ an ad hoc heating function with a prescribed heating profile α(r) and a target temperature T c . To derive the form of Eq. (8), we first seek the heating function Q iso which maintains a constant temperature T c everywhere, irrespectively of γ. This is done by inserting the corresponding isothermal equation of state into the MHD equations (1-4) and solving them analytically for Q, which yields Therefore, local heating (or cooling) at any heliocentric radius r h is conveniently achieved by choosing α(r h ) > 1 (or α(r h ) < 1) in Eq. (8). Test runs demonstrating the validity of this method have been carried out by Kleimann (2005).
3 Numerical implementation

Algorithm
In order to integrate Eqs.
(1-4) forward in time, we employ a 3-D variant (Kleimann et al., 2004) of a recent semi-discrete Central weighted essentially non-oscillatory (CWENO) scheme by Kurganov and Levy (2000) with thirdorder Runge-Kutta time stepping. Notable advantages of CWENO include its third-order accuracy in smooth regions (which automatically becomes second order near strong gradients to minimize spurious oscillations) and an easy generalization to multi-dimensional systems of equations due to the fact that no (exact or approximate) Riemann solver is needed. The CWENO scheme thus allows simultaneously to achieve high shock resolution comparable with the best shock capturing schemes and high order convergence in smooth regions dominated by plasma waves. Although this scheme is not strictly total variation diminishing (TVD), simulations by Levy et al. (2000) do indicate an upper bound for the total variation of their solutions. Moreover, Havlík and Liska (2006) use a set of astrophysically relevant test cases to compare the performance of several methods for ideal MHD and stress CWENO's superior accuracy.

Example test case: Alfvén wings
Various elementary tests of our implementation have been completed successfully (Kleimann et al., 2004), such as advection in one and two dimensions (also for propagation directions inclined at angles 0 < ν < π/2 to the coordinate surfaces), shock tubes with and without magnetic field inclusion, etc. While those standard tests will not be reproduced here, one rather advanced test setting, which is also of astrophysical relevance involving so-called "Alfvén wings" is worth being mentioned. While the the finite extent of the wave-generating obstacle does not allow for an exact analytical solution, the usefulness of this simple but meaningful test problem stems from the fact that it incorporates several types of MHD waves (Alfvén and slow/ fast magnetosonic), the expansion speed and characteristics of which can be verified quantitatively with theoretical expectations to ensure proper implementation of the relevant physics. As shown by Drell et al. (1965), the movement of a conductive obstacle (e.g. a satellite or small planet) through a homogeneous fluid with a perpendicular magnetic field will generate standing MHD waves in the (u, B) plane called Alfvén wings. This phenomenon plays a major role for the interaction of the moon Io with the Jovian magnetosphere (Neubauer, 1980;Linker et al., 1988), and also for artificial satellites in the magnetosphere of Earth, see e. g. Kopp and Schröer (1998) and references therein. The corresponding numerical test case, which works well both in 2-D and 3-D, involves an initially homogeneous flow u = u 0 e x of constant density, which is combined with a perpendicular, equally homogeneous magnetic field B = B 0 e z . A solid, spherical obstacle is then implemented by performing an artificial deceleration after each time step, such that for t ≥ 1, the flow will vanish within r ≤ 1. Figures 1 and 2 illustrate the emanating wing structure. Direction and speed of propagation agree well with their respective theoretical expectations.

Choice of coordinates
At first sight, the Sun's obviously spherical shape would suggest the use of spherical coordinates [r, ϑ, ϕ], especially since the radial convergence of lines of constant ϑ, ϕ entails the additional benefit of increased spatial resolution near the Sun's surface. On the other hand, the Courant-Friedrichs-Lewy (CFL) criterion of numerical validity and stability (Courant et al., 1928), which requires the "velocity" ∆x/∆t to be greater than the maximum physical propagation velocity, imposes a limit on the time step ∆t based on the cell size ∆x.
The very choice of a coordinate system with varying grid cell sizes, together with the requirement that the time step be uniform on the entire grid, thus implies that ∆t will be set by the ∆x of the grid's smallest cell. For spherical coordinates, this means that the increased resolution at small r, however welcomed for physical reasons, would force ∆t to be much lower than what the CFL criterion would require for most parts of the computational domain. This 'problem of small time steps' is avoided by using Cartesian coordinates, which have equal cell spacing everywhere and thus do not waste computing time on the larger cells. Even worse is the problem of coordinate singularities at the poles ϑ ∈ {0, π}, which require delicate numerical treatment. For these reasons, we opt for Cartesian coordinates [x, y, z], for which numerics are faster, simpler (esp. with respect to multi-dimensional extension), and more stable. This is especially true since our CWENO code is built within a framework that allows for Cartesian Adaptive Mesh Refinement (AMR, see (Kleimann et al., 2004)). This is of high interest for more detailed studies of, e.g., the inner structure of a CME. While AMR can, in principle, be used with spherical coordinates, its advantage is over-compensated by the fact that the convergence of grid spacing implies unacceptably low CFL numbers. (Note also that since CMEs generally do not exhibit any clear spatial symmetry, the use of non-Cartesian coordinates is not expected to entail any particular advantages for their description.)

Divergence cleaning
Like many other algorithms, CWENO does not exactly conserve the solenoidality condition ∇ · B = 0 for the magnetic field, and a correction scheme becomes mandatory to avoid unphysical artifacts. From the wealth of existing schemes (for an overview see, e. g., Tóth (2000)), we have evaluated the performance of the Generalized Lagrange Multiplier (GLM) approach by Dedner et al. (2002) against a classical projection scheme (see Sect. 3.4.2).

The GLM scheme
The GLM scheme solves an additional equation for a position-and time-dependent Lagrange multiplier Ψ, and adds a term −∇Ψ to the right hand side of Eq.
(3). This procedure causes Ψ (and hence ∇ · B) to be damped with decay constant τ d := λ/v f , while at the same time advection of Ψ towards the boundary of the computational volume occurs at the highest permissible speed v f (chosen to equal the global maximum of the fast magnetosonic speed in this case). Following Dedner et al. (2002), a value of 0.18 is used for the second constant λ.
The main advantage of this method is that Eq. (12) already possesses the correct conservative form, allowing for direct treatment with CWENO. In particular, physical conservation laws are not affected in any way. Figure 3 compares the performance of the two methods for a standard run. The obviously inferior performance of GLM can be explained by the fact that within a spherical layer L around the inner (solar) boundary, the boundary procedure described in section 5.1 entails an averaging of the inner boundary value B in and the newly computed outer solution B out via for some function f : L → [0, 1], which is bound to introduce a marked violation of the divergence constraint due to the first term of being clearly non-zero. This divergence-laden field is then advected outwards by the wind flow, thus causing the magnetic field to quickly become non-solenoidal in the outer region as well. (This behavior becomes particularly evident in the left plot of Fig. 3.) Since the resulting magnitude of ∇ · B in the non-solenoidal interface layer is inversely proportional to the layer's thickness, the problem cannot be avoided by choosing a different matching method (i. e. a different matching function r → f (r)). Note that this line of reasoning includes the case of doing no averaging at all: This simply corresponds to the limiting case f = f step , where Since this non-solenoidal layer is actively re-created every time the newly computed outer solution is connected to the inner boundary, we may conclude that a suitable divergence cleaning procedure must kill the divergence immediately afterwards in one step (as the projection scheme does), rather than only damping/ transporting it away on somewhat longer timescales (GLM). We must therefore conclude that for investigations of this kind, the presence of an inner inflow boundary is, at least, difficult and may, in some cases, even preclude the use of the GLM scheme for divergence cleaning. (Note however, that the applicability of GLM to other settings lacking such an internal boundary remains unimpeded by this finding.)

The projection scheme
The so-called 'projection method' was originally developed by Chorin (1967) for simulations of inviscid flow, and later applied in the context of MHD simulations by Brackbill and Barnes (1980). It solves the Poisson equation for Φ and then subtracts ∇Φ from B to ensure ∇ · B = 0. While numerically expensive, it is able to reduce divergence errors down to machine accuracy, and will therefore be used in all simulations presented here.

Types of boundaries
The computational volume consists of a brick-shaped region of space covering 100 × 70 × 50 cells in the x, y, and z direction, respectively. Each cell is a cube with a side length of typically ∆x = ∆y = ∆z = 0.1, implying a coverage of [5, 7, 10] R of real space. (We note that this relatively coarse resolution was chosen deliberately do demonstrate the excellent symmetry-maintaining properties of the employed scheme, see also Fig. 5 of Sect. 5.2. Higher spatial resolution, however desirable for the study of fine-scale structures, would tend to diminish the magnitude of numerical artifacts by which the scheme's performance could be judged, thereby hampering the usefulness of this demonstration.) The Sun's center is located at the origin, with the dipolar axis pointing into the positive z direction. The computational domain is surrounded by two layers of 'ghost cells', whose values are updated after each time step either from symmetry considerations (for 'mirror' boundaries intersecting the origin), or use of outflowing boundary conditions (at the actual 'outer' boundaries). The solar surface, which is represented by a sphere of unit radius located inside the computational volume, obviously does not coincide with any of the Cartesian coordinate surfaces, and therefore requires special treatment, which is discussed in detail in Sect. 5.1 and 5.2. This inner boundary is particularly delicate since it constitutes the surface from which the solar wind emanates, such that numerical artifacts imposed by an imperfect treatment of this boundary will be quickly advected through the entire domain.

Initial conditions
The generic setup for quiet-Sun solar wind simulations is as follows: At t = 0, the simulation is initialized with a radially symmetric wind flow u(r) = u(r) e r with such that a super-sonic value u m is reached at the innermost boundary point, r m . This is done to ensure that the initial velocity at the outer boundary is as small as possible to allow large time steps ∆t, while at the same time being large enough to prevent numerical boundary artifacts from being transported inwards. The density scales as ρ(r) ∝ 1/r 3 , and the temperature is equal to a constant T c . The initial magnetic field is implemented using where F (r) = P 0 /r yields a dipole of strength P 0 that is aligned with the z axis. Since the projection scheme described in Sect. 3.4.2 will operate on the entire grid, the singularity of Eq. (18) at r = 0 must be avoided. This is achieved by choosing a suitably matched polynomial for F (r) inside some small sphere around the origin. (Note that the radius of this sphere must be chosen at least several grid cell sizes smaller than unity to prevent the non-zero current density associated with F (r) = P 0 /r from causing unphysical Lorentz force accelerations just outside the r = 1 boundary.) 5 Numerical treatment of the solar surface boundary

The interpolation method
The inner (solar surface) boundary, which is just the sphere S := {r| r = 1}, obviously does not coincide with any of the Cartesian coordinate planes, which brings up the question of how these boundary conditions are best represented on the grid. Simple-minded attempts, such as keeping cell values inside the Sun fixed and integrating only those outside, have been tried but were shown to result in block-like artifacts at small radii (essentially tracing the envelope of the set of grid cells considered 'inside') where the problem's symmetry would stipulate spherical contours. While these artifacts plane for a standard solar wind run at time t = 1.5 without correction (left) and using GLM (right). The improvement is substantial but still insufficient due to massive divergence values introduced at the inner boundary (unit circle around the origin). The projection scheme achieves κ ∼ 10 −6 (not shown). Note the different color scales.
would of course diminish as spatial resolution is increased, it seems vital to obtain a high degree of symmetry-keeping already at this relatively coarse resolution, especially in view of the high numerical costs associated with increasing the number of grid cells in a 3-D simulation. After several possibilities have been tried, the following procedure was adopted: 1. At initialization, all grid points which are located outside S but have at least one of their 3 3 − 1 = 26 neighbors inside S are stored in a list I of 'interface points'. (The set neighbors of a cell r ijk is defined as the set of cells r i j k with |i − i |, |j − j |, |k − k | ∈ {0, 1} excluding r ijk itself.) 2. After each time step (which only advances grid points outside S in time), a weighted average for each variable w ∈ {ρ, ρu x , ρu y , ρu z , B x , B y , B z , e} is computed for each r I ∈ I viā and L Iα := r I − r α , where the sums in Eq. (19) are taken over all neighbors of r I . The choice of weights ∝ (L Iα ) −1 ensures that for r I → 1,w I smoothly tends to the appropriate boundary value. Figure 4 serves to illustrate the situation. When the above procedure is applied to the Cartesian components of vectors such as u and B, it will usually destroy any possibly existing symmetry of these vector fields (e. g., if u is purely radial, the averagedū vectors will slightly deviate from the radial direction). In order to preserve such symmetries, all Cartesian vector components entering the averaging process of Eq. (19) are first rotated until they are parallel to r I before the averaging takes place, thus ensuring that the symmetry is preserved.
3. In order to guarantee that the newly computed grid values for I are independent of the ordering within that list, all computed averages are first stored in a separate field. Only when all thew I are known will they be copied onto the actual grid.
Note that step 1 is executed only once, while steps 2 and 3 are called after each integration time step. The above procedure gives the best results when applied to a scalar field that varies approximately linear in space. Near the solar surface, however, strong radial gradients of density are present. For this reason, it has been found to be advantageous to artificially reduce the density gradient in Eq. (19) by multiplying ρ(r α ) with r α 3...4 before averaging, and consequently dividingρ I by r I 3...4 afterwards.

Velocity extrapolation versus fixed boundary
The averaging procedure of Sect. 5.1 keeps all quantities fixed on S. However, if solar wind configurations such as the Parker wind solution (Parker, 1958) are to be reproduced, it seems questionable to apply this procedure to the velocity, since the requirement that r → u(r) must pass through a critical (sonic) point completely determines the solution topology, and thus eliminates the freedom to prescribe a fixed (Dirichlet) boundary value at r = 1. Different possibilities are conceivable to handle this problem: 1. Allow the velocity near S to adjust freely by radial inward extrapolation of the time-advanced solution (Keppens and Goedbloed, 2000), or 2. enforce a fixed value for u on S in spite of the above problem, and accepting that (hopefully small) inaccuracies will be introduced at small radii (Manchester et al., 2004).
While the second alternative is just what the above averaging procedure does, the first option, while being straightforward in spherical coordinates, is clearly non-trivial to implement on the present Cartesian grid. In analogy to the averaging scheme used for the other variables, the adopted procedure (which replaces the procedure of Sect. 5.1 for u) is as follows: 1. Prior to initialization, a list A of all grid points r A with 1 ≤ r A ≤ 0.5 is set up and sorted by decreasing r A (such that the outermost points will be processed first).
2. For each r A ∈ A, a sub-list of grid points r A,i is created, such that (i) r A < r A,i and (ii) r A − r A,i < r 0 (with r 0 ≈ (2...3) ∆x).
In other words, the sub-list for r A contains grid points close to r A which are located at larger radii than r A itself. (Note again that steps 1 and 2 are executed only once.) 3. After each time step, the radial mass flux is computed from the sub-list at each r A , and a least-squares fit of the linear function g A : r → c 0,A + c 1,A r is used to find the mass flux at r A (which is then given by g A ( r A ) = c 0,A + c 1,A r A ). Finally, the corresponding radial momentum is immediately afterwards written to the grid, such that its value is available to the extrapolation at the next point in the list. Figure 5 shows a comparison of both methods for the (unmagnetized) Parker wind case. The interpolation method's superior performance is in the range of a few percent only and has to be gauged against its increased computational effort. Consequently, all of the simulations presented here use the Dirichlet method. (We note, however, that this may not always be appropriate when different parameter ranges are used. For instance, a higher base temperature T c will move the sonic point sunwards, leading to a higher flow velocity at the solar surface, and a presumably larger discrepancy to the zero-velocity condition.) 6 Solar wind and CME simulations

Creating equilibrium wind solutions
For the initialization of our CME expansion studies, we first seek a well-defined MHD equilibrium resembling a 'quiet' (i.e. stationary) setting during solar minimum. While this is of course not strictly required for such studies, it is nevertheless vital for the interpretation of the obtained results, since only then can structures like CMEs be clearly disentangled from the dynamics of the background flow. For magnetized, isothermal (γ = 1) winds, the system starts from the initial conditions of Sect. 4.2 and then quickly (within a few sound crossing times) settles into a stable equilibrium similar to the one depicted in the first frame of Fig. 6. At a distance of 5 R from the origin, the outflow velocity in x direction differs from that at the polar field line by a factor of about u(5, 0, 0) u(0, 0, 5) = 1.26 c s 3.13 c s ≈ 160 km/s 400 km/s = 0.4 , which is due to the retaining force of the closed magnetic field lines in the equatorial region, and reminiscent of the  = 1, green), and the situation after the stationary equilibrium has been reached (t = 8, red). Since all grid points are shown, the scatter at a given radius can be seen as a measure of the simulation's spurious departure from radial symmetry.
speed difference between the fast and slow solar wind. Examples of non-isothermal hydrodynamical runs integrating the full energy equation (4) with the heating source term (8) can be found in (Kleimann, 2005). It is noteworthy that essential features of the quiet inner heliosphere, such as a latitudinal dependence of outflow velocities resembling the fast and slow solar wind and the poleward transition from closed magnetic field lines (which span a static 'dead zone') below about 40 • of latitude to an open, more radial field, are self-consistently reproduced by our model. In particular, it was found to be unnecessary to invoke the method of latitude-dependent inner boundary conditions used by other authors (Keppens and Goedbloed, 2000;Manchester et al., 2004) to reproduce this dichotomy: The magnetic dipole strength P 0 proved fully sufficient to control the latitudinal extent of the closed-field helmet zone. As can be intuitively expected, a stronger B field at the surface will tend to conserve its arch-shaped closed structure, while in the limit P 0 → 0, all field lines will be stretched out radially by the flow, and spherical symmetry is recovered. The choice of P 0 = 4 results in the intermediate case with an open/ closed transition near ±40 • of solar latitude.

Initialization of CME onset
The present investigation focuses on the aspect of CME propagation, rather than on their actual nascency. Therefore, a simplifying approach similar to the one already employed by Groth et al. (2000) and Keppens and Goedbloed (2000) will be used. This approach is based on a time-dependent boundary condition at the solar surface, generating a transient, isothermal increase in density (and thus in pressure). If chosen sufficiently strong, this density excess is able to tear open the equatorial helmet streamer, causing the detachment of the excess matter as a rapidly expanding bubble. In order to initiate an eruption in the time interval an additional, localized mass flux ρ add u add with is released at a pre-defined location on the solar surface (implying r = 1 for the remainder of this section). Without loss of generality, let the center of the eruption region be in the plane ϕ = 0, such that its location is just For fixed time t, the value of ρ cme (r, t) should only depend on the angular distance α(r, r cme ) = arccos(r · r cme ) between r and r cme , such that ρ cme (r, t) possesses axial symmetry with respect to the r cme axis. Following Keppens and Goedbloed (2000), we employ the function with E(r, t) := sin 2 π t − t cme τ cme cos 2 π 2 α(r, r cme ) δ cme (27) which connects smoothly to the undisturbed state in both space and time. Here 2δ cme denotes the angular diameter of the circular eruption region Ω cme , which is defined as the region where ρ cme (r, t) gives a non-zero contribution according to Eq. (26), and which thus covers a total solid angle on the Sun's surface. The total mass released by the CME's eruption can be estimated as with ω being the solid angle. For δ cme = 30 • = π/6, this translates to physical units as a typical value for a strong CME.
6.3 CME expansion runs CME expansion runs have been carried out at various combinations of CME strength, heliographic latitude, dipolar field strength, etc. We first describe typical runs involving only a single CME, while the case of multiple events is deferred to the ensuing section. Unless indicated otherwise, the launch parameters f 0 = 16 and τ cme = 1 were used.

Single-event runs
The panel of Fig. 6 shows selected snapshots of a typical simulation run involving an isolated CME event. The first frame depicts the equilibrium situation of the pre-eruptive state. Following its initiation at the solar equator, the CME rapidly expands outwards, thereby quickly gaining both in size and speed. Note in particular the structures indicated by the kinked magnetic field lines that could lead to shocks in the non-isothermal case. In the third frame, while still continuing to accelerate, the CME reaches the volume boundary, thereby dragging the field lines outwards and deforming them almost radially. In the final frame, the CME has completely left the simulation volume and the system has relaxed into a state similar to the quiet initial situation.
Taking advantage of the fully 3-D nature of our simulations, we can also access the dynamics in the perpendicular planes. Figure 7 shows time frames from the same run, this time viewed as a contour of the sharp, wall-shaped B z signature which arises when the CME runs into the background magnetic field and forces it to pile up ahead of it. As can be expected, the magnetic front moves fastest in the x direction, thus forming an elongated shell around the CME's core. On the opposite side, the CME's wake shows a marked reduction in field strength, which even includes an expanding region of reversed field direction trailing the CME. The region's growing extent is particularly evident from the dotted wedge discernible in Fig. 8. Note that the steep outward slope of B (∝ r −3 for a dipole) makes it necessary to normalize the values appropriately.
To analyses the dynamics of the CME as a whole, a reliable tracer of its position is required. While the CME's density shows relatively large and irregular fluctuations which make it difficult to use it to monitor its location, we found the magnetic field signature of Fig. 7 to be more suitable for this purpose. Figure 8 may serve to illustrate this idea. From the resulting [t, x(t)] curves, we derive terminal velocities of 3.5 c s ≈ 450 km/s and 4.5 c s ≈ 580 km/s for CMEs launched with an initial velocity of u cme = 0 and 2, respectively.

Interacting CMEs
With the rate of CME occurrence reaching several events per day during solar maximum, it is not unusual to find more than one CME to be present in a given section of interplanetary space, a fact which motivates the numerical study of the interaction of CMEs. Simulations of this kind have been carried out by various authors (Vandas et al., 1997;Odstrčil et al., 2003;Schmidt and Cargill, 2004;Wang et al., 2005). Interacting CMEs have also been linked to the modulation of type II radio bursts (Nunes, 2007), and their importance for SEP generation has been investigated by Gopalswamy et al. (2005) and Vandas and Odstrčil (2004) using 2.5-D flux rope simulations. More recently, Lugaz (2008) has connected earlier simulations (Lugaz et al., 2005(Lugaz et al., , 2007 employing the BATS-R-US code (Manchester et al., 2004) to actual LASCO data by means of synthetic observations. Fig. 6. Time sequence of simulated CME expansion from pre-eruption (t = 10.0) to expansion and return to equilibrium (at t = 20.0), showing contours of u (top) and log 10 n (bottom) in the poloidal plane (y = 0), with magnetic field lines superimposed. The CME speed at onset was chosen to be ucme = 2. An MPEG movie of this simulation, which covers the entire simulation from initialization (t = 0) to convergence into steady-state (near t = 10), CME expansion and back to near-equilibrium (t ≈ 20), is available from the supplementary material at http://www.ann-geophys.net/0/1/2021/angeo-0-1-2021-supplement.mpg.
While it is clear that at this initial stage, our simulations cannot be expected to rival the existing work in either detail or scientific content, we can nevertheless demonstrate our code's general applicability to this important sub-class of CME phenomenae. Figure 9 shows selected snapshots of a corresponding simulation run: At t = 10, a slow (u cme,1 = 0) CME is initiated along the x axis, to be quickly followed by a faster one (u cme,2 = 2) launched at t = 11.0 into the same direction. (Note that this terminology is merely used to distinguish both entities from each other. We do not intend to relate these to the slow/ fast dichotomy known from actual CME observations. As was shown at the end of the preceding section, both simulated CMEs would qualify as 'slow' in this sense.) Both CMEs not only exhibit the individual effects of acceleration, expansion, and field line kinking already found and discussed in the previous case of an isolated CME, but there is apparently also a noticeable interaction between the two as the second CME gains speed and eventually collides with its predecessor. Note again the kinked magnetic field lines along with a corresponding density gradient, both related to discontinuities that would develop into shocks in the non-isothermal case. It is also interesting to observe that the prescribed density excess is sufficient to trigger a spontaneous, self-consistent outward acceleration, without the need to artificially 'push' the CME forward by enforcing a non-zero initial velocity at the instant of its launch.

Connecting to observations
Due to lack of in-situ data at small distances from the Sun, a direct comparison between simulation and actual CME data is currently not feasible. In order to at least qualitatively connect the simulations presented here to observations, five fixed locations at radii r b ∈ {2, 4, 6, 8, 10} were chosen along the CME's trajectory (i. e. the x axis). At every time step, the values of the non-vanishing variables [B z , n, u x ] at these locations were extracted from the simulation data and then combined in the panel of Fig. 10. Thus, a time profile of these quantities is generated, as it would be seen by a stationary observer while the CME moves past his location. (It should be noted that the profiles for particle density n and magnetic field B z at radius r b have been multiplied with (r b ) 3 and −(r b ) 2 , respectively, since otherwise the effect of radial dilution would not have allowed curves of various radii to be presented compactly in a single viewgraph. This obviously only changes the relative size of two profiles against one another but leaves the shapes of individual profiles un- The thick dashed line connects the respective maxima, thus forming an t → x(t) position curve for the magnetic peak leading the CME. The solid line shows the corresponding x(t) plot for the faster CME (ucme = 2 rather than 0). The corresponding contour stripes for this second CME are not shown. changed.) Using Fig. 11, these plots can be contrasted with a compilation of the temporal evolution of the solar wind's MHD properties, as measured by the Advanced Composition Explorer (ACE) for a magnetic cloud passing the probe's location, the inner Lagrange point at a heliocentric radius of 0.99 AU. A number of qualitative similarities between observation and our simulation can indeed be identified; especially the sharp rise and slow decay of the velocity's maxima is clearly discernible in both cases. The sharp, almost needleshaped peaks of the magnetic field profiles also exhibit a striking similarity. These pronounced field enhancements are induced by the CME's compression wave, and even seem to increase further as the driving CME accelerates outward. The notable differences in the total duration of passage (about one day for the magnetic cloud opposed to about one hour in the simulations) can easily be accounted for by the very different sites of observation. The cloud had much more time to extend from a presumably rather compact object to its full length of up to 1 AU. Also, the transit time cannot be expected to be totally independent of the duration of CME initiation (which in our case amounts to just 1.5 h real time). However, since observation and simulation stem from very different heliocentric radii, a direct, quantitative comparison between the respective profiles of Figs. 10 and 11 is of course not feasible. Our attempts to identify common features between them can therefore merely serve as a "reality check" on the general usefulness of these first simulation runs. Be-sides, they may serve to illustrate the type of comparison that are intended for future simulations covering the whole radial range up to Earth orbit.

Conclusions
We have reported on the creation of a 3-D MHD model of the near-Sun heliosphere, its numerical implementation and subsequent application to the propagation of CMEs. In order to adequately implement the Sun's spherical surface as an inner boundary on the Cartesian grid, a weighted averaging procedure was devised which is able to handle the huge gradients (most notably of mass density) present at this boundary. The use of this procedure also contributed to a reduction of spurious departures from the problem's underlying symmetry, which result from the fact that the Sun's spherical (boundary) surface cannot be mapped to a Cartesian grid of finite cell spacing. Comparing a Dirichlet boundary condition for the velocity against free inward extrapolation, the latter was found to yield slightly more accurate results, albeit requiring a more complex numerical implementation. To ensure a solenoidal magnetic field, the GLM scheme was found to be inappropriate due to the presence of internal boundaries, and was thus abandoned in favor of a classical projection method. After the model's CWENO-based numerical realization had satisfactorily passed various test cases, it was successfully employed to generate stable, self-consistent MHD equilibria of the quiet, magnetized solar wind. These were then themselves used as initial configurations to simulate the expansion of CME-like plasma bubbles. Since the modelling is fully three-dimensional, the CME's direction of expansion can be chosen independently of the system's axis of symmetry; in particular, it is possible to study expansion within the ecliptic plane. The extracted time profiles of density, flow velocity, and magnetic field strength show qualitative similarities to actual in-situ data obtained from satellites at much larger heliospheric distances. The fact that such similarities can be found lends support to the notion that the main physical processes which shape the structure of a CME occur shortly after onset, whereas the ensuing phase of interplanetary propagation is merely characterized by dilution and (almost) self-similar expansion, although a direct simulation covering the entire range up to Earth orbit will be needed to make unambiguous statements about the CME's IP evolution and its persistent self-similarity (or lack thereof). In a future extension of this work, we intend to merge heated (i. e. non-isothermal) scenarios with magnetized wind models, a step which, however desirable, could not yet be carried out due to remaining numerical difficulties. This direction seems even more promising since both aspects have been proven to yield satisfactory solutions individually. On the model side, we plan to include additional aspects (such as localized heating and changes in magnetic topology) into the CME's initialization to trigger its eruption. Since this will require a much higher grid resolution near the solar  Fig. 10. Simulated profiles of magnetic field, density, and fluid velocity as seen by static observers situated along the CME expansion direction at radii r/R ∈ {2, 4, 6, 8, 10} for the sequence of Fig. 6. The minus sign at Bz compensates for the magnetic field's north-south polarity (which has Bz < 0 at z = 0). Time is given in hours after CME onset. Temperature profiles are not shown due to γ = 1. Fig. 11. Actual in situ measurements of a magnetic cloud near 1 AU, adopted from Burlaga et al. (2001). The general shape of these profiles is to be compared with the simulation results depicted in Fig. 10. surface, a recourse to adaptive mesh refinement and/or parallelization becomes mandatory.