© N.M. Narechania et al., Published by EDP Sciences 2021

1 Introduction

Space weather (SW) and its effects on technology have become an important area of research over the past few decades. The importance of the SW field is recognized not only by the research community, but also by government and industry stakeholders. This interest is driven by the fact that SW can have adverse effects on space and ground-based technologies such as spacecrafts, satellites, navigation systems, communications, pipelines and electric power grids. Since human activities increasingly depend on such technology, SW poses a hazard to modern society. In order to predict and mitigate the hazards of SW, research and operational SW forecast activities have been undertaken. For example, a number of forecast centres, frequently established or supported by government organizations, monitor and forecast SW and its impacts.

Propagation of solar wind disturbances through interplanetary space and their subsequent arrival at Earth are certainly one of the focuses of the SW forecast community. High-speed solar wind streams and coronal mass ejections (CMEs) (Howard, 2011; Webb & Howard, 2012; Richardson, 2018) can cause geomagnetic storms and have a significant impact on critical infrastructure. As the typical Sun-to-Earth propagation time for CMEs is 1–4 days, there is generally sufficient time to provide an advanced or early SW warning based on observations of the solar corona. Coronal holes, which are acknowledged to be the source of high-speed solar wind, and CMEs can be captured in the images of the solar corona and provide information about possible adverse SW conditions.

The development of large-scale space plasma simulations has also been a focus of the scientific research community, particularly since the 1980s (see, e.g., Matsumoto & Sato, 1985). Various efforts have been made to describe the solar wind and propagation of transient phenomena, such as CMEs, using physics-based global magnetohydrodynamics (MHD) models. For example, Pizzo (1982) studied the structure of the solar wind between 35 R₀ and 1 AU using an idealized model of the inner heliosphere. Here, R₀ is the radius of the Sun. Usmanov (1993) developed a fully three-dimensional (3D) steady-state model of the solar corona and heliosphere using magnetic field observations for Carrington Rotation (CR) 1682 and made comparisons with spacecraft observations at 1 AU. Lionello et al. (1998) also solved the MHD equations to study the propagation of the solar wind in cylindrical geometry. Linker et al. (1999) modelled the solar corona during Whole Sun Month from 1 R₀ to 30 R₀ using photospheric field observations as boundary conditions and compared results with SOHO, Ulysses and WIND data. A parallel block-adaptive numerical framework was developed for the global MHD simulation of space weather in a series of studies by Powell et al. (1999), Groth et al. (2000), Roussev et al. (2003), Manchester et al. (2004, 2008), Tóth et al. (2007, 2012) and van der Holst et al. (2014). This global MHD model was first applied to the simulation of the interaction between the solar wind and a planetary magnetosphere by Powell et al. (1999). Moreover, Groth et al. (2000) subsequently studied the propagation of a coronal mass ejection (CME) in a steady background solar wind, using an octupole model for the Sun’s magnetic field, a pressure pulse for modelling the CME and source terms for the solar wind acceleration and heating effects. The latter represented one of the first global MHD simulation of a complete fully three-dimensional space weather event, spanning the initiation of a solar wind disturbance at the Sun’s surface to its interaction with the Earth’s magnetosphere. In other follow-on studies, Roussev et al. (2003) modelled the corona-heliosphere system by implementing a continuous variation in the polytropic index in a radially outward direction from the Sun and Manchester et al. (2004) used an idealized model of the steady state solar wind conditions near solar minimum and the global MHD model to study the propagation of a flux-rope-driven CME, making comparisons of the simulated results with coronagraph observations of CMEs. Additionally, Tóth et al. (2007) and Manchester et al. (2008) used synoptic magnetograms to model a CME event of October 28, 2003 and made comparisons with observations. The global MHD model of Powell et al. (1999), Groth et al. (2000), Roussev et al. (2003), Manchester et al. (2004, 2008) and Tóth et al. (2007) was also eventually extended to solve other forms of the MHD equations, including the Hall MHD, multi-fluid MHD, and radiative MHD models (Tóth et al., 2012) and van der Holst et al. (2014) implemented a two-temperature MHD model wherein low-frequency Alfvén wave turbulence was modelled to account for coronal heating and solar wind acceleration. The latter was used to simulate extreme ultraviolet (EUV) images of CR 2107 and comparisons were made to both SOHO and SDO observations.

In other global MHD modelling efforts, Riley et al. (2001) decoupled the MHD simulation in the corona and heliosphere using suitable polytropic indices for each model and subsequently used the coronal simulation as a driver for the heliospheric calculation. The combined model was used to generate the heliospheric structure during CRs 1913, 1892 and 1947. Similarly, Odstrčil et al. (2002) studied CME propagation using a resistive MHD model with a ratio of specific heats, γ = 1.05 for the coronal calculation up to 20 R₀ and an ideal MHD model with a ratio of specific heats, γ = 5/3 for the heliospheric calculation. Lionello et al. (2009) performed MHD simulations of the corona using various coronal heating models and reproduced observed multispectral properties of the corona. More recently, Feng et al. (2010) employed a six-component overset grid to study the background solar wind from Sun to Earth during CR 1911 using line-of-sight photospheric field observations and validated their MHD model using SOHO and WIND observations. Feng et al. (2012) later added isotropic block-based adaptive mesh refinement (AMR) capabilities to this MHD-based model. Lastly, Merkin et al. (2016) developed and applied an ideal MHD model to the heliosphere and a more realistic MHD model to the corona where the latter contained additional terms to account for radiative losses, coronal heating, thermal conduction and magnetic resistivity.

Finally, Keppens et al. (2012) and Porth et al. (2014) have developed the MPI-AMRVAC software, an open source toolkit for parallel, block-adaptive global MHD simulations of solar and non-relativistic astrophysical plasmas (a relativistic counterpart of this software, BHAC, was also developed by Porth et al. (2017) for solving the equations of ideal general-relativistic magnetohydrodynamics and studying astrophysical phenomena such as black holes). The MPI-AMRVAC software has been applied to the study of the formation of prominences in the solar corona using magnetic flux rope models of various levels of complexity (Xia et al., 2012, 2013, 2014; Xia & Keppens, 2016; Zhou et al., 2018). More recently, the MPI-AMRVAC software has also been used for the study of magnetic reconnection of solar flares and to simulate the trans-Alfvénic solar wind from the Sun to the Earth using a solar wind model replicating solar minimum conditions with artificial heating/cooling source terms (Xia et al., 2018), and to study the effect of background solar wind on breakout CMEs using solution-dependent AMR with an ultrahigh-resolution to capture current sheets and small-scale magnetic structures (Hosteaux et al., 2018).

Although significant advances have been made in understanding and predicting SW, from the solar phenomena which are the source of SW to the SW effects on Earth and its magnetic field, more research is still needed to fill the gaps in understanding. Additionally, while SW remains a subject of on-going research efforts, operational SW monitoring and forecasting requires the development of new numerical models and computational tools which can be used in SW operations. The development of the latter is not an easy task. Building a bridge between research and operations and effective research-to-operations transfer are identified as some of the key challenges (Araujo-Pradere, 2009; Merceret et al., 2013; Folini, 2018). Furthermore, the development of operationally oriented physics-based numerical models, in particular Sun-to-Earth simulations, requires significant computational effort. For this reason, a large number of existing models used in SW forecast operations are empirical and/or semi-empirical in nature.

It is well established that the interplay between the Sun’s magnetic field and coronal plasma is the source of solar disturbances and defines the structure of the solar corona. For this reason, SW forecasting using a physics-based numerical approach requires information about the Sun’s magnetic field. Due to the favorable signal-to-noise ratio, regular measurements of the Sun’s magnetic field are restricted to a region close to the Sun’s surface. Nevertheless, solar magnetograms observed from the photospheric/chromospheric region can be used in combination with numerical models to derive estimates of the magnetic field of the solar corona (Mackay & Yeates, 2012). Various models can be used to determine the coronal field. They range from simple models based on potential field theory (Altschuler & Newkirk, 1969; Schatten et al., 1969; Schatten, 1971; Altschuler et al., 1977; Zhao & Hoeksema, 1994), to more complex MHD models. The so-called potential field source surface (PFSS) model, proposed by Altschuler & Newkirk (1969), and Schatten et al. (1969), represents the most widely used model for the coronal magnetic field. Although relatively simple, the PFSS model offers some advantages in comparison to a MHD treatment of the solar corona. For example, PFSS numerical models are computationally significantly less costly than MHD descriptions. In coronal MHD simulations, the length and time scales are very much smaller as compared to those for interplanetary MHD simulations, beyond about 20 R₀ from the Sun, leading to rather high computational costs. Furthermore, as PFSS models neglect plasma dynamics, boundary conditions for plasma properties such as the plasma density, velocity and pressure are not required as in MHD models. Conversely, it is difficult to model phenomena such as coronal heating and solar wind acceleration in the corona using simple PFSS models. Despite the limitations of PFSS models (e.g., they do not include plasma dynamics), Riley et al. (2006) and Owens et al. (2008) have found that the models give similar predictions of the global topology of the coronal magnetic field as those of more complete MHD models.

To simplify Sun-to-Earth solar wind numerical modelling and to avoid the computationally intensive calculations associated with determining the coronal magnetic field, a combined two-model approach is often adopted in a number of SW prediction frameworks in which a semi-empirical data-driven model is used to describe the solar corona, and a global MHD model is used to describe the interplanetary space plasma of the heliosphere. The coupling of the two models is typically enforced at a radius of 20–30 R₀ from the solar surface. The inner boundary conditions for the global MHD model in this case are straightforward to implement as the solar wind attains its asymptotic super-Alfvénic speed by the time it reaches this boundary. Such data-driven solar wind models are usually based on the PFSS and Schatten current sheet (SCS) (Schatten, 1971) models of the solar corona. While such models do not include plasma dynamics and therefore cannot capture acceleration behaviour and regions where the solar wind has slow, Alfvén, and fast magnetosonic speeds, which depend on the magnetic field topology (see, e.g., Keppens & Goedbloed, 2000), the PFSS-SCS based solar wind models make use of empirical relations to associate solar wind plasma properties with the coronal magnetic field (Wang & Sheeley, 1990; Wang et al., 1997; Arge & Pizzo, 2000; Hakamada et al., 2002, 2005). The latter includes the widely used Wang–Sheeley–Arge (WSA) empirical description of the solar wind (Arge & Pizzo, 2000; Arge et al., 2003). For example, Detman et al. (2006) used a PFSS-SCS model for the coronal magnetic field to drive an heliospheric MHD model. A spherical coordinate grid was used in the MHD model, covering 360° in longitude and 90° in latitude. The ENLIL code (Odstrcil, 2003; Odstrcil et al., 2008) is another example of a data-driven solar wind model. The global MHD modelling in ENLIL for the simulation of solar wind structures uses a fixed, uniform or non-uniform mesh, which covers a latitudinal range of +60° to −60° and the full longitudinal range of 360°. The WSA model is used in ENLIL to determine the solar wind speeds, and the solar magnetic field is prescribed using the PFSS model in the region from the Sun’s surface to 2.5 R₀ and by the SCS model in the region beyond 2.5 R₀ out to 21.5 R₀. Baker et al. (2013) employed the WSA-ENLIL code to reproduce solar wind observations from the MESSENGER spacecraft in orbit around Mercury. ENLIL also has the capability to model coronal mass ejections (Taktakishvili et al., 2009, 2011). As another example of a data-driven solar wind prediction framework, Shiota et al. (2014) and Shiota & Kataoka (2016) used a PFSS model for the coronal magnetic field calculations and a WSA model for the solar wind speed to drive global MHD simulations in the heliosphere from 25–30 R₀ onwards using spherical Yin–Yang grids (Kageyama & Sato, 2004).

It should be noted that, while many advances in the understanding of SW phenomena have been made using the preceding large-scale data-driven simulation codes, their implementation for operational solar wind and CME forecasting has been rather slow and is still at a nascent stage. The ENLIL code which operates at the Space Weather Prediction Center (SWPC), National Oceanic and Atmospheric Administration, USA, was the first large-scale MHD simulation code to be transitioned from research to operations in 2011–2012 (Parsons et al., 2011; Steenburgh et al., 2014). This represented an important advancement in SW forecasting. Since then, the ENLIL code has also been implemented at the Korean Space Weather Center and Met Office Space Weather Operations Centre (UK). The ENLIL solar wind and CME simulations, provided by SWPC, have been widely used by the SW forecast community over the past several years.

