© J.J. González-Avilés et al., Published by EDP Sciences 2024

1 Introduction

Knowledge about the physical properties of the plasma from the solar atmosphere is crucial for understanding the solar wind (SW) dynamics and explosive events, such as coronal mass ejections (CMEs). These phenomena connect the Sun and the Earth, leading to an essential area of research known as Space Weather. For example, the CMEs cause geomagnetic storms that could affect the Earth’s magnetic field environment. These geomagnetic storms can disrupt the electric power supply, perturb navigation systems, and interrupt satellite functionality on Earth. Moreover, these events could impact the global economy, making it essential to develop novel scientific research dedicated to space weather forecasting (Schrijver et al., 2015).

In situ observations of the core properties of the SW and CMEs provide potentially useful scientific and operational information. In particular, e.g., the Solar and Heliospheric Observatory (SOHO), the Solar Terrestial Relations Observatory (STEREO), the Solar Dynamic Observatory (SDO), and the Global Oscillation Network Group (GONG) produce valuable information about the solar surface. In addition, missions such as the Advanced Composition Explorer (ACE), WIND, and the Deep Space Climate Observatory (DSCOVR) cover the near-Earth region. Finally, the Parker Solar Probe (PSP) and Solar Orbiter (SolO) observe the properties of the SW and CMEs in the inner heliosphere. These observatories provide in situ measurements of physical properties, such as velocity, proton density, temperature, and magnetic field of the SW and CMEs. However, they do not accurately estimate key forecasting parameters, such as arrival times of SW currents and CMEs near the Earth, i.e., at about 1 AU. Therefore, many current investigations employ numerical models to improve these limitations to develop a broader picture of space weather events, such as SW streams, CMEs, solar energetic particles (SEPs), Corotating Interaction Regions (CIRs), and stream interaction regions (SIRs).

There is a noticeable advance in models applied to study the propagation and dynamics of the SW and the interplanetary counterpart of the CMEs, commonly called Interplanetary Coronal Mass Ejections (ICMEs), in the inner heliosphere. Generally, we classify these models into three main categories: empirical, semi-empirical, and numerical. Empirical models use a probabilistic forecasting approach to analyze SW observations at Sun–Earth L1 Lagrangian point (see, e.g., Bussy-Virat & Ridley, 2014; Owens et al., 2017; Riley et al., 2017). Semi-empirical models employ semi-empirical relations of SW speed based on observations of coronal holes (see, e.g., Reiss et al., 2016; Arge & Pizzo, 2000). Numerical models are based on photospheric magnetograms to determine SW plasma properties in the inner heliosphere. These models numerically solve the MHD equations to describe the propagation and evolution of the SW streams and CMEs in the inner heliosphere, and among these models are, e.g., MAS (Riley et al., 2001), ENLIL (Odstrcil, 2003), Space Weather Modeling Framework (SWWF, Tóth et al., 2005), Solar–Interplanetary space–time conservation element and solution element (SIP-CESE, Feng et al., 2010),Space-weather-forecast-Usable System Anchored by Numerical Operations and Observations (SUSANOO, Shiota et al., 2014), European Heliospheric Forecasting Information Asset (EUHFORIA, Pomoell & Poedts, 2018), and more recently, SWASTi-SW (Mayank et al., 2022). However, despite all these models, there is still no model capable of capturing the fundamental parameters such as velocity, density, temperature, magnetic field polarity, and magnetic field strength of SW that fully match with in situ observations. Hence, there are still opportunities to continue exploring new models.

In this paper, we use sunRunner3D (SR3D) to interpret the global structure of the inner heliosphere in terms of SW from in situ measurements. SR3D is a community-developed open-source package and an improvement of sunRunner1D (Riley & Ben-Nun, 2022), which is a tool for exploring ICMEs evolution through the inner heliosphere, considering spherical symmetry in 1D. We plan to apply SR3D to specific scientific problems in heliophysics research, for example, to interpret the signatures of CMEs observed by various heliospheric spacecraft, including PSP and SolO. It can also be applied to explore the dynamic evolution of CME events out to 1 AU and to investigate whether CME events would have produced extreme space weather phenomena if they were directed toward the Earth. Besides, recently, in Aguilar-Rodriguez et al. (2024), SR3D has been successfully used to perform numerical global MHD simulations of SIRs and CIRS observed by PSP and STEREO-A.

In particular, for this paper, SR3D uses the boundary conditions generated by CORHEL (Riley et al., 2012; Linker et al., 2016) to drive the heliospheric model. In the inner heliosphere domain, that goes from 0.14 AU to 1.1 AU, SR3D employs the PLUTO code (Mignone et al., 2007) to compute the plasma properties of SW with the MHD approximation in the inertial and rotating frames of reference. One of the primary purposes of this paper is to simulate SW streams and produce background states in spherical coordinates in three dimensions in both inertial and rotating frames and compare with in situ measurements observed by Earth-based spacecraft (OMNI dataset) and STEREO-A (STA) and with the MAS code for a set of Carrington rotations that cover a period that lays in the declining phase of solar cycle 24.

We organize the paper as follows. First, in Section 2, we describe the CORHEL and PLUTO models. Then, in Section 3, we present the results of the numerical simulations of the relaxed, steady-state SW solutions, the comparisons of the model results with in situ measurements of OMNI and STA, the statistical analysis for a set of Carrington Rotations, and the comparison with similar existing model. Finally, we draw our conclusions in Section 4.

2 Model

SR3D is an MHD model to calculate SW stream steady-state solutions and evolving CMEs in a 3D spherical coordinate system in the inner heliosphere in the inertial and rotating frames of reference. Specifically, it contains two main parts: i) boundary conditions and ii) heliospheric model. The CORHEL model provides the boundary conditions described below, while the PLUTO code constitutes the MHD heliospheric model. The medium and long-term plans are to distribute SR3D as a community-developed open-source package.

