sunRunner1D: A Tool for Exploring ICME Evolution through the Inner Heliosphere

Abstract

Accurate forecasts of the properties of interplanetary coronal mass ejections (ICMEs) prior to their arrival at Earth are unquestionably a key goal for space weather. Currently, there are several promising techniques for accomplishing this, including the more focused but limited objective of predicting the time of arrival (ToA) of the ICME at Earth. In this study, we describe a new tool, sunRunner1D, with the initial goal of being able to reproduce the structure and evolution of four categories of CMEs as they propagate from the corona to 1 AU. We demonstrate that sunRunner1D can reproduce the essential properties of these ICMEs to varying degrees of success. We suggest that, ultimately, this tool could assist operational forecasters in predicting space weather events, and their associated geomagnetic consequences. In the nearer term, we anticipate that it could potentially provide useful forecasts for ToA.

1. Introduction

Coronal mass ejections (CMEs) are the manifestation of explosive releases of solar plasma and embedded magnetic fields from the solar atmosphere. Their interplanetary counterparts (ICMEs) are the primary cause of large geomagnetic storms on Earth (e.g., Gosling [1]). CMEs and ICMEs have been observed, described, and catalogued for more than 50 years, and their basic properties and evolution are reasonably well understood. A variety of numerical models, from 1 to dimensional (1-D) hydrodynamic to 3-D magnetohydrodynamic (MHD) have helped further interpret the broad variety of events that make up, at least, coherent CME structures (e.g., Hundhausen [2], Dryer et al. [3], Manchester et al. [4], Riley et al. [5], Odstrcil [6], Török et al. [7], Owens et al. [8]). A subset of these events, the so-called magnetic clouds, have proven even more amenable to modelling efforts, being able to be described by simpler force-free models (e.g., Riley et al. [9]).

Although CMEs and ICMEs continue to provoke interesting scientific questions, a major effort has shifted towards being able to predict, rather than understand, their properties. Currently, NOAA’s National Weather Service’s Space Weather Prediction Center, as well as many other national space weather programs, rely on several tools to make effective space weather predictions. These range from simple empirically derived rules to sophisticated 3-D MHD simulations (e.g., Pizzo et al. [10]), while the former approaches are easy to implement, they are often unreliable, and, similarly, while the latter techniques are often (although not always) more reliable, they can only be run by experts and, moreover, require substantial computational resources.

In this study, we present a new tool, sunRunner1D, which, we believe, bridges the gap between these two extremes. It is based on a sophisticated astrophysical MHD code (PLUTO) [11]; however, as the name suggests, sunRunner1D is run in a 1D, spherically symmetric mode. This, we argue, is sufficient for exploring the evolution of ICME profiles along the Sun-Earth line, which are driven primarily by radial dynamical interactions. By virtue of its simplicity, sunRunner1D can be easily run by novice users, requiring only modest computational resources (e.g., a laptop). We show how sunRunner1D can be used to reproduce the essential features of a range of observed ICMEs, and describe how it can be incorporated into a CME time-of-arrival (ToA) pipeline. Finally, we discuss its limitations and their potential solutions.

2. Materials and Methods

This study describes the development and application of a simple 1-D MHD model. Data from observed events are used to validate the model results and, in some cases, used to provide iterative feedback to improve the match with observations.

2.1. Events

For this study, we focus on reproducing the in situ properties of ICMEs at various locations in the inner heliosphere. These events are summarized in

Table 1. Summary of the four interplanetary coronal mass ejections (ICMEs) modeled in this study.

ICME No.	Date	Spacecraft	Type of CME	References
1	23 July 2012	STEREO-A	Extreme (fast)	Baker et al. [12], Riley et al. [13]
2	21 April 1994	Ulysses	Overexpanding	Gosling et al. [14], Gosling et al. [15]
3	14 July 2000	ACE/WIND	Multiple (3) interacting	Török et al. [7], Tokumaru et al. [16]
4	12 November 2018	PSP	Streamer-blowout (slow)	Nieves-Chinchilla et al. [17], Korreck et al. [18]

