© R. Biondo et al., Published by EDP Sciences 2021

1 Introduction

Due to interactions between the interplanetary magnetic field and the planetary magnetosphere, Sun and Earth form a complex, closely connected environment whose understanding is crucial to maintaining our industrial society, increasingly dependent on space-based technologies. For this reason, over the last decade the science community, governments and private partners have shown a growing interest in the topic of Space Weather: a broad discipline aimed to the forecasting of potential threats to space and ground-based systems and operations (Schrijver, 2015). Compared to our knowledge of terrestrial weather, the field is still in its infancy: the accuracy of long term space weather forecasting is, at best, modest. It is currently impossible to predict a solar event as a flare or a coronal mass ejection (CME): solar energetic particles (SEPs) can be detected at Earth nearly 10 minutes after their emission on the Sun, while CMEs, reaching speeds up to 2500 km/s, take at least a day to arrive (Cargill & Harra, 2007). Still, it is not easy to predict their geo-effectiveness (hence how the Earth’s magnetosphere will react), because their inner magnetic structure and their velocity is very often unknown. There are not unified models capable to predict a solar eruption’s evolution from the Sun up to the impact with the Earth’s magnetosphere and besides this, the amount of available in situ information is very low with respect to the size of volume in study (see review by Schwenn, 2006).

In order to provide an effective forecasting service and possibly to take protection measures against space-weather threats, it is essential to improve the current ability to predict arrival times and geo-effectiveness of major solar disturbances. To reach this goal it is first necessary to reconstruct the ambient background condition of the heliospheric plasma, at least because these conditions determine the interplanetary evolution of CMEs (see e.g. van der Holst et al., 2007; Maloney & Gallagher, 2010; Vršnak et al., 2013; Temmer & Nitta, 2015), and the interplanetary propagation of SEP streams (e.g. He et al., 2011). The complexity involved in modeling the outer corona, the solar wind and the rest of the heliosphere requires a numerical approach: unfortunately, the amount of information about plasma parameters available to correctly shape this phenomena (mostly provided by in situ spacecrafts located at 1 AU and remote-sensing measurements) is very limited. In addition to this, the dynamic processes that underlie the acceleration of solar wind, CMEs, and SEPs are not yet fully understood: what is known is that the configuration of the magnetic fields of their regions of origin, as well as of their internal structure, must play a key role (Forbes, 2000). However, the determination of magnetic field configurations is not a trivial task, because routine direct measurements of magnetic fields are currently available only for the solar photosphere (see review by Howard & Tappin, 2009).

Different methods and approaches have been developed to reconstruct the Parker spiral, that can be classified into two major categories: backward analytic reconstructions starting from the in situ data at 1 AU, and forward reconstructions starting from the remote sensing data on the Sun. The first category of methods was employed already in the first pioneering works by Schatten et al. (1968); Wilcox (1968), even taking into account multiple spacecraft data (e.g. Nolte & Roelof, 1973; Behannon, 1978) (we refer the reader to the reviews by Schatten, 1971; Ness & Burlaga, 2001, to have a broad historical perspective). These reconstructions, often employed to investigate the solar wind sources at the Sun (e.g. Neugebauer et al., 2002), usually assume a stationary flow from the Sun, and are complicated by the interplanetary evolution of plasma due to stream interaction regions (Burlaga, 1974), magnetic clouds (Klein & Burlaga, 1982), shocks (Dryer, 1974) and any other transient phenomena. The limits of this approach were discussed already decades ago (e.g. Pizzo, 1981; Burlaga, 1983). More recently, thanks also to more advanced computational capabilities currently available, these ballistic-mapping backward methods were reconsidered (Florens et al., 2007), and further developed also to take into account non-zero azimuthal field component at the source surface (Schulte in den Bäumen et al., 2011), angular momentum conservation (Tasnim & Cairns, 2016), and lack of corotation at the source surface (Tasnim et al., 2018).

On the other hand, works based on the second approach usually start from the magnetic fields measured in the photosphere, to reconstruct the magnetic fields in the inner and intermediate corona, and then to expand the derived plasma at a given altitude in the interplanetary medium. The reconstruction of plasma parameters in the whole 3D corona has been performed so far mainly in two different ways: (1) by running a full self-consistent MHD numerical simulation (e.g. Mikić et al., 1999; Yang et al., 2012; van der Holst et al., 2014), and (2) by extrapolating the photospheric fields and filling them with plasma with semi-empirical relationships (e.g. Riley et al., 2006; Pinto & Rouillard, 2017) or evolving them with time-relaxation methods (Feng et al., 2015). Results from many different coronal models have been compared to each other (e.g. Gressl et al., 2014; Asvestari et al., 2019), and constraints from the 3D coronal models expanded up to 1 AU have been also compared with in situ observations (e.g. Lee et al., 2009; Réville & Brun, 2017). Among different extrapolation methods, the Potential Field Source Surface (PFSS) extrapolation (Hundhausen, 1972) gives an approximate description of the coronal fields from the solar surface to the so-called source surface (a sphere with radius between 1.6 and 3.25 solar radii, with 2.5 often used as standard value (Altschuler & Newkirk, 1969; Hoeksema, 1984), by assuming the absence of currents in the inner corona (see Discussion by Lee et al., 2011). This is a strong assumption (solar eruptions occur only in coronal regions where non-potential field configurations are created), and recent results from Parker Solar Probe are now allowing to measure the non-sphericity of the Source Surface (Panasenco et al., 2020), but the PFSS is often used as a well established technique providing a quite good description of the overall coronal field configuration and the location of open and closed field regions (Nitta et al., 2006; Mandrini et al., 2014). Other extrapolation methods exist, but requires high-resolution vector magnetic fields measurements (De Rosa et al., 2009; Aschwanden, 2016), something that (alike LOS magnetograms) is also not currently available for the whole photosphere.

Once the coronal fields have been reconstructed in quasi real-time with similar methods up to a spherical surface concentric on the Sun, forecasting the interplanetary expansion of solar wind requires a method to convert these fields into solar wind plasma parameters. This can be done for instance by coupling a PFSS model with the Schatten Current Sheet (SCS) model (Schatten, 1971), and by applying the so-called Wang–Sheeley–Arge (WSA) technique (Wang & Sheeley, 1990, 1992; Sheeley & Wang, 1991; Arge & Pizzo, 2000), which is the most practical procedure currently in use. In this way, starting from PFSS extrapolation, it is possible to compute the global coronal field up to a surface usually set at 5 solar radii. Then, the WSA model deduces the solar wind speed from an empirical relationship between the divergence of the magnetic field and the proximity of a selected open field line footpoint to the nearest coronal hole.

Alternatively, the coronal magnetic fields and plasma conditions can be reconstructed by running a numerical simulation of the coronal plasma, starting again from photospheric field measurements. Among different models that have been developed by different groups, one of the most established models is the Magnetohydrodynamic Algorithm outside a Sphere (MAS) model (Mikić et al., 1999), that has been optimized over the last decades to provide the best agreement between observed and simulated full disk images of the Sun acquired in EUV bands by different instruments (Riley et al., 2011) as well as to predict the appearance of the solar corona during total solar eclipses (e.g. Mikić et al., 2018). Despite the significant level of agreement between observations and simulations in the distribution of coronal features reached by current numerical models (such as MAS and others), and the level of agreement between the extrapolated magnetic field lines and the orientation of coronal structures, all these methods have in common one limitation. Given the large scales involved in the propagation of solar disturbances, all these methods require a reconstruction of the whole corona, and this requires the knowledge of the magnetic fields over the entire photospheric surface at a given time.

Nevertheless, current real-time measurements of the photospheric fields over the full solar disk are limited only to the visible hemisphere, only to the line-of-sight component of these fields, and little is known about polar fields. To overcome these limitations, the observations of polar magnetic fields on the Sun and of photospheric fields on the far-side of the Sun are currently two of the main scientific objectives of the ESA Solar Obiter mission (Müller et al., 2013). At present, complex flux-transport and flux-balance models (one of which the most successful is the Air Force Data Assimilative Photospheric Flux Transport ADAPT, see Arge et al., 2010, 2013) are used to take into account the evolution of photospheric fields in the hemisphere not visible from the Earth during a half solar rotation (~2 weeks). In any case, these models cannot fully capture the complexity of photospheric evolution (in particular during the maximum phases of solar activity cycle), that is dictated in the end by the complexity of solar dynamo (responsible for emergence of magnetic flux), that is far from being fully understood, and cannot forecast the emergence of sunspots (a problem only partially mitigated by the availability of far-side active region observations derived with helioseismology; González Hernández et al., 2007).

At higher levels, the above coronal models are in general coupled with other models simulating the interplanetary propagation of solar plasma. At the present, the first (Sheeley Jr., 2017) and most used (MacNeice et al., 2018) heliospheric model routinely used for forecasting (ENLIL; Pizzo et al., 2011), works in the WSA frame. This model solves the three-dimensional equation of magnetohydrodynamics starting from inner boundary condition information derived by a WSA model. It uses a flux-corrected transport algorithm to deal with MHD shocks. More recently, the European Heliospheric Forecasting Information Asset (EUHFORIA; Pomoell & Poedts, 2018). As the acronym suggest, this model is entirely developed in European research centers, mainly in Belgium and Finland (Fig. 1). Based on a correction of the WSA approach coupled with an improved model for CME propagation, it still starts from a PFSS model, using synoptic magnetograms provided by the GONG (Harvey et al., 1996) to extrapolate coronal fields on a source surface located at 2.3 solar radii. Then a SCS model propagates these field up to 0.1 AU, where the WSA relation provides plasma speeds and densities, and a MHD code performs simulations of the inner heliosphere based on these data. A schematic of the steps involved in EUHFORIA for the construction of boundary conditions for the heliospheric MHD model is shown in Figure 2.

Fig. 1 Example of results from a EUHFORIA simulation run, the two panels show the distribution on the ecliptic plane (left) and at fixed longitude (right) of the solar wind speed (adapted from Pomoell & Poedts, 2018).

Fig. 2 A schematic of the steps involved in EUHFORIA for the construction of boundary conditions for the heliospheric MHD model.

While powerful, all these methods require considerable computational capabilities, mostly consumed by said extrapolations or reconstructions of the inner heliospheric conditions, even if with a significant reduction of time with respect to the much more complex full MHD 3D models mentioned above. To make matters worse, their forecasting ability is still very far from being optimal. These models are able, in an average sense, to predict the arrival time of major solar disturbances such as CMEs to within ±10 h, but their associated standard deviations often exceed 20 h (Riley et al., 2018). Furthermore, their ability to reconstruct the stationary conditions of the interplanetary plasma is still limited, as shown by the comparisons between model predictions and in situ measurements acquired by various spacecraft (Jian et al., 2016).

The above introduction was aimed at reviewing the current state of the art on this research topic, and also to provide the main motivations for this work. As mentioned, all the above works (and many others not cited here) belong to one of this two categories: 1) ballistic back-mapping starting from in situ measurements acquired at 1 AU (or other distances) and reconstructing the Parker spiral with analytical methods back to the Sun or the source surface, and 2) forward simulations starting from photospheric field measurements on the Sun and reconstructing the Parker spiral with MHD numerical methods out to 1 AU (or beyond). Nevertheless, on the one hand empirical and analytical methods are computationally simple and have a good accuracy for the reconstruction of large scale features, but the models are not time dependent, usually assume the same average speed at all longitudes for the back reconstruction to avoid the well known problem of the stream-line crossing, and need to assume stationarity. On the other hand, “cheaper models provide only partial information about the solar wind, while the MHD models offer insight into the underlying physics which we are trying to understand” (MacNeice et al., 2018), and more than that MHD simulations are able to provide the real-time evolution of the whole spiral at the same time, which is missing in analytical methods, but these models have limited knowledge of real plasma conditions at the inner boundaries.