While ENLIL simulations today represent an important part of current forecasting efforts, it is also important to undertake independent development, performance assessments, and model improvements of solar wind – CME propagation forecast models and codes. For example, the European Heliospheric Forecasting Information Asset (EUHFORIA) was recently developed by Pomoell & Poedts (2018). EUHFORIA also uses a semi-empirical data-driven approach based on the PFSS, SCS, and WSA models to describe the solar corona, and to provide boundary conditions at 0.1 AU to a global MHD description of the inner heliosphere. The computational domain in EUHFORIA, with a uniform spherical mesh in all directions, spans 360° in longitude and ±60° in latitude.

Sun-to-Earth solar wind modelling efforts have also been recently undertaken in Canada. Canada is amongst the countries most affected by SW effects, due to its high-latitude location (see, e.g., Lam, 2011; Boteler, 2019), and there are on-going efforts to both better understand and forecast SW. These efforts include activities at the Canadian Space Weather Forecast Centre (CSWFC), as well as research activities at various Canadian universities. This paper reports on the authors efforts to develop a new integrated data-driven solar wind – CME numerical framework for SW forecasting which couples the coronal and heliospheric subdomains. The primary aim of the study is to provide a flexible new testbed for SW research and advances, as well as to provide a potentially powerful new forecasting tool for solar wind disturbances. In the present configuration of the proposed framework, the computation of the coronal magnetic field is provided by a combination of the PFSS and SCS models (Nikolić, 2017), with associated empirically-based expressions for solar wind flow properties. For the global MHD description of the heliosphere, a second-order-accurate upwind finite-volume scheme (Ivan et al., 2011, 2013, 2015; Susanto et al., 2013) is used to solve the governing ideal MHD equations on a cubed-sphere mesh (Ivan et al., 2011, 2013, 2015) and this finite-volume scheme is combined with the highly-scalable and efficient parallel block-based anisotropic AMR technique developed previously by Williamschen & Groth (2013), Freret & Groth (2015) and Freret et al. (2019). As such, the proposed data-driven solar wind – CME numerical framework presented here includes combined capabilities which are not available in other similar models of the solar wind. Firstly, Global Oscillation Network Group (GONG) synoptic maps of the photospheric field are used to drive the PFSS model. This includes standard and zero-point corrected QuickReduce synoptic maps, as well as Air Force Data Assimilative Photospheric Flux Transport (ADAPT) GONG maps. Additionally, the combination of an advanced high-fidelity finite-volume scheme with parallel anisotropic block-based AMR in the global MHD modelling allows local solution-dependent adaptation of the mesh and affords significantly increased mesh resolution. Furthermore, as mentioned above, a cubed-sphere grid (Ivan et al., 2011, 2013, 2015) is used. The latter avoids the singularities of spherical grids at the poles and readily provides full coverage of the entire range of solid angles associated with 3D space, including high-latitude polar regions. Lastly, the solar wind simulations can be performed in either the inertial or the Sun’s co-rotating frame of reference, offering considerable flexibility for SW modelling.

The organization of the remainder of this paper is as follows. The proposed domains and data-driven mathematical modelling approaches adopted in the framework to describe the solar wind and its disturbances are described in Section 2. This includes the governing equations of the PFSS, SCS and global MHD models. The numerical solution methods, including a description of the GONG synoptic maps of the photospheric field which are used to drive the simulations, the second-order finite-volume spatial-discretization and time-marching schemes which are used to solve the 3D form of the ideal MHD equations, and the parallel block-based AMR approach which provides an efficient and flexible computational mesh and solution procedure for the global MHD model of the inner heliosphere, are all presented in Section 3. Several example numerical results are then discussed in Section 4 to illustrate the capabilities of the proposed numerical framework. This section includes simulation results obtained for a steady-state background solar wind and an unsteady simulation based on time-evolving synoptic maps performed to illustrate potential solar wind forecasting capabilities. Finally, a discussion of solar wind – CME simulations and challenges is provided in Section 5. This includes also a comparison of current predictions to previously reported simulation results. The paper concludes with a summary of the results and findings in Section 6.

2 Simulation subdomains and mathematical models of data-driven framework

As discussed in the introduction, the proposed numerical framework for forecasting the solar wind and its disturbances follows a scheme which is now commonly used in Sun-to-Earth solar wind simulations. In particular, the space is divided into coronal and inner heliosphere subdomains, with the boundary between the two typically set at a radius of 20–30 R₀ from the Sun, where the solar wind speed is both supersonic and super-Alfvénic (see, e.g., Usmanov, 1993; MacNeice et al., 2018). Figure 1 provides an illustration of these simulation subdomains and indicates the various mathematical models used to describe the solar wind in each region. A description of the models for each domain now follows.

Fig. 1 Schematic diagram showing the subdomains and mathematical models used in each region in the proposed Sun-to-Earth data-driven solar wind model. The PFSS model is used to describe the coronal magnetic field from the Sun’s surface up to a spherical surface where magnetic field lines are forced to open. From this surface, which is usually placed at ≈2.5 R0, the SCS model describes the coronal field up to 20–30 R0, where the inner MHD boundary conditions are applied. A global MHD model is then used to describe the solar wind plasma and its propagation through interplanetary space in the inner heliosphere.

2.1 Solar corona

While the solar corona subdomain has been modelled using MHD approaches (see, e.g., Riley et al., 2006), an established semi-empirical approach is adopted herein to describe the coronal magnetic field based on photospheric observations of the magnetic field and the PFSS and SCS models, along with additional empirical expressions to associate the solar wind properties with the magnetic field (see, e.g., Odstrcil et al., 2008; Shiota et al., 2014; Pomoell & Poedts, 2018).

2.1.1 Potential field source surface (PFSS) model

The PFSS model (Altschuler & Newkirk, 1969; Schatten et al., 1969) is used to estimate the global magnetic field, B, of the solar corona from photospheric field observations. The PFSS model is based on the assumption that there are no currents (j = 0) in the region between the Sun’s surface (r = R₀) and the so-called “source surface” (r = R_ss), which means that the magnetic field can be expressed as the gradient of a scalar potential, Ψ, as in

$$\nabla \times \mathbf{B}=0\Rightarrow \mathbf{B}=-\nabla \mathrm{\Psi } .$$(1)

Along with the solenoidal property of the magnetic field, ∇·B = 0, equation (1) leads to the Laplace equation for Ψ

$${\nabla }^2\mathrm{\Psi }=0 .$$(2)

Using a separation of variables with Ψ(r, θ, ϕ) = R(r)Θ(θ)Φ(ϕ) in spherical coordinates, where θ ∈ [0, π] and ϕ ∈ [0, 2π], and assuming that at the source surface, r = R_ss, the magnetic field is purely radial, i.e., Ψ(R_ss, θ, ϕ) = const., the solution of equation (2) in the region R₀ ≤ r ≤ R_ss can be expressed as