. We chose four of them to illustrate the diversity in the types of events that can be modelled by sunRunner1D. The first event was observed by the STEREO-A spacecraft and could best be described as an extreme ICME. The second event was observed during the early phase of the Ulysses mission, while the spacecraft was at high heliographic latitudes. It is a classic example of an over-expanding CME. The third example is the so-called “Bastille-day” event, which includes at least three interacting CMEs. The fourth example illustrates a streamer blow-out CME, which produced a clear flux rope signature at the location of Parker Solar Probe (PSP).

2.2. Model

The sunRunner1D model is based on a simplified form of the astrophysical MHD code, PLUTO [11]. Although the full code base provides many more features and capabilities than we leverage for sunRunner1D, we chose to leave the code “untouched”, relying only on user-defined boundary conditions to read in data, set the initial conditions, and generate the ICME pulses. By so doing, sunRunner1D should remain reasonably compatible with any future upgrades to PLUTO, and will not require the maintenance of a separate and ad hoc version of the code.

In addition to several user-defined routines that interface directly with the PLUTO code, a more user-friendly Python script is provided for initiating runs. This can be run from the command line. In the current “stable” release, the user specifies the following parameters: (1) density ($n_p$), radial velocity ($v_r$), transverse magnetic field ($B_{\varphi }$), and temperature ($T_p$) of the solar wind upstream of the ICME; (2) $n_p$, $v_r$, and $B_{\varphi }$ for the peak values within the ICME; and (3) the duration of the ICME (${\tau }_{duration}$). Several additional features are provided in the “unstable” version, including: (1) choice of ICME profile (e.g., bell, box, declining profile); (2) multiple ICMEs (currently up to three ICMEs are supported); and (3) support for multiple realizations allowing ensemble solutions to be computed. For this study, we use the “stable” version of sunRunner1D.

CMEs are modelled using simple pulses at the inner radial boundary of the simulation (30 R_S). This evades the difficult (but important) questions associated with their initiation, eruption, and early evolution; however, it allows us to explore the evolution of a wide range of pulses efficiently and using modest computational resources. The pulse can be specified by a variety of profiles, however, here, we limit them to sinusoidal variations. Each of the variables ($v_r$, $n_p$, and $B_{\varphi }$) can be modulated independently in an effort to best match the structure of the CME at 30 R_S. Finally, the duration of the pulse can be set arbitrarily, which impacts the resulting evolution of the structure. To the extent possible, the boundary conditions for the pulses are guided by the available remote observations. White-light observations, for example, can be used to infer the initial speed of the CME in the corona. However, other important amplitudes, like the peak density and field strength, must be indirectly inferred, say, from the observed signature at 1 AU. A rough estimate of the peak density can be made for events where a reliable mass of the CME can be determined. However, this cannot be done with any regularity. Additionally, in principle, estimates for the peak transverse field within the ejecta could be inferred from properties of the structure at the Sun producing the event (e.g., active region). However, again, here, we rely on inferences drawn from the in situ measurements to specify these values.

For completeness, in

Table 2. Summary of the input boundary conditions for the four simulations analysed in this study.

Event	${\mathit{v}}_0$ (km s−1)	${\mathsf{\Delta }}_{\mathit{v}}$ (km s−1)	${\mathit{B}}_0$ (nT)	${\mathsf{\Delta }}_{\mathit{B}}$ (nT)	${\mathit{n}}_0$ (cm−3)	${\mathsf{\Delta }}_{\mathit{n}}$ (cm−3)	${\mathit{T}}_0$ (K)	${\mathit{\tau }}_{\mathit{duration}}$ (h)
1	410	1977	15.6	126.6	69	1870	7.15 × 105	19.6
2	462	0.0	10.4	55.4	226	1208	5.0 × 106	7.0
3	616	238	48.7	976	377	0.0	1.43 × 106	17.7
4	250	0.0	59.8	1195	1276	1914	9.2 × 105	3.0