Hence, what we propose here is a different approach, where the above methods are mixed together combining for the first time the ballistic back-mapping technique with time-dependent MHD numerical simulations. In particular, here we propose to skip the complexity related with the reconstructions of the inner corona to determine the inner boundary conditions for the interplanetary plasma propagation, and to bind the interplanetary numerical simulations to direct measurements of the plasma parameters measured in situ at the distance of 1 AU. To do this, the ambient plasma parameters of the inner heliosphere are reconstructed, using the 1 AU data provided by different spacecrafts located in the Lagrangian point L1 such as DSCOVR and GGS Wind, by going backwards to 0.1 AU following the spiraling arms of the Interplanetary Parker Spiral. This back-reconstruction is possible by assuming, for the conditions of the interplanetary plasma, stationarity in the periods selected for data acquisition. The back-reconstructed plasma parameter at 0.1 AU are thus employed as inner boundary conditions for the solar wind expansion up to 1 AU based on a MHD simulation. This approach provides a new method to reconstruct the real-time conditions of interplanetary plasma from the Sun to 1 AU and beyond. Moreover, the approach described here will be also aimed at preserving the smaller-scale density and velocity inhomogeneities dragged by the solar wind for purposes and future applications that are discussed in details in the paper.

2 Theoretical framework for backward solar wind reconstruction

2.1 The Parker model

One of the most important steps in the work presented here is the back-reconstruction of solar wind plasma parameters in the inner heliosphere (0.1 AU) starting from the measurements acquired in situ at 1 AU. For this purpose, it is necessary to introduce a theoretical model describing the solar wind acceleration and propagation through the heliosphere. Historically, the first theoretical demonstration for the existence of the solar wind was provided in a classic paper by Parker (1958) neglecting the influence of the magnetic field B and of solar rotation, and by assuming that the governing forces leading to coronal expansion are only the pressure gradient and gravitation. In the same paper, Parker (1958) also predicted the general shape of the Interplanetary Magnetic Field (IMF). Including solar rotation and magnetic forces and considering the applicability of the frozen-in Alfvén theorem due to the high electrical conductivity of the solar wind plasma and assuming a high value for β, the magnetic field follows a path which is a velocity streamline in a frame corotating with the Sun and defined by

(1)

where Ω_⊙ = 2.8 × 10⁻⁶ rad s⁻¹ is the angular rotational velocity of the Sun’s equator. Treatment of differential rotation is neglected for the purposes of the following discussion, since the attention will be limited to the plane of the ecliptic and far from the solar surface.

The Parker model predicts that for distances greater than a certain critical radius, the wind speed is almost constant. So, considering v_r = v₀, equation (1) can be integrated from r = b (close to R_⊙) to r to give

(2)

where φ₀ is the initial angle at the surface of radius b. Hence, the overall effect of the solar rotation is the bending of the plasma trajectories in the arms of an Archimedean spiral, commonly called the Interplanetary Parker Spiral, defined by equation (2).

In this description it has been assumed that the wind velocity has a uniform value at all longitudes around the Sun, but in general this is not true, and this has important consequences. For instance, let us consider in the ecliptic plane a stream of slow wind with velocity v₁ issuing from (b, φ₁) accompanied by a stream of fast wind with v₂ > v₁ from (b, φ₂). The two paths collide in a certain (r, φ) where

(3)

This is the basis for the formation of Stream Interaction Regions (SIRs) and Corotating Interaction Regions (CIRs; Gosling & Pizzo, 1999; Vršnak et al., 2017). If v₁ = 400 km s⁻¹ and v₂ = 600 km s⁻¹, with φ₂ −φ₁ = π/2, the two streams intersect at r – b ≈ 6 × 10¹³ cm = 4 AU, where the spiral is more azimuthal than radial. The fast stream runs into the sunward side of the slow stream with a relative velocity of v₂ − v₁ = 200 km s⁻¹. Magnetic fields in both streams are compressed and there two shock waves are formed, both forward and backward, away from the contact surface, together with long-lasting shock discontinuity at the boundary tracing the actual (Gosling & Pizzo, 1999).

2.2 The Weber and Davies model

While the magnetic terms in the radial equation of motion are small compared to the gravitational and pressure ones, they are instead the dominant ones in the rotational equation of motion. Weber & Davis (1967) modelled the interaction between the solar wind and a magnetic field that is uniform and radial at the solar surface, restricted their discussion to equatorial winds, by assuming B_θ = 0 = v_θ, and that variations along φ may be negligible, either for all φ, or in small regions (flux tubes) where ∂_φ B_φ and ∂_φ v_φ can be neglected. Hence the null divergence of the magnetic field ∇ · B = 0 implies that

(4)

where B_b is the field strength at the distance b where the field is assumed to be purely radial. Moreover, for a stationary magnetic field the induction equation in ideal MHD states that ∇ × (v × B) = 0, of which the θ component is neglected for an equatorial wind. The azimuthal component can be integrated considering that the magnetic field is radial at the corotating surface r = b, so that the azimuthal component of the magnetic field is

