Issue 
J. Space Weather Space Clim.
Volume 11, 2021



Article Number  8  
Number of page(s)  25  
DOI  https://doi.org/10.1051/swsc/2020068  
Published online  28 January 2021 
Research Article
An integrated datadriven solar wind – CME numerical framework for space weather forecasting
^{1}
University of Toronto Institute for Aerospace Studies, Toronto, ON M3H 5T6, Canada
^{2}
Canadian Hazards Information Service, Natural Resources Canada, Ottawa, ON K1A 0E7 Canada
^{3}
Department of Applied Mathematics, University of Waterloo, Waterloo, ON N2L 3G1, Canada
^{*} Corresponding author: ljubomir.nikolic@canada.ca
Received:
19
August
2020
Accepted:
18
November
2020
The development of numerical models and tools which have operational space weather potential is an increasingly important area of research. This study presents recent Canadian efforts toward the development of a numerical framework for SuntoEarth simulations of solar wind disturbances. This modular threedimensional (3D) simulation framework is based on a semiempirical datadriven approach to describe the solar corona and an MHDbased description of the heliosphere. In the present configuration, the semiempirical component uses the potential field source surface (PFSS) and Schatten current sheet (SCS) models to derive the coronal magnetic field based on observed magnetogram data. Using empirical relations, solar wind properties are associated with this coronal magnetic field. Together with a coronal mass ejection (CME) model, this provides inner boundary conditions for a global MHD model which is used to describe interplanetary propagation of the solar wind and CMEs. The proposed MHD numerical approach makes use of advanced numerical techniques. The 3D MHD code employs a finitevolume discretization procedure with limited piecewise linear reconstruction to solve the governing partialdifferential equations. The equations are solved on a bodyfitted hexahedral multiblock cubedsphere mesh and an efficient iterative Newton method is used for timeinvariant simulations and an explicit timemarching scheme is applied for unsteady cases. Additionally, an efficient anisotropic blockbased refinement technique provides significant reductions in the size of the computational mesh by locally refining the grid in selected directions as dictated by the flow physics. The capabilities of the framework for accurately capturing solar wind structures and forecasting solar wind properties at Earth are demonstrated. Furthermore, a comparison with previously reported results and future space weather forecasting challenges are discussed.
Key words: solar wind / coronal mass ejections / space weather forecasting / MHD modelling
© N.M. Narechania et al., Published by EDP Sciences 2021
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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 groundbased 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. Highspeed 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 SuntoEarth 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 highspeed 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 largescale 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 physicsbased global magnetohydrodynamics (MHD) models. For example, Pizzo (1982) studied the structure of the solar wind between 35 R_{0} and 1 AU using an idealized model of the inner heliosphere. Here, R_{0} is the radius of the Sun. Usmanov (1993) developed a fully threedimensional (3D) steadystate 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_{0} to 30 R_{0} using photospheric field observations as boundary conditions and compared results with SOHO, Ulysses and WIND data. A parallel blockadaptive 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 threedimensional 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 followon studies, Roussev et al. (2003) modelled the coronaheliosphere 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 fluxropedriven 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, multifluid MHD, and radiative MHD models (Tóth et al., 2012) and van der Holst et al. (2014) implemented a twotemperature MHD model wherein lowfrequency 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_{0} 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 sixcomponent overset grid to study the background solar wind from Sun to Earth during CR 1911 using lineofsight photospheric field observations and validated their MHD model using SOHO and WIND observations. Feng et al. (2012) later added isotropic blockbased adaptive mesh refinement (AMR) capabilities to this MHDbased 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 MPIAMRVAC software, an open source toolkit for parallel, blockadaptive global MHD simulations of solar and nonrelativistic astrophysical plasmas (a relativistic counterpart of this software, BHAC, was also developed by Porth et al. (2017) for solving the equations of ideal generalrelativistic magnetohydrodynamics and studying astrophysical phenomena such as black holes). The MPIAMRVAC 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 MPIAMRVAC software has also been used for the study of magnetic reconnection of solar flares and to simulate the transAlfvé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 solutiondependent AMR with an ultrahighresolution to capture current sheets and smallscale 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 ongoing 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 researchtooperations transfer are identified as some of the key challenges (AraujoPradere, 2009; Merceret et al., 2013; Folini, 2018). Furthermore, the development of operationally oriented physicsbased numerical models, in particular SuntoEarth simulations, requires significant computational effort. For this reason, a large number of existing models used in SW forecast operations are empirical and/or semiempirical 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 physicsbased numerical approach requires information about the Sun’s magnetic field. Due to the favorable signaltonoise 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 socalled 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_{0} 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 SuntoEarth solar wind numerical modelling and to avoid the computationally intensive calculations associated with determining the coronal magnetic field, a combined twomodel approach is often adopted in a number of SW prediction frameworks in which a semiempirical datadriven 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_{0} 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 superAlfvénic speed by the time it reaches this boundary. Such datadriven 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 PFSSSCS 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 PFSSSCS 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 datadriven solar wind model. The global MHD modelling in ENLIL for the simulation of solar wind structures uses a fixed, uniform or nonuniform 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_{0} and by the SCS model in the region beyond 2.5 R_{0} out to 21.5 R_{0}. Baker et al. (2013) employed the WSAENLIL 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 datadriven 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_{0} 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 largescale datadriven 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 largescale 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 semiempirical datadriven 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.
SuntoEarth solar wind modelling efforts have also been recently undertaken in Canada. Canada is amongst the countries most affected by SW effects, due to its highlatitude location (see, e.g., Lam, 2011; Boteler, 2019), and there are ongoing 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 datadriven 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 empiricallybased expressions for solar wind flow properties. For the global MHD description of the heliosphere, a secondorderaccurate upwind finitevolume scheme (Ivan et al., 2011, 2013, 2015; Susanto et al., 2013) is used to solve the governing ideal MHD equations on a cubedsphere mesh (Ivan et al., 2011, 2013, 2015) and this finitevolume scheme is combined with the highlyscalable and efficient parallel blockbased anisotropic AMR technique developed previously by Williamschen & Groth (2013), Freret & Groth (2015) and Freret et al. (2019). As such, the proposed datadriven 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 zeropoint corrected QuickReduce synoptic maps, as well as Air Force Data Assimilative Photospheric Flux Transport (ADAPT) GONG maps. Additionally, the combination of an advanced highfidelity finitevolume scheme with parallel anisotropic blockbased AMR in the global MHD modelling allows local solutiondependent adaptation of the mesh and affords significantly increased mesh resolution. Furthermore, as mentioned above, a cubedsphere 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 highlatitude polar regions. Lastly, the solar wind simulations can be performed in either the inertial or the Sun’s corotating frame of reference, offering considerable flexibility for SW modelling.
The organization of the remainder of this paper is as follows. The proposed domains and datadriven 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 secondorder finitevolume spatialdiscretization and timemarching schemes which are used to solve the 3D form of the ideal MHD equations, and the parallel blockbased 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 steadystate background solar wind and an unsteady simulation based on timeevolving 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 datadriven 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 SuntoEarth 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_{0} from the Sun, where the solar wind speed is both supersonic and superAlfvé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 SuntoEarth datadriven 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 R_{0}, the SCS model describes the coronal field up to 20–30 R_{0}, 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 semiempirical 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_{0}) and the socalled “source surface” (r = R_{ss}), which means that the magnetic field can be expressed as the gradient of a scalar potential, Ψ, as in
Along with the solenoidal property of the magnetic field, ∇·B = 0, equation (1) leads to the Laplace equation for Ψ
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_{0} ≤ r ≤ R_{ss} can be expressed as
where 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
The coefficients g_{nm} and h_{nm}, can be obtained using equation (4) for the case r = R_{0}, and employing the orthogonality of Legendre polynomials
The following expressions for g_{nm} and h_{nm} can be obtained:
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_{0}, θ, ϕ), 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_{0} ≤ r ≤ R_{ss} can then be obtained using equations (4)–(6).
2.1.2 Schatten current sheet (SCS) model
The PFSS model uses the currentfree approximation of the solar corona in the region R_{0} ≤ 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 nonradial 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 reoriented 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 reoriented field at r = R_{cp} with the potential field solution for r ≥ R_{cp} (see, e.g., Schatten, 1971; Nikolić, 2017), and thus
The role of the PFSS field reorientation 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 reorientation 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 sphericallyshaped 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_{0} and 30R_{0} (see, e.g., Usmanov, 1993; MacNeice et al., 2018). In this range, the solar wind speed has already reached its asymptotic value and is superAlfvénic. For example, Odstrcil et al. (2004a) investigated the coupling of coronal and heliospheric simulation models with the boundary located at both 25 R_{0} and 50 R_{0} and their findings justifies the use of an interface boundary located at 25 R_{0}.
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
as well as the angular separation, θ_{b}, between an open magnetic field line footpoint 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 WSAtype models is given by
where a_{1}–a_{8} 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_{0} (Hayashi et al., 2003; Kataoka et al., 2009). Modified versions of these relations for a boundary R_{inner} = 25 R_{0} 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
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_{B}T, where k_{B} = 1.38 × 10^{−23} J K^{−1} is the wellknown 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 largescale 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 empiricallybased 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 socalled cone model of Xie et al. (2004) and Zhao et al. (2002) is an example of a datadriven 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 timedependent 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 spheromaktype 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^{1}, based on several other approximations.
In the proposed integrated datadriven 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 physicsbased 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, quasineutral, inviscid, ideal gas and, as noted in the introduction, are commonly used in the modelling of space plasmas. The equations of ideal MHD in nondimensional weak conservative form, for a rotating frame, are given by
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 noninertial 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
where p is the thermal pressure of the plasma and the total energy e is given by
with γ = C_{p}/C_{v} being the ratio of specific heats. Equations (16)–(19) can be expressed in matrixvector form as
where the vector of conserved variables, U, and the flux dyad, F, are given by
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
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 finitevolume 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 Galileaninvariant and also therefore hyperbolic with an eigenstructure that is not degenerate. The application of the proposed finitevolume 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 nondimensional solution variables of the ideal MHD equations can be related to their dimensional counterparts by , u = /a_{0}, , B = /B_{0}, , and r = /l_{0}. Here, ρ_{0} = m_{p}/cm^{3} where m_{p} = 1.672 × 10^{−27} kg is the mass of a proton. The length scale l_{0} is taken to be the radius of the Sun, l_{0} = R_{0} = 6.96 × 10^{8} m. The velocity scale is given by a_{0} = l_{0}/τ_{0} = 193.333 × 10^{3} m/s where the time scale is taken to be τ_{0} = 1, h = 3600 s. The value used to normalize the magnetic field is given by B_{0} = = 8.8642 × 10^{−9} T where μ_{0} = 4π × 10^{−7} H/m is the magnetic permeability of free space. The value used to normalize the pressure is given by p_{0} = ρ_{0} = 6.252 × 10^{−11} Pa. The angular velocity of the Sun is taken to be = (2π/27.27) rad/day = 2.67 × 10^{−6} rad/s .
3 Numerical solution methods
3.1 Synoptic maps of the photospheric field
Global Oscillation Network Group (GONG) magnetogram synoptic maps^{2} 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 realtime synoptic maps of the photospheric field (Hill, 2018). For this reason, GONG maps are frequently used in operationally oriented SW applications. This includes WSAENLIL 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 zeropoint 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^{−4} 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_{0}, θ, ϕ), which is used by the PFSS model to obtain the coefficients g_{nm} and h_{nm} by means of equation (8).
Both standard and zeropoint corrected maps can be used in the proposed datadriven 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 WSAENLIL 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 zeropoint corrected GONG maps promise to reduce these false variations. The zeropoint error is reduced using a software system that compares observations from different sites and in time, better control of the liquidcrystal 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 WSAENLIL simulation model was recently upgraded to allow for the use of both standard and zeropoint corrected maps. In this study, the differences between simulation results obtained using the standard and zeropoint 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 remeshing 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. Remeshed 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 remeshed B_{r}(R_{0}, θ_{i}, ϕ_{j}) synoptic map.
Following remeshing, the coefficients g_{nm} and h_{nm} are then determined by using a discretized form of equation (8) which can be expressed as
where ϵ_{1} = = 1/2, and ϵ_{i} = 1 for i ≠ 1, N_{θ}, and w_{i} are Clenshaw–Curtis weights given by
with H = (N_{θ} − 1)/2, = 1/2, and = 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 finitenumber 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 and in equations (4)–(6), the following relations are used:
along with the recursion relation
and (Altschuler et al., 1977; Nikolić, 2017).
3.3 Numerical solution of SCS model
A leastsquares 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 reorientation Schatten (1971) and thereby determine values for the coefficients and . Denoting the reoriented field at R_{cp} as (θ_{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 finitenumber 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, and , of equations (9)–(11), are then determined in the leastsquares approach by minimizing the sum of squared residuals, F, defined by
such that = 0 and = 0. In equation (30), α_{knm} and β_{knm} are given by
For each (m, n), the minimization of equation (30) results in a linear system of equations for and which can then be solved following some linear algebra. The solutions for and 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_{0} ≤ 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_{0} and the inner boundary of the global MHD model is placed at r = R_{inner} = 25 R_{0}. Equations (13)–(15) are then used to provide the solar wind plasma properties at R_{inner}. The WSA relation of equation (13) with a_{1} = 250, a_{2} = 875, a_{3} = 0.2, a_{4} = 1, a_{5} = 0.8, a_{6} = 2.6, a_{7} = 1.25 and a_{8} = 2.5 is used. For the properties needed for the open magnetic field lines, f_{s} and θ_{b}, a secondorder RungeKutta 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 PFSSSCS modelling on to a (θ, ϕ, t) grid representing the inner boundary of the MHD model. This can be done for both noninertial corotating and nonrotating inertial frames which account for the solar rotation. For calculations performed in the Sun’s corotating 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 corotating 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 nonrotating case to account for the timedependent nature of the solar wind plasma properties at the inner boundary (Hayashi, 2012).
3.5 Solutionadaptive upwind finitevolume scheme for ideal MHD equations
3.5.1 Spatial discretization and semidiscrete form
In the proposed Godunovtype upwind finitevolume scheme (Godunov, 1959; Toro, 2013) for solving the ideal MHD equations of equation (22), the spatial discretization is accomplished herein by applying a secondorder cellcentered finitevolume 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 semidiscrete form:
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, represents the number of faces for the cell and = 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 piecewiselinear leastsquares reconstruction is used for calculating the primitive flow variables at each cell face. A Riemannsolverbased function (Godunov, 1959; Toro, 2013), namely the socalled 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
Steadystate or timeinvariant solutions of the semidiscrete form of the governing equations given by equation (32) satisfy a large coupled system of nonlinear algebraic equations for all computational cells in the mesh given by
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
where U^{(n)} is the solution change associated with the nth iteration and successively improved estimates of the solution to the nonlinear equations satisfy
and where J = ∂R/∂U is the Jacobian of the solution residual. Given an initial estimate of the solution, U^{(n=0)}, 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 userdefined tolerance. A value of ϵ = 10^{−8} 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 nonsymmetric and sparse and also typically very large. As such, Krylov subspace iterative methods can be very effective in determining their solution. A “matrixfree” or “Jacobianfree” 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 lowerupper (BILU) local preconditioning (Northrup & Groth, 2013).
3.5.3 Explicit timemarching scheme for unsteady solar wind flows
For timedependent or unsteady solar wind flows, solutions of the coupled system of nonlinear ordinary differential equations (ODEs) represented by the semidiscrete form of the governing equations given by
are sought. A standard explicit, twostage, secondorder accurate, Runge–Kutta timemarching 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 blockbased AMR
The hexahedral computational cells of the global MHD grid are contained within multiblock bodyfitted structured grids that permit a general unstructured connectivity of the grid blocks thereby readily allowing the use of cubedsphere grids (Ivan et al., 2011, 2013, 2015). The proposed multiblock grid structure also readily facilitates automatic, solutiondependent, 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 highperformance 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 solutiondependent refinement criteria. Coarsening of the mesh is accomplished by reversing this process. The particular form of the anisotropic blockbased 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 gridblock 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 selfsimilar nature of grid blocks, leads to efficient and highly scalable implementations of the combined finitevolume scheme and AMR procedure on highperformance 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 gridblock 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 leastsquares 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 cellaveraged values from coarser to finer cells and restriction from finer to coarser cells, when exchanging information between blocks. When there are nonconforming cells or socalled “hanging nodes” where neighbouring blocks meet, the fluxes through the nonconforming 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 finitevolume 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 gridblock 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 datadriven solar wind – CME prediction framework are demonstrated by considering two representative simulations. In the first example, a steadystate 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 steadystate 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 secondorder, explicit, timemarching scheme described in Section 3.5.3 was used and the predictions of the solar wind speed, density, and magneticfield predictions are compared with available observations. Both simulations were performed in the corotating frame of the Sun and a cubedsphere 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, ∇(Br^{2}), was used as the criterion for the anisotropic mesh adaptation. The weight r^{2} 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^{2} 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 cubedsphere 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_{0}, are shown in Figure 6. Here, the central meridian and subEarth 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 highpressure 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 R_{0} for the steady solar wind simulation obtained using the standard GONG map for 12 February 2007 (11:54 UT). The central meridian and subEarth locations are denoted with dashdot and dashdash line, respectively. 
Using the inflow boundary conditions arising at r = 25 R_{0} shown above and after obtaining the timeinvariant 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 Br^{2} at the outer boundary r = 225 R_{0} (≈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 Br^{2} and the comparison between the structure of the current sheet on xz and yzslices 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 Br^{2} and magnetic field lines on a portion of the ecliptic plane (xyplane) 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 highquality ideal MHD solutions of the solar wind unless a much higher number of grid points were used.
Fig. 8 Predicted contours of Br^{2} at the outer boundary r = 225 R_{0} (≈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 Br^{2} in the xzplane (top) and yzplane (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 xyplane (ecliptic plane) showing the predicted contours of Br^{2} 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 greencoloured sphere, denoted with an E, represents the approximate location of the Earth. Finally, Figure 12 provides the predicted contours of Br^{2} in the ecliptic xyplane 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 greencoloured sphere, E, corresponds to the approximate location of Earth. 
Fig. 12 (a) Predicted contours of Br^{2} in the ecliptic xyplane 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 timedependent 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 timedependent 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_{0} 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 subEarth 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 highspeed 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 R_{0} 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 subEarth locations are denoted with dashdot and dashdash line, respectively. 
To initiate the unsteady simulation for the sequence of 64 synoptic maps, a steadystate 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 secondorder timemarching 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 timesteps and a Courant–Friedrichs–Lewy number of 0.3 was used throughout the simulation.
Figures 14a and 14b show the predicted contours of Br^{2} 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 cubedsphere 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 browncoloured sphere in the figure represents the position of the Earth at each instance in time. Note that, in the corotating reference frame, the Earth moves in a clockwise direction at the rate of (2π/27.27) rad/day.
Fig. 14 Predicted contours of Br^{2} in the xy, yz and xzplanes, 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 magneticfield 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 magneticfield 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 datadriven 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 retuning the empirical solar wind relations of equations (13)–(15). However, while readjusting 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 highspeed 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 11year solar activity cycle. For example, the GONG standard magnetogram synoptic maps are available from late 2006 and the WSAENLIL 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 highspeed wind, derived using the standard GONG maps exhibits significant deviations after 2013, particularity in highlatitude 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 zeropoint error, using the zeropoint 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 zeropoint 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 zeropoint 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 dashdash line while the dashdot 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_{0})) 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 zeropoint 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 zeropoint 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 zeropoint 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 zeropoint 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 datadriven solar wind – CME numerical framework and to further demonstrate the differences between simulations performed with the standard and zeropoint 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)^{3}, gives these CMEs on June 18 (20:00 UT) at (θ, ϕ) = (11°, −50°) with a speed of 1000 km/s (CME1), June 19 (14:59 UT) at (θ, ϕ) = (−33°, 9°) with a speed of 603 km/s (CME2), June 21 (05:01 UT) at (θ, ϕ) = (7°, 8°) with a speed of 1250 km/s (CME3), June 22 (21:10 UT) at (θ, ϕ) = (14°, 3°) with a speed of 1155 km/s (CME4), and June 25 (10:51 UT) at (θ, ϕ) = (23°, 46°) with a speed of 1450 km/s (CME5). The angular halfwidth 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^{4} were used. As such, the CME masses were set to 18.0 × 10^{12} kg, 6.1 × 10^{12} kg, 9.6 × 10^{12} kg, 4.4 × 10^{12} kg, and 31.0 × 10^{12} kg, respectively, and the radial width of the CMEs was taken to be 2 R_{0}. 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 zeropoint 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 xyplane, 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 WSAENLIL 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 highspeed 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 highlatitude 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 zeropoint corrected GONG synoptic maps in the PFSSSCS 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 (SMap) was used and, in the other case, the zeropoint corrected GONG map (ZPCMap) 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 (SMap MD and ZPCMap 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, CME1, 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, CME5, 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 (SMap results), particularly for CME2 and CME3, show better agreement with the observed CME arrivals than those of the previous EUHFORIA simulations, in general, the current SMap and ZPCMap results obtained here are quite similar to each other and those of EUHFORIA for all CMEs except for CME4. In the present simulations, CME4 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. CME4 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 semiempirical approach in the coronal domain to specify properties of the background solar wind in terms of the PFSSSCS derived coronal magnetic field. As a baseline configuration, the proposed framework uses WSAtype 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 semiempirical 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 zeropoint 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 (SMap) was used and, in the other case, the zeropoint corrected GONG map (ZPCMap) 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 (SMap MD and ZPCMap 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 SMap and ZPCMap, 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 zeropoint corrected GONG synoptic maps along with a modified plasma density, cases SMap MD and ZPCMap 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, CME4 and CME5 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 semiempirical 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 ongoing 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 ongoing 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 semiempirical 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 datadriven 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, remeshing 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, zeropoint 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 corotating frame of reference using a combination of advanced highfidelity numerical methods which are not available in other solar wind simulation codes and models and offer significantly enhanced capabilities. A secondorderaccurate upwind finitevolume scheme (Ivan et al., 2011, 2013, 2015) is used to solve the ideal MHD equations on cubedsphere meshes (Ivan et al., 2011, 2013, 2015) and this spatial discretization scheme is combined with a highlyscalable and efficient parallel blockbased anisotropic AMR technique (Williamschen & Groth, 2013; Freret & Groth, 2015; Freret et al., 2019). The use of cubedsphere grids avoids the singularities of spherical grids at the poles and thereby readily allows full coverage of 3D space, including highlatitude polar regions. The anisotropic blockbased AMR scheme permits dynamic and local refinement of the mesh while allowing efficient implementations of the MHD solution algorithm on highperformance 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 datadriven 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 zeropoint corrected GONG maps has been investigated. The results suggest that instead of using standard GONG maps, zeropoint corrected maps would appear to be more appropriate for coronal field and solar wind modelling.
The proposed new integrated datadriven solar wind – CME framework will provide a basis for future followon 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 zeropoint 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 zeropoint corrected and ADAPT maps. In particular, improvements in the semiempirical 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 MHDbased 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
 Alfvén H. 1957. On the theory of comet tails. Tellus 9(1): 92–96. https://doi.org/10.3402/tellusa.v9i1.9064. [Google Scholar]
 Altschuler MD, Newkirk G. 1969. Magnetic fields and the structure of the solar corona. Sol Phys 9(1): 131–149. https://doi.org/10.1007/BF00145734. [NASA ADS] [CrossRef] [Google Scholar]
 Altschuler M, Levine R, Stix M, Harvey J. 1977. High resolution mapping of the magnetic field of the solar corona. Sol Phys 51(2): 345–375. https://doi.org/10.1007/BF00216372. [NASA ADS] [CrossRef] [Google Scholar]
 AraujoPradere EA. 2009. Transitioning space weather models into operations: The basic building blocks. Space Weather 7(10): S10006. https://doi.org/10.1029/2009SW000524. [Google Scholar]
 Arge CN, Pizzo VJ. 2000. Improvement in the prediction of solar wind conditions using nearreal time solar magnetic field updates. J Geophys Res: Space Phys 105(A5): 10465–10479. https://doi.org/10.1029/1999JA000262. [CrossRef] [Google Scholar]
 Arge CN, Odstrcil D, Pizzo VJ, Mayer LR. 2003. Improved method for specifying solar wind speed near the Sun. In: Solar wind ten, Vol. 679 of American Institute of Physics Conference Series, pp. 190–193. https://doi.org/10.1063/1.1618574. [NASA ADS] [CrossRef] [Google Scholar]
 Arge CN, Henney CJ, Koller J, Compeau CR, Young S, MacKenzie D, Fay A, Harvey JW. 2010. Air force data assimilative photospheric flux transport (ADAPT) model. AIP Conf Proc 1216(1): 343–346. https://doi.org/10.1063/1.3395870. [CrossRef] [Google Scholar]
 Baker DN, Poh G, Odstrcil D, Arge CN, Benna M, et al. 2013. Solar wind forcing at Mercury: WSAENLIL model results. J Geophys Res: Space Phys 118(1): 45–57. https://doi.org/10.1029/2012JA018064. [NASA ADS] [CrossRef] [Google Scholar]
 Boteler DH. 2019. A 21st century view of the March 1989 magnetic storm. Space Weather 17(10): 1427–1441. https://doi.org/10.1029/2019SW002278. [CrossRef] [Google Scholar]
 Butcher J. 1996. A history of RungeKutta methods. Appl Numer Math 20(3): 247–260. https://doi.org/10.1016/01689274(95)001085. [CrossRef] [Google Scholar]
 Detman T, Smith Z, Dryer M, Fry CD, Arge CN, Pizzo V. 2006. A hybrid heliospheric modeling system: Background solar wind. J Geophys Res: Space Phys 111(A7): A07102. https://doi.org/10.1029/2005JA011430. [NASA ADS] [CrossRef] [Google Scholar]
 Einfeldt B. 1988. On Godunovtype methods for gas dynamics. SIAM J Numer Anal 25: 294–318. https://doi.org/10.1137/0725021. [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Feng X, Yang L, Xiang C, Wu ST, Zhou Y, Zhong D. 2010. Threedimensional solar wind modeling from the Sun to Earth by a SIPCESE MHD model with a sixcomponent grid. Astrophys J 723(1): 300. https://doi.org/10.1088/0004637X/723/1/300. [NASA ADS] [CrossRef] [Google Scholar]
 Feng X, Yang L, Xiang C, Jiang C, Ma X, Wu ST, Zhong D, Zhou Y. 2012. Validation of the 3D AMR SIPCESE solar wind model for four Carrington rotations. Sol Phys 279(1): 207–229. https://doi.org/10.1007/s1120701299699. [CrossRef] [Google Scholar]
 Folini D. 2018. Climate, weather, space weather: Model development in an operational context. J Space Weather Space Clim 8: A32. https://doi.org/10.1051/swsc/2018021. [CrossRef] [Google Scholar]
 Freret L, Groth CPT. 2015. Anisotropic nonuniform blockbased adaptive mesh refinement for threedimensional inviscid and viscous flows. In: 22nd AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, AIAA 20152613. https://doi.org/10.2514/6.20152613. [Google Scholar]
 Freret L, Ivan L, De Sterck H, Groth CPT. 2019. Highorder finitevolume method with blockbased AMR for magnetohydrodynamics flows. J Sci Comput 79(1): 176–208. https://doi.org/10.1007/s1091501808441. [CrossRef] [Google Scholar]
 Gao X, Groth CPT. 2010. A parallel solutionadaptive method for threedimensional turbulent nonpremixed combusting flows. J Comput Phys 229(5): 3250–3275. https://doi.org/10.1016/j.jcp.2010.01.001. [CrossRef] [Google Scholar]
 Gao X, Northrup SA, Groth CPT. 2011. Parallel solutionadaptive method for twodimensional nonpremixed combusting flows. Prog Comput Fluid Dyn 11(2): 76–95. https://doi.org/10.1504/PCFD.2011.038834. [CrossRef] [Google Scholar]
 Godunov SK. 1959. Finitedifference method for numerical computations of discontinuous solutions of the equations of fluid dynamics. Mat Sb 47: 271–306. https://hal.archivesouvertes.fr/hal01620642. [Google Scholar]
 Groth CPT, De Zeeuw DL, Gombosi TI, Powell KG. 2000. Global threedimensional MHD simulation of a space weather event: CME formation, interplanetary propagation, and interaction with the magnetosphere. J Geophys Res 105(A11): 25053–25078. https://doi.org/10.1029/2000JA900093. [NASA ADS] [CrossRef] [Google Scholar]
 Hakamada K, Kojima M, Tokumaru M, Ohmi T, Yokobe A, Fujiki K. 2002. Solar wind speed and expansion rate of the coronal magnetic field in solar maximum and minimum phases. Sol Phys 207(1): 173–185. https://doi.org/10.1023/A:1015511805863. [CrossRef] [Google Scholar]
 Hakamada K, Kojima M, Ohmi T, Tokumaru M, Fujiki K. 2005. Correlation between expansion rate of the coronal magnetic field and solar wind speed in a solar activity cycle. Sol Phys 227(2): 387–399. https://doi.org/10.1007/s1120700533047. [CrossRef] [Google Scholar]
 Harten A, Lax PD, van Leer B. 1983. On upstream differencing and Godunovtype schemes for hyperbolic conservation laws. SIAMRev 25(1): 35–61. https://doi.org/10.1137/1025002. [CrossRef] [MathSciNet] [Google Scholar]
 Hayashi K. 2012. An MHD simulation model of timedependent corotating solar wind. J Geophys Res: Space Phys 117(A8): A08105. https://doi.org/10.1029/2011JA017490. [CrossRef] [Google Scholar]
 Hayashi K, Kojima M, Tokumaru M, Fujiki K. 2003. MHD tomography using interplanetary scintillation measurement. J Geophys Res: Space Phys 108(A3): 1102. https://doi.org/10.1029/2002JA009567. [CrossRef] [Google Scholar]
 Hill F. 2018. The Global Oscillation Network Group facility: An example of research to operations in space weather. Space Weather 16(10): 1488–1497. https://doi.org/10.1029/2018SW002001. [CrossRef] [Google Scholar]
 Hinterreiter J, Magdalenic J, Temmer M, Verbeke C, Jebaraj IC, et al. 2019. Assessing the performance of EUHFORIA modeling the background solar wind. Sol Phys 294(12): 170. https://doi.org/10.1007/s1120701915588. [CrossRef] [Google Scholar]
 Hosteaux S, Chané E, Decraemer B, Talpeanu DC, Poedts S. 2018. Ultrahighresolution model of a breakout CME embedded in the solar wind. A&A 620: A57. https://doi.org/10.1051/00046361/201832976. [Google Scholar]
 Howard T. 2011. Coronal mass ejections, SpringerVerlag, New York. https://doi.org/10.1007/9781441987891. [CrossRef] [Google Scholar]
 Ivan L, De Sterck H, Northrup SA, Groth CPT. 2011. Threedimensional MHD on cubedsphere grids: Parallel solutionadaptive simulation framework. In: 20th AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, AIAA 20113382. https://doi.org/10.2514/6.20113382. [Google Scholar]
 Ivan L, De Sterck H, Northrup SA, Groth CPT. 2013. Multidimensional finitevolume scheme for hyperbolic conservation laws on threedimensional solutionadaptive cubedsphere grids. J Comput Phys 255: 205–227. https://doi.org/10.1016/j.jcp.2013.08.008. [CrossRef] [Google Scholar]
 Ivan L, De Sterck H, Susanto A, Groth CPT. 2015. Highorder central ENO finitevolume scheme for hyperbolic conservation laws on threedimensional cubedsphere grids. J Comput Phys 282: 157–182. https://doi.org/10.1016/j.jcp.2014.11.002. [CrossRef] [Google Scholar]
 Kageyama A, Sato T. 2004. “YinYang grid”: An overset grid in spherical geometry. Geochem Geophys Geosyst 5(9): Q09005. https://doi.org/10.1029/2004GC000734. [NASA ADS] [CrossRef] [Google Scholar]
 Kataoka R, Ebisuzaki T, Kusano K, Shiota D, Inoue S, Yamamoto TT, Tokumaru M. 2009. Threedimensional MHD modeling of the solar wind structures associated with 13 December 2006 coronal mass ejection. J Geophys Res: Space Phys 114(A10): A10102. https://doi.org/10.1029/2009JA014167. [Google Scholar]
 Keppens R, Goedbloed JP. 2000. Stellar winds, dead zones, and coronal mass ejections. Astrophys J 530(2): 1036–1048. https://doi.org/10.1086%2F308395. [NASA ADS] [CrossRef] [Google Scholar]
 Keppens R, Meliani Z, van Marle AJ, Delmont P, Vlasis A, van der Holst B. 2012. Parallel, gridadaptive approaches for relativistic hydro and magnetohydrodynamics. J Comput Phys 231(3): 718–744. https://doi.org/10.1016/j.jcp.2011.01.020. [NASA ADS] [CrossRef] [Google Scholar]
 Lam HL. 2011. From early exploration to space weather forecasts: Canada’s geomagnetic odyssey. Space Weather 9(5): S05004. https://doi.org/10.1029/2011SW000664. [Google Scholar]
 Lang MS, Browne PA, van Leeuwen PJ, Owens M. 2017. Data assimilation in the solar wind: Challenges and first results. Space Weather 15: 1490–1510. https://doi.org/10.1002/2017SW001681. [CrossRef] [Google Scholar]
 Linker JA, Mikic Z, Biesecker DA, Forsyth RJ, Gibson SE, Lazarus AJ, Lecinski A, Riley P, Szabo A, Thompson BJ. 1999. Magnetohydrodynamic modeling of the solar corona during Whole Sun Month. J Geophys Res: Space Phys 104(A5): 9809–9830. https://doi.org/10.1029/1998JA900159. [CrossRef] [Google Scholar]
 Lionello R, Mikić Z, Schnack DD. 1998. Magnetohydrodynamics of solar coronal plasmas in cylindrical geometry. J Comput Phys 140: 172–201. https://doi.org/10.1006/jcph.1998.5841. [CrossRef] [Google Scholar]
 Lionello R, Linker JA, Mikic Z. 2009. Multispectral emission of the Sun during the first Whole Sun Month: Magnetohydrodynamic simulations. Astrophys J 690(1): 902–912. https://doi.org/10.1088/0004637X/690/1/902. [NASA ADS] [CrossRef] [Google Scholar]
 Lomax H, Pulliam TH, Zingg DW. 2013. Fundamentals of computational fluid dynamics, SpringerVerlag, New York. https://doi.org/10.1007/9783662046548. [Google Scholar]
 Mackay DH, Yeates AR. 2012. The Sun’s global photospheric and coronal magnetic fields: Observations and models. Living Rev Sol Phys 9: 6. https://doi.org/10.12942/lrsp20126. [CrossRef] [Google Scholar]
 MacNeice P. 2009. Validation of community models: 2. Development of a baseline using the WangSheeleyArge model. Space Weather 7(12): S12002. https://doi.org/10.1029/2009SW000489. [Google Scholar]
 MacNeice P, Jian LK, Antiochos SK, Arge CN, BussyVirat CD, et al. 2018. Assessing the quality of models of the ambient solar wind. Space Weather 16(11): 1644–1667. https://doi.org/10.1029/2018SW002040. [CrossRef] [Google Scholar]
 Manchester WB, Gombosi TI, Roussev I, De Zeeuw DL, Sokolov IV, Powell KG, Toth G, Opher M. 2004. Threedimensional MHD simulation of a flux rope driven CME. J Geophys Res: Space Phys 109(A1): A01102. https://doi.org/10.1029/2002JA009672. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Vourlidas A, Tóth G, Lugaz N, Roussev II, Sokolov IV, Gombosi TI, De Zeeuw DL, Opher M. 2008. Threedimensional MHD simulation of the 2003 October 28 coronal mass ejection: Comparison with LASCO coronagraph observations. Astrophys J 684(2): 1448–1460. https://doi.org/10.1086/590231. [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto H, Sato T (Eds.). 1985. Computer simulation of space plasmas, Terra Scientific Publishing, Tokyo/Dordrecht. [CrossRef] [Google Scholar]
 McGregor SL, Hughes WJ, Arge CN, Owens MJ. 2008. Analysis of the magnetic field discontinuity at the potential field source surface and Schatten Current Sheet interface in the WangSheeleyArge model. J Geophys Res: Space Phys 113(A8): A08112. https://doi.org/10.1029/2007JA012330. [CrossRef] [Google Scholar]
 Merceret FJ, O’Brien TP, Roeder WP, Huddleston LL, Bauman WH III, Jedlovec GJ. 2013. Transitioning research to operations: Transforming the “Valley of Death” into a “Valley of Opportunity”. Space Weather 11(11): 637–640. https://doi.org/10.1002/swe.20099. [CrossRef] [Google Scholar]
 Merkin VG, Lionello R, Lyon JG, Linker J, Török T, Downs C. 2016. Coupling of coronal and heliospheric magnetohydrodynamic models: Solution comparisons and verification. Astrophys J 831(1): 23–36. https://doi.org/10.3847/0004637X/831/1/23. [CrossRef] [Google Scholar]
 Mulder WA, van Leer B. 1985. Experiments with implicit upwind methods for the Euler equations. J Comput Phys 59: 232–246. https://doi.org/10.1016/00219991(85)901445. [CrossRef] [Google Scholar]
 Nikolić L. 2017. Modelling the magnetic field of the solar corona with potentialfield sourcesurface and Schatten current sheet models. Geol Surv Canada Open File 8007. https://doi.org/10.4095/300826. [Google Scholar]
 Nikolić L. 2019. On solutions of the PFSS model with GONG synoptic maps for 2006–2018. Space Weather 17(8): 1293–1311. https://doi.org/10.1029/2019SW002205. [CrossRef] [Google Scholar]
 Northrup SA, Groth CPT. 2013. Parallel implicit adaptive mesh refinement scheme for unsteady fullycompressible reactive flows. In: 21st AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, AIAA 20132433. https://doi.org/10.2514/6.20132433. [Google Scholar]
 Odstrcil D. 2003. Modeling 3D solar wind structure. Adv Space Res 32(4): 497–506. https://doi.org/10.1016/S02731177(03)003326. [NASA ADS] [CrossRef] [Google Scholar]
 Odstrčil D, Linker JA, Lionello R, Mikić Z, Riley P, Pizzo VJ, Luhmann JG. 2002. Merging of coronal and heliospheric numerical twodimensional MHD models. J Geophys Res: Space Phys 107(A12): 1493–1–1493–11. https://doi.org/10.1029/2002JA009334. [Google Scholar]
 Odstrcil D, Pizzo VJ, Linker JA, Riley P, Lionello R, Mikic Z. 2004a. Initial coupling of coronal and heliospheric numerical magnetohydrodynamic codes. J Atmos SolTerr Phys 66(15): 1311–1320. https://doi.org/10.1016/j.jastp.2004.04.007. [NASA ADS] [CrossRef] [Google Scholar]
 Odstrcil D, Riley P, Zhao XP. 2004b. Numerical simulation of the 12 May 1997 interplanetary CME event. J Geophys Res: Space Phys 109(A2): A02116. https://doi.org/10.1029/2003JA010135. [NASA ADS] [CrossRef] [Google Scholar]
 Odstrcil D, Pizzo VJ, Arge CN, Bissi MM, Hick PP, et al. 2008. Numerical simulations of solar wind disturbances by coupled models. In: Numerical modeling of space plasma flows, Pogorelov NV, Audit E, Zank GP (Eds.), Vol. 385 of Astronomical Society of the Pacific Conference Series, 167 p. [Google Scholar]
 Owens MJ, Spence HE, McGregor S, Hughes WJ, Quinn JM, Arge CN, Riley P, Linker J, Odstrcil D. 2008. Metrics for solar wind prediction models: Comparison of empirical, hybrid, and physicsbased schemes with 8 years of L1 observations. Space Weather 6(8): S08001. https://doi.org/10.1029/2007SW000380. [CrossRef] [Google Scholar]
 Parker EN. 1958. Dynamics of the interplanetary gas and magnetic fields. Astrophys J 128(3): 664–676. https://doi.org/10.1086/146579. [NASA ADS] [CrossRef] [Google Scholar]
 Parsons A, Biesecker D, Odstrcil D, Millward G, Hill S, Pizzo V. 2011. WangSheeleyArgeEnlil cone model transitions to operations. Space Weather 9(3): S03004. https://doi.org/10.1029/2011SW000663. [CrossRef] [Google Scholar]
 Pizzo VJ. 1982. A threedimensional model of corotating streams in the solar wind: 3. Magnetohydrodynamic streams. J Geophys Res: Space Phys 87(A6): 4374–4394. https://doi.org/10.1029/JA087iA06p04374. [CrossRef] [Google Scholar]
 Pomoell J, Poedts S. 2018. EUHFORIA: European heliospheric forecasting information asset. J Space Weather Space Clim 8: A35. https://doi.org/10.1051/swsc/2018020. [CrossRef] [Google Scholar]
 Porth O, Xia C, Hendrix T, Moschou SP, Keppens R. 2014. MPIAMRVAC for solar and astrophysics. Astrophys J Suppl Ser 214(1): 4–29. https://doi.org/10.1088/00670049/214/1/4. [NASA ADS] [CrossRef] [Google Scholar]
 Porth O, Olivares H, Mizuno Y, Younsi Z, Rezzolla L, Moscibrodzka M, Falcke H, Kramer M. 2017. The black hole accretion code. Comput Astrophys Cosmol 4(1): 1–42. https://doi.org/10.1186/s4066801700202. [NASA ADS] [CrossRef] [Google Scholar]
 Powell KG. 1997. An approximate Riemann solver for magnetohydrodynamics (that works in more than one dimension). In: Upwind and highresolution schemes. Hussaini MY, van Leer B, Van Rosendale J (Eds.), Springer, Berlin, Heidelberg. https://doi.org/10.1007/9783642605437_23. [Google Scholar]
 Powell KG, Roe PL, Linde TJ, Gombosi TI, De Zeeuw DL. 1999. A solutionadaptive upwind scheme for ideal magnetohydrodynamics. J Comput Phys 154: 284–309. https://doi.org/10.1006/jcph.1999.6299. [NASA ADS] [CrossRef] [Google Scholar]
 Reiss MA, Temmer M, Veronig AM, Nikolic L, Vennerstrom S, Schöngassner F, Hofmeister SJ. 2016. Verification of highspeed solar wind stream forecasts using operational solar wind models. Space Weather 14(7): 495–510. https://doi.org/10.1002/2016SW001390. [CrossRef] [Google Scholar]
 Richardson IG. 2018. Solar wind stream interaction regions throughout the heliosphere. Living Rev Sol Phys 15: 1. https://doi.org/10.1007/s411160170011z. [CrossRef] [Google Scholar]
 Riley P, Linker JA, Mikic Z. 2001. An empiricallydriven global MHD model of the solar corona and inner heliosphere. J Geophys Res: Space Phys 106(A8): 15889–15901. https://doi.org/10.1029/2000JA000121. [NASA ADS] [CrossRef] [Google Scholar]
 Riley P, Linker JA, Mikić Z, Lionello R, Ledvina SA, Luhmann JG. 2006. A comparison between global solar magnetohydrodynamic and potential field source surface model results. Astrophys J 653(2): 1510–1516. https://doi.org/10.1086/508565. [NASA ADS] [CrossRef] [Google Scholar]
 Roussev II, Gombosi TI, Sokolov IV, Velli M, Manchester W, DeZeeuw DL, Liewer P, Toth G, Luhmann J. 2003. A threedimensional model of the solar wind incorporating solar magnetogram observations. Astrophys J Lett 595(1): L57. https://doi.org/10.1086/378878. [NASA ADS] [CrossRef] [Google Scholar]
 Saad Y, Schultz MH. 1986. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear equations. SIAM J Sci Stat Comput 7(3): 856–869. https://doi.org/10.1137/0907058. [CrossRef] [MathSciNet] [Google Scholar]
 Schatten KH. 1971. Current sheet magnetic model for the solar corona. Cos Electrodyn 2: 232. [Google Scholar]
 Schatten K, Wilcox J, Ness N. 1969. A model of interplanetary and coronal magnetic fields. Sol Phys 6(3): 442–455. https://doi.org/10.1007/BF00146478. [NASA ADS] [CrossRef] [Google Scholar]
 Sheeley NR. 2017. Origin of the WangSheeleyArge solar wind model. Hist Geo Space Sci 8(1): 21–28. https://doi.org/10.5194/hgss8212017. [Google Scholar]
 Shiota D, Kataoka R. 2016. Magnetohydrodynamic simulation of interplanetary propagation of multiple coronal mass ejections with internal magnetic flux rope (SUSANOOCME). Space Weather 14(2): 56–75. https://doi.org/10.1002/2015SW001308. [NASA ADS] [CrossRef] [Google Scholar]
 Shiota D, Kataoka R, Miyoshi Y, Hara T, Tao C, Masunaga K, Futaana Y, Terada N. 2014. Inner heliosphere MHD modeling system applicable to space weather forecasting for the other planets. Space Weather 12(4): 187–204. https://doi.org/10.1002/2013SW000989. [CrossRef] [Google Scholar]
 Steenburgh RA, Biesecker DA, Millward GH. 2014. From predicting solar activity to forecasting space weather: Practical examples of researchtooperations and operationstoresearch. Sol Phys 289(2): 675–690. https://doi.org/10.1007/s1120701303086. [CrossRef] [Google Scholar]
 Susanto A, Ivan L, Sterck HD, Groth CPT. 2013. Highorder central ENO finitevolume scheme for ideal MHD. J Comput Phys 250: 141–164. https://doi.org/10.1016/j.jcp.2013.04.040. [NASA ADS] [CrossRef] [Google Scholar]
 Taktakishvili A, Kuznetsova M, MacNeice P, Hesse M, Rastätter L, Pulkkinen A, Chulaki A, Odstrcil D. 2009. Validation of the coronal mass ejection predictions at the Earth orbit estimated by ENLIL heliosphere cone model. Space Weather 7(3): S03004. https://doi.org/10.1029/2008SW000448. [CrossRef] [Google Scholar]
 Taktakishvili A, Pulkkinen A, MacNeice P, Kuznetsova M, Hesse M, Odstrcil D. 2011. Modeling of coronal mass ejections that caused particularly large geomagnetic storms using ENLIL heliosphere cone model. Space Weather 9(6): S06002. https://doi.org/10.1029/2010SW000642. [CrossRef] [Google Scholar]
 Toro EF. 2013. Riemann solvers and numerical methods for fluid dynamics: A practical introduction, SpringerVerlag, New York. https://doi.org/10.1007/9783662034903. [Google Scholar]
 Tóth G, De Zeeuw DL, Gombosi TI, Manchester WB, Ridley AJ, Sokolov IV, Roussev II. 2007. Suntothermosphere simulation of the 28–30 October 2003 storm with the Space Weather Modeling Framework. Space Weather 5(6): S06003. https://doi.org/10.1029/2006SW000272. [Google Scholar]
 Tóth G, van der Holst B, Huang Z. 2011. Obtaining potential field solutions with spherical harmonics and finite differences. Astrophys J 732(2): 102. https://doi.org/10.1088/0004637X/732/2/102. [NASA ADS] [CrossRef] [Google Scholar]
 Tóth G, van der Holst B, Sokolov IV, De Zeeuw DL, Gombosi TI, et al. 2012. Adaptive numerical algorithms in space weather modeling. J Comput Phys 231(3): 870–903. https://doi.org/10.1016/j.jcp.2011.02.006. [NASA ADS] [CrossRef] [Google Scholar]
 Usmanov A. 1993. A global numerical 3D MHD model of the solar wind. Sol Phys 146(2): 377–396. https://doi.org/10.1007/BF00662021. [CrossRef] [Google Scholar]
 van der Holst B, Sokolov IV, Meng X, Jin M, Manchester WB, Tóth G, Gombosi TI. 2014. Alfvén wave solar model (AWSoM): Coronal heating. Astrophys J 782(2): 81–95. https://doi.org/10.1088/0004637X/782/2/81. [NASA ADS] [CrossRef] [Google Scholar]
 Verbeke C, Pomoell J, Poedts S. 2019. The evolution of coronal mass ejections in the inner heliosphere: Implementing the spheromak model with EUHFORIA. A&A 627: A111. https://doi.org/10.1051/00046361/201834702. [CrossRef] [EDP Sciences] [Google Scholar]
 Vourlidas A, Subramanian P, Dere KP, Howard RA. 2000. Largeangle spectrometric coronagraph measurements of the energetics of coronal mass ejections. Astrophys J 534(1): 456–467. https://doi.org/10.1086%2F308747. [NASA ADS] [CrossRef] [Google Scholar]
 Wang YM, Sheeley NR. 1990. Solar wind speed and coronal fluxtube expansion. Astrophys J 355: 726–732. https://doi.org/10.1086/168805. [NASA ADS] [CrossRef] [Google Scholar]
 Wang Y, Sheeley NR, Phillips JL. 1995. Coronal fluxtube expansion and the polar wind. Adv Space Res 16(9): 365–365. https://doi.org/10.1016/02731177(95)00371K. [CrossRef] [Google Scholar]
 Wang YM, Sheeley NR, Phillips JL, Goldstein BE. 1997. Solar wind stream interactions and the wind speedexpansion factor relationship. Astrophys J 488(1): L51–L54. https://doi.org/10.1086/310918. [CrossRef] [Google Scholar]
 Webb DF, Howard TA. 2012. Coronal mass ejections: Observations. Living Rev Sol Phys 9: 1. https://doi.org/10.12942/lrsp20123. [Google Scholar]
 Williamschen MJ, Groth CPT. 2013. Parallel anisotropic blockbased adaptive mesh refinement algorithm for threedimensional flows. In: 21st AIAA Computational Fluid Dynamics Conference, American Institute of Aeronautics and Astronautics, AIAA 20132442. https://doi.org/10.2514/6.20132442. [Google Scholar]
 Wold AM, Mays ML, Taktakishvili A, Jian LK, Odstrcil D, MacNeice P. 2018. Verification of realtime WSAENLIL+Cone simulations of CME arrivaltime at the CCMC from 2010 to 2016. J Space Weather Space Clim 8: A17. https://doi.org/10.1051/swsc/2018005. [CrossRef] [Google Scholar]
 Xia C, Keppens R. 2016. Formation and plasma circulation of solar prominences. Astrophys J 823(1): 22–40. https://doi.org/10.3847/0004637X/823/1/22. [NASA ADS] [CrossRef] [Google Scholar]
 Xia C, Chen PF, Keppens R. 2012. Simulations of prominence formation in the magnetized solar corona by chromospheric heating. Astrophys J Lett 748(2): L26. https://doi.org/10.1088/20418205/748/2/L26. [NASA ADS] [CrossRef] [Google Scholar]
 Xia C, Keppens R, Guo Y. 2013. Threedimensional prominencehosting magnetic configurations: Creating a helical magnetic flux rope. Astrophys J 780(2): 130–154. https://doi.org/10.1088/0004637X/780/2/130. [NASA ADS] [CrossRef] [Google Scholar]
 Xia C, Keppens R, Antolin P, Porth O. 2014. Simulating the in situ condensation process of solar prominences. Astrophys J Lett 792(2): L38. https://doi.org/10.1088/20418205/792/2/L38. [NASA ADS] [CrossRef] [Google Scholar]
 Xia C, Teunissen J, El Mellah I, Chané E, Keppens R. 2018. MPIAMRVAC 2.0 for solar and astrophysical applications. Astrophys J Suppl Ser 234(2): 30–55. https://doi.org/10.3847/15384365/aaa6c8. [NASA ADS] [CrossRef] [Google Scholar]
 Xie H, Ofman L, Lawrence G. 2004. Cone model for halo CMEs: Application to space weather forecasting. J Geophys Res: Space Phys 109(A3): A03109. https://doi.org/10.1029/2003JA010226. [Google Scholar]
 Zhao X, Hoeksema J. 1994. A coronal magnetic field with horizontal volume and sheet currents. Sol Phys 151(1): 91–105. https://doi.org/10.1007/BF00654084. [NASA ADS] [CrossRef] [Google Scholar]
 Zhao XP, Plunkett SP, Liu W. 2002. Determination of geometrical and kinematical properties of halo coronal mass ejections using the cone model. J Geophys Res: Space Phys 107(A8): SSH 13–1–SSH 13–9. https://doi.org/10.1029/2001JA009143. [NASA ADS] [CrossRef] [Google Scholar]
 Zhou Y, Xia C, Keppens R, Fang C, Chen PF. 2018. Threedimensional MHD simulations of solar prominence oscillations in a magnetic flux rope. Astrophys J 856(2): 179–192. https://doi.org/10.3847/15384357/aab614. [NASA ADS] [CrossRef] [Google Scholar]
Cite this article as: Narechania NM, Nikolić L, Freret L, De Sterck H & Groth CPT. 2021. An integrated datadriven solar wind – CME numerical framework for space weather forecasting. J. Space Weather Space Clim. 11, 8. https://doi.org/10.1051/swsc/2020068.
All Figures
Fig. 1 Schematic diagram showing the subdomains and mathematical models used in each region in the proposed SuntoEarth datadriven 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 R_{0}, the SCS model describes the coronal field up to 20–30 R_{0}, 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 gridblock structure of computational mesh after several levels of refinement. 

In the text 
Fig. 4 Schematic diagram of heterogeneous gridblock 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 cubedsphere 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 R_{0} for the steady solar wind simulation obtained using the standard GONG map for 12 February 2007 (11:54 UT). The central meridian and subEarth locations are denoted with dashdot and dashdash 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 Br^{2} at the outer boundary r = 225 R_{0} (≈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 Br^{2} in the xzplane (top) and yzplane (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 xyplane (ecliptic plane) showing the predicted contours of Br^{2} 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 greencoloured sphere, E, corresponds to the approximate location of Earth. 

In the text 
Fig. 12 (a) Predicted contours of Br^{2} in the ecliptic xyplane 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 R_{0} 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 subEarth locations are denoted with dashdot and dashdash line, respectively. 

In the text 
Fig. 14 Predicted contours of Br^{2} in the xy, yz and xzplanes, 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 zeropoint 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 (SMap) was used and, in the other case, the zeropoint corrected GONG map (ZPCMap) 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 (SMap MD and ZPCMap MD). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.