we summarise the parameters used to drive each of the four simulations. The values with subscripts ‘0’ correspond to the properties of the ambient solar wind ahead of the CME, and into which the ICME is launched. These were obtained from the in situ measurements shortly before the ICME’s arrival at the spacecraft. Within the constraints of a predictive environment, to allow a longer lead-time, these could be inferred from observations of the background solar wind made one solar rotation earlier, or from forecasts made by global MHD models, such as WSA+ENLIL e.g., [19], EUHFORIA, e.g., [20], or CORHEL, e.g., [21]. These values are set at the inner radial boundary of the calculation and the simulation is allowed to run forward in time until a steady-state equilibrium is achieved. The values denoted by Δ describe the ICME pulse. These will be discussed in more detail below, when we present the results.

2.3. Data

In-situ measurements at Earth were obtained from NASA’s Space Physics Data Facility (SPDF) using SPEDAS (Space Physics Environment Data Analysis System). Specifically, we retrieved: (1) OMNI data, which is an aggregation of data from multiple near-Earth spacecraft (primarily IMP 8, Geotail, Wind and ACE), providing solar wind measurements upstream of the Earth [22]; (2) Ulysses data; and (3) PSP data. In all cases, we focused on the bulk solar wind parameters, including speed, number density, temperature, and magnetic field strength.

As noted earlier, where possible, initial properties of the CMEs in the corona were estimated from white-light observations from Earth-based spacecraft as well as from STEREO A/B. If the analysis had been performed and published before, we were guided by these estimates.

3. Results

In the sections that follow, we introduce each ICME modelled by summarizing the main features of the event, describing the boundary conditions used to generate the model solutions, and comparing the model output with the available observations. We focus on event 1 in more detail, omitting some of the supporting visualizations for the later cases in the interests of brevity. However, all visualizations were created for all four events and can be recreated using the scripts provided with the accompanying GitHub repository https://github.com/predsci/sunRunner1D-manuscript (accessed on 20 August 2022).

3.1. Event 1

Our first event is the well-studied extreme CME of 23 July 2012. With a peak magnetic field strength of more than 100 nT, a peak velocity in excess of 1300 km s⁻¹ and a travel time from Sun to Earth of ∼20 h, it ranks as one of the most extreme events in recorded history and certainly of the space era [13]. However, and fortunately, its trajectory took it past STEREO-A, which was well separated from Earth’s longitude. In fact, no perturbations at Earth were obviously associated with the event. Several studies have analyzed the solar and interplanetary aspects of this event [12,13,23,24,25,26], including the potential formation of a “cosmic-ray-modified” (CRM) shock [25] and the possible geomagnetic consequences that would have occurred had the event impacted the Earth’s magnetosphere [12]. Here, we summarize the main features of the observations as they relate to modelling the event. First, it is likely that there were two prominence eruptions, separated by 15 min, which led to a compound ICME at 1 AU [27], which, for simplicity, we will model as a single event. EUVI observations by the SECCHI instrument onboard STEREO-A suggested that the CME was launched at 0208 UT, while in situ measurements localized the ToA of the event to 2255 UT on the same day. These provided a relative start time and travel time for the model to reproduce. Based on this, the average speed of the ICME would be 2000 km s⁻¹. Additionally, time-elongation maps from STEREO-B’s COR 2 instrument suggested two CMEs being launched with speed estimates ranging from 2500 ± 500 km s⁻¹ [26] to 3050 km s⁻¹ [27].

The in situ measurements are shown in Figure 1, together with model profiles. At this point, we restrict our discussion to the observations, returning to the model results later in the discussion. Unfortunately, the Plasma and Suprathermal Ion Composition (PLASTIC) investigation onboard STEREO-A was not able to provide reliable estimates for the plasma density or temperature. The values shown here were inferred from the >45 eV electron density measurements for plasma density (by multiplying this by a factor of five) and velocity to infer a temperature estimate [27]. Thus, the “observed” temperature is an extremely poor proxy for what was likely observed, whereas the density can probably provide some value, albeit with caution.