(5)

The main difference in the solar wind behaviour between the Weber & Davis (1967) model and the simpler isothermal, spherically symmetric, unmagnetized Parker (1958) treatment is the introduction of two rotational regimes: while close to the Sun the rotation is almost rigid due to the stronger magnetic pressure, far from it the angular momentum is conserved with azimuthal speed dropping as 1/r. The behaviour of radial speed also reflects the influence of the magnetic field, with the flow attaining the speed of slow, Alfvénic and fast-mode waves as it moves away from the Sun. The theoretical model described above has been used in this work to start from in situ measurements acquired at the heliocentric distance of 1 AU, and trace different solar wind flowlines back to the distance of 0.1 AU. The in situ data that have been employed here are described in the next section.

3 Input in situ data

We start here with a quick recap of major in situ observatories dedicated to the study of the Sun and the interplanetary medium (thus excluding spacecraft orbiting inside the Earth’s magnetosphere, and excluding just for brevity other existing spacecraft orbiting other planets).

3.1 Available in situ observatories

Due to its proximity, the Lagrangian point L1 is a privileged site of observation for the solar wind that will immediately impact Earth’s magnetosphere. Currently, L1 hosts four operational spacecrafts dedicated to the investigation of the solar plasma: the Advanced Composition Explorer (ACE; Stone et al., 1998), the Solar and Heliospheric Observatory (SOHO; Domingo et al., 1995), the Global Geospace Science Wind (GGS Wind; Harten & Clark, 1995) and the Deep Space Climate Observatory (DSCOVR; Marshak et al., 2018). Data from the latter two missions were used in this work.

Unlike the spacecrafts previously mentioned, the NASA operated Solar Terrestrial Relations Observatory (STEREO) is not orbiting the Sun at L1. It is made of two twin probes, STEREO-A and B, which departed from Earth in 2006 and now are in heliocentric orbits around the Sun. STEREO-A’s orbit is inside Earth’s one (~347 days for a revolution), while STEREO-B orbit is outside the planet’s one (one revolution in ~387 days). Since October 1, 2014 communications with STEREO-B were lost. All STEREO scientific instrumentation has been specifically designed for space-weather purposes.

The above summary of major in situ observatories placed along the Earth’s orbit is important also to better outline that, even if in the first implementation of RIMAP presented here we used for each simulation run only the data acquired by a single spacecraft. In future implementations of RIMAP we plan to base the simulations on data acquired at the same time by multiple spacecraft; this will allow us to relax the hypothesis made here of stationarity of interplanetary plasma during one full solar rotation.

3.2 Preparing data for ingestion in the model

The density and speed measurements of the solar wind, as well as measurements of the components of the magnetic field carried by it, are collected by the above spacecrafts during their heliocentric orbit at regular time intervals. For the purposes of the work presented here, using the equatorial synodic period of the Sun T_⊙ = 26.5 days, a simple conversion is carried out on these time intervals t transforming them into longitudinal positions φ in a chosen heliocentric coordinate system, so as to transform the temporal profiles of the quantities q(t) into azimuthal profiles q(t). The above solar period T_⊙ was obtained by averaging estimates provided by Snodgrass & Ulrich (1990) and Wöhl et al. (2010). The first work gives an equatorial sidereal speed of 14.71°/day, while the second one estimates it as 14.499°/day: thus, a value of 14.6°/day was assumed.

The selected reference system was the Heliocentric Earth Ecliptic (HEE) system, in which the x-axis coincides with the Sun-Earth line, the z-axis is perpendicular to the ecliptic plane and the y-axis completes the triad. Another possible choice was the Heliocentric Earth Equatorial (HEEQ) system, where e_x is normal to the intersection of solar equator and solar central meridian as seen from Earth while e_z is the solar rotational axis. The choice fell on the first because of its greater immediacy in converting the measurements into data usable in the reconstruction of the spiral. The 0.99 AU and/or the 1 AU data for different selected periods were normalized as arrays {q_i} of N = 361 elements, and tabulated as

(6)

After the appropriate conversions to the cgs units system, and the subsequent normalizations, these data are ready to act as input for the backward reconstruction of in situ plasma parameters, and to run the forward MHD simulation with PLUTO (Mignone et al., 2007, 2012). The next Section 4 describes in details the setup for the simulations with PLUTO, while the back-reconstruction is described in Section 5.

4 Numerical method

We employ the PLUTO code (Mignone et al., 2007, 2012) for the solution of the time-dependent MHD equations in multiple spatial dimensions. The equations are solved in their usual conservative form,

(7)

where ρ denotes the gas density, v is the fluid velocity, p = p_g +B²/2 represents the total (gas + magnetic pressure) while

(8)

is the total energy density for a gas obeying the ideal gas law with Γ = 5/3 the specific heat ratio.

PLUTO employs a conservative finite-volume discretization in which inter-cell fluxes, computed by means of a Riemann solver, are used to update conservative variables in time. For the problem at hand, the system of equations (7) is solved using a 3D spherical coordinate system, identifying (r, θ, φ) with the respective converted coordinates in the HEE system. For conveniency, we employ a system corotating with the Sun so that the force F can be written as

(9)

where .

The code employs a secon-order Runge Kutta time stepping, linear reconstruction in primitive variables, V = {ρ, v, p, B} and the Riemann solver of Roe. More details on the numerical method in curvilinear coordinates may be found in Mignone (2014). The divergence-free constraint of magnetic field is controlled using the divergence cleaning method, see Dedner et al. (2002) (and also Mignone & Tzeferacos, 2010, for the actual implementation in the PLUTO code).

The simulation time step is subject to the Courant-Friedrichs-Lewy (CFL) condition (Courant et al., 1928),

(10)

with λ_d and Δx_d being, respectively the largest characteristic velocity of the system and the grid spacing in the direction given by d, while the maximum is taken over the whole computational domain. In the equation above, N_d represents the number of spatial dimensions while C_a(1/N_d) is the Courant number.

We point out that our approach solves the full time-dependent MHD equations in 3D dimensions and, as such, it is not limited to the restrictions imposed by the 1D upwind method of Riley et al. (2011), which seeks for stationary solutions neglect pressure gradient, gravity terms and magnetic fields.

4.1 Grid layout and boundary conditions

Our computational grid extends from r_e = 21.5 R_⊙ (≈0.1 AU) up to 236.5 R_⊙ (≈1.1 AU) in the radial direction with N_r = 384 uniformly spaced zones and covers the entire azimuthal interval, φ ∈ [−π, π].

Due to the previously mentioned lack of information outside the ecliptic, the latitudinal aperture is restricted to 2° centered around θ =π/2. In order to have grid zones with aspect ratio ~1 in the central region of the domain (r ≈ 108 R_⊙), we adopt a uniform grid spacing on r and choose the number of zones in the φ direction as N_φ = 4N_r. Correspondigly, only few zones are used in the latitudinal coordinate, N_θ = 9 (an odd number was chosen in order to have a set of cells centered on the ecliptic plane).

The boundary conditions are prescribed as follows. At the inner radial boundary we impose physical parameters of the Parker spiral (reconstructing them from the in situ data) while outflow (zero-gradient) boundary conditions are imposed at the outer boundary. At the boundaries in θ we impose a reflecting boundary condition while periodicity holds in the azimuthal direction.

5 Back-reconstruction of the inner boundary

5.1 Going from 1 to 0.1 AU

The solar wind plasma measured at ~1 AU = r_E by the spacecrafts described in Section 3.1 arrived there by propagating almost radially and forming the interplanetary fieldlines of the Parker spiral. Assuming a certain stationarity for the overall plasma flow in the ecliptic plane, thus neglecting impulsive solar disturbances, turbulence, resistive effects and interactions between adjacent flow tubes, it is possible to trace back the spiraling paths followed by the plasma flow lines. Because in the back-reconstruction it will not be possible to automatically separate signatures of transient phenomena from the steady background wind on the measured time-series, the model will also back-project transients as steady features (and then propagate them forward eventually with MHD simulation). To limit the possible consequences of this effect, simulations presented here in the first implementation of RIMAP have been restricted to periods of solar activity minima. In the future developments of this work, more active periods will be analyzed and the consequences of transient phenomena will be further investigated.