2.1 The CORHEL model

To generate the boundary conditions, we used the CORHEL framework, which is capable of modeling the ambient solar corona and the inner heliosphere for a specific period of interest. Mainly, it derives the boundary conditions using maps of the Sun’s photospheric magnetic field derived from magnetograms. These magnetograms are obtained principally from the Helioseismic and Magnetic Imager (HMI) instrument of the SDO. Then, it runs the coronal model using the MAS code (Mikić & Linker, 1996; Mikić et al., 1999) until obtaining a relaxed state, which serves to produce the boundary conditions for the heliospheric models. CORHEL solutions are available to the community at the Community Coordinated Modelling Center (CMMC, https://ccmc.gsfc.nasa.gov) and the Predictive Science website (https://www.predsci.com).

In this paper, we take advantage of the fact that boundary conditions can be downloaded directly from https://www.predsci.com/data/runs/. There, the boundary conditions are already in a readable format for the PLUTO code, so we have all the variables, such as number density, radial velocity, radial magnetic field, and temperature, to set at the inner boundary (R _b = 0.14 AU) and the whole domain as initial conditions. We use the relaxed steady-state solution of the variables to drive the inner heliospheric MHD model, described in the following section.

This paper uses the MAS polytropic solutions obtained by solving the set of MHD equations in spherical coordinates, employing an adiabatic energy equation. To reproduce the structure of the magnetic field and generate solutions with sufficient variation in SW speed or densities, the model employs the empirically based approach called the “Distance from the Coronal Hole Boundary” (DCHB), which helps to specify the SW speed at the inner boundary of the chosen heliospheric code (Riley et al., 2001). The DCHB model assumes that the flow is fast within coronal holes (i.e., away from the boundary between open and closed magnetic field lines). The flow is slow at the boundary between open and closed field lines. Over a relatively short distance, the rise of the flow speed is smoothed to match the fast coronal hole flow (Riley et al., 2001). In summary, the MAS polytropic solutions, which incorporate the empirical DCHB model and are relaxed, serve as boundary conditions for SR3D. These polytropic solutions ensure consistency with the inner heliospheric MHD model, which utilizes the polytropic approximation with γ = 5/3 to characterize the plasma properties of the solar wind.

2.2 The inner heliosphere MHD model

The inner heliosphere model employs the PLUTO code from R _b = 0.14 AU outwards to solve the three-dimensional time-dependent MHD equations in spherical coordinates. Particularly, we adopt the ideal MHD equations written in the following dimensionless conservative form,

(1)

(2)

(3)

(4)

(5)

where ϱ is the plasma density, v represents the fluid velocity, E is the total energy density, B is the magnetic field, p _t is the total pressure (thermal+magnetic), and I is the unit matrix. In the equations, p _t = p + B ²/2, where by the ideal gas law . Here, T is the temperature of the plasma, is the particle mass specified by a mean molecular weight value μ = 0.6 for a fully ionized gas, m _H is the mass of the hydrogen atom, and k _B is the Boltzmann constant. In addition, E = p/(γ − 1) + ϱ v ²/2 + B ²/2, being γ = 5/3 the polytropic index. The source terms include the gravitational forces defined in terms of the gravitational acceleration , with G representing the gravitational constant and M _⊙ the solar mass. Besides, F contains the Coriolis and centrifugal forces that in PLUTO code are treated conservatively (Kley, 1998).

Equations (1)–(5) are solved in a 3D spherical coordinate system (r, θ, ϕ) for both rotating and inertial frames of reference. In the case of the solution in the rotating frame, we consider F and let the system rotate with constant angular velocity Ω_c = 2.8653 × 10⁻⁶ Hz, which represents the rotation rate of the solar equator. In the rotating frame, we define the azimuth component of the velocity V _ϕ = −rΩ_c, in the whole domain and at the boundary R _in = 0.14 AU. This definition helps us better match between solutions with PLUTO in the inertial and rotation frames, as we will see in the simulation results. On the other hand, in the inertial frame, we set F = 0 and account for the solar rotation by rotating the boundary conditions in ϕ (longitude) at a rate equal to Ω_c, which is similar to the method employed in CORHEL, ENLIL, and EUHFORIA. Finally, to ensure a zero electric field (v × B = 0) in the rotating frame, we introduce an azimuthal component of the magnetic field B _ϕ = −B _r sinθ(v _rot/v _r ), where v _rot = Ωr, representing the rotating speed of the inner boundary. The latter definition of B _ϕ is consistent with the requirement for a steady-state solution in the rotating frame (see, e.g., Pomoell & Poedts, 2018; Mayank et al., 2022).

The MHD domain extends from 0.14 to 1.1 AU along the radial coordinate, 0° to 180° in latitude, and 0° to 360° in longitude, with a grid resolution of 141 × 111 × 128, respectively. We chose the RK2 time-stepping algorithm for the time integration, a second-order linear reconstruction scheme in combination with minmod limiter, and the Harten–Lax–van Leer–Contact (HLLC) Riemann solver (Li, 2005). To ensure the divergence-free condition (Eq. (5)), we selected Powell’s eight-wave formulation (Powell, 1997). At the inner-radial boundary, we specify the values of the radial velocity v _r , the number density n, the temperature, and the radial component of the magnetic field B _r , given by the CORHEL coronal solution at R _in = 0.14 AU. At the outer radial boundary, we set outflow boundary conditions, which allow the fluid to leave the computational domain. Furthermore, we employ polar axis boundary conditions at both poles, which works in conjunction with the ring average technique that helps to remove time step restrictions near a singular axis, i.e., at θ = 0, π (see, e.g., Zhang et al., 2019). Finally, we set periodic boundary conditions in longitude.

3 Simulation results

We select the following CRs to validate SR3D: CR2190, CR2199, CR2202, CR2203, CR2205, CR2208, CR2209, CR2210, CR2211, CR2214, CR2215, and CR2221. These CRs are between April 2017 and August 2019, coinciding with the declining phase of cycle 24. In Figure 1, we show the synoptic maps of B _r with a resolution of 180 × 360 at about r = 30.1045 R_☉ ∼ 0.14 AU obtained with CORHEL for some representative CRs. These maps are interpolated to a coarser 3D spherical grid of PLUTO code and represent the input of the MHD model at the inner boundary (r ∼ 0.14 AU). Throughout the manuscript, we name the corresponding solutions for PLUTO in the inertial and rotating frames as SR3D-I and SR3D-R, respectively.

Figure 1 Synoptic maps of the radial magnetic field Br (nT) at the inner radial boundary (∼0.14 AU) for CR2190 (top-left), CR2202 (top-right), CR2210 (bottom-left), and CR2221 (bottom-right). For each panel, the horizontal axis is the longitude, and the vertical axis is the latitude in the Carrington coordinate system. The black line in each panel is the contour of the radial magnetic field drawn at Br = 0, which marks the location of the heliospheric current sheet.

In Figure 2, we show the results for the relaxation of the SW conditions corresponding to CR2210 using the solution of SR3D-R. For this simulation, we ran out PLUTO long enough to convect any/all transient phenomena created at t = 0 past the outer radial boundary. In Figure 2, we normalize N multiplying by the factor r ², which helps to illustrate the results better. In particular, we display equatorial and meridional cuts, observing standard features of interplanetary solutions expected for the steady state SW. We recognize the fast and slow SW streams in the radial velocity, which are discernible as low and high-density SW streams, respectively, in the normalized number density maps. Notably, in the plot of the radial velocity V _r, it is also visible that the high-speed (>700 km s⁻¹) SW dominates in the north and south poles, as shown in the meridional plane. We also observe a mix of slow and fast wind at all latitudes, which the coronal structure during the period of CR2210 could produce. In addition, the temperature cuts are also a mix of low and high-temperature regions, which are related to the presence of high and low-density SW streams there. Finally, in the radial magnetic field cuts, we identify the formation of the Parker spiral represented by the magnetic field lines colored in black.

Figure 2 Maps of simulated SW parameters in the heliospheric equatorial and meridional planes for CR2210 in the rotating frame. From the upper left to the lower right, the panels show velocity Vr (km s−1), scaled number density r2 N (cm−3), temperature T (K), and the scaled radial magnetic field r2 Br (nT), overlay with magnetic field lines (black).

In Figure 3, we present snapshots of the output variables for the relaxation, i.e., the steady state of the SW corresponding to a representative set of CRs: CR2208, CR2210, and CR2221. In Figures 3a–3c, we display the radial velocity in km s⁻¹ in the equatorial and meridional planes and a slice of constant radial distance at 1 AU, correspondingly. In Figures 3d–3f, we show the radial magnetic field in nT, the number density in cm⁻³, and the temperature in Kelvin, respectively. The number density and temperature are shown on a logarithmic scale to discern their structure more clearly. The typical SW structures are visible in the radial velocity cuts, as shown in Figures 3a–3c. Mainly, we identify steady-state slow and fast SW streams in the equatorial cuts in Figure 3a, typically denominated as CIRs. In addition, in Figure 3d, the radial magnetic field component shows a polarity change, identified with the heliospheric current sheet (HCS). Furthermore, in all panels in Figure 3e, we visualize high-density SW structures near the HCS while low-density SW streams develop in surrounding regions. This behavior is consistent with the radial velocity as schematized in all panels in Figure 3c, where slow SW streams are located along the HCS. In contrast, fast SW streams dominate at higher latitudes. Also, the fast SW streams have a low density, as expected due to the polytropic thermodynamic relation between the two variables, e.g., Riley et al. (2001). The radial slices of the temperature displayed on panels in Figure 3c show regions of lower temperatures near the HCS and higher temperatures in the surrounding regions. Furthermore, in the surrounding regions of the HCS, in the slices at 1 AU, we observe some compression and rarefaction regions that are likely formed due to the interaction of the high-speed streams (HSSs) and the low-speed streams (Gosling et al., 1972; Gosling & Pizzo, 1999).

Figure 3 Snapshots of the outputs of the inner heliospheric model for CR2208 (A), CR2210 (B), and CR2221 (C). In all the plots, panels (a), (b), and (c) show the radial velocity Vr, (d) the radial magnetic field Br, (e) the proton number density Np, and (e) the proton temperature in logarithm scale. Top-left panels are in the r–ϕ equatorial plane, top-middle panels are in the r–θ meridional plane at 0° longitude and bottom panels are in the θ–ϕ plane at r = 1 AU. In all panels (a) and (b), the black and red circles represent the projected location of Earth and STA.

3.1 Comparisons with in situ measurements of OMNI (Earth) and STA

We compare the steady-state SW solutions obtained with SR3D-I and SR3D-R with data from OMNI (Earth-based spacecraft) and STA and the MAS results in the inertial frame for all the CRs listed above. To do so, we use the PsiPy tool (Stansby & Riley, 2022), which is helpful for reading, processing, and visualizing MHD models developed by Predictive Science Inc.

In Figure 4, we show the results for comparisons between the steady-state SW solutions of SR3D-I (green curves) and SR3D-R (brown curves) and MAS (blue curves) with the OMNI in situ measurements 12 h averaged (red curves) for three representative CRs: CR2210, CR2214, and CR2221. The 12-hour average means that we make the average 6 h behind and 6 h ahead of OMNI in situ measurements. We include the results of MAS since it has been extensively compared with observations (Riley et al., 2021); therefore, MAS serves as a reference model to the SR3D solutions. For example, the comparisons between the models and the radial velocity V _r for CR2210 show regions of slow SW (∼300 km s⁻¹) in the first ten days observed by OMNI, which the three models capture. Notably, the three models capture the rise in speed observed at about 2018-11-07; however, there is a slight phase difference. In particular, the three models match the variations from slow to fast wind observed from 2018-11-10 to 2018-11-22. The solutions of SR3D-I and MAS behave similarly, but the SR3D-R results show a slight phase difference compared with the latter models. Regarding number density, the models overestimated the observed values. It is also evident that SR3D-I and MAS overlap, but SR3D-R solutions are sharper than MAS. SR3D-R captures lower densities at the sharp regions, and the phase shift is visible. For the comparisons with temperature, we see that the three models capture the regions of low temperatures; however, the models overestimate the regions with high temperatures. This behavior is consistent with the results shown in the number density. In the case of the radial magnetic field, we note that the three models only capture global variations, i.e., the changes of sign, but underestimate their amplitude. This behavior has already been reported by, for example, Riley et al. (2012), Linker et al. (2017), and it can be improved if we multiply the numerical results by a factor of about three. Finally, for CR2210, we display comparisons with the magnetic field magnitude, where it is discernible that the three models underestimate its strength.

Figure 4 Comparison of model results of SR3D-I (green curves), SR3D-R (brown curves), and MAS (blue curves) for the radial velocity Vr (km s−1), number density N (N/cm3), temperature T (MK) in logarithmic scale, radial magnetic field Br (nT), and the magnetic field magnitude |B| (nT) with OMNI in situ 12-hours averaged measurements (black curves) for CR2210 (a), CR2214 (b), and CR2221 (c).

In Figure 4b, we show the results for CR2214. In this case, the slow and fast SW streams are well captured by the three models, as shown in the radial velocity time series. However, the three models overestimate the two rises in velocity occurring at about 2019-02-21 and 2019-03-01, respectively. Again, the slight phase shift between SR3D-I and SR3D-R is visible. In addition, OMNI observations show streams of low density (∼20 cm⁻³) in most of CR2214; however, the three models overestimated it, and instead, they obtained higher density streams. The comparisons with temperature show the three models capture the global temperature variations. However, they obtain solutions with hotter regions where the observed temperature is lower and colder regions where the in situ measurements indicate high temperatures. This result is related to the density behavior since the three models overestimated low-density SW streams observed by the in-situ measurements. Generally, low-density regions are hot, so the models struggle to match them, as noticed in the temperature plots. The radial magnetic field and the magnetic field magnitude time series obtained with the three models underestimate the global variations observed by OMNI, similar to the results of CR2210.

In Figure 4c, we display the results for CR2221. In this CR, the three models captured the HSSs observed on about 2019-08-30. However, the low density that characterizes these regions was considerably overestimated by SR3D-I and MAS, while SR3D-R solutions describe lower-density HSSs, which are close to the observed values by OMNI. Furthermore, the temperature of the HSSs is colder than the one estimated by the models. Despite this mismatch, the three models can describe the increase in temperature related to the HSSs. Again, the three models underestimate the radial and magnetic fields’ magnitude.

In Figure 5, we show the results for comparisons between the model solutions and STA in situ measurements 12-hour averaged for the case of CR2190, CR2203, and CR2210. In the case of CR2190 (Fig. 5a), the three models matched the slow SW streams but underestimated the high-speed SW streams. For the number density, we note that the three models overestimated the low-density streams; however, they reproduce the rise. Again, we identify that three models underestimated the magnitude of the magnetic field. In Figure 5b, we show the comparisons for the CR2203. Here, the three models match the radial velocity of the SW streams. Again, they overestimated the low-density SW streams for the number density, but they captured the large-scale increases and decreases. In the case of the radial magnetic field and the magnetic field magnitude, the three models underestimated the strength; however, they describe the global behavior of the magnetic field. In the case of CR2210 (Fig. 5c), STA observed three CIRs as reported by Allen et al. (2021). The first CR was observed between 2018-10-26 at 06:30 and 2018-10-27 at 20:50; the second started at about 06:00 of 2018-11-02 and finalized at 18:10 of 2018-11-02; while the third one was observed between 2018-11-22 at 20:40 and 2018-11-23 at 00:35. From the three CIRs, the models match reasonably well the radial velocity and temperature for the first CIR. The match to the second CIR is acceptable; however, the three models did not match the third CIR. Notably, the three models underestimated the radial velocity and overestimated the density of the CIR.

Figure 5 Comparison of model results of SR3D-I (green curves), SR3D-R (brown curves) and MAS (blue curves) for the radial velocity Vr (km s−1), number density N (N/cm3), temperature T (MK) in logarithmic scale, radial magnetic field Br (nT), and the magnetic field magnitude |B| (nT) with STA in situ 12-hour averaged measurements (black curves) for CR2190 (a), CR2203 (b), and CR2210 (c).

In all the time series of Figures 4 and 5, we identify slight differences regarding the results obtained with SR3D-R, SR3D-I, and MAS. In particular, we note a difference in the sharpness of features and a phase shift between both solutions. These differences could be related to how we solve the MHD equations in each frame of reference. However, the differences could be related to other reasons, as discussed in more detail in Appendix A.

Furthermore, from the three model comparisons with OMNI and STA in situ measurements as described in the above paragraphs and Figures 4 and 5, we note that MAS results differ from SR3D-I and SR3D-R results. It is appropriate to compare MAS and SR3D-I solutions since both models solve the same MHD equations in the inertial frame and employ identical boundary conditions. At the same time, SR3D-R solves the MHD equations in the rotating frame, accounting for centrifugal and Coriolis forces. This introduces additional numerical diffusion, as elaborated further in Appendix A. Focusing on the comparison between SR3D-I and MAS, both models employ different numerical algorithms; for instance, MAS uses the finite differences method to discretize the equations, and SR3D-I employs the finite volume discretization and high-resolution-capturing methods to solve the MHD equations, which are the algorithms used by PLUTO. Besides, MAS evolves the magnetic field through a magnetic vector potential, which maintains the ∇ · B condition to a round-off value. In contrast, SR3D-I uses Powell’s eight-wave formulation that maintains ∇ · B to a truncation order error. Based on the differences above, the main reason for the disparities between MAS and SR3D-I solutions is mainly related to the numerical algorithms used to solve the MHD equations. Commonly, finite difference methods produce more diffusion than finite volume solutions for MHD equations. So, based on the explanation above, the differences between MAS and SR3D-I solutions relate to numerical artifacts.

3.1.1 Statistical analysis

We performed a statistical analysis to validate the quality of the three models, SR3D-I, SR3D-R, and MAS, compared to in situ measurements. In particular, we compared measurements and model solutions by estimating the Root Mean Square Error (RMSE) and the Pearson Correlation Coefficient (PCC). The RMSE represents the mean squared difference between measurements and models and is less sensitive to atypical values (see, e.g., Reiss et al., 2022). At the same time, the PCC measures the linear correlation between two data sets. These two statistical measures are defined as follows:

(6)

(7)where f _k and O _k are the kth elements of N total model values and in situ observations, respectively. The and in equation (7) represent the mean values. Despite the small shifts in time between simulation and data that could lead to misleadingly poor performance using to point-to-point comparisons, more sophisticated techniques, such as dynamic time warping (see, e.g., Owens & Nichols, 2021; Samara et al., 2022), are beyond the scope of this paper.

In Tables B1 and B2 of Appendix B, we show all the statistical results in terms of RMSE and PCC for OMNI and STA, correspondingly. To avoid over-information, we use Taylor diagrams (Taylor, 2001), combining correlation, centered RMSE difference, and variances to summarize all the statistics and discern the best results for all the SW variables and all CRs mentioned above. Note that there is a difference between RMSE shown in Tables B1 and B2 of Appendix B with the centered RMSE difference shown in the Taylor diagrams. In particular, the only difference is that before computing the RMSE, the average values of the data and model are first subtracted from the former. These diagrams are helpful to validate multiple aspects in a single diagram, and they have already been applied to the validation of ambient solar wind numerical models against in situ measurements (see, e.g., Riley et al., 2013; Owens, 2018; Kumar et al., 2020; Reiss et al., 2020). The main objective of Taylor diagrams is to recognize the geometric relationship between the correlation coefficient, the RMSE, and the amplitude of variations in the modeled and reference time series. In Figures 6a–6c, the azimuthal position indicates the PCC, the radial distance from the circle at x-axis is proportional to the centered RMSE difference, and the distance from the origin is proportional to the normalized standard deviation. Therefore, models with good agreement with in situ measurements should be close to the red star on the x-axis, indicating a similar standard deviation, a high correlation, and a low centered RMSE difference. In particular, we normalize the standard deviation in terms of the lowest deviation estimated for observations in each CR; therefore, the standard deviation values of the models can be less or greater than one in all the cases. Also, we set the observation standard value for all cases equal to one, representing the reference value. This normalization helps tighten all the values with similar scales and compares all the model results more quickly. Furthermore, we set points with three different colors corresponding to each model (SR3D-R, SR3D-I, and MAS) and labeled each point with a number from 1 to 12; these labels correspond to the 12 CRs analyzed. This way, it is easier to identify which of the three models for each CRs is closest to the observation values.

Figure 6 Taylor diagrams to models with OMNI 12-hr avg measurements for all CRs corresponding to Vr (a), Np (b), Br (c), and |B| (d).

For example, Figure 6a, we show the Taylor diagram for the radial velocity V _r corresponding to the comparisons with OMNI 12-hr averaged data. In this panel, it is discernible that the best performance of models is for CR2199 and CR2205, whose values are close to the observation data, while CR2215 achieves the worst match of the three models. Globally, from the diagram, MAS results are better than SR3D-I and SR3D-R results for matching the radial velocity. In Figure 6b, the distribution of the points for proton density N _p shows that three models barely match with in situ measurements for most of the CRs; however, MAS matches acceptable to the observation values for CR2199 and CR2205. In summary, the diagram for proton density shows a clustering of the models for most of the CRs; therefore, they are comparable. The distribution of the model results for the radial magnetic field B _r, as displayed in Figure 6c, shows that the three models are close to each other and are close to observation values; MAS shows an appropriate performance for CR2211. Finally, in Figure 6d, for the magnetic field magnitude |B|, we identify that the results are generally poor compared to the other variables. Specifically, in this panel, we see SR3D-R has a negative PCC for CR2190 and small values for CR2199, CR2208, and CR2211. Despite the latter results in Figure 6d, the three models reach comparable results; mainly, SR3D-I and MAS perform better than SR3D-R.

Similarly to the analysis performed for OMNI, we also show Taylor diagrams to summarize the statistical comparisons between models and STA-12 hrs averaged in situ measurements. In Figure 7a, we display the results for the radial velocity V _r, where the three models show an acceptable performance for most of the CRs, except for CR2209, where three models have negative PCC and high-centered RMSE difference values. At the same time, SR3D-R is closest to the observed radial velocity for CR2190. The proton number density N _p displayed in Figure 7b indicates that the three models are not close to observations. However, SR3D-R for CR2205 matches reasonably well with STA in-situ measurements. For the magnetic radial field, shown in Figure 7c, the three models cluster similarly, and they lie between 0.75 and 1.00 centered RMSE difference and between 0.6 and 0.7 PCC, which are consistent with the results obtained in Riley et al. (2015). Lastly, in Figure 7d, the Taylor diagram for the magnetic field strength shows that SR3D-R has negative PCC values for CR2199, CR2205, CR2208, CR2215, and CR2221. However, the three models reach values reasonably close to the reference magnetic field strength observed by STA.

Figure 7 Taylor diagrams to compare models with STA 12-hr avg measurements for all CRs corresponding to Vr (a), Np (b), Br (c), and |B| (d).

3.2 Comparison with similar existing model

In this work, we have also conducted a comparative study with another existing model that utilizes the PLUTO code, mainly focusing on the impact of initial boundary conditions on the final results at 1 AU. This feature involved the application of various coronal models including GONG (https://gong.nso.edu/data/magmap/) and Air Force Data Assimilative Photospheric Flux Transport (ADAPT: https://nso.edu/data/nisp-data/adapt-maps/) based Wang-Sheeley-Arge (WSA), as well as HMI-based MAS, integrated with the SWASTi solar wind module (Mayank et al., 2022). The SWASTi framework, recently developed to forecast and assess the solar wind and CME characteristics within the inner heliosphere, operates on the ideal MHD module of the PLUTO code. The framework also incorporates gravity and adopts an ideal equation of state with γ = 1.5. Notably, SWASTi offers a high degree of flexibility regarding time cadence, allowing for intervals as brief as 5 minutes or less, depending on user specifications. This flexibility is facilitated through a virtual spacecraft approach, which stores plasma properties in real-time simulations. Such a generic approach is crucial to understanding the dynamic behavior of heliospheric transients, particularly CMEs (see, e.g., Mayank et al., 2024).

In this study, we have employed the SWASTi’s default numerical setup, wherein the form of the WSA speed relation is the follows:

(8)where v ₁, v ₂, and β are independent parameters whose values are taken to be 250 km/s, 675 km/s and 1.25, respectively (Narechania et al., 2021). f_s and d are the flux tube’s areal expansion factor and the foot-point’s minimum angular separation from the coronal hole boundary. In contrast to the original WSA relation, the parameter w is not independent; it is determined by the median of the d values corresponding to the field lines that reach Earth’s location. This means that the adapted speed relation operates with fewer independent parameters, offering a more adaptable approach than the initial WSA relation. Additionally, the MHD domain range spans from 0.1 AU to 2.1 AU radially, −60° to 60° latitudinally, and 0° to 360° longitudinally, with respective grid resolutions of 150 × 120 × 360. For further details, readers may refer to Mayank et al. (2022).

To clearly understand the differences between SR3D and SWASTi, we have encapsulated the unique features of both models in Table 1. We specify the supported magnetograms and the models employed in the coronal region. Then, it is discernible that both models use different semi-empirical coronal models and magnetograms but employ PLUTO as the heliospheric model. Hence, the purpose of having both models (SWASTi and SR3D) in this section is to compare each other and explore their capabilities of capturing solar wind properties, for example, at about 1 AU.

Figure 8 illustrates the time series graphs depicting solar wind characteristics at 1 AU, utilizing different coronal models across three CR periods: CR2210, CR2214, and CR2221. We compared three distinct solutions – SR3D-I, SWASTi-GONG, and SWASTi-ADAPT at L1. These solutions differ primarily in their coronal domains, which serve as the initial boundaries for the heliospheric domain. As indicated by their names, the first two solutions derive from GONG and ADAPT magnetograms, respectively, utilizing the modified WSA relation to estimate solar wind speed at the boundary. Conversely, SR3D-I employs a physics-based MAS coronal solution grounded on HMI magnetogram data to establish the initial boundary for the heliospheric domain. We note that the three solutions significantly differ in Figure 8a–8c. Specifically, the semi-empirical (SWASTi-GONG and SWASTi-ADAPT) obtain higher solar wind speeds than the physics-based initial conditions of SR3D-I, which is similar to CR2210 (Fig. 8a) and CR2214 (Fig. 8b). However, for CR2221 (Fig. 8c), all approaches demonstrated comparable results for the solar wind speed profiles. Regarding proton number density analysis, SR3D-I obtains higher density peaks than SWASTi-GONG and SWASTi-ADAPT for the three CRs. Conversely, the temperature profiles revealed that SWASTi-ADAPT and SRD3D-I estimate higher temperatures than SWASTi-GONG, as shown in the three panels. For the radial magnetic field and the magnetic field magnitude, it is discernible that SWASTi-GONG and SWASTi-ADAPT obtain larger magnitudes than SR3D-I for the three CRs.

Figure 8 SR3D-I (green lines) and SWASTi results using default coronal models based on GONG (brown lines) and ADAPT (blue lines) magnetograms. All plots of CR2210 (a), CR2214 (b), and CR2221 (c) are of 12-hour averaged measurements at L1.

Figure A1 The picture demonstrates the differences between the solar wind solutions of the inertial and rotating frames. The differences become evident at larger radial distances.

While showing similar behavior between the SWASTi and SRD3D-I results in the level of matching number density, temperature, and magnetic field for all CR periods, there are differences between them in the setups, particularly in the location of the inner radial boundary in the heliospheric domain, the extent of boundaries in radial and latitudinal directions, grid resolution, and the value of the polytropic index. A significant distinction was noted in the relative flatness of the profiles for the magnetic field obtained by SR3D-I compared to those from SWASTi. This difference might be due to the lower numerical grid resolution in the SR3D-I solution and the underestimation of the magnetic field magnitude in the boundary conditions provided by CORHEL.

4 Conclusions

In this paper, we have used SR3D to interpret the global structure of the heliosphere from in situ measurements. This model uses boundary conditions generated by CORHEL which uses coronal conditions to generate steady-state solutions that extend up to 0.14 AU. In the inner heliosphere domain, the PLUTO code is employed to compute the plasma properties of SW with the MHD approximation of up to 1.1 AU in the inertial and rotating frames of reference. According to the solutions for the steady-state SW for various CRs, we found that SR3D can capture the global properties of the SW streams in the inner heliosphere, including the CIRs and the Parker Spiral magnetic field structure.

SR3D can capture global solar wind structures observed by OMNI and STA. The statistical analysis shows that SR3D produces acceptable values for the RMSE and PCC for the selected CRs. For most CRs, SR3D matched the number density better than MAS. The results for the radial magnetic field showed an underestimation of the strength; however, the three models showed reasonable values of RMSE and PCC. Overall, the statistical analysis indicates that SR3D is comparable to MAS in most CRs, which gives us the confidence to apply it to simulations of background SW solutions in the heliosphere. Importantly, in the Taylor diagrams, it is discernible that the cluster patterns show that different solutions for a particular CR are closer than the types of runs over different CRs, which implies consistency in the solutions of the three models. In addition, it is notable that PLUTO gets the solutions using a medium grid resolution running on a local workstation, which represents an advantage compared to the models that need supercomputers to exploit their full capabilities.

In our comparative analysis of SR3D with SWASTi, we observed that both models offer comparable accuracy in their results, even though these models have the following primary differences: SWASTi employs GONG and ADAPT magnetic field maps based on the WSA model and also incorporates gravity and adopts an ideal equation of state with γ = 1.5. While SR3D employs HMI magnetic field maps, the distance from the DCHB model ignores gravity and uses γ = 5/3 in the ideal equation of state. Nevertheless, there is room for SR3D to enhance the precision in computing the magnetic field magnitude. Our study indicates that physics-based initial boundary conditions might yield accurate results. However, semi-empirical boundary conditions can offer consistent results during specific periods. It is crucial to highlight that the computational time required for solving physics-based coronal domains is more significant than the semi-empirical approach, even at a considerably lower resolution. Therefore, the ultimate choice depends on the user’s preferences and requirements. A sensible strategy might be to cultivate an ensemble of solutions, incorporating both approaches, to optimize the potential for accurate solar wind property forecasting. Furthermore, a relevant potential improvement of our models could be creating seamless links with different coronal models to integrate a modular code capable of running WSA, DCHB, or fully time-dependent output from coronal models. This time-dependent output could propagate, for example, the CME directly from one region to another.

Additionally, SR3D describes the global structure of the SW streams observed near the Earth environment, i.e., at about 1 AU, and also features observed by STA. Hence, it is a helpful model for interpreting the global structure of the heliosphere from background SW conditions, which include CIRs and SIRs. Our medium- and long-term plans aim to share SR3D as a community-developed open-source model, allowing us to receive constant feedback from, for example, the international heliospheric community dedicated to modeling. Also, given the capabilities, SR3D could be part of CCMC as a heliospheric model in the future. Finally, we envisage SR3D as a tool for scientific understanding and forecasting of ICMEs in the inner heliosphere.

Acknowledgments

We appreciate the constructive comments and suggestions of the two anonymous referees that helped to increase the quality of the results presented in this paper. JJGA acknowledges the “Programa de Apoyo a Proyectos de Investigación e Innovación Tecnológica” (PAPIIT) IN115423 project and “Consejo Nacional de Humanidades Ciencias y Tecnologías (CONAHCYT)” 319216 project “Modalidad: Paradigmas y Controversias de la Ciencia 2022,” for the financial support of this work. PR and MBN gratefully acknowledge the support from NASA (80NSSC20K0695, 80NSSC20K1285, 80NSSC23K0258, and the Parker Solar Probe WISPR contract NNG11EK11I to NRL (under subcontract N00173-19-C-2003 to PSI). We thank Goddard Space Flight Center’s Space Physics Data Facility (https://omniweb.gsfc.nasa.gov/) for providing the solar wind data. In addition, we produce the synoptic maps of the radial magnetic field and the time series of the comparisons between models and in situ measurements with the PsiPy package (https://psipy.readthedocs.io/en/stable/), while we generate Figure 2 using VisIt (https://visit-dav.github.io/visit-website/index.html) and Figure 3 with pyPLUTO (https://github.com/coolastro/pyPLUTO) packages. We are also grateful to the developers of the PLUTO software, which provides a general and sophisticated interface for the numerical solution of mixed hyperbolic/parabolic systems of partial differential Equations (conservation laws) targeting high Mach number flows in astrophysical fluid dynamics. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Appendix A

Rotating vs. inertial frame analysis

The differences briefly mentioned in Section 3 are related to the results in the time series comparisons, showing that SR3D-R and SR3D-I are sharper than MAS in regions with higher variations in density, which is a consequence of the less diffusive numerical algorithms employed by PLUTO compared to MAS. Moreover, we highlight two primary characteristics in analyzing the simulation of solar wind in both inertial and rotating frames: the diffusion phenomenon and the presence of a phase shift. Notably, the outcomes within the rotating frame exhibited a higher diffusion level when contrasted with those in the inertial frame. Additionally, the time series plots demonstrated phase shift between rotating and inertial frames.

The peaks in solar wind properties exhibited greater sharpness in the inertial frame. This discrepancy in behavior might stem from numerical diffusion resulting from the rotation of the solar wind within the rotating frame. In the inertial frame, solar wind properties are adjusted with the appropriate angle at each time step, independently of the simulation’s numerical calculations. Conversely, the rotating frame produces the solar wind rotation due to the source terms defined by the Centrifugal and Coriolis forces. These source terms introduce extra diffusion, in contrast to rotating the inner boundary conditions in the inertial frame. The impact of this added diffusion becomes increasingly noticeable with greater radial distance from the simulation’s inner boundary. As a result, at a distance of 1 AU, discernible shifts in the profile of solar wind properties are observed, while these differences remain minimal closer to the simulation’s inner boundary.

Figure A1 illustrates the profile of solar wind density for different resolutions (rows) at varying radial distances (columns) in the equatorial plane. At 0.10 AU, the inner boundary, there is no difference between the inertial and rotating plots. However, the differences are apparent at 0.55 and 1.0 AU. The peaks in the inertial frame results are sharp, and the profile exhibits smaller patterns than the rotating frame results. Furthermore, Figure A1 also illustrates a phase shift between the two profiles. When we initially assign V _ϕ = 0 at the internal boundary R _b in the rotating frame, we achieve solutions with a noticeable phase shift relative to those in the inertial frame. This condition was mitigated by setting V _ϕ = −rΩ in the rotating frame, bringing the solutions closer to those in the inertial frame with only a minor phase shift. The latter means that if we employ a purely radial flow at the inner boundary in the inertial frame, it is necessary to define V _ϕ = −rΩ in the rotating frame to make the two models consistent. This aspect of the results has yet to be documented in global heliospheric modeling and, we believe, merits further detailed investigation.

We simulated the same solar wind condition for different resolutions to explore the reason behind the abovementioned differences further. In Figure A1, the discrepancies between the two rotating frames diminished as the resolution increased. In the rotating frame, the density peaks became sharper. For the simulation at the 4-degree resolution, the density profile appeared more diffuse than that in the inertial frame, especially in regions of low density (blue-shaded regions). At 2-degree resolution, this disparity slightly lessened, and the features became more pronounced. The plots at 1-degree resolution are markedly more similar, albeit with some differences in the phi-direction shift.

The highest degree of shift is also highlighted in Figure A1, showing 11, 9, and 8 degrees for 4, 2, and 1-degree resolutions, respectively. To further quantify the correlation between the solutions in both frames, we computed the Pearson correlation coefficient (cc) for both. They were calculated after adjusting the rotating frame by its maximum shift value to lend more significance to the cc values. We observed that the cc value increased from 0.67 to 0.76 to 0.8 with an increase in resolution. This gradual convergence of density plots in rotating and inertial frames strongly suggests that numerical diffusion is the primary factor behind the observed discrepancies. With further increased resolutions, the alignment between these plots is expected to improve even more.

The computational time required to model the same solar wind conditions in each frame differs significantly. Rotating frames demands considerably more time than inertial frames. For instance, in a simulation with 1-degree resolution, the inertial frame required approximately 4 hours, whereas the rotating frame took about 44 hours to complete on a 48-core CPU. Given that higher resolution results in the rotating frame increasingly resembling those from the inertial frame, opting for the inertial frame in medium-resolution simulations seems a prudent choice.

Appendix B

Full statistical results of the comparisons of models with in situ measurements

The full statistical results obtained in comparing SR3D-I, SR3D-R, and MAS with OMNI and STA 12-hour averaged data for all the CRs considered in this paper. Table B1 contains the statistical results represented by RMSE/PCC between SR3D-R, SR3D-I, and MAS with OMNI’s 12-hour averaged data for the radial velocity V _r in km s⁻¹, proton number density N _p in cm⁻³, proton temperature T _p in MK, radial magnetic field B _r in nT, and magnetic field magnitude |B| in nT. Table B2 contains the STA 12-hour averaged data comparisons statistical results.

References

All Tables

All Figures

Figure 1 Synoptic maps of the radial magnetic field Br (nT) at the inner radial boundary (∼0.14 AU) for CR2190 (top-left), CR2202 (top-right), CR2210 (bottom-left), and CR2221 (bottom-right). For each panel, the horizontal axis is the longitude, and the vertical axis is the latitude in the Carrington coordinate system. The black line in each panel is the contour of the radial magnetic field drawn at Br = 0, which marks the location of the heliospheric current sheet.

In the text

	Figure 2 Maps of simulated SW parameters in the heliospheric equatorial and meridional planes for CR2210 in the rotating frame. From the upper left to the lower right, the panels show velocity Vr (km s−1), scaled number density r2 N (cm−3), temperature T (K), and the scaled radial magnetic field r2 Br (nT), overlay with magnetic field lines (black).
In the text

Figure 3 Snapshots of the outputs of the inner heliospheric model for CR2208 (A), CR2210 (B), and CR2221 (C). In all the plots, panels (a), (b), and (c) show the radial velocity Vr, (d) the radial magnetic field Br, (e) the proton number density Np, and (e) the proton temperature in logarithm scale. Top-left panels are in the r–ϕ equatorial plane, top-middle panels are in the r–θ meridional plane at 0° longitude and bottom panels are in the θ–ϕ plane at r = 1 AU. In all panels (a) and (b), the black and red circles represent the projected location of Earth and STA.

In the text

	Figure 4 Comparison of model results of SR3D-I (green curves), SR3D-R (brown curves), and MAS (blue curves) for the radial velocity Vr (km s−1), number density N (N/cm3), temperature T (MK) in logarithmic scale, radial magnetic field Br (nT), and the magnetic field magnitude |B| (nT) with OMNI in situ 12-hours averaged measurements (black curves) for CR2210 (a), CR2214 (b), and CR2221 (c).
In the text

	Figure 5 Comparison of model results of SR3D-I (green curves), SR3D-R (brown curves) and MAS (blue curves) for the radial velocity Vr (km s−1), number density N (N/cm3), temperature T (MK) in logarithmic scale, radial magnetic field Br (nT), and the magnetic field magnitude |B| (nT) with STA in situ 12-hour averaged measurements (black curves) for CR2190 (a), CR2203 (b), and CR2210 (c).
In the text

	Figure 6 Taylor diagrams to models with OMNI 12-hr avg measurements for all CRs corresponding to Vr (a), Np (b), Br (c), and |B| (d).
In the text

	Figure 7 Taylor diagrams to compare models with STA 12-hr avg measurements for all CRs corresponding to Vr (a), Np (b), Br (c), and |B| (d).
In the text

	Figure 8 SR3D-I (green lines) and SWASTi results using default coronal models based on GONG (brown lines) and ADAPT (blue lines) magnetograms. All plots of CR2210 (a), CR2214 (b), and CR2221 (c) are of 12-hour averaged measurements at L1.
In the text

	Figure A1 The picture demonstrates the differences between the solar wind solutions of the inertial and rotating frames. The differences become evident at larger radial distances.
In the text