The measurements show that a strong F shock arrived at the spacecraft approximately 20 h after the eruption. This corresponds with an increase in density and magnetic field (the temperature increase merely mirrors the fluctuations in speed). The event from the peak in speed (shock front) until the end of the ICME (as determined from the individual components in the magnetic field) lasts ∼28 h or ∼24 h from the start of the ICME to its end. Although the two CMEs merged within the field of view of COR 2, they can still be identified as distinct magnetic structures within the compound ICME at 1 AU (blue dashed lines). However, we will model only the first ICME, recognizing that some of the structure contributed by the second CME, will not be reproduced. The reason this is possible is that the first event is traveling faster than the second one. Thus, the presence of the second one added relatively little to the dynamics of the first one, and, in particular, its interaction with the ambient solar wind ahead. On the other hand, had the dynamics of the second CME been the objective, this would certainly have required consideration of the earlier event.

To simulate this ICME, we developed an ambient solar wind solution that matched the in situ measurement at the STEREO-A spacecraft approximately 12 h prior to the arrival of the first event. To mimic the dynamic properties of the CME at the inner boundary of 30 R_S, we smoothly raised the speed profile by 1977.4 km s⁻¹, the density profile by a factor of 30.0, and the magnetic field strength by a factor of 9.1. All of this was over a duration of 19.6 h. These boundary conditions are summarized in Figure 2. The ambient solar wind conditions were first established during the interval from $t=0$ h to $t=332$ h, at which point the pulse was inserted. At the beginning and end of the pulse, we introduced a massless tracer particle to track the location of the boundaries of the modelled ICME (dashed vertical lines). This also provides a method to calculate the radial width of the ICME (see below).

In Figure 1, we directly compare the model results and the in situ measurements. In general, the speed profiles match one another well and the results are consistent with those of earlier simulations (e.g., Riley et al. [13], Ngwira et al. [26], Desai et al. [28]). The ramp-up prior to the shock is not reproduced in the model. This is believed to be related to the presence of a large population of energetic particles creating a foreshock, and is beyond the capabilities of the MHD code. The general form of the density and transverse field are well matched, at least for the first CME. The second CME had a significant flux rope embedded within it, providing the sustained $B_{\varphi }$ values after the initial peak. Our modelled CME attempted to treat both events as a single CME and, thus, without allowing for a multi-peak profile, would not be able to mimic this structure. However, it did match the general leading edge of CME 1 and the trailing edge of CME2 (dashed vertical lines), which were obtained by advecting the tracer particles through the heliosphere and identifying the time that they crossed 1 AU. As noted earlier, the “observed” temperatures before and during the event were only indirectly estimated. Nevertheless, we used the upstream temperature value to drive the ambient solar wind solution, and this likely contributed to the large and unrealistic peak in the computed temperature profile within the ICME. Finally, we emphasize that the ToA of the modelled ICME is a free parameter in the sense that we can adjust the model time to best match the observations. However, as we will show below, the ToA for this event was very close to that observed.

In Figure 3, we show the profiles of the ICME at three different points in time, roughly corresponding to its arrival at 0.75 AU, 1.0 AU, and 1.5 AU. We note that, although to a first approximation, the profiles shown here (and particularly the red profile corresponding to values near ∼1 AU) can be inferred (at least qualitatively) by reversing the profiles in Figure 1, they do differ in some key respects. This presentation shows the intrinsic spatial structure at three points in time, whereas in Figure 1, we see the structure as it passes 1 AU. Thus, the trailing edge of the ICME in the earlier plot has evolved more by the time it reaches 1 AU. In this case, because the ICME was travelling so quickly, these differences are small; however, the effect can be more pronounced for other events, particularly those “coasting” along with the slow solar wind. A second point worth noting is that the ejecta itself has expanded considerably, moving from 150 R_S to almost 300 R_S. However, it has retained its overall two-step velocity profile. Finally, we note that by visually interpolating between the F shock location of the blue (20 h) and red (26 h) curves, we see that the disturbance arrived at 1 AU in under 23 h following launch, which matches well with the ∼20 h transit time estimated from the observations.