Within these assumptions, one can back-reconstruct in the inner heliosphere a map of the plasma parameters distribution that led to the profiles recorded by the instruments on board the probes.

To this end, the paths traced by Weber and Davis (Sect. 2.2) are followed, by assuming that longitudinal invariance holds for each ith flow tube, and considering the section of these centered around the angular positions φ_i. The radial velocity trend (theoretically described in Sect. 2.2) is then exploited as follows. Since for r ≫ r_c, r_A (the critical radius from Parker model and the Alfvén radius) both the solar gravity and magnetic field are too tenuous to influence plasma motion, its speed remains constant. Typical values for the Parker and the Alfvén critical points may vary depending on solar temperatures and coronal fields strength, but they are always estimated way below 20 solar radii. Therefore the reconstruction can safely be extended backwards to 0.1 AU ≈ 21.5 R_⊙ = r_e, considering within this domain (at least for the ecliptic plane) constant radial speed in each streamline:

(11)

and since v_yi,v_zi ≪ v_xi ∀i it is possible to neglect them and to consider the total radial speed equal to the speed on the x axis of the HEE system.

In situ data of equation (6) can be thought of as a ring of N cells marked with the label φ_i that contain the respective values for the parameters of speed, density and magnetic field. Since the longitudinal profiles of this ring of radius r_E = 1 AU are known, those of the inner rings remain to be determined back to the radius r_e = 0.1 AU. This is easily achieved by inverting the Parker Spiral (Eq. (2)) for each cell (r_E; φ_i), so as to obtain their respective cells of origin (r_e; φ′_i), where

(12)

Since the back-traced footpoints given by the previous expression will not be evenly spaced, we use linear interpolation in order to remap the corresponding values to the code azimuthal grid. The radial speed is assumed here to be constant v′_r,i r_e = v_r,i r_E, hence the values for numerical densities are directly obtained from the continuity equation

(13)

The pressure is prescribed directly from the inner density using the adiabatic equation of state. Further details about the temperature values will be discussed in Section 6.

The treatment of the magnetic field is similar.

Radial and azimuthal components are scaled, according to equations (4) and (5) for the equatorial magnetized wind, as

(14)

However, according to the Weber & Davis (1967) model for r much greater than the Alfvén critical radius the longitudinal speed v_φ decreases as 1/r, since the field is not intense enough to make the solar wind co-rotate. Therefore, at 1 AU the azimuthal component of the magnetic field B_φ,i should already be

(15)

while at the Alfvénic radius there must be co-rotation, v′_φ,i (r_A) = Ω_⊙r_A, and a purely radial field B′_φ,i ≡ 0. In the theoretical model there is no longitudinal variability, so r_A is the same for each streamline. Here, on the other hand, there are different values of radial field and density at each φ_i, and therefore different Alfvén velocities v_A,i and radii r_A,i should be defined. To avoid this and to use the in situ data of azimuthal component, it was chosen the following formulation of B_φ,i(r_e) that includes B_φ,i(r_E):

(16)

where b ≈ 2.5 R_⊙ is the source surface radius at which B_φ = 0.

For simplicity reasons, the latitudinal component B_θ was set to zero in the internal boundary in order to leave the two-dimensional Weber & Davis (1967) treatment as unaltered as possible, as well as to not have to deal with the complete and still mostly unknown 3D trend of the heliospheric magnetic field. As will be illustrated in the Results section, this choice has clearly consequences on the output latitudinal fields at 1 AU. In this way, by following the streamlines backwards, the plasma conditions at the inner boundary of the simulation (0.1 AU) have been reconstructed as schematically shown in Figure 3.

Fig. 3 Cells of the outer ring are mapped into cells of the inner ring according to equation (12). It may be necessary to interpolate between the mapped φ′i and the equidistant grid to have cell-centered parameter values. Then the quantities q′i of the inner cells are scaled according to the theoretical models.

This is the general roadmap to follow in the reconstruction of the internal boundary. However, in the procedure described so far a caveat has not been mentioned, which will be illustrated in the following subsection.

5.2 Removal of crossing wind streams

As described by equation (2), the curvature of an arm of the spiral constituting a streamline is inversely proportional to the speed of its flow. Because of this, two flows departed from (r; φ_i) and (r; φ_j) with speed v_i and v_j will necessarily cross in (R_x; φ_x), where the radius R_X is

(17)

as already described by equation (3a). This phenomenon is responsible for the formation of Corotating Interaction Regions (CIRs), solar wind compression and decompression regions in the Parker spiral caused by the fast streams impacting on the slower ones. It is clear that, by interacting, beyond R_X the two streams will have different properties than those of departure. It would then be unrealistic to trace back the arms of the spiral even where they intersect with each other within the radial domain [r_e; r_E], as shown in Figure 4.

Fig. 4 Having two different values of speed that do not satisfy condition given by equation (18), the streamlines of cells i and j cross at (Rx; φx) before being mapped in i′ and j′. One of them must be discarded, then the azimuthal profiles are refilled by interpolating the adiacent non-crossing streams.

It was thus decided here to remedy to this incongruity by calculating the crossing radius for each pair of streamlines (i, j) and by verifying that

(18)

for each pair. In the pairs of streamlines that do not satisfy the above relation the data relative to the streamline with the lowest speed is discarded (being the one that most likely would cross other streamlines). The corresponding “emptied” cells of the inner boundary are then filled by interpolating between adjacent non-crossing cells, until the condition given by equation (18) is satisfied by all pairs. In principle, the opposite selection criterion (i.e. to remove among the crossing streamlines the one with the highest speed) would also be acceptable. Nevertheless, we decided to keep the information about the high-speed streams because these are more of interest for the possible space weather effects when impacting on planetary magnetospheres. Moreover, this criterion allows to keep the information about any possible transient phenomenon propagating faster than the background solar wind speed.

The number of removed streamlines clearly varies depending on the specific configuration analyzed. Periods associated with solar activity minima generally have more uniform radial velocity profiles and will be subject to fewer removals (about 5–10% of total streamlines). On the other hand, solar maxima and periods containing eruptive events present more complex configurations, with a much greater number of possible crossings (and consequently even 40% of the cells may have to be removed). Finally the inner boundary under the stationary condition of the ambient solar wind parameters is ready to be used for simulations.

6 RIMAP general pipeline description

As described in Section 3.1, the in situ data were first converted into azimuthal profiles of ρ, v and B by assuming constant angular speed. For the first implementation described here of the model we selected three time intervals in the data, corresponding to the solar activity minimum of March 2009, to the minimum of January 2018 and to the solar activity maximum of March 2000.

A simulation run then proceeds as follows. After setting up the PLUTO code as previously described in Section 4, the code accesses the in situ data files, reading them and storing them in a 361 × 8 elements matrix where each row represents a point of the 1 AU azimuthal profiles. At that point, the Parker spiral back-reconstruction is performed, as described in Section 5, to map back to 0.1 AU the values of each ith point, and by eliminating the pair of intersecting streamlines as explained in Section 5.2.

After that, a data conversion is applied: all the physical quantities are first converted into cgs units, then normalized, and finally scaled following the assumptions on the conservation laws described in Section 5. At the last, the new reconstructed profiles are made periodic to facilitate the subsequent interpolation between them and the grid of the internal radial boundary of PLUTO, thus making each data point to fall at the center of the corresponding ghost cell.

Following these preliminary steps, the code sets in motion the apparatus of numerical techniques described in Section 4 to make the input data evolve through the magnetohydrodynamic equations. Different configurations have been tested for the initial condition of the computational domain without finding particular variations either in the final results of the simulations or in the time in which they reach steady state. In the runs presented in Section 7, we selected spherically symmetric profiles for density, pressure and magnetic field components (taken from Perrone et al., 2019) and a constant value of radial speed (400 km s⁻¹) for each cell.

At the beginning of the time evolution, each simulation crosses a sequence of transient states, characterized by the natural expansion of the solar wind plasma filling the available volume and ejected at the inner boundary with the parameters provided by the in situ data. In this first implementation described here, ideal MHD conditions were assumed hence by neglecting turbulence and resistive effects. Albeit our computations rely on the numerical solution of the ideal MHD equations, some level of numerical resistivity is inherited from the discretization process. The amount of numerical resistivity, however, does not seem to affect our results which are consistently reproduced also for smaller grid sizes.