$$\mathrm{\Psi }=\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\sum _{m=0}^n\mathrm{}{P}_n^m\left(\theta \right)(g_{\mathrm{nm}}\cos \mathrm{m\varphi }+h_{\mathrm{nm}}\sin \mathrm{m\varphi })\times \left[R_0{\left(\frac{R_0}{r}\right)}^{n+1}-R_{\mathrm{s}\mathrm{s}}{\left(\frac{R_0}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n+2}{\left(\frac{r}{R_{\mathrm{s}\mathrm{s}}}\right)}^n\right],$$(3)

where $P_n^m\left(\theta \right)$ represent the associated Legendre polynomials with Schmidt normalization (see, e.g., Nikolić, 2017). Using equations (1) and (3), the components of the magnetic field can be written as

$$B_r=-\frac{\mathrm{\partial }\mathrm{\Psi }}{\mathrm{\partial }r}=\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\left[(n+1){\left(\frac{R_0}{r}\right)}^{n+2}+n{\left(\frac{R_0}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n+2}{\left(\frac{r}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n-1}\right]\times \sum _{m=0}^n\mathrm{}{P}_n^m\left(\theta \right)(g_{\mathrm{nm}}\cos \mathrm{m\varphi }+h_{\mathrm{nm}}\sin \mathrm{m\varphi }),$$(4)

$$B_{\theta }=-\frac{1}{r}\frac{\mathrm{\partial }\mathrm{\Psi }}{\mathrm{\partial }\theta }=-\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\left[{\left(\frac{R_0}{r}\right)}^{n+2}-{\left(\frac{R_0}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n+2}{\left(\frac{r}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n-1}\right]\times \sum _{m=0}^n\mathrm{}\frac{\mathrm{d}{P}_n^m\left(\theta \right)}{\mathrm{d}\theta }(g_{\mathrm{nm}}\cos \mathrm{m\varphi }+h_{\mathrm{nm}}\sin \mathrm{m\varphi }),$$(5)

$$B_{\varphi }=-\frac{1}{r\sin \theta }\frac{\mathrm{\partial }\mathrm{\Psi }}{\mathrm{\partial }\varphi }=\sum _{n=1}^{\mathrm{\infty }}\mathrm{}\left[{\left(\frac{R_0}{r}\right)}^{n+2}-{\left(\frac{R_0}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n+2}{\left(\frac{r}{R_{\mathrm{s}\mathrm{s}}}\right)}^{n-1}\right]\times \sum _{m=0}^n\mathrm{}{P}_n^m\left(\theta \right)\frac{m}{\sin \theta }(g_{\mathrm{nm}}\sin \mathrm{m\varphi }-h_{\mathrm{nm}}\cos \mathrm{m\varphi }) .$$(6)

The coefficients g_nm and h_nm, can be obtained using equation (4) for the case r = R₀, and employing the orthogonality of Legendre polynomials

$$\frac{1}{4\pi }{\int }_0^{\pi }\mathrm{}{\int }_0^{2\pi }\mathrm{}{P}_n^m\left(\theta \right)\left\{\begin{array}{l}\cos \mathrm{m\varphi }\\ \sin \mathrm{m\varphi }\end{array}\right\}{P}_{n^{\mathrm{\prime }}}^{m^{\mathrm{\prime }}}\left(\theta \right)\left\{\begin{array}{l}\cos m^{\mathrm{\prime }}\varphi \\ \sin m^{\mathrm{\prime }}\varphi \end{array}\right\}\sin \theta \mathrm{d}\theta \mathrm{d}\varphi =\frac{1}{2n+1}{\delta }_n^{n\mathrm{\prime }}{\delta }_m^{m\mathrm{\prime }} .$$(7)

The following expressions for g_nm and h_nm can be obtained:

$$\left\{\begin{array}{c}{g}_{\mathrm{nm}}\\ h_{\mathrm{nm}}\end{array}\right\}=\frac{2n+1}{4\pi \left(n+1+n{\left(\frac{R_0}{R_{\mathrm{s}\mathrm{s}}}\right)}^{2n+1}\right)}{\int }_0^{\pi }\mathrm{}\mathrm{d}\theta \sin \theta P_n^m\left(\theta \right){\int }_0^{2\pi }\mathrm{}\mathrm{d}\varphi B_r(R_0,\theta ,\varphi )\left\{\begin{array}{c}\cos \mathrm{m\varphi }\\ \sin \mathrm{m\varphi }\end{array}\right\} .$$(8)

Note that in equations (3)–(6), the term n = 0 is omitted, since the condition that ∇·B = 0 requires that this term vanish. As the radial component of the photospheric magnetic field, B_r(R₀, θ, ϕ), of equation (8) can be derived from the observed solar magnetograms, the coefficients g_nm and h_nm can be evaluated from the magnetogram data, and the coronal field in the region R₀ ≤ r ≤ R_ss can then be obtained using equations (4)–(6).

2.1.2 Schatten current sheet (SCS) model

The PFSS model uses the current-free approximation of the solar corona in the region R₀ ≤ r ≤ R_ss, and forces the magnetic field lines to be radial at r = R_ss. In order to include effects of plasma currents and describe a non-radial coronal field structure, Schatten (1971) proposed the introduction of a new spherical source, at r = R_cp. From this surface, transverse currents are allowed between regions of opposite polarity of the magnetic field where the Lorentz force, j × B, is small. In the coupling of the PFSS and SCS models, this surface can be at the same location as the PFSS source surface, i.e. R_cp = R_ss, or R_cp can be set below R_ss. A benefit of using R_cp < R_ss is the removal of kinks in the magnetic field lines at the interface of the models (McGregor et al., 2008).

To derive the SCS coronal field for r > R_cp, the magnetic field obtained by the PFSS model is first re-oriented at r = R_cp to point outwards. This means that if B_r(R_cp) ≥ 0, no changes are needed to the field, but if B_r(R_cp) < 0, the signs of magnetic field components B_r(R_cp), B_θ(R_cp) and B_ϕ(R_cp) are reversed. The coronal magnetic field beyond R_cp is obtained by matching this re-oriented field at r = R_cp with the potential field solution for r ≥ R_cp (see, e.g., Schatten, 1971; Nikolić, 2017), and thus

$$B_r=\sum _{n=0}^{\mathrm{\infty }}\mathrm{}(n+1){\left(\frac{R_{\mathrm{c}\mathrm{p}}}{r}\right)}^{n+2}\sum _{m=0}^n\mathrm{}{P}_n^m\left(\theta \right)(g_{\mathrm{nm}}^{\mathrm{\prime }}\cos m\mathrm{\varphi }+h_{\mathrm{nm}}^{\mathrm{\prime }}\sin \mathrm{m\varphi }),$$(9)

$$B_{\theta }=-\sum _{n=0}^{\mathrm{\infty }}\mathrm{}{\left(\frac{R_{\mathrm{c}\mathrm{p}}}{r}\right)}^{n+2}\sum _{m=0}^n\mathrm{}\frac{\mathrm{d}{P}_n^m\left(\theta \right)}{\mathrm{d}\theta }(g_{\mathrm{nm}}^{\mathrm{\prime }}\cos \mathrm{m\varphi }+h_{\mathrm{nm}}^{\mathrm{\prime }}\sin \mathrm{m\varphi }),$$(10)

$$B_{\varphi }=\sum _{n=0}^{\mathrm{\infty }}\mathrm{}{\left(\frac{R_{\mathrm{c}\mathrm{p}}}{r}\right)}^{n+2}\sum _{m=0}^n\mathrm{}{P}_n^m\left(\theta \right)\frac{m}{\sin \theta }(g_{\mathrm{nm}}^{\mathrm{\prime }}\sin \mathrm{m\varphi }-h_{\mathrm{nm}}^{\mathrm{\prime }}\cos \mathrm{m\varphi }) .$$(11)

The role of the PFSS field re-orientation at R_cp is to provide conditions, so that the derived coronal field beyond R_cp using equations (9)–(11), consists strictly of open magnetic field lines. After this field is derived, the final step is to assign proper polarity to the field lines in the region r ≥ R_cp using the polarity obtained prior to the field re-orientation at R_cp. This sign restoration of the magnetic field lines ensures that ∇·B = 0 is not violated. The resulting structure of the field implies that current sheets are introduced between the magnetic fields of opposite polarity.

2.1.3 Solar wind parameters at the inner MHD boundary

The boundary between the inner subdomain (solar corona) and the outer subdomain (inner heliosphere) associated with the global MHD modelling is defined here by a spherically-shaped surface of radius r = R_inner where R_inner ≥ R_ss ≥ R_cp. This surface represents the inner boundary of the global MHD model and, on this spherical shell, boundary conditions are required for the magnetic field of the solar wind as well as the solar wind plasma properties such as velocity, density, and pressure. This boundary is typically defined to be in the range between 20R₀ and 30R₀ (see, e.g., Usmanov, 1993; MacNeice et al., 2018). In this range, the solar wind speed has already reached its asymptotic value and is super-Alfvénic. For example, Odstrcil et al. (2004a) investigated the coupling of coronal and heliospheric simulation models with the boundary located at both 25 R₀ and 50 R₀ and their findings justifies the use of an interface boundary located at 25 R₀.

The coronal magnetic field derived from the PFSS and SCS models is used to directly specify the solar wind magnetic field at the interface between the inner and outer subdomains at r = R_inner. Additionally, the boundary values for the solar wind velocity, density, and pressure are associated with this prescribed magnetic field at r = R_inner using empirical correlations based on solar wind observations. As mentioned before, in solar wind numerical modelling, the WSA model (Wang & Sheeley, 1990; Wang et al., 1995; Arge & Pizzo, 2000; Arge et al., 2003; Sheeley, 2017) is frequently used to correlate the solar wind speed, V_SW, with the magnetic field and is used here. In the WSA model, V_SW is taken to depend on the flux tube expansion factor, f_s, given by

$$f_s=\frac{\left|\mathbf{B}\right(R_0\left)\right|}{\left|\mathbf{B}\right(R_{\mathrm{s}\mathrm{s}}\left)\right|}\frac{R_0^2}{R_{\mathrm{s}\mathrm{s}}^2},$$(12)

as well as the angular separation, θ_b, between an open magnetic field line foot-point and the coronal hole boundary at the photosphere. A general form for the empirical relationship between V_SW and f_s and θ_b used in WSA-type models is given by

$$V_{\mathrm{S}W}=a_1+\frac{a_2}{(1+f_s{)}^{a_3}}{\left[a_4-a_5\exp \left\{-{\left(\frac{{\theta }_b}{a_6}\right)}^{a_7}\right\}\right]}^{a_8} \mathrm{km}/\mathrm{s},$$(13)

where a₁–a₈ are empirical numerical coefficients (MacNeice, 2009). The latter are tunable parameters that depend on the magnetogram source used to derive the coronal magnetic field and calculation parameters.

In order to provide values for the plasma density and temperature at R_inner, additional empirical relations as derived by Hayashi et al. (2003) are also used. These additional empirical correlations are based on Helios observations of the solar wind and were used in previous MHD simulations in which the inner boundary of the global MHD subdomain was specified at R_inner = 50 R₀ (Hayashi et al., 2003; Kataoka et al., 2009). Modified versions of these relations for a boundary R_inner = 25 R₀ were subsequently derived by Shiota et al. (2014). These modified relations for the particle number density and temperature, which depend on the solar wind speed, V_SW, are given by

$$n\left(V_{SW}\right)=4\left\{62.98+866.4{\left(\frac{V_{SW}}{100}-1.549\right)}^{-3.402}\right\} \mathrm{cm}{ }^{-3},$$(14)

$$T\left(V_{SW}\right)=4^{\gamma -1}\left\{-0.455+0.1943\frac{V_{SW}}{100}\right\}10^6 \mathrm{K},$$(15)

respectively. The plasma mass density, ρ, can subsequently be obtained from the particle number density by multiplying with the proton mass m_p, and the plasma thermal pressure p can be obtained from the number density and temperature using the ideal gas equation of state p = nk_BT, where k_B = 1.38 × 10⁻²³ J K⁻¹ is the well-known Boltzmann constant.

2.1.4 CME initialization at the inner MHD boundary

The preceding boundary conditions for the solar wind at the inner MHD boundary are appropriate for representing the “quiescent” background solar wind as a function of the derived coronal magnetic field; however, they cannot accurately represent large-scale transient solar wind features such as CMEs. For the latter, additional boundary treatments are required. CMEs are however extremely complex phenomena and their numerical modelling is not an easy task. Fortunately, images of the solar corona, such as coronagraphs from the LASCO instrument on board the SOHO satellite, can provide an insight into key properties of the CMEs and then simplified empirically-based theoretical models can be used in solar wind – CME simulations to initiate these solar wind disturbances based on the actual measured and estimated parameters from the coronagraph data.

The so-called cone model of Xie et al. (2004) and Zhao et al. (2002) is an example of a data-driven theoretical model for CME initiation. This model has been implemented and used in both the ENLIL (Odstrcil et al., 2004b, 2008) and EUHFORIA (Pomoell & Poedts, 2018) simulation codes. In the cone model, the CME is launched in the interplanetary simulation subdomain at R_inner as a time-dependent plasma cloud. The onset time, location, angular width (i.e. diameter), and speed of the plasma cloud are derived from CME observations, while estimations of the plasma density and temperature are based on additional assumptions. A weakness of the cone model is that the plasma cloud introduced at the inflow boundary does not carry a magnetic field. To include the magnetic field as part of the CME modelling, Shiota & Kataoka (2016) proposed an alternative model. The latter uses a spheromak-type flux rope description for the magnetic field which is deformed into a “pancake” shape. A model of this type was also recently implemented in the EUHFORIA simulation code by Verbeke et al. (2019). Similar to the cone model, the spheromak pancake CME model requires input parameters such as the CME onset time, source location (latitude and longitude), and propagation speed. Additional parameters are used to define the CME shape, including the angular and radial widths of the CME. The magnetic properties of the CME are described using the toroidal magnetic flux, chirality, and two parameters for the spheromak orientation: the tilt and inclination (Shiota & Kataoka, 2016). As is the case with the cone model, some of these CME parameters can be derived directly from coronagraphs; however, the CME magnetic parameters, density, and temperature are generally difficult to infer from observations. Shiota & Kataoka (2016) propose estimating the CME magnetic flux using the class of the solar flare which is associated with the CME. Additionally, estimates of the CME density can be obtained from the SOHO LASCO CME catalog¹, based on several other approximations.

In the proposed integrated data-driven solar wind – CME simulation framework described herein, the spheromak pancake model of Shiota & Kataoka (2016) is used for CME initiation at the inner boundary of the global MHD subdomain. In the present modelling, it is assumed that the CME has a uniform density and pressure. The CME pressure is evaluated using the density from the LASCO catalog and assuming adiabatic expansion of the CME plasma from an initial temperature of 0.8 MK.

2.2 Inner heliosphere subdomain

To propagate the solar wind plasma and CMEs from the inner boundary at r = R_inner and through the computational subdomain extending outward from r = R_inner representing interplanetary space in the inner heliosphere, a global 3D MHD model is used based on the equations of ideal MHD. This global MHD model provides a physics-based description of plasma processes in the solar wind and is now briefly reviewed.

2.2.1 Ideal 3D MHD model

The ideal MHD equations describe the behaviour of a compressible, perfectly electrically conducting, fully ionized, quasi-neutral, inviscid, ideal gas and, as noted in the introduction, are commonly used in the modelling of space plasmas. The equations of ideal MHD in non-dimensional weak conservative form, for a rotating frame, are given by

$$\frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}+\nabla ·\left(\rho \mathbf{u}\right)=0,$$(16)

$$\frac{\mathrm{\partial } }{\mathrm{\partial }t}\left(\rho \mathbf{u}\right)+\nabla ·\left(\rho \mathbf{uu}-\mathbf{BB}+p_T\mathbf{I}\right)=-\rho \left[\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathbf{r}\right)\right]-2\rho \left(\mathbf{\Omega }\times \mathbf{u}\right)-\left(\nabla ·\mathbf{B}\right)\mathbf{B},$$(17)

$$\frac{\mathrm{\partial }e}{\mathrm{\partial }t}+\nabla ·\left((e+p_T)\mathbf{u}-(\mathbf{u}·\mathbf{B})\mathbf{B}\right)=-\rho \mathbf{u}·\left[\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathbf{r}\right)\right]-\left(\nabla ·\mathbf{B}\right)\mathbf{u}·\mathbf{B},$$(18)

$$\frac{\mathrm{\partial }\mathbf{B} }{\mathrm{\partial }t}+\nabla ·\left(\mathbf{Bu}-\mathbf{uB}\right)=-\left(\nabla ·\mathbf{B}\right)\mathbf{u},$$(19)

where ρ, u, e and B are the plasma density, velocity, specific total energy and magnetic field, respectively, t is time and r is the position vector. For simulations performed in a non-inertial rotating frame, Ω is the angular velocity of the reference frame (i.e., the angular rate of rotation of the Sun). The total pressure, p_T, is given by

$$p_T=p+\frac{{\mathit{B}}^2}{2},$$(20)

where p is the thermal pressure of the plasma and the total energy e is given by

$$e=\frac{{\rho \mathbf{u}}^2}{2}+\frac{p}{\gamma -1}+\frac{{\mathbf{B}}^2}{2},$$(21)

with γ = C_p/C_v being the ratio of specific heats. Equations (16)–(19) can be expressed in matrix-vector form as

$$\frac{\mathrm{\partial }\mathbf{U} }{\mathrm{\partial }t}+\nabla ·\mathbf{F}=\mathbf{Q}+\mathbf{S},$$(22)

where the vector of conserved variables, U, and the flux dyad, F, are given by

$$\mathbf{U}=\left[\begin{array}{l}\rho \\ \rho \mathbf{u} \\ e\\ \mathbf{B} \end{array}\right], \hspace{1em}F=\left[\begin{array}{c}\rho \mathbf{u}\\ \rho \mathbf{uu} - \mathbf{BB} +p_T\mathbf{I}\\ \begin{array}{c}(e+p_T) \mathbf{u}-(\mathbf{u}·\mathbf{B})\mathbf{B}\\ \mathbf{Bu}-\mathbf{uB}\end{array}\end{array}\right],$$(23)

and the source term for rotational effects, Q, and the vector containing terms arising from expressing Faraday’s law in divergence form, S, can be expressed as

$$\mathbf{Q}=\left[\begin{array}{c}0\\ -\rho \left[\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathbf{r}\right)\right]-2\rho \left(\mathbf{\Omega }\times \mathbf{u}\right)\\ \begin{array}{c}-\rho \mathbf{u} ·\left[\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathbf{r}\right)\right]\\ \mathbf{0}\end{array}\end{array}\right], \hspace{1em}\mathbf{S}=-\left[\begin{array}{c}0\\ \mathbf{B}\\ \begin{array}{c}\mathbf{u}·\mathbf{B}\\ \mathbf{u}\end{array}\end{array}\right]\nabla ·\mathbf{B},$$(24)

for which Ω × (Ω × r) is the centrifugal force and 2Ω × u is the Coriolis force. The source term S is associated with the enforcement of the divergence constraint on the magnetic field using the approach as proposed by Powell (1994). Theoretically, as in the solar corona subdomain, the magnetic field of the global MHD solutions would be expected to satisfy the solenoidality condition ∇·B = 0 such that S = 0. However, in the proposed finite-volume scheme used for the solution of equation (22) as described in Section 3, this additional constraint is not strictly enforced. By retaining the terms contained in S, the resulting system of conservation laws can be shown to be formally both symmetrizable and Galilean-invariant and also therefore hyperbolic with an eigenstructure that is not degenerate. The application of the proposed finite-volume scheme to equation (22) then ensures that the discrete version of the divergence constraint for B is satisfied to the order of the truncation error of the scheme and these divergence errors are generally convected out of the computational domain (Powell, 1994; Powell et al., 1999).

The preceding non-dimensional solution variables of the ideal MHD equations can be related to their dimensional counterparts by $\rho =\tilde{\rho }/{\rho }_0$, u = $\tilde{\mathbf{u}}$/a₀, $p=\tilde{p}/p_0$, B = $\tilde{\mathbf{B}}$/B₀, $\mathbf{\Omega }=\tilde{\mathbf{\Omega }}(l_0/a_0)$, $t=\tilde{t}/{\tau }_0$ and r = $\tilde{\mathbf{r}}$/l₀. Here, ρ₀ = m_p/cm³ where m_p = 1.672 × 10⁻²⁷ kg is the mass of a proton. The length scale l₀ is taken to be the radius of the Sun, l₀ = R₀ = 6.96 × 10⁸ m. The velocity scale is given by a₀ = l₀/τ₀ = 193.333 × 10³ m/s where the time scale is taken to be τ₀ = 1, h = 3600 s. The value used to normalize the magnetic field is given by B₀ = $\sqrt{{\mu }_0 {\rho }_0 {a_0}^2}$ = 8.8642 × 10⁻⁹ T where μ₀ = 4π × 10⁻⁷ H/m is the magnetic permeability of free space. The value used to normalize the pressure is given by p₀ = ρ₀${a_0}^2$ = 6.252 × 10⁻¹¹ Pa. The angular velocity of the Sun is taken to be $\tilde{\mathbf{\Omega }}$ = (2π/27.27) rad/day $\widehat{\mathbf{k}}$ = 2.67 × 10⁻⁶ rad/s $\widehat{\mathbf{k}}$.

3 Numerical solution methods

3.1 Synoptic maps of the photospheric field

Global Oscillation Network Group (GONG) magnetogram synoptic maps² are used herein to drive the PFSS model and solar wind simulations. GONG provides 24 h coverage of the Sun and is a reliable source of near real-time synoptic maps of the photospheric field (Hill, 2018). For this reason, GONG maps are frequently used in operationally oriented SW applications. This includes WSA-ENLIL at SWPC (Parsons et al., 2011; Steenburgh et al., 2014), as well as EUHFORIA (Pomoell & Poedts, 2018). The GONG magnetic products include standard QuickReduce magnetogram synoptic maps (“mrbqs” in the GONG file name) and newer zero-point corrected QuickReduce maps (“mrzqs” in the GONG file name).

The GONG maps are available in the Flexible Image Transport System (FITS) data format and cover time starting from September 2006. The inferred magnetic field is in units of Gauss (1G = 10⁻⁴ T), and the maps are given on a uniform (sinθ_i, ϕ_j) mesh with i ∈ [1, 180] and j ∈ [1, 360]. The synoptic maps provide the radial component of the photospheric magnetic field B_r(R₀, θ, ϕ), which is used by the PFSS model to obtain the coefficients g_nm and h_nm by means of equation (8).

Both standard and zero-point corrected maps can be used in the proposed data-driven solar wind framework. Furthermore, the framework can also be driven with ADAPT maps (Arge et al., 2010). The latter are based on flux transport modelling which promises better estimates of the global photospheric field distribution. The standard GONG maps have been widely used in other previous research (see, e.g., Shiota et al., 2014; Steenburgh et al., 2014; Pomoell & Poedts, 2018). They have also been used in operational WSA-ENLIL simulations. However, these standard maps, which have been available since 2006, show degradation of their accuracy over time. For example, it was shown recently by Nikolić (2019) that the PFSS solutions obtained using the standard GONG maps exhibit issues, particularly in the polar regions since 2013. The issue with the standard GONG maps can be attributed to background (i.e., zero point) magnetic field variations. The zero-point corrected GONG maps promise to reduce these false variations. The zero-point error is reduced using a software system that compares observations from different sites and in time, better control of the liquid-crystal variable retarders, and a modified data acquisition system (Hill, 2018).

The issues with the standard GONG maps have been recognized by others in the community and the WSA-ENLIL simulation model was recently upgraded to allow for the use of both standard and zero-point corrected maps. In this study, the differences between simulation results obtained using the standard and zero-point corrected maps are highlighted and discussed.

3.2 Numerical solution of PFSS model

As is mentioned, the GONG maps are provided on a uniform sin(θ) − ϕ mesh. It has been shown by Tóth et al. (2011) that PFSS model accuracy in polar regions can be improved by re-meshing the maps on to a uniform (θ_i, ϕ_j) grid, where i ∈ [1, N_θ] and j ∈ [1, N_ϕ]. The new mesh includes the poles and contains an odd number of θ_i points. Re-meshed GONG maps with N_θ = 181 and N_ϕ = 360 are used in this work. Linear interpolation is used to assign magnetic field values from the original GONG map to the re-meshed B_r(R₀, θ_i, ϕ_j) synoptic map.

Following re-meshing, the coefficients g_nm and h_nm are then determined by using a discretized form of equation (8) which can be expressed as

$$\left\{\begin{array}{c}{g}_{\mathrm{nm}}\\ h_{\mathrm{nm}}\end{array}\right\}=\frac{2n+1}{4\pi \left(n+1+n{\left(\frac{R_0}{R_s}\right)}^{2n+1}\right)}\times \frac{2\pi }{N_{\varphi }}\sum _{i=1}^{N_{\theta }}\mathrm{}\sum _{j=1}^{N_{\varphi }}\mathrm{}{ϵ}_i{w}_i{P}_n^m\left({\theta }_i\right)B_r(R_0,{\theta }_i,{\mathrm{\varphi }}_j)\left\{\begin{array}{c}\cos m{\varphi }_j\\ \sin m{\varphi }_j\end{array}\right\},$$(25)

where ϵ₁ = ${ϵ}_{N_{\theta }}$ = 1/2, and ϵ_i = 1 for i ≠ 1, N_θ, and w_i are Clenshaw–Curtis weights given by

$$w_i=-\frac{2}{H}\sum _{k=0}^H\mathrm{}\frac{{ϵ}_k^{\mathrm{\prime }}}{4k^2-1}\cos \left(\frac{\pi k(i-1)}{H}\right),$$(26)

with H = (N_θ − 1)/2, ${ϵ}_0^{\mathrm{\prime }}={ϵ}_H^{\mathrm{\prime }}$ = 1/2, and ${ϵ}_k^{\text{'}}$ = 1 for k ≠ 1, H (see, e.g., Tóth et al., 2011; Nikolić, 2017).

Using equation (25) along with equations (4)–(6), the magnetic field components can then be evaluated. Instead of using an infinite sum (i.e., n → ∞) in equations (4)–(6), a finite-number of harmonic coefficients are considered with the maximum number of terms limited to N. In the present study, a value of N = 120 is used which has been shown to ensure that the numerical errors/artifacts in the resulting PFSS magnetic field solutions are small (see, e.g., Nikolić, 2017). To generate $P_n^m\left(\theta \right)$ and $\mathrm{d}{P}_n^m\left(\theta \right)/\mathrm{d}\theta$ in equations (4)–(6), the following relations are used:

$$R_m^m\left(\theta \right)=\frac{{\left((2-{\delta }_m^0)(2m+1)\left(2m\right)!\right)}^{1/2}}{2^{m}m!}{\sin }^m \theta ,$$(27)

$$R_{m+1}^m\left(\theta \right)=(2m+3{)}^{\frac{1}{2}}\cos \theta R_m^m(\theta ),$$(28)

along with the recursion relation

$$R_n^m\left(\theta \right)={\left(\frac{2n+1}{n^2-m^2}\right)}^{1/2}[(2n-1{)}^{\frac{1}{2}}\cos \theta R_{n-1}^m(\theta )-{\left(\frac{(n-1{)}^2-m^2}{2n-3}\right)}^{1/2}{R}_{n-2}^m\left(\theta \right)] \mathrm{for} n\ge m+2,$$(29)

and $P_n^m\left(\theta \right)=R_n^m\left(\theta \right)/(2n+1{)}^{1/2}$ (Altschuler et al., 1977; Nikolić, 2017).

3.3 Numerical solution of SCS model

A least-squares procedure is used to fit the components of the magnetic field given by equations (9)–(11) at R_cp in the SCS model to the values arising from the PFSS model with the field re-orientation Schatten (1971) and thereby determine values for the coefficients $g_{\mathrm{nm}}^{\text{'}}$ and $h_{\mathrm{nm}}^{\text{'}}$. Denoting the re-oriented field at R_cp as $B_k^{\mathrm{cp}}$(θ_i, ϕ_j), where k = 1, 2, 3 corresponds to the r, θ, ϕ, coordinates, respectively, the magnetic field of the SCS model is defined discretely on a uniform mesh with grid points θ_i and ϕ_j (i ∈ [1, I], j ∈ [1, J]) and using a finite-number of harmonics with maximum degree n = N_s. In the current study, values of I = 181, J = 360 and N_s = 10 are used. The coefficients, $g_{\mathrm{nm}}^{\text{'}}$ and $h_{\mathrm{nm}}^{\text{'}}$, of equations (9)–(11), are then determined in the least-squares approach by minimizing the sum of squared residuals, F, defined by

$$F=\sum _{k=1}^3\mathrm{}\sum _{i=1}^I\mathrm{}\sum _{j=1}^J\mathrm{}{\left[B_k^{\mathrm{c}p}({\theta }_i,{\mathrm{\varphi }}_j)-\sum _{n=0}^{N_s}\mathrm{}\sum _{m=0}^n\mathrm{}\left(g_{\mathrm{nm}}^{\mathrm{\prime }}{\alpha }_{\mathrm{knm}}({\theta }_i,{\varphi }_j)+h_{\mathrm{nm}}^{\mathrm{\prime }}{\beta }_{\mathrm{knm}}({\theta }_i,{\varphi }_j)\right)\right]}^2,$$(30)

such that $\mathrm{\partial }F/\mathrm{\partial }{g}_{\mathrm{nm}}^{\text{'}}$ = 0 and $\mathrm{\partial }F/\mathrm{\partial }{h}_{\mathrm{nm}}^{\text{'}}$ = 0. In equation (30), α_knm and β_knm are given by

$$\begin{array}{c}{\alpha }_{1\mathrm{nm}}({\theta }_i,{\varphi }_j)=\left(n+1\right)P_n^m\left({\theta }_i\right)\cos m{\varphi }_j,\\ {\alpha }_{2\mathrm{nm}}({\theta }_i,{\varphi }_j)=-{\frac{\mathrm{d}{P}_n^m\left(\theta \right)}{\mathrm{d}\theta }|}_{\theta ={\theta }_i}\cos m{\varphi }_j,\\ {\alpha }_{3\mathrm{nm}}({\theta }_i,{\varphi }_j)=\frac{m}{\sin {\theta }_i}{P}_n^m\left({\theta }_i\right)\sin m{\varphi }_j,\end{array}$$(31)

$$\begin{array}{c}{\beta }_{1\mathrm{nm}}({\theta }_i,{\varphi }_j)=\left(n+1\right)P_n^m\left({\theta }_i\right)\sin m{\varphi }_j,\\ {\beta }_{2\mathrm{nm}}({\theta }_i,{\varphi }_j)=-{\frac{\mathrm{d}{P}_n^m\left(\theta \right)}{\mathrm{d}\theta }|}_{\theta ={\theta }_i}\sin m{\varphi }_j,\\ {\beta }_{3\mathrm{nm}}({\theta }_i,{\varphi }_j)=\frac{m}{\sin {\theta }_i}{P}_n^m\left({\theta }_i\right)\cos m{\varphi }_j .\end{array}$$

For each (m, n), the minimization of equation (30) results in a linear system of equations for $g_{\mathrm{nm}}^{\text{'}}$ and $h_{\mathrm{nm}}^{\text{'}}$ which can then be solved following some linear algebra. The solutions for $g_{\mathrm{nm}}^{\text{'}}$ and $h_{\mathrm{nm}}^{\text{'}}$ can subsequently be used to define the SCS magnetic field for r ≥ R_cp. As a final step, the polarity of magnetic field lines is restored in the r ≥ R_cp region to match the polarity of the R₀ ≤ r ≤ R_cp coronal magnetic field solution of the PFSS model.

3.4 Coupling of PFSS, SCS, and global MHD models

The locations of the coupling surfaces of the PFSS, SCS, and global MHD models are in general free parameters in the proposed integrated solar wind – CME framework. In the present study, the PFSS and SCS models are coupled on a surface of radius r = R_cp = R_ss = 2.5 R₀ and the inner boundary of the global MHD model is placed at r = R_inner = 25 R₀. Equations (13)–(15) are then used to provide the solar wind plasma properties at R_inner. The WSA relation of equation (13) with a₁ = 250, a₂ = 875, a₃ = 0.2, a₄ = 1, a₅ = 0.8, a₆ = 2.6, a₇ = 1.25 and a₈ = 2.5 is used. For the properties needed for the open magnetic field lines, f_s and θ_b, a second-order Runge-Kutta method to trace the field lines is used.

For the global MHD model, a suitable coordinate transformation and linear interpolation is used to map the solution obtained from the combined PFSS-SCS modelling on to a (θ, ϕ, t) grid representing the inner boundary of the MHD model. This can be done for both non-inertial co-rotating and non-rotating inertial frames which account for the solar rotation. For calculations performed in the Sun’s co-rotating frame, the variables at the inner boundary of the MHD model are taken to be fixed for a given synoptic map. The latitudinal component of the solar wind velocity is assumed to be zero, i.e. V_θ = 0, while the radial component of the solar wind velocity V_r is obtained directly from equation (13). Additionally, in this co-rotating frame, the longitudinal component is also taken to be zero, i.e. V_ϕ = 0. The radial component of the magnetic field at the inner boundary of the MHD model is obtained from the SCS model. The latitudinal and longitudinal components are assumed to be zero, i.e., B_θ = 0 and B_ϕ = 0. A modified procedure with interpolation is adopted for the non-rotating case to account for the time-dependent nature of the solar wind plasma properties at the inner boundary (Hayashi, 2012).

3.5 Solution-adaptive upwind finite-volume scheme for ideal MHD equations

3.5.1 Spatial discretization and semi-discrete form

In the proposed Godunov-type upwind finite-volume scheme (Godunov, 1959; Toro, 2013) for solving the ideal MHD equations of equation (22), the spatial discretization is accomplished herein by applying a second-order cell-centered finite-volume scheme (Ivan et al., 2011, 2013, 2015; Susanto et al., 2013). Application of this discretization procedure to the integral form of the governing MHD equations for a control volume as defined by a hexahedral computational cell, (i, j, k), of a structured 3D grid (see Fig. 2) results in the following semi-discrete form:

$$\frac{{\mathrm{d}\mathbf{U}}_{i,j,k}}{\mathrm{d}t}=- {\mathbf{R}}_{i,j,k}\left(\mathbf{U}\right)=-\frac{1}{V_{i,j,k}}\sum _{f=1}^{N_{f_{i,j,k}}}\mathrm{}{\left({\mathbf{F}}_f. {\mathbf{n}}_f\Delta A_f\right)}_{i,j,k}+\frac{ {\mathbf{S}}_{i,j,k}}{V_{i,j,k}}\sum _{f=1}^{N_{f_{i,j,k}}}\mathrm{}{\left({\mathbf{B}}_f. {\mathbf{n}}_f\Delta A_f\right)}_{i,j,k}+{\mathbf{Q}}_{i,j,k},$$(32)

where U_i,j,k is the average value of the conserved solution vector for cell (i, j, k), and R_i,j,k is called the discrete residual representing the sum of the face fluxes for cell (i, j, k), as well as the contributions of the volume sources. Here, $N_{f_{i,j,k}}$ represents the number of faces for the cell and $N_{f_{i,j,k}}$ = 6 for hexahedral computational cells. The variables V_i,j,k, S_i,j,k, Q_i,j,k denote the cell volume, Powell source term and rotational effect source term, F_f, B_f, n_f and ΔA_f are the flux vector, magnetic field vector, outward pointing unit normal vector and the area of the cell face, f, respectively. For the evaluation of the inviscid fluxes, limited piecewise-linear least-squares reconstruction is used for calculating the primitive flow variables at each cell face. A Riemann-solver-based function (Godunov, 1959; Toro, 2013), namely the so-called HLLE approximate Riemann solver proposed by Einfeldt (1988) based on the HLL flux function of Harten et al. (1983), is then used for evaluation of the inviscid fluxes at the cell faces.

Fig. 2 Hexahedral cell at grid location i, j, k showing face normals.

3.5.2 Newton method for steady solar wind flows

Steady-state or time-invariant solutions of the semi-discrete form of the governing equations given by equation (32) satisfy a large coupled system of non-linear algebraic equations for all computational cells in the mesh given by

$$\frac{\mathrm{d}\mathbf{U} }{\mathrm{d}t}=- \mathbf{R}\left(\mathbf{U}\right)=0,$$(33)

and an iterative technique based on Newton’s method can be effective in their solution. In particular, the parallel inexact Newton’s method developed by Northrup & Groth (2013) is used. Application of Newton’s method to the solution of equation (33) yields a system of linear equations of the form

$$\left[\frac{\mathrm{\partial }\mathbf{R} }{\mathrm{\partial }\mathbf{U} }\right]{\Delta \mathbf{U}}^{\left(n\right)}={\mathbf{J}\Delta \mathbf{U}}^{\left(n\right)}=-\mathbf{R}\left({\mathbf{U}}^{\left(n\right)}\right),$$(34)

where U⁽ⁿ⁾ is the solution change associated with the nth iteration and successively improved estimates of the solution to the non-linear equations satisfy

$${\mathbf{U}}^{(n+1)}= {\mathbf{U}}^{\left(n\right)}+{\Delta \mathbf{U}}^{\left(n\right)},$$(35)

and where J = ∂R/∂U is the Jacobian of the solution residual. Given an initial estimate of the solution, U⁽ⁿ⁼⁰⁾, the linear system given by equation (34) is solved repeatedly and the solution updated at each step, n, of Newton’s method until the solution residual is deemed to be sufficiently small, i.e., ||R(U)|| < ϵ where ϵ is a user-defined tolerance. A value of ϵ = 10⁻⁸ is used for this work. A slightly modified version of the scheme outlined above, based on the switched evolution/relaxation (SER) approach as proposed by Mulder & van Leer (1985), can be used here to aid in the global convergence of the Newton method.

The system of linear equations represented by equation (34) is both non-symmetric and sparse and also typically very large. As such, Krylov subspace iterative methods can be very effective in determining their solution. A “matrix-free” or “Jacobian-free” version of the generalized minimum residual (GMRES) method originally developed by Saad & Schultz (1986) is used here to obtain solutions to equation (34). For the GMRES method to be effective, preconditioning is required and an additive Schwarz global preconditioner is used in conjunction with block incomplete lower-upper (BILU) local preconditioning (Northrup & Groth, 2013).

3.5.3 Explicit time-marching scheme for unsteady solar wind flows

For time-dependent or unsteady solar wind flows, solutions of the coupled system of non-linear ordinary differential equations (ODEs) represented by the semi-discrete form of the governing equations given by

$$\frac{\mathrm{d}\mathbf{U} }{\mathrm{d}t}=-\mathbf{R}\left(\mathbf{U}\right),$$(36)

are sought. A standard explicit, two-stage, second-order accurate, Runge–Kutta time-marching scheme (Butcher, 1996; Lomax et al., 2013) is used here to evolve the solutions of the ideal MHD equations forward in time.

3.5.4 Anisotropic block-based AMR

The hexahedral computational cells of the global MHD grid are contained within multi-block body-fitted structured grids that permit a general unstructured connectivity of the grid blocks thereby readily allowing the use of cubed-sphere grids (Ivan et al., 2011, 2013, 2015). The proposed multi-block grid structure also readily facilitates automatic, solution-dependent, local adaptation of the mesh (Williamschen & Groth, 2013; Freret & Groth, 2015; Freret et al., 2019), and leads to efficient and scalable parallel implementations of the solution algorithm on distributed memory high-performance computing systems. In the proposed approach, the computational domain is divided into subdomains or “grid blocks”, where each block contains a predefined number of cells. During anisotropic refinement, the blocks can be divided into two, four, or eight blocks, thereby doubling the mesh resolution in preferred directions as dictated by solution-dependent refinement criteria. Coarsening of the mesh is accomplished by reversing this process. The particular form of the anisotropic block-based AMR scheme adopted here (Williamschen & Groth, 2013; Freret & Groth, 2015; Freret et al., 2019) has been shown to be highly efficient in reducing the overall mesh size for a given flow problem.

The data structure used to store the grid-block connectivity is a hierarchical flexible binary tree. This data structure not only provides the connectivity between individual blocks but also the level and sequence of block refinements associated with the computational mesh. Figure 3 depicts the grid block structure and the resulting binary tree after several refinements of an initial mesh consisting of a single grid block. Note that the binary data structure is relatively light and compact in terms of memory, as it only has to account for the connectivity between blocks and not the individual cells. This light storage requirement for connectivity readily allows for dynamic and local refinement of the mesh and, combined with the self-similar nature of grid blocks, leads to efficient and highly scalable implementations of the combined finite-volume scheme and AMR procedure on high-performance parallel computing systems using domain decomposition (Gao & Groth, 2010; Gao et al., 2011; Freret & Groth, 2015).

Fig. 3 Schematic diagram of binary tree and corresponding 3D grid-block structure of computational mesh after several levels of refinement.

Each grid block of the computational mesh is surrounded by two layers of ghost cells which carry solution information from the neighboring blocks. This information is exchanged via message passing procedures between neighbouring blocks and allows for the simulation and update of the solution to proceed on each block in an independent fashion. In this implementation, which allows for the use of heterogeneous neighbouring blocks, the ghost cells for a given block are stored at the same refinement level as the neighbouring blocks and no restriction and prolongation procedures are required since the least-squares reconstruction that is used automatically handles irregular stencils with varying grid resolution (Freret & Groth, 2015; Freret et al., 2019). This heterogeneous approach is illustrated in Figure 4. The domain in Figure 4 consists of 8 blocks adjacent to each other as shown in Figure 4a. The ghost cells for the block of interest at the centre, ID, are shown in Figure 4b. The ghost cells and their solution values are directly provided by the neighboring blocks. This eliminates the need for prolongation of cell-averaged values from coarser to finer cells and restriction from finer to coarser cells, when exchanging information between blocks. When there are non-conforming cells or so-called “hanging nodes” where neighbouring blocks meet, the fluxes through the non-conforming faces are calculated in a systematic way, in a manner that is the same for neighbouring blocks, so as to maintain the conservation properties of the finite-volume scheme (Freret & Groth, 2015; Freret et al., 2019). In this way, the need for flux correction strategies to ensure the conservation properties of the scheme at block interfaces with resolution changes is also completely eliminated.

Fig. 4 Schematic diagram of heterogeneous grid-block structure showing (a) the grid block of interest, marked ID, and the associated neighbouring blocks; and (b) the ghost cells of block ID taken directly from the neighbouring blocks.

4 Simulation capabilities of the numerical framework

In this section, the capabilities of the proposed data-driven solar wind – CME prediction framework are demonstrated by considering two representative simulations. In the first example, a steady-state simulation representing a relatively quiescent solar wind is examined based on a single standard GONG synoptic map. The inexact Newton’s method discussed in Section 3.5.2 was used to obtain the steady-state solution and this example is used to illustrate both the effectiveness of the fully converged steady Newton solutions and the efficiency of the anisotropic AMR method in capturing the current sheet. In the second example, multiple synoptic maps were used to produce an example of an unsteady solar wind simulation and provide a forecast of the solar wind conditions at Earth. For this unsteady simulation, the second-order, explicit, time-marching scheme described in Section 3.5.3 was used and the predictions of the solar wind speed, density, and magnetic-field predictions are compared with available observations. Both simulations were performed in the co-rotating frame of the Sun and a cubed-sphere computational mesh was used (Ivan et al., 2011, 2013, 2015) which initially consisted of 144 blocks, with each block containing 8 × 8 × 8 = 512 cells. A portion of the initial mesh is shown in Figure 5. To accurately capture the current sheet, the gradient of the magnetic field weighted with the square of the radial distance from the center of the Sun, ∇(|B|r²), was used as the criterion for the anisotropic mesh adaptation. The weight r² was introduced to account for the fact that in a supersonic spherical outflow with constant velocity, the magnetic field strength rapidly decreases in proportion to 1/r² in the radially outward direction. Blocks having gradients higher than a predefined refinement threshold were flagged for refinement while those with gradients below a predefined coarsening threshold were flagged for coarsening. A maximum refinement level of seven was imposed for both simulations, with zero being the initial mesh refinement level, for a total of eight levels.

Fig. 5 Portion of the initial cubed-sphere mesh used in the simulation of both the steady and unsteady solar wind. The inner boundary of the global MHD subdomain is shaded in green. The thicker black lines are the grid block edges and the thinner red lines are the edges of the computational cells. Each block contains 8 × 8 × 8 = 512 cells and the total number of grid blocks is 144 for a total number of 144 × 512 = 73,728 computational cells.

4.1 Steady solar wind simulation

For the steady solar wind case, a single standard GONG synoptic map was used as the input for the coronal model calculations and PFSS and SCS models. In particular, the GONG map for 12 February 2007 (11:54 UT) was used. The resulting surface plots of the predicted radial component of the magnetic field, solar wind speed, plasma density, and thermodynamic pressure at the inner MHD boundary, r = 25 R₀, are shown in Figure 6. Here, the central meridian and sub-Earth locations are denoted with dash–dot line and dash–dash line, respectively. The current sheet across which the magnetic field changes sign can be distinctively seen in Figure 6a. The higher density plasma originates from the region of the current sheet as can be seen in Figure 6c. The fast high-pressure solar wind originates mainly from the higher latitudes closer to the poles as can be seen in Figures 6b and 6d.

Fig. 6 Surface plots of the estimated solar wind properties at r = 25 R0 for the steady solar wind simulation obtained using the standard GONG map for 12 February 2007 (11:54 UT). The central meridian and sub-Earth locations are denoted with dash-dot and dash-dash line, respectively.

Using the inflow boundary conditions arising at r = 25 R₀ shown above and after obtaining the time-invariant ideal MHD solution on the initial mesh, seven successive refinements of the global MHD mesh were performed and a fully converged solution was obtained on each mesh using the proposed inexact Newton’s method. The final mesh after 7 refinements contained a total of 45,799 blocks, that is 45,799 × 512 = 23,449,088 cells, ≈23.5 million cells. A total of 2048 intel processor cores were used over 3 h to perform the entire simulation. A convergence plot of the global MHD computations on the 8 successive meshes showing the density residual as a function of the number of Newton steps is given in Figure 7. It can be seen that a fully converged steady state solution was obtained on the final mesh comprising 23.5 million cells in a total of about 1200 iterations, with an average of about 150 iterations per mesh level.

Fig. 7 Convergence history of density residual for steady solar wind simulation showing convergence of Newton method on the sequence of anisotropic AMR meshes.

Figure 8 shows the predicted contours of |B|r² at the outer boundary r = 225 R₀ (≈1.05 AU) of the computational domain on meshes obtained after each successive adaptive mesh refinement. With each refinement, there is an increase in the number of blocks in the vicinity of the current sheet. The proposed anisotropic AMR proves to be both very effective and efficient in capturing the discontinuous nature and complex topology of the current sheet by refining extensively only across the current sheet while keeping the grid resolution low along the current sheet surface. Figure 9 shows contours of |B|r² and the comparison between the structure of the current sheet on xz- and yz-slices of the computational domain on the initial and final meshes. This comparison qualitatively shows how the current sheet is captured not only at the outer boundary but throughout the entire computational domain. Finally, Figure 10 shows a comparison between contours of |B|r² and magnetic field lines on a portion of the ecliptic plane (xy-plane) on the initial and final meshes. On the initial mesh, the oppositely directed lines across the current sheet artificially reconnect due to numerical diffusion. On the final mesh, these oppositely directed field lines do not exhibit reconnection and are predicted to lie very close together on either side of the narrow current sheet, accurately capturing the magnetic field discontinuity across the sheet. Without the anisotropic AMR approach, it would be extremely computationally expensive to obtain such high-quality ideal MHD solutions of the solar wind unless a much higher number of grid points were used.

Fig. 8 Predicted contours of |B|r2 at the outer boundary r = 225 R0 (≈1.05 AU) of the computational domain after 0, 2, 4, and 7 refinements. The black lines are the block boundaries.

Fig. 9 Predicted contours of |B|r2 in the xz-plane (top) and yz-plane (bottom) for the initial (73,728 cells) (left) and final (23,449,088 cells) (right) meshes showing capturing of the current sheet. The black lines are the block boundaries.

Fig. 10 Portion of the xy-plane (ecliptic plane) showing the predicted contours of |B|r2 and magnetic field lines on initial (73,728 cells) (left) and final (23,449,088 cells) (right) meshes.

The predicted distributions of the magnitude of the solar wind velocity obtained on the final mesh for each of the three Cartesian planes are shown in Figure 11. The slow and fast streams of the simulated solar wind can be clearly identified here. The fast solar wind emerges mainly from the polar regions of the Sun while the slow solar wind emerges mainly from the region near the solar equator for periods of low solar activity. The green-coloured sphere, denoted with an E, represents the approximate location of the Earth. Finally, Figure 12 provides the predicted contours of |B|r² in the ecliptic xy-plane of the Sun–Earth system obtained on the final mesh along the isosurface of the current sheet lying on the northern side of the plane. The Parker spiral (Alfvén, 1957; Parker, 1958) can be clearly seen in the ecliptic plane of both Figures 11 and 12.

Fig. 11 Predicted distributions of the magnitude of the solar wind velocity obtained on the final mesh (23,449,088 cells) after 7 successive anisotropic adaptive mesh refinements. The black lines are the block boundaries. The green-coloured sphere, E, corresponds to the approximate location of Earth.

Fig. 12 (a) Predicted contours of |B|r2 in the ecliptic xy-plane obtained on the final mesh (23,449,088 cells) along with (b) the isosurface of the current sheet lying on the northern side of the plane. The black lines correspond to the grid block boundaries.

4.2 Unsteady solar wind simulation

In the example of an unsteady solar wind simulation with time-dependent inner boundary conditions for the global MHD model considered here, a sequence of 64 consecutive standard GONG magnetogram synoptic maps was used to estimate the coronal magnetic field and solar wind parameters at the inner boundary of the MHD subdomain. These synoptic maps were obtained over a span of 16 days starting from 18 January to 2 February, 2007. Four magnetograms per day were used with an interval of typically 6 hours between any two consecutive maps. The first synoptic map was dated 18 January (5:54 UT), while the last magnetogram was dated 2 February (23:54 UT). To illustrate time-dependent behavior of the solar wind properties at the inner boundary of the global MHD subdomain, surface plots of the estimated radial component of the magnetic field obtained at r = 25R₀ using the PFSS and SCS models, along with the corresponding solar wind speed obtained using the WSA relation are given in Figure 13. The results shown in the figure correspond to the GONG maps for 22 and 26 January 2007 (5:54 UT). The central meridian and sub-Earth positions are indicated in Figure 13 with dash–dot and dash–dash lines, respectively. It is evident from the surface plots of Figure 13 that, due to the Sun’s rotation, there is a lateral shift in the location of the current sheet between 22 and 26 January. Furthermore, while the location of the high-speed solar wind sources evolves with time, looking at the equatorial sources, it is also evident that there are changes in the source area itself.

Fig. 13 Surface plots of the estimated radial component of the magnetic field obtained at the inner boundary r = 25 R0 obtained using the PFSS and SCS models with synoptic maps for 22 and 26 January 2007 (5:54 UT), and the corresponding solar wind speed obtained using the WSA model. The central meridian and sub-Earth locations are denoted with dash-dot and dash-dash line, respectively.

To initiate the unsteady simulation for the sequence of 64 synoptic maps, a steady-state solution corresponding to the first synoptic map was first obtained on the initial mesh and 7 adaptive mesh refinements were subsequently applied to obtain the starting initial solar wind solution. The proposed Newton method was used to obtain this initial solution. After this, the simulation was advanced in time until the date of the last magnetogram using the explicit second-order time-marching scheme. The inner boundary conditions for the global MHD subdomain were continuously updated using linear interpolation between successive synoptic maps. Adaptive mesh refinement and coarsening was periodically applied to the solution after every 60 time-steps and a Courant–Friedrichs–Lewy number of 0.3 was used throughout the simulation.

Figures 14a and 14b show the predicted contours of |B|r² and the isosurfaces of the current sheet obtained on intermediate adapted meshes for 26 January (5 pm UT) and 31 January (11 pm UT), 2007, respectively. Figure 15 provides the predicted distributions of the magnitude of the solar wind velocity in the ecliptic plane at various instances in time during the unsteady simulation. The dates and approximate times of these snapshots are: (a) 23 January (12 pm UT); (b) 25 January (1 pm UT); (c) 27 January (10 am UT); and (d) 29 January (10 am UT). The sizes of the cubed-sphere AMR mesh at these times are: 5,569,024 cells as shown in Figure 15a; 13,434,880 cells as depicted in Figure 15b; 9,246,720 cells as given in Figure 15c; and 8,822,272 cells as shown in Figure 15d. The brown-coloured sphere in the figure represents the position of the Earth at each instance in time. Note that, in the co-rotating reference frame, the Earth moves in a clockwise direction at the rate of (2π/27.27) rad/day.

Fig. 14 Predicted contours of |B|r2 in the xy-, yz- and xz-planes, and current sheet isosurfaces on adapted meshes showing the grid block boundaries (not individual cells) on 26 January (5 pm UT) (17,671,680 cells) and 31 January (11 pm UT) (31,553,536 cells) 2007, respectively, for the unsteady solar wind simulation. The black lines are the block boundaries.

Fig. 15 Solar wind speed on the ecliptic plane at various time instances. The dates and approximate times are (a) 23 January (12 pm UT) (5,569,024 cells), (b) 25 January (1 pm UT) (13,434,880 cells), (c) 27 January (10 am UT) (9,246,720 cells), and (d) 29 January (10 am UT) (8,822,272 cells). The brown spot indicates the approximate location of the Earth.

The time evolution of the solar wind velocity, plasma density, and magnetic-field strength at the Earth can be obtained from the unsteady simulation. In Figure 16, these predicted quantities are compared to available ACE satellite observations. The predicted sharp increase in the solar wind speed depicted in Figure 16a agrees qualitatively well with the ACE satellite observations. This increase in solar wind speed is also reflected in the snapshots provided in Figure 15 where the Earth is observed to move into the faster solar wind at later times during the simulation, i.e., sometime just after the snapshot of Figure 15c. It is also evident that the rapid increases and peaks in the predicted temporal variations of the plasma number density and magnetic-field strength of Figure 16 are also in rather good agreement with the other corresponding ACE satellite observations.

Fig. 16 Comparison of (a) solar wind speed, (b) plasma density, and (c) magnetic field strength at Earth predicted by the proposed SW framework to ACE satellite observations for the unsteady solar wind simulation.

5 Solar wind – CME simulations and forecasting

The previous section illustrates some of the capabilities of the proposed data-driven solar wind – CME framework for SW forecasting. In the examples given, synoptic maps from a relatively quiet period of solar activity were used to directly drive the simulations and, as illustrated by the comparisons of Figure 16, good agreement between the simulation results and observations can be obtained. The agreement between predicted solar wind properties and observations, such as those shown in Figure 16, can be further improved, for example, by re-tuning the empirical solar wind relations of equations (13)–(15). However, while re-adjusting empirical relations can help a particular simulation case, the operationally oriented solar wind – CME framework should perform well during periods of low solar activity, when high-speed streams dominate, as well as during an active period when CMEs play a significant role. Frequently, forecasting models are tuned and tested using datasets from a relatively small time window of solar activity, in comparison to the longer 11-year solar activity cycle. For example, the GONG standard magnetogram synoptic maps are available from late 2006 and the WSA-ENLIL code, which uses these maps, was transitioned to SWPC operations in 2011 (Parsons et al., 2011).

The forecasting accuracy of numerical models can change as the solar cycle progresses. The changes in the accuracy can not only be due to the fact that the solar activity changes, but also due to changes in the input data quality. Recently, for example, solutions of the PFSS and SCS models based on standard GONG maps have been investigated by Nikolić (2019) for the period from 2006–2018. It was shown that the area of the coronal holes, which are known to be sources of the high-speed wind, derived using the standard GONG maps exhibits significant deviations after 2013, particularity in high-latitude regions. As already mentioned in Section 3.1, the standard GONG maps have been widely used in research and operational solar wind models and their performance assessments (see, e.g., Hinterreiter et al., 2019). Since these synoptic maps exhibit significant zero-point error, using the zero-point corrected GONG maps to derive the coronal magnetic field in the solar wind models should be considered.

To illustrate the differences between standard and newer zero-point corrected GONG maps, the coronal holes arising from the PFSS and SCS models and the corresponding solar wind speed maps obtained using the WSA relation (13) are shown in Figure 17 as obtained using (a) the standard GONG map and (b) the zero-point corrected GONG map. Both maps are for 25 June 2015 (01:04 UT). The map time is associated with the central meridian which is denoted in the figures by the dash-dash line while the dash-dot line denotes projections of the Earth onto the map at different times. The estimated coronal holes represent foot points of open magnetic field lines at the photosphere. Here, 1/F_s is associated with coronal holes and the predicted contours of this quantity are shown in the figure. The quantity, F_s, represents the flux tube expansion factor of a magnetic field line multiplied with its magnetic field polarity, i.e. F_s = sign(B_r(R₀)) f_s. In Figure 17, 1/F_s is saturated at ±0.1, and the red and blue color represent the magnetic field lines which are pointed away from and toward the Sun, respectively.

Fig. 17 Predicted coronal hole (left) and solar wind speed (right) maps of the PFSS and SCS models obtained using (a) the standard GONG map and (b) the zero-point corrected GONG map for 25 June 2015 01:04 UT.

The results shown in Figure 17a were computed using the same standard GONG synoptic map as in Pomoell & Poedts (2018). A comparison of the current results with Figure 4 of this previous study, reveals a very good agreement between the estimated coronal holes of the standard map. It is evident in this case that the southern coronal hole does not cover all of the south polar region. However, as can be seen in Figure 17b, use of the zero-point corrected map in the PFSS model produces a southern coronal hole that occupies the entire polar region. Furthermore, there are also some differences in the low latitude coronal holes predicted by the two maps and some portions of the northern coronal holes, present in the results for the standard map, are missing from the results of the zero-point corrected map. As can be seen from the corresponding solar wind speed maps of Figure 17, the differences in the coronal holes representing the open magnetic field lines obtained using standard and zero-point corrected GONG maps can also result in significantly different inputs to the inner boundary of global MHD interplanetary simulations. Note that, in Figure 17, the contour plots of the solar wind speed is saturated below 300 km/s and above 800 km/s.

As further illustration of the capabilities of the proposed data-driven solar wind – CME numerical framework and to further demonstrate the differences between simulations performed with the standard and zero-point corrected maps, an unsteady simulation with multiple CMEs is now considered. In particular, simulations were performed in the inertial frame for the same series of CMEs as examined in the previous study by Pomoell & Poedts (2018). This series includes five CMEs which erupted between 18th and 25th June, 2015. The CME analysis, from the Space Weather Database Of Notifications, Knowledge, Information (DONKI)³, gives these CMEs on June 18 (20:00 UT) at (θ, ϕ) = (11°, −50°) with a speed of 1000 km/s (CME-1), June 19 (14:59 UT) at (θ, ϕ) = (−33°, 9°) with a speed of 603 km/s (CME-2), June 21 (05:01 UT) at (θ, ϕ) = (7°, 8°) with a speed of 1250 km/s (CME-3), June 22 (21:10 UT) at (θ, ϕ) = (14°, 3°) with a speed of 1155 km/s (CME-4), and June 25 (10:51 UT) at (θ, ϕ) = (23°, 46°) with a speed of 1450 km/s (CME-5). The angular half-width of the five CMEs were 45°, 54°, 47°, 45°, and 41°, respectively. In the present simulations, the estimated masses of the CMEs from the SOHO LASCO CME catalog⁴ were used. As such, the CME masses were set to 18.0 × 10¹² kg, 6.1 × 10¹² kg, 9.6 × 10¹² kg, 4.4 × 10¹² kg, and 31.0 × 10¹² kg, respectively, and the radial width of the CMEs was taken to be 2 R₀. As in Pomoell & Poedts (2018), the simulated CMEs did not include the magnetic field. Furthermore, to make simulations more comparable with the previous simulation results of Pomoell & Poedts (2018), AMR was not used in the present CME simulations. The use of AMR would have allowed better capture of CME structures.

Figure 18 shows the predicted distributions of the solar wind speed for 27 June 2015 (11:05 UT) obtained using (a) the standard GONG map and (b) the zero-point corrected map as discussed above and shown in Figure 17. The results of Figure 18a corresponds to Figure 8 from Pomoell & Poedts (2018) (note, the distance in their simulations is 2 AU). In Figure 18a, contours of the solar wind speed are given in the ecliptic xy-plane, as well as the plane perpendicular to the ecliptic, with the Earth at x = 1 AU, y = 0, x = 0. The CME, which is delivering a glancing blow to the Earth, erupted on 25 June 2015 (10:51 UT), and was the last event in the CME series.

Fig. 18 Solar wind speed obtained with (a) standard GONG map and (b) zero point corrected GONG map for 27 June 2015 (11:05 UT).

Figure 18 reveals one of the main differences between the proposed new SW simulation framework and the EUHFORIA and WSA-ENLIL solar wind – CME simulation codes. While the latter cover latitudes of θ = ±60°, the proposed framework is capable of describing all latitudes, including the poles. A comparison of the current simulation results with those of Pomoell & Poedts (2018) (their Fig. 8) shows that, for this particular CME simulation, there is no significant difference in the predicted arrival times and the CME impact on Earth. In both cases, the CME delivers only a small glancing blow to Earth. However, in the current simulations of this case, it can be seen that the CME propagation is influenced by the high-speed stream on the west side (i.e., the solar wind stream that emanates from the Sun around x = 0 and y > 0). This stream is predicted to have a lower speed in the results of Pomoell & Poedts (2018). Furthermore, the high-latitude northern solar wind stream, depicted in the region x > 0 and z > 0 of Figure 18a, is not visible in the previous simulations of Pomoell & Poedts (2018), possibly due to the ±60° restrictions in latitude of their computational domain. Finally, as can be seen from the results of Figure 18, the use of the zero-point corrected GONG synoptic maps in the PFSS-SCS modelling of the coronal field produces noticeable changes in the coronal hole properties and solar wind speed distributions. While the differences in the polar regions are obvious, it should be noted that there are also differences in the predicted solar wind speed in the equatorial plane.

A final set of CME simulation results are given in Figure 19. In the figure, the observed and forecasted (a) solar wind speed, and (b) particle number density at Earth are compared for the period 19–28 June, 2005. The time period shown captures the signatures of all five CMEs. The results for the solar wind speed and density can also be compared with those of Figures 9 and 10 from Pomoell & Poedts (2018), respectively. Four sets of simulated results are shown in Figure 19. The first two sets of predicted results were obtained with identical parameters for the background solar wind as those used for the results shown in Figure 18 except that, in one case, the standard GONG map (S-Map) was used and, in the other case, the zero-point corrected GONG map (ZPC-Map) was used. These results are denoted with solid lines in Figure 19. The third and fourth set of results were again obtained identical parameters and the same two synoptic maps, but with a modified expression for the plasma density (S-Map MD and ZPC-Map MD). The latter are denoted with dashed lines in Figure 19. Concentrating on the first set of two simulation results (solid lines), it can be seen from Figure 19a that there is very good agreement between observed and forecasted solar wind speed for the first three CMEs. As in Pomoell & Poedts (2018), the results for the first CME, CME-1, show a gradual increase in the plasma speed at Earth, while the observed solar wind speed is indicative of a shock. The results for the fifth CME, CME-5, shown also in Figure 18, are also in agreement with the results of Pomoell & Poedts (2018). However, both set of simulated results, those obtained using the proposed numerical framework described herein and those from Pomoell & Poedts (2018), show late arrival of this CME at Earth and lower wind speed than that observed. While the forecasted CME arrivals obtained with the standard GONG map (S-Map results), particularly for CME-2 and CME-3, show better agreement with the observed CME arrivals than those of the previous EUHFORIA simulations, in general, the current S-Map and ZPC-Map results obtained here are quite similar to each other and those of EUHFORIA for all CMEs except for CME-4. In the present simulations, CME-4 arrives late at Earth with a much slower speed than was observed. In contrast, the speed associated with this CME is significantly overestimated in the predictions of Pomoell & Poedts (2018). This disagreement is not entirely surprising. For example, one of the main challenges in CME forecasting is appropriate choice of CME parameters. Some of these parameters, which are estimated from coronagraphs, such as the CME mass, are based on a number of assumptions and these can have large uncertainties (see, e.g., Vourlidas et al., 2000). As already mentioned, the simulations shown herein have used estimated CME masses from the SOHO LASCO CME catalog, while Pomoell & Poedts (2018) assumed the same value of the density for all five CMEs. CME-4 has the lowest mass in the CME series considered here and increasing the mass of this CME in the simulations (not shown here), can improve the agreement between observed and modelled solar wind speeds. Furthermore, while adjusting the CME parameters can improve simulation results, one of the main factors which limits forecast accuracy is associated with the use of a semi-empirical approach in the coronal domain to specify properties of the background solar wind in terms of the PFSS-SCS derived coronal magnetic field. As a baseline configuration, the proposed framework uses WSA-type solar wind speed relation as given in equation (13). The empirical coefficients of this relation are the same as those used to validate a recently proposed semi-empirical solar wind speed forecast model which utilize kinematic propagation of solar wind streams Reiss et al. (2016). It should also be mentioned that these empirical coefficients have been tuned based on standard GONG magnetograms. The tuning of the empirical coefficients and performance of the solar wind speed forecasts with zero-point corrected maps deserves further investigation, particularly since it is felt that these synoptic maps should be adopted in solar wind modelling. Additionally, the plasma density and temperature are evaluated herein using equations (14) and (15), respectively, which are based on a model derived from Helios solar wind observations (Hayashi et al., 2003; Kataoka et al., 2009; Shiota et al., 2014). Finally, while it is noted that that Pomoell & Poedts (2018) have used the same form of the WSA relation considered herein, different empirical coefficients were used. Furthermore, they assumed constant kinetic energy density (i.e., the plasma number density is inversely proportional to the square of the solar wind speed) and plasma thermal pressure at the inner boundary of the global MHD model.

Fig. 19 Comparisons of observed and simulated (a) solar wind speed and (b) particle number density at Earth for 19–28 June 2015. The time period shown includes the signatures of all five simulated CMEs. Four sets of simulated results are shown. The first two sets of predicted results were obtained with identical parameters for the background solar wind as those used for the results shown in Figure 18 except that, in one case, the standard GONG map (S-Map) was used and, in the other case, the zero-point corrected GONG map (ZPC-Map) was used. The third and fourth set of results were again obtained identical parameters and the same two synoptic maps, but with a modified expression for the plasma density (S-Map MD and ZPC-Map MD).

As the background solar wind plasma properties can directly affect CME propagation, improvements to the solar corona model and solar wind description will be the subject of future research efforts. For example, from Figure 19b, it is evident that the predicted plasma density of the background solar wind is overestimated compared to the observed density (see, e.g., 19–20 June) for simulations S-Map and ZPC-Map, which suggests that the value of the number density given by equation (14) should be reduced. To illustrate changes in the simulation results introduced by altering this empirical relation for the plasma density, additional simulations were performed using both the standard and zero-point corrected GONG synoptic maps along with a modified plasma density, cases S-Map MD and ZPC-Map MD, respectively, where 30 km/s was added to the solar wind speed in the expression of equation (14), (i.e., V_SW → V_SW + 30) so as to reduce the speficied value of the background solar wind number density. This second set of simulation results are also shown in Figure 19. As can be seen, the predicted arrival times of the five CMEs with this modified background density, which is now closer to the observed density, are now all shorter corresponding to faster CME propagation speeds. In particular, CME-4 and CME-5 are now predicted to arrive more than 6 h earlier.

The results presented in this section illustrate the forecasting potential of the proposed solar wind – CME numerical framework. Nevertheless, there are aspects, particularly related to the semi-empirical model of the solar corona and the solar wind, which deserve further investigation and development. Future studies are planned to investigate and improve the physics of the baseline models and numerical techniques described herein, as well as to assess the capabilities of the framework as a forecasting tool at the Canadian Space Weather Forecast Centre.

6 Summary and conclusions

As discussed in the introduction, SW and its effects on technology are now widely recognized as a hazard to modern society and there are many on-going research activities in this field. The adverse effects of space weather are recognized by both industry stakeholders as well as government organizations. For this reason, a number of operational forecast centers have been established to provide information about SW conditions. Furthermore, the development and further improvement of numerical models with enhanced forecasting capabilities for use in SW operations represent an important component of the on-going research effort.

Due to computational constraints, common approaches used in the current generation of solar wind – CME forecast models are data driven and make use of simplified and semi-empirical modelling to describe the solar corona combined with a global MHD model for the interplanetary propagation of the solar wind and CMEs (Odstrcil et al., 2008; Shiota & Kataoka, 2016; Pomoell & Poedts, 2018). In the current study, such an approach is considered and a new integrated data-driven solar wind – CME numerical framework for potential use in SW forecasting has been developed and described. In the proposed framework, the PFSS (Altschuler & Newkirk, 1969; Schatten et al., 1969) and SCS (Schatten, 1971) models are used to derive baseline predictions of the coronal magnetic field from GONG magnetogram observations and additional empirical relations, including the WSA solar wind speed relation (MacNeice, 2009), and are subsequently applied to associate solar wind properties with open magnetic field lines. To improve the accuracy of the PFSS model, particularly in the polar regions, re-meshing of the magnetogram synoptic maps is utilized (Tóth et al., 2011). This baseline coronal component of the proposed framework is augmented with an approximate CME model based on a spheromak pancake flux rope model (Shiota & Kataoka, 2016) so as to provide the inner boundary conditions to the global MHD model of the interplanetary subdomain. The proposed coronal model presented here can be driven by a single magnetogram synoptic map or using a series of maps and accepts standard, zero-point corrected, as well as ADAPT GONG synoptic maps.

The global MHD model used in the proposed solar wind – CME framework is based on the solution of the ideal MHD equations. The latter can be solved in either the inertial or the Sun’s co-rotating frame of reference using a combination of advanced high-fidelity numerical methods which are not available in other solar wind simulation codes and models and offer significantly enhanced capabilities. A second-order-accurate upwind finite-volume scheme (Ivan et al., 2011, 2013, 2015) is used to solve the ideal MHD equations on cubed-sphere meshes (Ivan et al., 2011, 2013, 2015) and this spatial discretization scheme is combined with a highly-scalable and efficient parallel block-based anisotropic AMR technique (Williamschen & Groth, 2013; Freret & Groth, 2015; Freret et al., 2019). The use of cubed-sphere grids avoids the singularities of spherical grids at the poles and thereby readily allows full coverage of 3D space, including high-latitude polar regions. The anisotropic block-based AMR scheme permits dynamic and local refinement of the mesh while allowing efficient implementations of the MHD solution algorithm on high-performance computers. For this reason, the proposed framework offers significantly increased grid resolution while maintaining the computational cost at reasonable levels.

Example simulation results for both steady and unsteady solar winds corresponding to a relatively quiet period of solar activity, as well as additional comparisons to the previous results of Pomoell & Poedts (2018) for a CME event, have also been described herein in order to illustrate the potential SW forecasting capabilities of the proposed data-driven solar wind – CME computational framework. The simulation results demonstrate that the anisotropic AMR does an efficient job in accurately capturing the solar wind current sheet by refining extensively only across the sheet while keeping the grid resolution low along the current sheet surface. The refinement normal to the current sheet minimizes the numerical diffusion which can cause unphysical reconnection of the oppositely directed magnetic field lines associated with the current sheet. The comparisons with Pomoell & Poedts (2018) have also shown the similarities of the current solar wind – CME simulation results with those of the previous study for the considered event. Although there are some differences in the predicted background solar wind, the CME propagation characteristics appear to be very similar. Differences in the solar wind predictions can be attributed to differently tuned empirical relations for the solar wind properties and to the fact that the simulations of Pomoell & Poedts (2018) do not cover higher latitudes. Finally, differences in the coronal hole and solar wind solutions predicted by the PFSS and SCS models using standard and zero-point corrected GONG maps has been investigated. The results suggest that instead of using standard GONG maps, zero-point corrected maps would appear to be more appropriate for coronal field and solar wind modelling.

The proposed new integrated data-driven solar wind – CME framework will provide a basis for future follow-on research. Reliable data and empirical relations that perform well over time are important, particularly for operational forecast models where the updates are usually less frequent. Although it is expected that many solar wind and CME cases, particularly before 2013 (Nikolić, 2019), can be well described using standard GONG maps, to obtain more consistent results over time, the use of zero-point corrected as well as ADAPT synoptic maps in SW simulations should be explored further. As the predicted background solar wind depends intimately on the synoptic maps, studies similar to those of Hinterreiter et al. (2019), MacNeice (2009) and Wold et al. (2018), should be undertaken to assess performances of solar wind models with the GONG zero-point corrected and ADAPT maps. In particular, improvements in the semi-empirical solar corona component of the proposed framework, which provides the boundary conditions to the global MHD model, will be explored, including an assessment of the influence of the various forms of synoptic maps on solar wind predictions as well as the more optimal tuning of the empirical relations used to relate the coronal magnetic field and solar wind parameters. Additionally, the application of data assimilation techniques (Lang et al., 2017) will be explored for use with the proposed framework.

The solar wind – CME numerical framework will also provide a platform for testing new research developments in SW forecasting. It will also be more thoroughly tested as a forecasting tool at the Canadian Space Weather Forecast Centre. Since MHD-based CME forecasting models of probable arrival and impact use parameters which are derived from CME observations, such as the speed, direction and size, as well as other parameters which are more difficult to estimate, such as the CME mass and magnetic field, several runs with different parameters are typically performed as part of the forecasting (see, e.g., Steenburgh et al., 2014). Using different forecasting codes, including the solar wind – CME simulation framework presented here, which run independently with different set of parameters, will help efforts to improve SW forecasting and will benefit the global SW forecasting community.

Acknowledgments

This research was partially supported by the Canadian Space Agency (CSA) through a grant under the GO Canada, Science and Applications Program and was performed as part of Natural Resources Canada’s Public Safety Geoscience program. Additionally, the General Purpose Science Cluster (Government of Canada) and the Compute/Calcul Canada and SciNet High Performance Computing Consortium at the University of Toronto, funded by the Canada Foundation for Innovation (CFI) and the Province of Ontario, Canada, are both acknowledged for providing computing resources. The research utilizes data obtained by the Global Oscillation Network Group (GONG) Program (http://gong.nso.edu/data/magmap/), managed by the National Solar Observatory, which is operated by AURA Inc. under a cooperative agreement with the U. S. National Science Foundation. The data were acquired by instruments operated by the Big Bear Solar Observatory, High Altitude Observatory, Learmonth Solar Observatory, Udaipur Solar Observatory, Instituto de Astrofísica de Canarias, and Cerro Tololo Interamerican Observatory. Data from the Advanced Composition Explorer (ACE, http://www.srl.caltech.edu/ACE/) were also used and are acknowledged. Finally, the authors would like to thank the reviewers for their time and efforts reviewing the manuscript as well as their useful comments and suggestions, which have helped to improve the final version of the paper. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Figures

Fig. 1 Schematic diagram showing the subdomains and mathematical models used in each region in the proposed Sun-to-Earth data-driven solar wind model. The PFSS model is used to describe the coronal magnetic field from the Sun’s surface up to a spherical surface where magnetic field lines are forced to open. From this surface, which is usually placed at ≈2.5 R0, the SCS model describes the coronal field up to 20–30 R0, where the inner MHD boundary conditions are applied. A global MHD model is then used to describe the solar wind plasma and its propagation through interplanetary space in the inner heliosphere.

In the text

	Fig. 2 Hexahedral cell at grid location i, j, k showing face normals.
In the text

	Fig. 3 Schematic diagram of binary tree and corresponding 3D grid-block structure of computational mesh after several levels of refinement.
In the text

	Fig. 4 Schematic diagram of heterogeneous grid-block structure showing (a) the grid block of interest, marked ID, and the associated neighbouring blocks; and (b) the ghost cells of block ID taken directly from the neighbouring blocks.
In the text

Fig. 5 Portion of the initial cubed-sphere mesh used in the simulation of both the steady and unsteady solar wind. The inner boundary of the global MHD subdomain is shaded in green. The thicker black lines are the grid block edges and the thinner red lines are the edges of the computational cells. Each block contains 8 × 8 × 8 = 512 cells and the total number of grid blocks is 144 for a total number of 144 × 512 = 73,728 computational cells.

In the text

	Fig. 6 Surface plots of the estimated solar wind properties at r = 25 R0 for the steady solar wind simulation obtained using the standard GONG map for 12 February 2007 (11:54 UT). The central meridian and sub-Earth locations are denoted with dash-dot and dash-dash line, respectively.
In the text

	Fig. 7 Convergence history of density residual for steady solar wind simulation showing convergence of Newton method on the sequence of anisotropic AMR meshes.
In the text

	Fig. 8 Predicted contours of |B|r2 at the outer boundary r = 225 R0 (≈1.05 AU) of the computational domain after 0, 2, 4, and 7 refinements. The black lines are the block boundaries.
In the text

	Fig. 9 Predicted contours of |B|r2 in the xz-plane (top) and yz-plane (bottom) for the initial (73,728 cells) (left) and final (23,449,088 cells) (right) meshes showing capturing of the current sheet. The black lines are the block boundaries.
In the text

	Fig. 10 Portion of the xy-plane (ecliptic plane) showing the predicted contours of |B|r2 and magnetic field lines on initial (73,728 cells) (left) and final (23,449,088 cells) (right) meshes.
In the text

	Fig. 11 Predicted distributions of the magnitude of the solar wind velocity obtained on the final mesh (23,449,088 cells) after 7 successive anisotropic adaptive mesh refinements. The black lines are the block boundaries. The green-coloured sphere, E, corresponds to the approximate location of Earth.
In the text

	Fig. 12 (a) Predicted contours of |B|r2 in the ecliptic xy-plane obtained on the final mesh (23,449,088 cells) along with (b) the isosurface of the current sheet lying on the northern side of the plane. The black lines correspond to the grid block boundaries.
In the text

	Fig. 13 Surface plots of the estimated radial component of the magnetic field obtained at the inner boundary r = 25 R0 obtained using the PFSS and SCS models with synoptic maps for 22 and 26 January 2007 (5:54 UT), and the corresponding solar wind speed obtained using the WSA model. The central meridian and sub-Earth locations are denoted with dash-dot and dash-dash line, respectively.
In the text

	Fig. 14 Predicted contours of |B|r2 in the xy-, yz- and xz-planes, and current sheet isosurfaces on adapted meshes showing the grid block boundaries (not individual cells) on 26 January (5 pm UT) (17,671,680 cells) and 31 January (11 pm UT) (31,553,536 cells) 2007, respectively, for the unsteady solar wind simulation. The black lines are the block boundaries.
In the text

	Fig. 15 Solar wind speed on the ecliptic plane at various time instances. The dates and approximate times are (a) 23 January (12 pm UT) (5,569,024 cells), (b) 25 January (1 pm UT) (13,434,880 cells), (c) 27 January (10 am UT) (9,246,720 cells), and (d) 29 January (10 am UT) (8,822,272 cells). The brown spot indicates the approximate location of the Earth.
In the text

	Fig. 16 Comparison of (a) solar wind speed, (b) plasma density, and (c) magnetic field strength at Earth predicted by the proposed SW framework to ACE satellite observations for the unsteady solar wind simulation.
In the text

	Fig. 17 Predicted coronal hole (left) and solar wind speed (right) maps of the PFSS and SCS models obtained using (a) the standard GONG map and (b) the zero-point corrected GONG map for 25 June 2015 01:04 UT.
In the text

	Fig. 18 Solar wind speed obtained with (a) standard GONG map and (b) zero point corrected GONG map for 27 June 2015 (11:05 UT).
In the text

Fig. 19 Comparisons of observed and simulated (a) solar wind speed and (b) particle number density at Earth for 19–28 June 2015. The time period shown includes the signatures of all five simulated CMEs. Four sets of simulated results are shown. The first two sets of predicted results were obtained with identical parameters for the background solar wind as those used for the results shown in Figure 18 except that, in one case, the standard GONG map (S-Map) was used and, in the other case, the zero-point corrected GONG map (ZPC-Map) was used. The third and fourth set of results were again obtained identical parameters and the same two synoptic maps, but with a modified expression for the plasma density (S-Map MD and ZPC-Map MD).

In the text