We can investigate the expansion of the ICME by following the difference in location of the tracer particles at the leading and trailing edges of the ICME. This is shown in Figure 4. Interestingly, in this case, the ICME maintains a constant expansion with respect to distance from the Sun, at least in the range from ∼90 R_S to ∼130 R_S. Since the simulation’s inner/outer boundary was set at 30 R_S/1.2 AU, only the interval where both boundaries are within the simulation domain are shown. This profile is not surprising for such an extreme event, which moves rapidly through the solar wind. Thus, it arrives at 1 AU remarkably quickly and is thus relatively “unevolved”, at least in terms of the pressure gradients acting on the ejecta itself.

As a final illustration, we can combine the visualizations in Figure 1 and Figure 3 into a single plot. This is shown in Figure 5. From this perspective, we can see how an initially super-fast velocity pulse bifurcates into two distinct waves; a fast forward (F) wave to the left and a fast reverse (R) wave further to the right. With increasing distance from the Sun, the two separate. We can also use this visualization to directly estimate the speed of the modelled fast-mode shock (as described in more detail by Riley et al. [13]). A tangent drawn along the shock front at 1 AU has a slope of almost ∼3000 km s⁻¹, which is roughly consistent with the speeds estimated using the in situ measurements [13].

In summary, the simulated profiles for event 1 appear to have captured the essential features of the observations, with the notable exceptions of: (1) the foreshock, which is believed to have been mediated by energetic particles, and, thus, obviously beyond the capabilities of an MHD-only algorithm; and (2) the structure of the second embedded flux rope.

3.2. Event 2

Our second case study comes from a class of CMEs known as “over-expanding” CMEs [15]. In general, they are CMEs observed in the high-speed solar wind (and, thus, generally observed by the Ulysses spacecraft) that maintain speeds comparable to that of the ambient, surrounding solar wind ( 750 km s⁻¹). Additionally, they are bounded by a pair of fast-mode shocks; a F shock propagating ahead of the ejecta and a R shock propagating behind [29]. The shock pair is produced from an initial high central pressure CME that drives F and R waves away from it and which then steepen into shocks.

Here, we use sunRunner1D to better interpret one particular event occurring in April 1994. At the time, Ulysses was located at 3.2 AU from the Sun, at a southern latitude of ${61}^{\circ }$. The event has been described in more detail by Gosling et al. [14]. The CME pulse was created as a pure density pulse, with the CME speed and magnetic pressure being set equal to the ambient solar wind values. The density was increased from 226 to 1434 cm⁻³ over a period of 3.5 h. In reality, it is quite possible that some or the majority of this additional pressure was provided by the magnetic field. However, given no guidance from remote solar observations, we elected to use plasma density to drive the evolution. The extent to which this might be incorrect can be inferred from comparisons with the data at 3.2 AU.

Figure 6 compares the plasma and magnetic field variables with the in situ measurements in the same format as Figure 1. We note the following points regarding the observations. First, the speed of the CME was approximately the same as the ambient wind both ahead and behind it. Thus, all evolution is driven by gradients in pressure. Second, ejecta material was inferred based on the presence of counterstreaming suprathermal electrons, which are usually a good indication that the field lines are connected back to the Sun at both ends [30]. Third, the speed profile within the ejecta was monotonically declining, suggesting that this structure was expanding. This is consistent with the minima in density and temperature within the event, suggesting that the material had expanded. Fourth, the R shock lies farther away from the centre of the CME than the F shock. The radial extent of the CME was ∼0.5 AU, although the shock-to-shock distance was 1.3 AU [14].

Comparing the simulation results with the observations reveals several notable similarities, but also some differences. Overall, the profiles in speed, field amplitude, density and temperature match well. This is particularly true with respect to the duration of the event and the size of the fluctuations. The F wave, in particular, seems to have been well captured by the model. The R wave, while modelled well in terms of amplitude, arrives some 4 h after the observed event. With several more iterations, perhaps reducing the duration of the pulse modestly, but possibly also adjusting the peak in density, this could be “fixed” in the sense of bringing the curves closer to each other. However, it would add little value to interpreting the observations and would represent more of a curve-fitting exercise.