A stationary state is reached after an initial transient period. In particular, the simulations stop at t = 500 t₀ ≈ 10 days, when the system has already reached a steady state (visually determinable when the ecliptic maps, as well as the individual quantity profiles, have stopped time evolution). A scheme summarizing the main steps in the RIMAP pipeline is shown in Figure 5, also for direct comparison with the scheme relative to EUHFORIA shown in Figure 2. The scheme also includes the final verification (described later on) based on a comparison between the input in situ data at 1 AU and the output plasma parameters from RIMAP. The results from RIMAP simulations are presented in details in the following sections.

Fig. 5 A schematic of the steps involved in RIMAP for the reconstruction of heliospheric plasma conditions starting from in situ measurements.

7 Results from RIMAP

In this section the simulation results are presented in two different sub-section, corresponding to input data taken during two different periods: at the minimum of solar activity cycle 23 (March 2009) and the minimum of solar activity cycle 24 (January 2018).

7.1 March 2009 solar minimum

The input data used in this subsection were collected in 2009 by GGS Wind from L1 point (at heliocentric distance of 0.99 AU = 212.85 R_⊙) between March 3, 08:02:09.05 UT and March 29, 20:02:09.11 UT.

This interval corresponds to a total time of about 26.5 days. Due to the conversion described in Section 3, time-direction in the azimuthal comparisons of the following subsections goes from right to left.

7.1.1 Distributions on the ecliptic plane

Figure 6 shows the resulting 2D distributions on the ecliptic plane of plasma density (top panel) and radial velocity (bottom panel). The first important verification for the model results is the reproduction of the expected physical trends of the quantities as a function of heliocentric distance. By averaging over all longitudes φ we found that the resulting plasma density drops like 1/r², as one would expect from simple mass flux conservation.

Fig. 6 Pseudocolor plots of the plasma number density (left, Log color scale) and radial speed (right, linear color scale) on the ecliptic plane as obtained with RIMAP from the March 2009 in situ observations. In the flow below r ~ 50 R⊙ it is possible to notice some deviations: this could be due to inaccuracies in the interpolation carried out by the visualization software between the spherical grid of RIMAP and the Cartesian one of the image.

7.1.2 Radial trends

The corresponding averaged radial trend for the flow speed v_r is shown in Figure 7 normalized to the average value v₀. Since a radial gradient of density, and therefore of pressure, is present, following from the continuity equation one should observe a progressive acceleration as the solar wind is blown away, as observed. In any case, at these distances the wind almost reached its final speed, and the variations shown in Figure 7 are small with respect to the average wind speed v₀, in general less that ~10%.

Fig. 7 The mean trend for the radial speed, as derived by averaging at each altitude over all longitudes; a slight acceleration is still present, mainly due to the thermal pressure gradient.

A very interesting parameter is the radial trend of the average plasma (proton) temperature T_p. In the literature, different trends (and therefore different values) of temperature were tried for the inner boundary of this work: the one that seems to have given the best overall agreement with the azimuthal profiles described in the next subsection was taken from Perrone et al. (2019). In that work, the total proton temperature is estimated to decrease with distance r more slowly than simple adiabatic trend, being

(19)

and thus giving at the inner radius of 21.5 R_⊙ = r_e a value of T_p ≈ 1.5 MK. However, as it is shown in Figure 8, the radial temperature trend resulting with RIMAP drops faster, approximately as r^−1.2. Lowering the polytropic index could lead in a slower trend, which in turn would increase the wind acceleration and worsen the match of the 1 AU plasma parameters described in Sections 7.1.3 and 7.2.3. This has been avoided, since the focus of this work is precisely the reconstruction of those parameters. Nevertheless, considering the simplified description of plasma assumed in the first implementation of our model, we consider our results in satisfactory agreement with observations.

Fig. 8 In solid red, the mean radial behaviour for the plasma temperature, as derived by averaging at each altitude over all longitudes. Overplotted in dashed black is the reference trend TE · (r/rE)−1.2, while in dashed-dotted black TE · (r/rE)−0.9, where TE = 1.9 × 105 K from Perrone et al. (2019).

7.1.3 Comparison between input and output data

One way to verify the implementation of the scheme briefly illustrated in Figure 5, is to perform a final comparison between the input in situ measurements acquired at 1 AU, and the output plasma parameters at the same distance obtained from the simulation with RIMAP. A good agreement between input and output parameters would indicate that this reconstruction as a viable approach in determining the ambient conditions of the heliospheric plasma. Moreover, because the reconstruction is performed here by using data acquired by a single spacecraft, this corresponds to assuming stationarity over the time interval required for one full solar rotation. The reliability of this assumption can be tested by making a comparison between in situ data acquired by spacecraft located at different longitudes along the Earth orbit, as it is done here in this section. Figure 9 shows the comparisons between input and output plasma densities and speeds. Although the structures carried by the solar wind show discrepancies with the input measurements at high frequencies (small scales), at low frequencies (larger scales) a greater agreement can be noted, particularly for the radial velocity. The large amount of short fluctuations clearly shows how strong the longitudinal variability along the ecliptic plane is. Nevertheless, the general azimuthal behaviour of these quantities is well reproduced, even with a relatively small number of points on φ. In particular, the root mean square error (RMSE) between input in situ measurements and output plasma parameters is about Δn = 4 cm⁻³ for the number density and Δv = 50 km s⁻¹ for the wind speed.

Fig. 9 March 2009: the output azimuthal profiles (red) of plasma density (top) and radial speed (bottom) at 1 AU plotted against input in-situ measurements (black). Longitude intervals corresponding to SIRs/CIRs identified with STEREO data (Jian et al., 2019) are superimposed as vertical dashed lines. Time direction runs from right to left.

The larger disagreements between the input and output curves are mainly related with the crossing-streamlines removal procedure discussed in Section 5. In fact, as it is shown by the v_r(φ) profile in Figure 9 (bottom panel), the slower wind streams (removed for the sake of consistency in the inner boundary reconstruction) are not well reproduced. In correspondence of longitudes where this happens, the density curves in Figure 9 (top panel) also show the worse agreement between input and output quantities. As expected, the main disagreements are observed in correspondence with the transit of SIRs/CIRs. In order to show this point, plots in Figure 9 show also the longitudinal locations of these features as identified based on the automated analysis of STEREO data¹ for the considered period, with the numbering provided by the on-line catalog. The comparison shows that the larger discrepancies are located at longitudinal angles where the slow wind streams are preceding (in the sense of Parker spiral rotation) the arrival of SIRs/CIRs. It is worth noting that the passage of a SIR/CIR can correspond to a sudden increase in speed: therefore, there is a connection between the passage of these structures and the need to eliminate intersecting streams described in Section 5.2.

Figure 10 shows the comparisons between input and output radial and azimuthal magnetic field components. The third component B_θ (φ), remains at zero coherently with the value set in the internal boundary. The latitudinal field is an important source of space-weather activity for the planetary magnetosphere. However, a more complete treatment of B_θ would have been beyond the scope of this work, since it would have required the determination of a different internal boundary for each θ value. In light of these considerations it is therefore positive that (as shown in Fig. 10) the average azimuthal behaviours of the radial and azimuthal fields are reproduced with such agreement, having used only the simple, two-dimensional Weber &s Davis (1967) model in the reconstruction of the internal boundary. The resulting RMSE for the magnetic field is ΔB_r = 2.6 nT and ΔB_φ = 2.9 nT for the radial and longitudinal components, respectively.

Fig. 10 March 2009: azimuthal profiles of the radial (top) and azimuthal (bottom) components of B at 1 AU resulting from RIMAP (red) in comparison with the input in situ measurements (black). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived automatically with STEREO data (Jian et al., 2019). Time direction runs from right to left.

As mentioned above for the density and radial velocity trends, also the longitudinal profiles of the magnetic field components exhibit a strong variability at small scales (high frequencies). Although it is possible that part of these structures can be reproduced by increasing the spatial resolution of the simulation, many of them are instead attributable to intense transient phenomena (see for example the peak of B near 1 rad), whose reproduction would require an extension of the model to the not-ideal MHD regime. This clearly goes beyond the scope of this work.

7.1.4 Comparison between in situ data from different spacecraft

For the purposes of this work, it is important also to compare the ACE measurements acquired in L1 with measurements acquired by the twin STEREO probes at the same time intervals, but from different longitudes along the Earth’s orbit. As explained in Section 3.1, the two probes STEREO-A and STEREO-B are always located, respectively, ahead of and behind the Earth along its orbit around the Sun. Therefore, a direct comparison between the quantities measured by these two probes and by the spacecraft that provided the model inputs, GGS Wind, can be also a test for the stationarity hypothesis of the solar wind plasma conditions. As a reference for the approximate position of STEREO spacecrafts during the March 2009 time interval considered in this analysis, we can assume the position measured 13.25 days after the beginning of observations. In particular, considering the reference date of March 16, 14:02:09.05 UT, the corresponding longitudinal angles between different spacecraft were at that time Δφ_STA-GGS = 44.639° and Δφ_STB-GGS = 47.512°, respectively between STEREO-A and GGS, and between STEREO-B and GGS spacecraft. These angles can be converted into time shifts Δt, by assuming again the same rotational period of the Sun by 26.50 days. The resulting time shifts are Δt_STA-GGS = 3.286 days and Δt_STA-GGS = 3.498 days, respectively between STEREO-A and GGS, and between STEREO-B and GGS spacecraft. In particular, the data acquired by GGS Wind in L1 have been compared with data acquired by STEREO-A 3.286 days later (i.e. between March 6, 14:53:59.05 UT and April 2, 02:53:59.11 UT), and by STEREO-B 3.498 days before (i.e. between February 27, 20:55:34.05 UT and March 26, 08:55:34.11 UT).

This comparison is shown in Figure 11 for the evolution of total density (top) and outflow speed (bottom). In particular, the comparison between the longitudinal distributions of solar wind speeds v_r measured by STEREO-A and -B and by the GGS Wind spacecraft (Fig. 11, bottom panel) shows a quite good agreement, with large scale fast and slow wind regions observed approximately at the same longitudinal location (as dragged by the solar rotation) by the three spacecraft at three different times. Again, while at small scales there are significant discrepancies due to transient phenomena (like the two density peaks of ~40 cm⁻³ associated with the transit of SIRs/CIRs), at large scales there is overall a good agreement between the measured profiles. This means that the fast and slow wind streams at solar minimum conditions are surviving for many days, making the assumption of stationarity performed to run our simulation quite realistic.

Fig. 11 March 2009: comparison between the azimuthal profiles of number density (top) and radial speed (bottom) as measured by the twin STEREO probes (in orange STEREO-A, in blue STEREO-B) and by GGS Wind (black line). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time runs from right to left.

The comparison between the longitudinal distributions of densities measured by the three spacecraft (Fig. 11, top panel) shows in general a quite good agreement as well, considering the large scale (low frequency) structures dragged by the solar wind expansion. On the other hand, at higher frequencies the three curves show significant disagreements, in particular with the formation of a narrow density peak observed by STEREO-B and not observed 3.498 days later (dragged by solar rotation) by GGS Wind. A similar peak is observed also by STEREO-A, but located in a totally different longitude, and not observed by GGS Wind 3.286 days before. This suggests that these two major density peaks (as well many as many other peaks) are short-duration transient phenomena with lifetime smaller than three days: similar transient phenomena cannot be reproduced by our model.

7.2 January 2018 solar minimum

January 2018 offers another quiet, ambient configuration of streamlines to test the Parker Spiral reconstruction procedures presented here. The measurements relative to 2018 analyzed here were collected by GGS Wind in L1 between January 3, 06:10:45.63 UT and January 29, 18:10:45.69 UT, to cover again an entire solar rotation.

7.2.1 Distributions on the ecliptic plane

The simulated map values of number density and radial speed on the ecliptic plane, obtained by evolving the in situ data using PLUTO, are shown in Figure 12. Unlike the more uniform configuration reconstructed from the March 2009 dataset (Fig. 6), here more higher-speed streams (at ~550 km s⁻¹) are visible, corresponding to the four peaks of Figure 15.

Fig. 12 Pseudocolor plots of the plasma number density (left, Log color scale) and radial speed (right, linear color scale) on the ecliptic plane as obtained with RIMAP from the January 2018 in situ observations.

7.2.2 Radial trends

The radial trends averaged over all longitudes φ for the solar wind radial outflow speed v_r and plasma temperature T are shown in Figures 13 and 14, respectively. None of them shows significant differences with the behaviours obtained from the March 2009 dataset (Figs. 7 and 8). The average radial velocity trend exhibits between 20 and 50 solar radii the characteristic steeper acceleration due (mainly) to the thermal pressure gradient, before settling down to an average nearly constant speed of 440 km s⁻¹. This is even in better agreement with the expected Parker-like behaviour with respect to results obtained for the March 2009 simulation (Fig. 7).

Fig. 13 Radial trend of vr averaged on all longitudes φ as obtained from the January 2018 dataset. The radial speed shows again a slight acceleration mainly due to the thermal pressure gradient, as in the March 2009 case.

Fig. 14 Radial trends of T averaged on all longitudes φ as obtained from the January 2018 dataset. The plasma temperature goes like r−1.2, as in the March 2009 case. In dashed-dotted black line, the Perrone r−0.9 trend is plotted.

Following the trend obtained in the previous case, here the plasma temperature averaged over all the wind streamlines decays again in a way similar to the adiabatic plasma, as 1/r^−1.2. For comments on this behaviour, see the previous section.

7.2.3 Comparison between input and output data

Figures 15 and 16 show a comparison between the azimuthal variations at 1 AU of the output (solid red) and input (dashed black) plasma parameters for the 2018 case. The density longitudinal profile measured by GGS Wind (Fig. 15, top panel) exhibits two relatively high peaks near φ = 0.25 rad and φ = 2 rad, respectively at (approximately) 30 cm⁻³ and 50 cm⁻³. By looking at other quantity profiles it appears that the highest peak (at 2 rad) is associated with an intense transient event similar to those recorded by STEREO in March 2009. Given the disturbance caused by this peak, the low-frequencies plasma density behaviour is nevertheless reproduced by the simulation. Also for this period we identified the SIRs/CIRs provided by the STEREO catalog (Jian et al., 2019) and plotted their longitudinal locations in Figures 15 and 16. The comparison shows that in this case larger discrepancies are located at longitudinal angles where the slow wind streams are following (in the sense of Parker spiral rotation) the arrival of SIRs/CIRs. The RMSE between the two curves is around Δn = 7.0 cm⁻³.

Fig. 15 January 2018: azimuthal profiles of the plasma density (top) and radial speed (bottom) at 1 AU resulting from RIMAP (red) in comparison with the input in situ measurements (black). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time direction goes from right to left.

Fig. 16 January 2018: azimuthal profiles of the radial (top) and azimuthal (bottom) components of B at 1 AU resulting from RIMAP (in red) in comparison with the input in situ measurements (black line). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time direction goes from right to left.

The bottom panel of Figure 15 shows the comparison between the input and output radial velocities, demonstrating how the simulation is capable of reproducing the four high speed wind streams. Note how, due again to the crossing streamlines removal, lowest speed values are not well reproduced by the model, as already noticed for 2009 simulations. This is particularly noticeable between 2 and 2.5 rad, just before the 53rd SIR of 2017: streamlines in this region have much lower speeds than those in the surrounding areas, so they intersect them in the back-projection and must therefore be removed. The result of the subsequent interpolation is a bump above 400 km s⁻¹. Although this can and should be improved, we remind that, for the purposes of Space Weather forecasting it is more important to correctly reproduce the high-speed solar wind streams, more than the low speed ones that are associated with calmer space-weather conditions (Kamide & Maltsev, 2007; Badruddin & Singh, 2009; Gerontidou et al., 2018), even if a reliable prediction of the formation of Corotating Interaction Regions requires a correct determination of the longitudinal distribution of both high- and low-speed streams. RMSE is here around Δv_r = 70 km s⁻¹.

The φ-profiles of radial and longitudinal components of the interplanetary magnetic field at 1 AU are shown in Figure 16, while (as for the March 2009 case) B_θ component remains zero across the entire computational domain. The agreement between input and output B_r and B_φ components is slightly worse than the 2009 case, with the RMSE being ΔB_r = 3.5 nT and ΔB_φ = 3.64 for the radial and azimuthal components, respectively. This is probably due to the greater intensity of the variations of their azimuthal profiles in this dataset, with deviations of even 14 nT in the space of a few degrees. Given the irregularities of the longitudinal profiles and the fact that PLUTO tries to preserve the solenoidality of a field of which a component (B_θ) has been suppressed, the general reproduction of the trend at large scales is to be considered a good result. Further progresses on this front will be achieved by using a staggered formulation of the field components in the internal radial boundary, allowing the use of solenoidality control via the constrained transport algorithm (Balsara & Spicer, 1999; Gardiner & Stone, 2005).