Model variations in temperature show a trough that is considerably lower than the observations suggest. This is likely driven by the strong expansion of the ejecta. Had we included a small positive variation in the temperature, this may have mitigated this difference. On the other hand, though, we have no reliable observations of the core CME temperature. Thus, again, we did not attempt to chase down this difference. Of more importance, perhaps is that while the plasma density and magnetic field amplitude are broadly consistent with the observations, the former contains a small positive excursion in the centre of the CME, the latter contains a small valley. This perhaps suggests that the density gradients for this event should have been more equally shared between the two parameters. This would not be inconsistent with our intuition that CMEs are magnetic structures and, at least close to the Sun, are likely to have stronger than average fields within them.

Focusing next on the width of the ICME, in Figure 7 we show its evolution from approximately 1.9 to 2.9 AU. Although the ejecta continues to expand, the rate is slowly and smoothly decreasing. These results are consistent with those of Gosling et al. [14], which relied on a hydrodynamic 1D code to study a selection of over-expanding CMEs observed by Ulysses.

Finally, as before, we can visualise the evolution of this over-expanding ICME in two dimensions. Figure 8 summarizes speed, magnetic field, density, and temperature. Unlike event 1, the profiles for event 2 are much more symmetric as the F and R waves break away from the ejecta, propagating further ahead and behind and into the undisturbed ambient wind. We note that the R wave lies proportionately farther behind the F wave, consistent with the observations. At any particular distance (horizontal slice) it remains substantially weaker than the F wave; however, at any specific point in time (vertical line), the difference is not as significant. We note, parenthetically, the small triangular region in the upper-left of each panel. This reflects the fact that the ambient solar wind had not quite reached the outer boundary of the simulation before the ICME was launched at the inner boundary. This plasma represents transient material created at $t=0$ and not yet “flushed out”; however, it is merely a visual artifact and had no impact on the evolution of the ICME and its associated disturbance.

In summary, for over-expanding ICMEs, sunRunner1D seems very capable of reproducing the observed in situ signatures. Moreover, given their simple initial profiles, the model appears to provide all the salient physical processes to understand their evolution.

3.3. Event 3

Event 3 is the so-called “Bastille Day” CME, named because it occurred on the 14 July 2000. It was and remains one of the most energetic ICMEs to arrive at Earth [7]. In fact, it was not a single event but a sequence of four, all occurring in rapid succession. For the purposes of modelling, however, the step-wise increase in the speed of each subsequent event allows us to model the final and fastest event as if it were a single CME [31]. All that we require is an accurate estimate of the solar wind conditions upstream of it, which is most salient for reasonably capturing the dynamics of the ejecta and its associated disturbance. To model this event, we inserted a combined velocity/field pulse, where the speed was smoothly increased to 854 km s⁻¹ and the field amplitude was increased by a factor of 21.0 (or from a background value of 49 to a maximum of 1024). The total duration of the event was 17.7 h.

Figure 9 compares speed, density, transverse field and temperature for this event. Again, the main features are reasonably well reproduced, including (1) the shock speed, (2) the declining speed profile, (3) the low density within the ejecta, (4) the enhanced field, and (5) the modest decrease in temperature. The most glaring disagreements are: (1) the elevated model field throughout the interval and (2) the strength of the R wave. Intuitively, we posit that the magnetic structure within the ejecta could explain the difference in the transverse fields. In turn, this may have impacted the evolution of the R wave. Without observations to support this, however, any explanation remains speculative.

Figure 10 summarizes the ICME’s evolution in (r,t) phase space. We note several points. First, the initially bell-shaped speed profile rapidly evolves into an asymmetric saw-tooth by 50 R_S. By 70 R_S the highest speed is at the shock front, with the speed smoothly declining behind it until the R shock is encountered. Second, a remarkable feature appears at roughly 70 R_S, corresponding to a distinct boundary in the speed within the disturbance that arcs to the right and upward. This, we suggest, is a second R wave created from the strong compression of plasma at the leading edge of the ejecta. Propagating back toward the Sun (although also being convected out by the solar wind flow) in the rarefied region of the expanding ejecta, it is able to catch up to and merge with the main R wave/shock before 100 R_S. Remarkably, the wave’s speed back toward the Sun is sufficiently high that it balances the outward speed of the entire structure, and it appears stationary (horizontal) in an inertial reference frame for a significant amount of time. Only as it slows down does it begin moving away from the Sun. The processing of the plasma can also be seen in faint wisps in the density (and even fainter in the temperature) along the same arc traced out by the radial speed. This second wave also modulates the structure of the transverse field, but in a more complicated interaction, since the CME began as a high-field region, it produces the lenticular shape at the trailing edge of the CME, which pinches off at roughly the location that the second R wave reaches the first R wave/shock.

Finally, we note that although there is not a direct relationship between CME width and disturbance width, since the F and R shocks move out along straight lines in this phase space, one might infer that the width of the ICME increases linearly with time, which it does (plot not shown).

3.4. Event 4

On November 11, 2018, PSP encountered a streamer blowout (SBO) CME at ∼57 R_S from the Sun [18]. At the time, PSP was already moving away from its first perihelion. It lay just behind the west limb as viewed from Earth, and approximately at quadrature, with respect to STEREO-A. At 0200 UT a small CME was observed to lift off in the COR 1 field of view (FOV) and then continue through the COR 2 FOV. Using a Gaussian fit to the brightness and tracking the maxima, it was estimated that the duration of the CME was 2.5 h and that the speed of the CME at 6 R_S was ∼250 km s⁻¹ and accelerating [18]. Slightly before 12 November 2018, the ICME arrived at PSP’s location. Its leading and trailing edges were readily identified from various features in the magnetic field and plasma parameters, including the presence of counterstreaming suprathermal electrons (suggesting field lines that connected back to the Sun at both ends), lower plasma-β, Mach number, and temperature, increased magnetic field amplitude, and coherent rotation of the magnetic field, suggesting the presence of a complex magnetic flux rope structure [17].

Based on the observations for this event (Figure 11), we attempted to reconstruct the profile using a combination of density and magnetic field fluctuations. Given the large internal field amplitude, we were forced to provide a relatively large peak value. However, since the density did not show any significant deficit, we also had to boost that, although to a smaller degree. As can be seen from the comparison, the model results do not agree well. The field/density perturbations drove a F wave for which the observations had no obvious counterpart, either in the velocity or density. Reducing the field amplitude, while lessening the strength of the wave would also have reduced the amplitude of the model field. Additionally, the observed temperature was substantially lower in the ejecta than the model suggested. Since there was no comparable decrease in the density, we do not believe this lower temperature was caused by the expansion, but likely reflects an intrinsically cooler temperature within the ejecta. However, no obvious remote observations would allow us to infer this. On a more positive note, the modest decrease in density within the event is relatively well matched by the model. Perhaps the most notable comparison lies with the transverse field component. Both the overall shape as well as the relative position of the compression and peak within the CME are captured in the model results. However, this is not surprising since this was the feature that we focused on when developing reasonable boundary conditions for the pulse.

Although it might be possible to iterate sequentially over a range of possible scenarios to infer the likely initial conditions for this event, we present it as a “failed” case to illustrate some basic limitations of the model, particularly when the model is applied within the forecasting domain. SBO CMEs might be expected to be the most difficult to model correctly since their dynamical structure is primarily defined by their initial intrinsic density, field, and temperature structure, none of which can be adequately described by the simple pulses we are using.

4. Discussion

In this study, we have developed an easy-to-use tool for exploring the dynamical evolution of ICME profiles from near the Sun (e.g., 20–30 R_S) to up to 5 AU. We demonstrated its application using four examples representing different classes of ICMEs. We noted where it worked well, and where it likely failed. This tool represents a compromise between simple analytic formulations for propagating CMEs from the Sun to Earth and complex 3-D MHD models, requiring more substantial investment in time and computational resources. We suggest that, in the future, sunRunner1D could be an effective tool for potentially predicting the ToA of ICMEs at Earth, as well as other ICME properties, such as final speed, plasma density, and even magnetic field amplitude.