7.2.4 Comparison with the ENLIL forecast

Since it is possible to access the archive of forecasts provided by the ENLIL code (earliest available simulations date back to November 19, 2013), it may be interesting also to compare the ENLIL with RIMAP reconstructions of the inner Interplanetary Parker Spiral. To perform this comparison, the density and velocity maps obtained with RIMAP for 2018 (Fig.12) have been replotted by using the same color scales used for ENLIL simulations. Figure 17 (left panels) shows the forecasts for January 16, 2018 as it appeared on the NOAA/NWS website, as well as the results from RIMAP. For the comparison, please notice that while the ENLIL computational domain extends up to the heliocentric distance of 2 AU, the RIMAP domain stops at 1.1 AU. Figure 18 shows the azimuthal profiles of solar wind density (top) and radial speed (bottom) at 1 AU simulated by ENLIL and RIMAP, overplotted with the measurements carried by GSS Wind.

Fig. 17 Visual comparison between the ENLIL forecast of the plasma conditions on the ecliptic plane of January 16, 2018 (on the left) and the simulation made with RIMAP starting from data acquired at 1 AU (on the right). Notice that the ENLIL domain extends nearly up to 2.0 AU, while the RIMAP one extends to 1.1 AU (≃236.5 R⊙).

Fig. 18 Comparison between the azimuthal profiles of number density (top) and radial speed (bottom) as measured at 1 AU by GGS Wind (black line) and simulated by WSA-ENLIL (blue) and RIMAP (red). Time direction goes from right to left. The RMSE between ENLIL forecasts and Wind measurements are Δn ≃ 6 cm−3 and Δv = 110 km s−1.

The comparison shows that, as one could expect, the general configuration of the plasma streamlines is similar at large scales, but not entirely identical. As also shown in Figure 2, the ENLIL simulations (similar to EUHFORIA) take their internal boundary from complex extrapolations starting from ground-based observation of the photospheric fields, followed by multiple assumptions and semi-empirical relationships, making the prediction at 1 AU challenging. In contrast, the reconstruction carried out by RIMAP is entirely based on direct observations of the plasma parameters at 1 AU. This means that errors are due only to transients that could temporarily weaken the validity of the stationary hypothesis (which we have seen in the previous section to be generally valid within 1 week, at least in solar minima) and to the necessary suppression of some crossing streamlines (as explained in Sect. 5).

The second thing notable is the greater number of individual streamlines discernible in the reconstruction calculated with RIMAP: this is due to the fact that, with the exception of the above-mentioned removals, each point of the in situ data profile is mapped into the internal boundary and from there propagated outward. Therefore, the code naturally maintains the azimuthal variability present in the parameters of the interplanetary plasma, which allows to explain the higher number of high-density streams found in our reconstruction, a variability that the WSA models struggle more to reproduce. On the other hand, ENLIL and EUHFORIA models have information to construct an internal boundary that varies with solar latitude θ through [0.4π; 0.6π] and therefore, as was shown in Figure 17, are able to provide forecasts about v_θ and B_θ, while the model described here focuses only on the ecliptic plane components.

8 Conclusions and future perspectives

In Space Weather research a very important step is the reconstruction of the Parker spiral of interplanetary magnetic field, because physical conditions of the interplanetary plasma dictate the propagation of various solar disturbances (Coronal Mass Ejections, Solar Energetic Particles, Corotating Interaction Regions) possibly affecting the Earth and other planets. In the literature, the usual approach to face this problem (e.g. ENLIL or EUHFORIA models) is to start from photospheric magnetic field measurements, and then apply different methods to reconstruct the entire solar corona and heliosphere, and determine the distribution of different solar wind streams (Pizzo et al., 2011; Pomoell & Poedts, 2018).

In this work we presented a different approach to this problem, that we called Reverse In situ data and MHD APproach – RIMAP. With this approach, the starting point are directly the in situ data acquired by spacecraft running around in the heliosphere, instead of photospheric field measurements. The in situ data acquired at 1 AU are employed here to remap, by applying the Weber & Davis (1967) solar wind theoretical description, the plasma parameters at the inner boundary located at 0.1 AU. These parameters are then used to reconstruct the 2D distribution of interplanetary plasma up to 1 AU by employing the PLUTO MHD numerical code (Mignone et al., 2007, 2012). Solar and space-weather disturbances are inherently 3D structures in which all components of the magnetic field contribute to the plasma motion through the heliosphere. Therefore, adopting a 2D description is a non-trivial simplification, albeit a necessary one, since in-situ measurements outside the ecliptic plane are scarce and will remain so in the near future.

The aim of the RIMAP model is to reconstruct the plasma ambient conditions on the ecliptic plane from the Sun to 1 AU and beyond, trying to maintain as faithfully as possible the solar wind smaller-scale streams observed at 1 AU which photospheric measurements-based models struggle to predict. Although the results presented here do not technically constitute a forecasts (since the data have to be first collected at 1 AU to perform the reconstruction), and the reconstructed Parker spiral cannot be considered as a snapshot of physical conditions of the interplanetary plasma at a given moment (because we assumed stationarity over one full solar rotation), the work presented here have many other interesting aspects and future developments that we summarize here.

First of all, the developed pipeline could be employed to provide forecast services, for instance once the ESA Lagrange mission will be launched to acquire the first ever measurements by a probe in orbit around L5 (Hapgood, 2017). Being located at a longitudinal distance of about 60° from the Earth, the Parker spiral reconstructed with RIMAP with L5 data will be almost the same hitting the Earth’s magnetosphere approximately 4.5 days later. A detailed reconstruction of the interplanetary plasma conditions between 0.1 and 1 AU are mandatory to have a good knowledge of the magnetic connectivity between the Sun and the Earth to forecast the propagation of SEPs, and also to predict the arrival time on Earth of interplanetary CMEs. Hence, a more detailed description of interplanetary plasma physical properties (such as the one provided by RIMAP) will improve our capabilities to forecast the arrival of disturbances originating on the Sun. Moreover, the model, now extending out to 1.1 AU, will be easily extended to farther distances from the Sun, thus allowing to predict the physical conditions of the interplanetary medium hitting the magnetospheres of outer planets in the solar system, or any other spacecraft orbiting close to the ecliptic plane. The RIMAP model would also be employed to validate future model data at L5 or L1 in real-time, thus testing and improving our forecasting capabilities.

Second, a comparison between the input in situ measurements and the output plasma parameters derived at 1 AU with the reconstruction model presented here have a relative agreement (based on the measured RMSE) comparable to what obtained with current analytical methods, and better than what currently provided by MHD numerical simulations. More in details, the current empirical models have RMSE for the wind velocity measured at 1 AU on the order of Δv = 50–100 km s⁻¹ (MacNeice et al., 2018), while MHD models reveal on average RMSE of around Δv = 100–150 km s⁻¹ in the wind speed and time shifts in the arrival of the peak speed of about one day and up to three days (Hinterreiter et al., 2019). On average the modeled solar-wind speed underestimate the in-situ measured velocity, while the modeled solar-wind density is considerably higher than the observed one (Hinterreiter et al., 2019). The disagreement is much larger for proton temperatures at 1 AU as derived with MHD codes, being systematically too small by about an order of magnitude (Gressl et al., 2014). On the other hand, even if RIMAP is not able to predict the interplanetary conditions as MHD models starting from photospheric measurements, the RMSE for the model presented here are lower (on the order of Δv = 50–70 km ⁻¹ and Δn = 4–7 cm⁻³). The radial profiles of plasma parameters (density, speed, temperature) from 0.1 to 1 AU are also in very good agreement with the current knowledge of interplanetary plasma conditions (Cranmer, 2002; Hellinger et al., 2011; Perrone et al., 2019); in future developments of this work we will investigate how modification of the polytropic index and/or addition of an energy and/or momentum source term could provide also more realistic temperature radial trends.