Our approach comes with several important caveats. First, the assumption that the evolution of ICMEs can be accurately captured in 1-D requires consideration. It is well known that in moving from 1 to D to 3-D, dynamical evolution can proceed more slowly since shear can allow flow to be deflected [32]. To anticipate this effect, one might, for example, specify slower initial CME peak speeds, or, as we have done here, specify a smooth ramp up to the peak and similar profile back down. However, more generally, the 1-D approach also requires that the precise perturbation is chosen along the radial trajectory that intercepts the Earth (or another location of interest). ICMEs modeled in 3-D naturally provide a complete range of profiles capturing evolution across the transverse cross-section of the CME. The 1-D approximation may also have played a role in the relative strength of the simulated F and R waves/shocks. In general, while the model was able to reproduce the properties of the F shock, it often suggested the presence of a stronger R wave/shock that could be inferred from the observations. Generalizations of heliospheric modeling from 1 to to 3-D suggest that R waves/shocks generally form sooner in the 1-D case [33]. Thus, this may be an intrinsic limitation of the 1-D approximation.

Second, although we provide a basic formulation for incorporating the magnetic field within the ejecta, it is certainly limited. At best, it provides a way to include the effects of magnetic pressure into the overall evolution; however, magnetic tension is obviously not present. Moreover, as a single component of what is a vector field, it is difficult to relate this to, say, IMF $B_z$, at least directly.

Third, by setting the inner boundary at 30 R_S, we conveniently evade the more difficult topic of CME initiation and early evolution in the corona, where non-radial effects are likely to be more significant. Additionally, we must specify the properties of the ICME at 30 R_S, while reasonable estimates from remote solar observations for velocity can be made, mass (or density) estimates are more difficult. Magnetic field estimates must currently be guessed, at least in part by comparing with 1 AU observations and then iterating. Hopefully, this exercise will inform future techniques that may be able to obtain magnetic field amplitude estimates from solar observations, such as the properties of the active regions involved in the eruption (e.g., Falconer et al. [34], Linker et al. [35]).

In spite of these caveats, our initial study, we believe, has demonstrated the potential power of sunRunner1D. One obvious application would be to estimate the ToA of ICMEs already observed through the CCMC scoreboard [36]. Initial tests suggest that sunRunner1D can provide potentially accurate estimates of ToA; however, such preliminary calculations require further and more substantial validation before any quantitative statements can be made.

sunRunner1D is a progenitor for a more sophisticated 3-D ICME tool, sunRunner3D, which is currently nearing completion. We are, however, releasing the 1D code as a stand-alone tool since it is likely that some users will either not need the additional capabilities of a global solution or will not have access to the computational requirements that such runs will need. Given the distinct use cases for the two tools, it seems reasonable to maintain them as separate and distinct repositories. Should there be a need to merge the codes together, this can be achieved relatively easily. Additionally, if there is interest in developing a 2-D version, for fitting to heliographic imager data, for example, or to multiple in situ spacecraft (e.g., STEREO A/B and Earth), this can be accomplished by simplifying the 3-D model. In the former case, the dimensions would likely be radial distance and latitude at some specified longitude, whereas in the latter, they would be radial distance and longitude in the ecliptic plane.

The solutions for the four ICMEs studied here were generally reasonable, but caution should be taken in over-interpreting the uniqueness of the boundary conditions that created them. In particular, it is quite possible that different combinations of the pulse’s speed, density and field amplitude can produce profiles at 1 AU that are, at least qualitatively, similar. Moreover, the initial duration of the pulse, in addition to modulating the extent of the ICME at 1 AU, can also lead to differences in the shape of the profile by affecting the initial gradients in the pulse’s other parameters. A wider pulse, in particular, will generate weaker F and R waves, which will take longer to steepen into fast-mode F and R shocks, respectively. This issue of degeneracy can only be broken by accurate inferences of the initial properties of the ICME in the corona. Speed and density (inferred from mass) can, in principle, be captured from white light observations, and temperature variations can be approximated or assumed to match ambient conditions without too much impact. However, at least presently, inferring the amplitude of the magnetic field variations cannot be made easily from available observations.