Moreover, different from previous similar works in the literature, the method described here for the back-mapping of the solar wind to the inner boundary was also aimed at preserving as much as possible the smaller scale features sampled at 1 AU with the in situ data. This implies that globally the 2D reconstruction performed with RIMAP is much more representative of the real physical conditions of interplanetary plasma, at least over the period of one solar rotation. The role of smaller scale features is usually neglected by current analytical and numerical reconstructions of the interplanetary medium, that provides only information on the larger scale features dragged by the solar wind. Smaller scale features can play also a role in the propagation of SEP streams (Ruffolo et al., 2003) that are affected by stochastic magnetic field components (e.g. Pei et al., 2006), and smaller scale structures may also play a role in the evolution of CMEs (see e.g. Hosteaux et al., 2018), topics which are far from being fully understood.

For these reasons, the model presented here could be also considered in the future as a test-bench to better understand for example the propagation of SEPs and/or CMEs across an interplanetary plasma which is more representative of real conditions that these disturbances encounter, in particular for what concern the short-duration density and speed anisotropies. This suggests that the model presented here could be used to correctly characterize the interplanetary plasma medium from the Sun to 1 AU, which is the starting point to provide reliable predictions not only on the arrival of wind streams on the Earth magnetosphere, but also on the propagation of solar eruptions. In a future development of this work (now in progress) a CME will be inserted at the inner boundary of the simulation, starting from kinematic and dynamic parameters derived from coronagraphic observations, to explore how the interplanetary plasma parameters derived with RIMAP (and smaller scale inhomogeneities) will affect the propagation of the disturbance.

More than that, as it is shown by Figures 2 and 5 making a comparison between the classical approach and the new approach presented here, the Parker spiral reconstruction performed with RIMAP is conceptually much more simple. This means that this method can be more easily implemented on local machines, without the need for complex and expensive computer clusters, making it a viable method that can be potentially applied and tested by any research group, without the need of significant resources. Also, the simpler conceptual approach of RIMAP allows to have a much more direct control on the results, and to test more easily the possible causes of disagreements between model plasma parameters and physical measurements.

It is also very important to consider that, in the first implementation of RIMAP Parker spiral reconstruction presented here, only the data acquired by a single spacecraft have been employed. Nevertheless, in a future development of this work, we plan to implement a remapping of the solar wind parameters at 0.1 AU by employing measurements acquired at the same time by different spacecraft mentioned in Section 3, as well as other spacecraft such as Parker Solar Probe, Solar Orbiter, Bepi Colombo, and any other spacecraft measuring density, speed, temperature, and magnetic field at different locations in the interplanetary space. This will allow to relax significantly the hypothesis of stationarity of heliospheric plasma performed in this first implementation. The use of ensemble modeling (Riley et al., 2013) or data-assimilative approach (Lang et al., 2017; Lang and Owens, 2019; Owens et al., 2019) will be also considered. In the future developments of this work, it will be relatively easy to run the RIMAP model directly with the input data acquired in situ at the inner boundary by the Parker Solar Probe spacecraft, thus skipping the problem of back-mapping form 1 AU to the Sun. Moreover, the inclusion of plasma parameters measured with remote sensing data (e.g. space based coronagraphs or heliospheric imagers) will be considered to extend the model to the third dimension.

In summary, the RIMAP method presented here will improve our theoretical knowledge of solar disturbance propagation, and will provide new methods to reconstruct and forecast the conditions in the interplanetary medium.

Acknowledgments

The authors acknowledge R. Susino for his help in the preparation of the in situ data during the first phases of this work, and in the comparison with the ENLIL model results. The authors acknowledge the use of data from the NASA STEREO and WIND spacecraft. The authors would also like to thank the reviewers of the manuscript for constructive comments and suggestions. R. Biondo acknowledges support by ASI-INAF contract nbr. 2018-30-HH.0 “Scientific support to the development of the instruments Metis, SWA DPU and STIX Phases D/E”. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Figures

	Fig. 1 Example of results from a EUHFORIA simulation run, the two panels show the distribution on the ecliptic plane (left) and at fixed longitude (right) of the solar wind speed (adapted from Pomoell & Poedts, 2018).
In the text

	Fig. 2 A schematic of the steps involved in EUHFORIA for the construction of boundary conditions for the heliospheric MHD model.
In the text

	Fig. 3 Cells of the outer ring are mapped into cells of the inner ring according to equation (12). It may be necessary to interpolate between the mapped φ′i and the equidistant grid to have cell-centered parameter values. Then the quantities q′i of the inner cells are scaled according to the theoretical models.
In the text

	Fig. 4 Having two different values of speed that do not satisfy condition given by equation (18), the streamlines of cells i and j cross at (Rx; φx) before being mapped in i′ and j′. One of them must be discarded, then the azimuthal profiles are refilled by interpolating the adiacent non-crossing streams.
In the text

	Fig. 5 A schematic of the steps involved in RIMAP for the reconstruction of heliospheric plasma conditions starting from in situ measurements.
In the text

Fig. 6 Pseudocolor plots of the plasma number density (left, Log color scale) and radial speed (right, linear color scale) on the ecliptic plane as obtained with RIMAP from the March 2009 in situ observations. In the flow below r ~ 50 R⊙ it is possible to notice some deviations: this could be due to inaccuracies in the interpolation carried out by the visualization software between the spherical grid of RIMAP and the Cartesian one of the image.

In the text

	Fig. 7 The mean trend for the radial speed, as derived by averaging at each altitude over all longitudes; a slight acceleration is still present, mainly due to the thermal pressure gradient.
In the text

	Fig. 8 In solid red, the mean radial behaviour for the plasma temperature, as derived by averaging at each altitude over all longitudes. Overplotted in dashed black is the reference trend TE · (r/rE)−1.2, while in dashed-dotted black TE · (r/rE)−0.9, where TE = 1.9 × 105 K from Perrone et al. (2019).
In the text

	Fig. 9 March 2009: the output azimuthal profiles (red) of plasma density (top) and radial speed (bottom) at 1 AU plotted against input in-situ measurements (black). Longitude intervals corresponding to SIRs/CIRs identified with STEREO data (Jian et al., 2019) are superimposed as vertical dashed lines. Time direction runs from right to left.
In the text

	Fig. 10 March 2009: azimuthal profiles of the radial (top) and azimuthal (bottom) components of B at 1 AU resulting from RIMAP (red) in comparison with the input in situ measurements (black). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived automatically with STEREO data (Jian et al., 2019). Time direction runs from right to left.
In the text

	Fig. 11 March 2009: comparison between the azimuthal profiles of number density (top) and radial speed (bottom) as measured by the twin STEREO probes (in orange STEREO-A, in blue STEREO-B) and by GGS Wind (black line). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time runs from right to left.
In the text

	Fig. 12 Pseudocolor plots of the plasma number density (left, Log color scale) and radial speed (right, linear color scale) on the ecliptic plane as obtained with RIMAP from the January 2018 in situ observations.
In the text

	Fig. 13 Radial trend of vr averaged on all longitudes φ as obtained from the January 2018 dataset. The radial speed shows again a slight acceleration mainly due to the thermal pressure gradient, as in the March 2009 case.
In the text

	Fig. 14 Radial trends of T averaged on all longitudes φ as obtained from the January 2018 dataset. The plasma temperature goes like r−1.2, as in the March 2009 case. In dashed-dotted black line, the Perrone r−0.9 trend is plotted.
In the text

	Fig. 15 January 2018: azimuthal profiles of the plasma density (top) and radial speed (bottom) at 1 AU resulting from RIMAP (red) in comparison with the input in situ measurements (black). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time direction goes from right to left.
In the text

	Fig. 16 January 2018: azimuthal profiles of the radial (top) and azimuthal (bottom) components of B at 1 AU resulting from RIMAP (in red) in comparison with the input in situ measurements (black line). This figure shows also the longitudinal locations of SIRs/CIRs (vertical dashed lines) as derived with STEREO data (Jian et al., 2019). Time direction goes from right to left.
In the text

	Fig. 17 Visual comparison between the ENLIL forecast of the plasma conditions on the ecliptic plane of January 16, 2018 (on the left) and the simulation made with RIMAP starting from data acquired at 1 AU (on the right). Notice that the ENLIL domain extends nearly up to 2.0 AU, while the RIMAP one extends to 1.1 AU (≃236.5 R⊙).
In the text

	Fig. 18 Comparison between the azimuthal profiles of number density (top) and radial speed (bottom) as measured at 1 AU by GGS Wind (black line) and simulated by WSA-ENLIL (blue) and RIMAP (red). Time direction goes from right to left. The RMSE between ENLIL forecasts and Wind measurements are Δn ≃ 6 cm−3 and Δv = 110 km s−1.
In the text