Issue 
J. Space Weather Space Clim.
Volume 11, 2021
Topical Issue  10 years of JSWSC



Article Number  42  
Number of page(s)  55  
Section  Agora  
DOI  https://doi.org/10.1051/swsc/2021020  
Published online  13 August 2021 
Agora – Project report
What sustained multidisciplinary research can achieve: The space weather modeling framework
^{1}
Department of Climate and Space Sciences and Engineering (CLaSP), University of Michigan, Ann Arbor, MI 48109, USA
^{2}
Geospace Physics Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD 20771, USA
^{3}
Department of Radio Science and Engineering, Aalto University, Espoo, Finland
^{4}
Department of Physics, University of Texas at Arlington, TX 76019, USA
^{*} Corresponding author: tamas@umich.edu
Received:
22
March
2021
Accepted:
20
May
2021
Magnetohydrodynamics (MHD)based global space weather models have mostly been developed and maintained at academic institutions. While the “free spirit” approach of academia enables the rapid emergence and testing of new ideas and methods, the lack of longterm stability and support makes this arrangement very challenging. This paper describes a successful example of a universitybased group, the Center of Space Environment Modeling (CSEM) at the University of Michigan, that developed and maintained the Space Weather Modeling Framework (SWMF) and its core element, the BATSRUS extended MHD code. It took a quarter of a century to develop this capability and reach its present level of maturity that makes it suitable for research use by the space physics community through the Community Coordinated Modeling Center (CCMC) as well as operational use by the NOAA Space Weather Prediction Center (SWPC).
Key words: Space weather / Solar flares and CMEs / Scientific computing / Space plasma physics / MHD
© T.I. Gombosi 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.
Special Issue  10 years of JSWSC
1 Introduction
Over the past few decades there has been an increasing awareness of the potentially devastating impact that the dynamic space environment can have on human assets. Extreme “space weather” events, driven by eruptive solar events such as Coronal Mass Ejections (CMEs), are widely recognized as critical hazards whose consequences cannot be ignored.
Conservation laws in BATSRUS.
Because of society’s reliance on the electrical grid, the internet, highfrequency communication, GPS (Global Positioning System) navigation signals and an increasing array of digital electronic devices, space weather events – such as severe solar storms – can wreak havoc on technological systems and trigger losses from business interruption and damaged physical assets (cf., Baker et al., 2009). While power outages from space weather are lowfrequency events, they have the potential to cause crippling longterm damage. In fact, the risk of high impact damages due to space weather fits the profile of a marketchanging catastrophe such as hurricane Katrina, the 9/11 attack, or the Japanese earthquake and tsunami (cf. FEMA, 2019). All were unprecedented and believed to be highly unlikely – and yet they occurred.
There is an additional, less publicized reason that policymakers care about space weather: its association to electromagnetic pulses (EMPs) (cf., Gombosi et al., 2017). An EMP is a natural or anthropogenic burst of electromagnetic energy that can damage all kinds of electronic and even physical objects. Understanding and mitigating space weather effects also have national defense implications.
Space weather involves a vast domain extending from the Sun to beyond Earth’s orbit, with regions governed by very different physics at different spatial and temporal scales. Simulating and predicting space weather with firstprinciples models requires space physics expertise for the various subdomains and advanced numerical algorithms. Since the subdomain models keep changing and evolving, they need to be coupled in a flexible manner using proper software engineering. Finally, the simulation needs to run faster than real time, which means that a deep understanding of highperformance computing is required. Clearly, developing a firstprinciples space weather model requires sustained multidisciplinary collaboration of space physicists, applied mathematicians, computer scientists and software engineers.
Presently there are only a couple of physicsbased space weather models that are capable of spanning the entire region from the low solar corona to the edge of the heliosphere. One is the European Space Agency’s Virtual Space Weather Modelling Centre (VSWMC, Poedts et al., 2020) and the other one is the Space Weather Modeling Framework (SWMF, Toth et al., 2005, 2012). In this paper we describe the evolution and current capabilities of the SWMF and its unique capabilities to address the myriad of processes involved in studying and predicting space weather. In the main text we focus on the the broad range of space weather simulations made possible by the advanced capabilities of BATSRUS (Block AdaptiveTree Solarwind Roetype Upwind Scheme) and SWMF. The fundamentals of the BATSRUS and SWMF codes are described in detail in Appendix A. The extended physics and algorithmic advances incorporated in these codes are important and we present a concise summary of these advances in Appendices C (physics) and D (algorithms). Finally, Appendix E describes our most advanced simulation capability that embeds fully kinetic domains inside extended MHD models.
2 Evolution of space weather models
Models capable of predicting space weather can be loosely divided into three broad categories: Empirical models, black box (mainly machine learned) models, and physicsbased models.
2.1 Empirical models
Empirical models aggregate data in different ways to make specific predictions of the current and future state of the system based on how the system has responded historically. Such models are mostly data driven and typically make limited or no assumptions of the underlying physics. The quality of the models is heavily dependent on the data coverage in space and over different geomagnetic conditions. Widely used examples are the MSIS (Mass Spectrometer and Incoherent Scatter) model (Hedin, 1987, 1991) of the upper atmosphere and the Tsyganenko (1989, 1995, 2002a, b) models of the terrestrial magnetic environment.
The MSIS model (Hedin, 1987, 1991) brings together mass spectrometer and incoherent scatter data to build an empirical model of the thermosphere. The model provides estimates of temperature and the densities of atmospheric constituents such as N_{2}, O, O_{2}, He, Ar, and H. Loworder spherical harmonics expansion is used to describe spatial (latitude, local time), and temporal (annual, semiannual) variations. The model is often used for data comparisons and theoretical calculations requiring a background atmosphere, for example in calculations of satellite orbital decay caused by atmospheric drag.
The extension of the geodipole field to the magnetosphere is sustained by currents flowing in the geospace. The magnetic field variations from these currents can be deduced from spaceborne magnetic field measurements, and have been collected into a large database. The Tsyganenko models (Tsyganenko, 1989, 1995, 2002a, b; Tsyganenko & Stern, 1996; Tsyganenko & Sitnov, 2005), and describe the largescale current systems with parametrized empirical functions, and the parameter values are found through leastsquares fitting to the large observational database. The models have been extensively used e.g., to connect magnetospheric substorm and storm dynamic processes to their ionospheric signatures (Pulkkinen et al., 1992, 2006; Baker et al., 1996).
2.2 Blackbox models
Linear prediction filters have been used to build models for a variety of space weather parameters, including the auroral electrojet (AE) indices and the ring current Dst (Disturbance stormtime) index. Predictions of magnetospheric storm conditions have been done using neural networks to construct nonlinear models to forecast the AL and/or Dst index using various solar wind driver parameters (Lundstedt & Wintoft, 1994; Weigel et al., 2003).
Recent machine learning models have been quite successful in predicting geomagnetic indices (see Leka & Barnes, 2018; Camporeale, 2019). To support their use in space weather research requires openaccess, robust, and effective software tools. Typically, the models are custommade and making use of a stack of standard computational frameworks for learning. Machine learning methods have also been employed for prediction of the ionospheric total electron content (TEC) (cf. Liu et al., 2000) and solar flares (cf. Chen et al., 2019b; Jiao et al., 2020; Wang et al., 2020). However, as most machine learning models are not interpretable, they typically do not help us to understand the underlying physics.
2.3 Physicsbased models
Physics based models directly solve equations representing the underlying physical processes in the system, often with observations based inputs, in order to study the evolution and dynamics of the space environment. Physicsbased space weather models have been found to be particularly valuable for predicting both the rare extreme events as well as more commonly observed space weather.
Extreme space weather events with the most severe implications for human assets and activities are lowfrequency events that create challenges for forecasting and prediction. Since the dawn of the space age, there have been a handful of events with major space weather impacts, as well as other events with more modest effects. For example, the 13 March 1989 event was a particularly strong case with a minimium Dst of −589 nT that induced currents in the power grid leading to the ultimate collapse of the HydroQuebec power system (Bolduc, 2002). There is a great deal of interest in both being able to predict such events in advance, as well as quantifying how strong events could result in widespread disruptions. The low frequency of such events is particularly challenging for empirical or machine learning models, which struggle with out of sample predictions.
Global magnetohydrodynamics (MHD) models for space science applications were first published in the early 1980s (Brecht et al., 1981, 1982; LeBoeuf et al., 1981; Wu et al., 1981). Later models applied more advanced algorithms to solve the MHD equations. These models include the Lyon–Fedder–Mobarry (LFM) (Lyon et al., 1986, 2004), the OpenGGCM (Open Geospace General Circulation Model, Raeder et al., 1996, 1995), the Watanabe–Sato (Watanabe & Sato, 1990; Usadi et al., 1993), the GUMICS (Grand Unified Magnetosphere Ionosphere Coupling Simulation model Janhunen, 1996), and the Integrated Space Weather Prediction Model (ISM) (White et al., 1998; Siscoe et al., 2000), models of the Earth’s magnetosphere. The solar codes include models for the solar corona (Magnetohydrodynamics Around a Sphere (MAS), Linker et al., 1994; Linker et al., 1999), (Hayashi, 2013), the heliosphere (Usmanov, 1993; Usmanov et al., 2000), the inner heliosphere ENLIL (Odstrčil, 2003; Odstrčil & Pizzo, 2009), as well as combined models of the corona and inner heliosphere (Solarinterplanetary adaptive mesh refinement spacetime conservation element and solution element MHD model (SIPAMRCESE MHD Model), Feng et al., 2014a, b). More generaluse models include Ogino’s planetary magnetosphere code (Ogino, 1986), Tanaka’s 3D global MHD model and, Winglee’s multifluid Hall MHD code (Washimi & Tanaka, 1996), Winglee’s multifluid HallMHD code (Winglee, 1998; Winglee et al., 2005), Toth’s general MHD Versatile Advection Code (VAC) (Toth, 1996) and its modern version, MPIAMRVAC (Keppens et al., 2021), KU Leuven’s European heliospheric forecasting information asset (EUHFORIA, Pomoell & Poedts, 2018) and the University of Michigan’s BATSRUS (Powell et al., 1999; Toth et al., 2012), model.
3 The origins of BATSRUS and SWMF
Advanced space plasma simulation codes became possible when leading applied mathematicians and computer scientists became integral parts of the teams developing models to solve physical systems. In the early 1990s, two pioneers of highorder Godunov (1959) schemes that revolutionized computational fluid dynamics (CFD), Bram van Leer (cf. van Leer, 1973, 1974, 1977a, 1977b, 1979) and Philip Roe (cf. Roe, 1981), became interested in space physics problems.This interest resulted in the extension of modern CFD methods to rarefied magnetized plasma flows and the development of the first modern, high performance MHD code, BATSRUS (Powell et al., 1999). Figure 1 summarizes the present capabilities of BATSRUS; the algorithms are discussed in detail in Appendix D.
Fig. 1 Overview of the BATSRUS multiphysics code. 
The BATSRUS (Powell et al., 1999; Toth et al., 2012) is a versatile, highperformance, generalized magnetohydrodynamic code with adaptive mesh refinement (AMR) that can be configured to solve the governing equations of ideal and resistive MHD (Powell et al., 1999), semirelativistic (Gombosi et al., 2002), anisotropic (Meng et al., 2012), Hall (Toth et al., 2008), multispecies (Ma et al., 2002), and multifluid (Glocer et al., 2009c), extended magnetofluid equations (XMHD) and, most recently, nonneutral multifluid plasmas (Huang et al., 2019). BATSRUS is used to model several physics domains (see Fig. 2). The efficiency of BATSRUS is crucial to reach faster than realtime performance with the SWMF while maintaining high resolution in the domains of interest.
Fig. 2 Schematic diagram of the Space Weather Modeling Framework. The SWMF and its core models are open source (https://github.com/MSTEMQUDA), while the full SWMF is available via registration under a user license (http://csem.engin.umich.edu/tools/smmf). 
In a number of fields in which computerbased modeling of complex, multiscale, multiphysics problems plays an important role, software frameworks have been developed. In the area of computational space physics there are only two operational software frameworks, the Space Weather Modeling Framework (SWMF, Toth et al., 2005, 2012) and the Virtual Space Weather Modelling Centre (VSWMC, Poedts et al., 2020). Other frameworks are either under development (Zhang et al., 2019a), abandoned (Luhmann et al., 2004), or are rarely used for space weather applications (Hill et al., 2004).
The SWMF (Toth et al., 2005, 2012) is a fully functional, documented software that provides a highperformance computational capability to simulate the spaceweather environment from the upper solar chromosphere to the Earth’s upper atmosphere and/or the outer heliosphere. The SWMF tackles the wide range of temporal and spatial scales as well as the different physical processes governing the different heliophysics domains through a modular approach. Each physics domain is covered by a numerical model developed particularly for that purpose. The framework couples several of these components together to execute the simulation in a setup best suited for the problem at hand.
4 The SWMF today
In 2021, the Space Weather Modeling Framework (SWMF) (Toth et al., 2012), consists of a dozen physics domains and a dozen different models that provide a flexible highperformance computational capability to simulate the spaceweather environment from the upper solar chromosphere to the Earth’s upper atmosphere and/or the outer heliosphere. It contains over 1 million lines of Fortran 2008 and C++ code, dozens of Perl, Python and Julia scripts, IDL visualization tools and XML descriptions of the input parameters. Figure 2 summarizes the main features and capabilities of the current SWMF.
The full SWMF suite, developed and maintained at the University of Michigan, has been openly available for a long time via registration under a user license (http://csem.engin.umich.edu/tools/swmf). Recently, a major part of the SWMF has been released on Github under a noncommercial opensource license (https://github.com/MSTEMQUDA). Figure 2 shows the open source and registration controlled components of the SWMF.
In addition, SWMF runs can be requested via the Community Coordinated Modeling Center (CCMC) at the NASA Goddard Space Flight Center (https://ccmc.gsfc.nasa.gov/index.php), where people even with little experience in advanced computer simulations can request specific runs through a userfriendly web interface. The user specifies the domains and the driving input parameters, and the CCMC runsonrequest system carries out the simulation. Once the CCMC completes the run, the output files and standard visualization images are made available through the web interface (https://ccmc.gsfc.nasa.gov/index.php).
For space weather related simulations, the SWMF is typically used in two basic configurations: The Alfvén Wave Solaratmosphere Model (AWSoM/AWSoMR) and the SWMF/Geospace Model.
AWSoM/AWSoMR (van der Holst et al., 2010, 2014; Sokolov et al., 2013, 2021; Gombosi et al., 2018), describes the solar corona (SC) from the low transition region where the plasma temperature is about 5 × 10^{4} K and goes out to about 20 R_{⊙}. This is the region where the hot, supersonic solar wind is generated. It also simulates the 3D inner heliosphere (IH) out to Neptune’s orbit. The outer boundary can be varied depending on the region of interest.
The SWMF/Geospace Model (cf. Haiducek et al., 2017; Welling et al., 2020) describes the tightly coupled basic elements of the magnetosphereionosphere system: the global magnetosphere (GM), the inner magnetosphere (IM), the ionospheric electrodynamics (IE). An operational version of the SWMF/Geospace model has been running 24/7 at SWPC since 2016.
4.1 AWSoM/AWSoMR configuration
It is commonly accepted that the gradient of the Alfvén wave pressure is the key driver for solar wind acceleration. Damping of Alfvén wave turbulence due to reflection from sharp pressure gradients in the solar wind is a critical component of coronal heating. For this reason, many numerical models explore the generation of reflected counterpropagating waves as the underlying cause of the turbulence energy cascade (e.g., Cranmer & Van Ballegooijen, 2010), which transports the energy of turbulence from the largescale motions across the inertial range of the turbulence spatial scale to shortwavelength perturbations. The latter can be efficiently damped due to waveparticle interaction. In this way, the turbulence energy is converted to random (thermal) energy (cf. Sokolov et al., 2013).
4.1.1 AWSoM
AWSoM (Sokolov et al., 2013, 2021; Gombosi et al., 2018; van der Holst et al., 2014), is a 3D global solar corona/solar wind model that selfconsistently incorporates lowfrequency Alfvén wave turbulence. The Alfvén waves are represented by the energy density distribution of two discrete populations propagating parallel and antiparallel to the magnetic field at the local Alfvén speed. The wave energy densities are imposed at the inner boundary with a Poynting flux of the outbound Alfvén waves assumed to be proportional to the magnetic field strength. In this model, outward propagating waves experience partial reflection on fieldaligned Alfvén speed gradients and the vorticity of the background. In addition, the two populations counterstream along closed field lines. The nonlinear interaction between oppositely propagating Alfvén waves results in an energy cascade from the large outer scale through the inertial range to the smaller perpendicular gyroradius scales, where the dissipation takes place. These processes are handled with analytic formulas that provide the resulting ion and electron heating. The solar wind is accelerated by the gradient of the Alfvén wave pressure. The main physics elements of the AWSoM model are shown in Figure 3.
Fig. 3 Overview of the AWSoM and AWSoMR physics. They solve XMHD equations with separate ion and electron temperatures. The energy densities of parallel and antiparallel propagating turbulence that are selfconsistently coupled to each other and to the plasma are solved together with the XMHD equations. Heat conduction and radiative cooling are also taken into account. The turbulence is powered by the Poynting flux leaving the solar photosphere. 
The boundary conditions for the MHD quantities are obtained from the synoptic or synchronic photospheric magnetograms. The outward propagating Poynting flux at the solar surface (S_{A}) is measured in units of “W/m^{2}” and it is taken to be proportional to the magnetic field magnitude, B (measured in units of Tesla or T). The proportionality constant α is measured in units of “MW/m^{2}/Tesla”, and its value varies between 0 and 1. The actual value of α depends on the phase of the solar cycle and on the choice of magnetogram (to account for the calibration differences between magnetograms).
The inner heliosphere (IH) component extends from about 20 R_{⊙} to anywhere between the orbits of the Earth and Neptune. It uses the BATSRUS and it solves the same equations as the solar corona model, but on a Cartesian grid in either corotating or inertial frame. The IH model can propagate interplanetary CMEs (ICMEs) from the Sun to the planets. Adaptive mesh refinement is used to increase the grid resolution along the path of the CME (cf. Roussev et al., 2004; van der Holst et al., 2009; Manchester et al., 2014a; Manchester & van der Holst, 2017).
4.1.2 ThreadedFieldLine Model and AWSoMR
In the transition region the plasma temperature increases some two orders of magnitude over ~10^{2} km, resulting in a temperature gradient of ~10^{4} K/km. To resolve this gradient, 3D numerical simulations require subkilometer grid spacing, making these simulations computationally very expensive. AWSoM uses an artificial broadening of the transition region (Lionello et al., 2009; Sokolov et al., 2013).
An alternative approach is to reformulate the mathematical problem in the region between the chromosphere and the corona in a way that decreases the computational cost. Instead of solving a computationally expensive 3D problem on a very fine grid, one can reformulate it in terms of a multitude of much simpler 1D problems along threads that allows us to map the boundary conditions from the the solar photosphere to the corona. This approach is called the ThreadedFieldLine Model (TFLM) (Gombosi et al., 2018; Sokolov et al., 2021).
The physics behind the reformulated problem is the assumption that between the solar surface and the top of the transition region (R_{⊙} ≤ r ≤ R_{b}) the magnetic field is potential and varies slowly in time. Each thread represents a field line and one can solve a 1D problem that describes evolution of the plasma in a magnetic flux tube around a given thread. The algorithm uses an implicit scheme to allow for large time steps. Using the TFLM methodology results in a significant speedup for timedependent simulations. The AWSoM model with TFLM inner boundary conditions is called AWSoMR, where the letter “R” implies that this version can run faster than real time on ~200 cores at a moderate grid resolution (about 2° near the Sun).
4.2 SWMF/Geospace configuration
While the BATSRUS can model many of the dynamical plasma processes in the solar wind and magnetosphere, it is widely accepted that MHD alone cannot sufficiently describe the coupled solar wind – magnetosphere – ionosphere system. The ionosphere and space close to the Earth is not suited for MHD, and is beyond the numerical capabilities due to the high magnetic field intensity, which increases the wave speeds, thus requiring very small time steps and high spatial resolution. Furthermore, the inner magnetosphere ring current, which is an integral part of the storm dynamics, cannot be described by a temperature of a Maxwellian plasma population, which calls for separate treatment of the dynamics in the quasidipolar region. To that end, the SWMF/Geospace couples three different models describing these three domains. Furthermore, additional models can be coupled to tackle multiple plasma populations, kinetic physics, or other phenomena and processes (see Sect. 5).
The base SWMF/Geospace configuration is illustrated in Figure 4. Under this setup, the global magnetosphere model BATSRUS is coupled to the Ridley ionosphere electrodynamics model (RIM) (Ridley et al., 2004), and the inner magnetosphere Rice Convection Model (RCM) (Harel et al., 1981). BATSRUS supplies nearbody fieldaligned currents (FACs) to the RIM, which, using an empirical specification of conductance (Ridley et al., 2004; Mukhopadhyay et al., 2020), solves for the electric potential. This electric potential is returned to BATSRUS to set the plasma tangential velocity at the inner boundary. The RCM receives its initial and boundary field and plasma conditions from BATSRUS as well as electric field from RIM. It returns total plasma pressure and density to BATSRUS inside the closed field line region, significantly improving the inner magnetosphere results of the MHD solution (De Zeeuw et al., 2004), especially during geomagnetic storm times (Liemohn et al., 2018). In addition, the current configuration can include the Radiation Belt Environment (RBE) model (Fok et al., 2008), that receives information from BATSRUS and RIM and solves for the energetic electron population in the radiation belts.
Fig. 4 Illustration of the models (components within SWMF) and couplings in the SWMF/Geospace configuration. Arrows denote the information that is passed between the components (adapted from Haiducek et al., 2017). 
The couplings default to 5second (GMIE) and 10second (all other) frequency; faster coupling frequencies are required under extreme driving or when highfrequency output is produced (Welling et al., 2020). While the explicit couplings are shown, the selfconsistent nature of multimodel SWMF simulations produces implicit couplings. For example, while region2 Birkeland currents are not explicitly passed from IM to GM physics modules, the improved pressure gradients in BATSRUS due to the pressure coupling from RCM drives region2 Birkeland currents (Welling et al., 2018). Under this model configuration, only upstream solar wind and IMF conditions, as well as F10.7 solar radio flux, are needed as inputs to the model.
The Geospace model is initialized by iterating GM and IE toward an approximate steady state solution using the initial solar wind, IMF and F10.7 values for boundary conditions. Using a local time stepping mode, this is done very efficiently. Next the IM component is switched on and the Geospace model is run in timedependent mode using the time varying boundary conditions. It takes about 5 h for the ring current to build up to a realistic strength. After this point the model can be used for simulation and prediction. In operational use, the Geospace model is run continuously. The model is only reinitialized from scratch if there is a long (an hour or more) gap in the solar wind observations.
In addition to the physics models and couplings, spatial resolution of the included models strongly affects the simulation results. RIM defaults to 2° × 2° grid spacing in geomagnetic longitude and latitude. BATSRUS has no default grid, but the base SWMF/Geospace configurations are illustrated in Figure 5 for version 1 and the more recent version 2. These configurations result in ~1 million grid cells with a nearbody resolution of 1/4 R_{E} and ~2 million grid cells with 1/8 R_{E} maximum resolution, respectively.
Fig. 5 Grid configurations for BATSRUS within the SWMF Geospace. The left and right hand panels illustrate the grid configuration of the operational Geospace model versions 1 and 2, respectively (from Haiducek et al., 2017). 
While capable of running faster than real time on a modest number (about 100) of CPU cores, the operational SWMF/Geospace models can well reproduce largescale features such as Cross Polar Cap Potential (CPCP) and Dst (Haiducek et al., 2017; Mukhopadhyay et al., 2021), and can also predict local ground magnetic perturbations with skill scores of practical value (Pulkkinen et al., 2013; Toth et al., 2014).
4.2.1 Virtual Magnetic Observatories
The coupledmodel approach of SWMF/Geospace allows for the production of virtual observatory simulations during code execution. The most widely used of these are virtual magnetometers. We use BiotSavart integrals to find the total surface magnetic perturbation at an arbitrary point about the globe due to the simulated magnetospheric and ionospheric current systems (Yu & Ridley, 2008; Yu et al., 2010; Welling, 2019). For a detailed description of the methodology see Appendix C.5. While tools exist to create such outputs as part of postprocessing (Rastätter et al., 2014), the SWMF/Geospace combines information from the IE and GM models onthefly to provide continuous output during the simulation. A recently developed mathematical reformulation of the problem replacing the volume integrals with surface integrals speeds up the calculation by an order of magnitude (see Appendix C.5).
In a similar fashion, advanced virtual satellite observations are created by mapping kinetic distributions from the IM and optional RB modules along selfconsistent global magnetic field lines obtained from GM. The net result is the ability to extract ring current and radiation belt flux distributions at arbitrary points about the inner magnetosphere. Virtual satellites have also been used to assess the simulation results through comparisons with insitu spacecraft observations (cf. Welling & Zaharia, 2012; Glocer et al., 2013).
4.2.2 Operational use at NOAA/SWPC and the CCMC
In 2015, NOAA’s Space Weather Prediction Center (NOAA/SWPC) decided to transition a research model to operational space weather prediction. As part of this effort, a systematic study was undertaken to evaluate the performance of various physicsbased and empirical models to predict ground magnetic perturbations (Pulkkinen et al., 2013; Glocer et al., 2016). The physicsbased SWMF/Geospace model in particular was found to systematically be a top performing model using the selected metrics. That code has since been used for routine space weather prediction at NOAA/SWPC and at the Community Coordinated Modeling Center (CCMC) located at NASA GSFC. The operational codes run in the configuration illustrated in Figure 4. In 2020, the NOAA/SWPC upgraded to version 2 of the SWMF/Geospace model, which has a higher grid resolution near the Earth and better ionospheric conductance.
5 Growing number of space weather applications
Space weather simulations using the SWMF have been carried out in multiple configurations and contexts, demonstrating that SWMF and its components are able to successfully simulate globalscale, mesoscale and microscale processes in a selfconsistent manner, and integrate these processes to form a truly multiscale space weather simulation capability. In addition, significant validation efforts have been made by a variety of comparisons with both insitu and remotesensing observations.
5.1 Ambient solar wind
CMEs and ICMEs do not propagate and evolve in vacuum. They travel through the ambient interplanetary medium and interact with its plasma and magnetic field. Therefore, in order to simulate real space weather events, it is critical to have a validated ambient corona/solar wind model in which the CME/ICME will propagate and cause significant distortions. These distortions can include plasma pileup, shock fronts, magnetic field line distortion and many other phenomena (cf. Manchester et al., 2004b, 2005, 2008, 2012, 2014a). The situation can be even more complicated when several CMEs are generated in rapid succession (cf. Lugaz et al., 2005b, 2008, 2009).
Sachdeva et al. (2019) performed a detailed validation study of the AWSoM for the quiettime solar wind for Carrington Rotations (CR) representative of the solar minimum conditions (CR2208 and CR2209). They compared simulation results with a comprehensive suite of observations extending from the solar corona to the heliosphere up to Earth’s orbit. In the low corona (r < 1.25 R_{⊙}), extreme ultraviolet (EUV) images from both the STEREOA (Solar TErrestrial RElations Observatory Ahead) EUVI (extreme ultraviolet imaging) instrument and the SDO (Solar Dynamics Observatory) AIA (atmospheric imaging assembly) were compared to 3D tomographic reconstructions of the simulated electron temperature and density. Model results were also compared to tomographic reconstructions of the electron density from the SOHO (Solar and Heliospheric Observatory) LASCO (Large Angle and Spectrometric Coronagraph) observations in the 2.55 R_{⊙} < r < 6.0 R_{⊙} region. In the heliosphere, model predictions of solar wind speed were compared to velocity reconstructions from interplanetary scintillation observations. Simulation results at the first Lagrange point between the Sun and Earth (L1) were compared to OMNI data. The results of Sachdeva et al. (2019) show that the AWSoM performs well in quantitative agreement with the observations between the inner corona and 1 AU.
Recently AWSoM/AWSoMR was also validated for solar maximum conditions. Using S_{⊙}/B_{⊙} (S_{⊙} is the Poynting flux of outward propagating Alfvén waves at the solar surface) as an adjustable parameter, good agreement was found for CR2123 that characterizes solar maximum conditions for solar cycle 24 (see Fig. 6). Figure 6 shows the comparisons of AWSoMR simulation results for CR2123 and CR2209 with AIA images and solar wind parameters at 1 AU. For both the rotations the AIA comparisons include six wavelengths (94, 171, 193, 131, 211 and 335 Å). The L1 parameters include radial speed (U_{r} in km/s), proton number density (N_{p} in cm^{−3}) and temperature (T in K) and magnetic field magnitude (B in nT). Inspection of Figure 6 reveals that we are able to match the observed slow/fast solar wind structure at 1 AU and, simultaneously, reproduce a number of optically thin coronal spectral observations. For AWSoM model results of CR2209 the reader is referred to Sachdeva et al. (2019).
Fig. 6 Background corona and solar wind solutions with the AWSoMR model for solar minimum and maximum conditions. The background solar wind is driven by an outward going photospheric turbulent energy flux per unit magnetic field of 1 MJ m^{−2} s^{−1} Tesla^{−1} (CR2209) and by 0.45 MJ m^{−2} s^{−1} Tesla^{−1} (CR2123). 
5.2 CME generation
The Eruptive Event (EE) generator algorithm of the SWMF is responsible for creating the initial conditions within the corona, which produces a CME eruption. This be done by inserting an unstable (or force imbalanced) flux rope into the steady solar corona solution, or inserting an arcade and applying shearing motion at the lower boundary of the corona model (Antiochos et al., 1999; van der Holst et al., 2009). This approach offers a relatively simple, and inexpensive model for CME initiation based on empirical understanding of preevent conditions. We also have a SWMF component (EE), which is a physicsbased extended MHD model (BATSRUS) of the convection zone (Fang et al., 2012b, a), where the domain is a localized wedge extending 30 Mm below the photosphere and hundreds of Mm into the corona. The wedge extends hundreds of Mm at the photosphere, sufficient to contain a large active region. The model includes optically thin radiative loss terms appropriate for the corona and empirical cooling terms to approximate optically thick radiative transfer near the photosphere, which drives cellular convection (Abbett & Fisher, 2003; Abbett et al., 2004). In the environment, a CME may be initiated by the emergence of a flux rope from the convection (Manchester et al., 2004a). Currently, the physicsbased EE model only works in a standalone mode (Fang et al., 2012b, a) and we use empirical models to generate CMEs in the SWMF (Borovikov et al., 2017a; Jin et al., 2017a).
Magneticallydriven CMEs were first modeled with the SWMF suite in the early 2000s. First, the distorted spheromactype Gibson & Low (1998) (GL) unstable fluxrope model was implemented (Manchester et al., 2004a, 2014a, b; Lugaz et al., 2005a, b). Later, the Titov & Démoulin (1999) (TD) twisted eruptive flux rope model was also added to the SWMF tool box as a CME initiation option (Roussev et al., 2003, 2007; Roussev & Sokolov, 2006). The TD eruption model was used in the first physicsbased SuntoEarth space weather simulation of two consecutive CMEs during the 2003 Halloween event (Toth et al., 2007; Manchester et al., 2008), showing quantitative agreement with several observations including insitu observations at 1 AU and coronagraph images from LASCO C2 and C3. An automated tool, the Eruptive Event Generator using GibsonLow (EEGGL) configuration was developed (Borovikov et al., 2017a; Jin et al., 2017a), and added to the SWMF suite to make CME simulations more widely available to the heliophysics community. In 2016, EEGGL was made available interactively through the CCMC’s runsonrequest service to provide CME simulations.
Representative results from EEGGLdriven CME simulations are shown in Figure 7 (Jin et al., 2017a), using a combination of two flux rope sizes and two magnetic field strength parameters. The left panel shows the initial configuration of the flux ropes with two density isosurfaces. The middle panel depicts the resulting CME evolution at 20 min. The background color shows the density ratio between the CME solution steady background solar wind. The right panel shows the synthesized (modelderived) SOHO/LASCO white light images. The color scale shows the white light total brightness divided by that of the preevent background solar wind. Comparing panels (a) and (d), we can see that with a higher magnetic field strength parameter, more plasma is added at the bottom of the flux rope (red isosurface). The second and third cases have the same magnetic field strength parameter but with different flux rope sizes. In this case, we can see the flux rope is considerably smaller at the beginning. With this smaller flux rope, the resulting CME speed is reduced and the morphology of CME in the synthesized white light image is quite different with narrower CME width angle.
Fig. 7 Three examples of Gibson & Low (1998) flux ropes with different size and magnetic strength parameters. Panels (a)–(f) and (g)–(i) show, respectively, flux ropes specified with radii of 0.8 and 0.6 R_{s}. Strength parameters are set to 0.6 for model run (a)–(c) and 2.25 for (d)–(i). The left column shows the initial configuration of the flux ropes with blue and red isosurfaces showing, respectively, the ratios of 0.3 and 2.5 of the mass density of the CME model divided by that of the preevent corona. The middle column shows the resulting CME evolution at t = 20 min. Here, magnetic field lines are colored red, grayshaded and green to illustrate the flux rope, largescale helmet streamers, and magnetic fields surrounding active regions and open flux. Color contour images show the ratio of the mass density of the CME divided by that of the preevent corona. The right column shows modelproduced SOHO/LASCO white light images, where the total brightness is normalized by dividing by that of the preevent background solar wind (from Jin et al., 2017a). 
5.3 ICME Simulation
The evolution of CMEs in the solar corona and interplanetary medium has been extensively simulated with the SWMF (Manchester et al., 2004a, 2004b, 2012; Roussev et al., 2004, 2004b, 2005, 2008; van der Holst et al., 2007, 2009; Roussev, 2008; Manchester & van der Holst, 2017). Current models (since 2014) start from the upper chromosphere with fixed temperature T = 5 × 10^{4} K and density n = 2 × 10^{17} m^{−3}. The Alfvén wave turbulence is launched at the inner boundary, with the Poynting flux scaling with the surface magnetic field. The electron and proton temperatures are solved separately. The smallest radial cell size is 10^{−3} R_{⊙} near the Sun to resolve the steep density and temperature gradients in the upper chromosphere. The initial condition for the radial magnetic field at the inner boundary is provided by synoptic/synchronic maps of the photospheric magnetic field using the Potential Field Source Surface (PFSS) model.
The inclusion of the lower corona in our model allows us to produce synthesized extreme ultraviolet (EUV) images, which are then compared with the EUV observations from SDO/AIA (Lemen et al., 2012) and STEREO/EUVI (Howard et al., 2008). Figure 8 shows an example of model results compared with observations of the 7 March 2011 CME event, which demonstrates enhanced emission from regions of the lower atmosphere compressed and heated by CMEdriven shocks and compressional waves.
Fig. 8 CMEdriven EUV waves in the simulation (left) and in the corresponding SDO/AIA observation (right). Both the simulation and observation images are produced by a triratio running difference method. The tricolor channels are AIA 211 Å (red), AIA 193 Å (green), and AIA 171 Å (blue). The ratio in each channel is identically scaled to 1 ± 0.2 for both observation and simulation (from Jin et al., 2017b). 
In addition to EUV images, our model also allows us to make synthetic Thomsonscattered white light images of the corona. Figure 9 shows a comparison between the observed white light images and the model synthesized images for the 7 March 2011 CME event (Jin et al., 2017b). The synthesized runningdifference images are able to reproduce the observed typical threepart CME structure comprising the bright core that represents the filament material, the dark cavity that corresponds to the flux rope, and the bright front that is due to the mass pileup in front of the flux rope (Illing & Hundhausen, 1985). Moreover, the model is also able to resolve the observed second faint front that is the outermost part of the increased intensity region associated with the CMEdriven shock, as was first quantitatively demonstrated in Manchester et al. (2008). The white light comparison from three points of view confirms that the simulated CME propagates in the observed direction. The model results in Figures 8 and 9 are produced by running the AWSoM with the magnetic field specified by GONG synoptic magnetograms for CR2107 and synchronous magnetograms for the month of September 2014, respectively.
Fig. 9 Comparison showing a general agreement between the whitelight observations from SOHO LASCO C2 (top left) and STEREOB COR1 (top right) and the respective synthesized whitelight images from the simulation (bottom). The color contours show the relative total brightness changes compared to the preevent background level (from Jin et al., 2018). 
EEGGL was designed to provide datadrive CME simulations that are capable of reproducing the solar wind disturbances at 1 AU that generate geomagnetic storms. To achieve this goal, the model must capture the bulk plasma properties, in particular the plasma velocity, mass density and magnetic field. An example of this capability is shown in Figure 10, where we show the simulated (shown with dashed lines) and L1observed plasma conditions (shown with solid lines) resulting from the Earthdirected CME that occurred on 12 July 2012. Here, timeseries data are shown (top to bottom) for the Cartesian components of the magnetic field, mass density and Earthdirected velocity. We shift the simulated time by roughly 10 h to provide a better comparison with observations. We find that the magnetic x and y components appear to be missmatched while the z component very well matches the observed magnitude and time profile of the observations. The velocity roughly matches the increase from the ambient background to the shocked value found in the sheath region, but then increases above observed values in the relaxation region. The model delivers mass density, early velocity and stormdriving B_{z}, which allows the model to successfully drive a magnetospheric simulation, while issues with flux rope rotation and streaminteraction remain to be addressed. This EEGGLdriven simulation was performed on demand at the CCMC where the model outputs are available to the public.
Fig. 10 1 AU results of the EEGGL simulation of the 12 July 2012 CME event simulated at CCMC. Shown are the simulated and observed plasma quantities plotted with dashed and solid lines, respectively. From top to bottom are the magnetic field component B_{x}, B_{y} and B_{z}, the mass density, and the Earthdirected velocity V_{x}. Simulation results are shifted 10 h to match the shock arrival. We find good agreement, with the exception of the B_{x}–B_{y} rotation and the excessive trailing velocity. 
5.4 Solar energetic particle simulations
The acceleration of energetic particles in a CMEdriven shock and the subsequent transport processes are modeled using the MFLAMPA module in SWMF (Sokolov et al., 2004; Borovikov et al., 2018). The distribution function of energetic particles are solved on a multitude of extracted magnetic field lines advecting with the background plasma (Lagrangian grids) (Sokolov et al., 2004). MFLAMPA is fully coupled with the solar corona (SC), inner heliosphere (IH), and the outer heliosphere (OH) components. The plasma and turbulence parameters along the magnetic field lines are extracted dynamically from the the BATSRUS simulations.
Figure 11 shows the application of MFLAMPA to model the acceleration and transport processes of energetic particles in an SEP event that occurred on 23 January 2012 (Borovikov et al., 2018). The ambient solar corona and interplanetary steadystate solar wind background are obtained as discussed in Section 5.1 and the CME, which is the source of this SEP event, is simulated by inserting a fluxrope into the active region on the Sun using the EEGGL model (see Sect. 5.2). In Figure 11, the green isosurface represents the leading edge of the CME. Hundreds of magnetic field lines whose footpoints on the solar surface are close to the active region are extracted using the coupled AWSoMR, EEGGL, and MFLAMPA modules. Left and right panels are at 10 min and 20 min after the CME eruption, respectively. The colors on the magnetic field lines represent the flux, in the unit of particle flux unit (pfu, particles/cm^{2}/s/sr) of the energetic protons, whose energies are greater than 10 MeV. Along single field lines, the proton flux is larger in the region close to the CME, where the acceleration takes place. And the flux decreases away from the CME when the accelerated protons stream into interplanetary space. The proton’s flux is higher at the center of the CME than at the flank, indicating a stronger acceleration at the center where the compression is larger.
Fig. 11 Distribution of the energetic particles (> 10 MeV) along the extracted magnetic field lines at 10 min (left panel) and 20 min (right panel) after the eruption of CME. The flux is in the unit of particle flux unit (pfu, particles/cm^{2}/s/sr). The green isosurface represents the leading edge of the CME. 
Figure 11 demonstrates the capability of using the selfconsistent physicsbased modules in SWMF to calculate the flux of the energetic particles at any location in the heliosphere, showing it to be a powerful tool to study the acceleration and transport processes of SEP events.
5.5 Rigidity cutoff simulations
Overall, the Earth’s radiation environment is very dynamic. Such fluxes of the energetic ions (above 1 MeV per nucleon) can be enhanced by several orders of magnitude during SEP events, which can last from a few hours to a week (Baker & Kanekal, 2008). SEPs are energetic particles ejected by the Sun in events that are correlated with coronal mass ejections (CMEs) and solar flares (Reames, 1999). The occurrence of SEPs is in positive correlation with ongoing solar activity.
The most stable component of the Earth’s radiation environment, galactic cosmic rays (GCRs), varies by an order of magnitude at energies below a few hundred MeV per nucleon due to heliospheric modulation (cf. Vainio et al., 2008). Variability of GCRs observed in the Earth’s magnetosphere is due to a combined effect of the IMF in the heliosphere and the geomagnetic field inside the magnetosphere on the GCR transport.
The Earth’s magnetosphere presents a shield against GCRs and SEPs. Those particles with energies below 100 MeV/n are effectively blocked by the Earth’s magnetosphere (Badavi et al., 2011). Usually, the geomagnetic interaction of SEPs and GCRs is described in terms of rigidity, R (momentum/unit charge) rather than energy. Transport of SEPs and GCRs in the geospace is a kinetic process due to a significant value of particles’ gyroradius that can reach the value of tens of Earth’s radii. An example of GCR’s proton gyroradius calculated for quiet geomagnetic conditions is presented in Figure 12. One can see that even for particles that are on the lower end of the energetic spectrum of SEPs and GCRs penetrating in the geospace, the gyroradius can be as large as tens of Earth’s radii, meaning that in practical calculations, kinetic methods that account for the gyromotion of the energetic particles must be employed. The effect of the gyromotion on the topology of the SEPs’ population in geospace is illustrated in Figure 13, which shows the density of SEPs in the plane orthogonal to the equatorial plane and the SEP density’s isosurface in geospace.
Fig. 12 Example of gyroradii of particles with 1 MeV < E < 16 MeV during quiet geomagnetic conditions. Left: Gyroradius map in the equatorial plane. Right: Gyroradius map in the meridional plane (X = 0). The gyroradii of these particles can be as large as tens of R_{E}. Here, Xaxis is directed toward the Sun, and Yaxis is in the equatorial plane, and Zaxis is such that the frame of reference is righthanded. The freespace energy spectrum of the simulated energetic particles is taken from Badavi et al. (2011). 
Fig. 13 Example of the calculated density of energetic protons with energies 1 MeV < E < 100 MeV) in geospace. Both the SEP’s energy spectrum and geomagnetic parameters are taken for quiet conditions. The figure demonstrates that the topology of the SEPs population in the geospace is affected by the particles’ gyromotion. Here, Xaxis is directed toward the Sun, Yaxis in the equatorial plane, and Zaxis is such that the coordinate frame is righthanded. 
An example of calculating cutoff rigidity detailed by Tenishev et al. (2021) is presented in Figure 14. The calculation is done using the Adaptive Mesh Particle Simulator (AMPS) employing particle timebackward tracing starting from the altitude of 500 km. The calculations presented in the figure were performed for quiet geomagnetic conditions (p_{SW} = 2 nPa, Dst = 1 nT, B_{y} = − 0.08 nT, and B_{z} = 2 nT) and for the conditions during the geomagnetic storm on 17 March 2015 (p_{SW} = 10 nPa, Dst = −200 nT, B_{y} = −7 nT, and B_{z} = −10 nT). The left panel of Figure 14 shows the rigidity cutoff map before the storm. The right panel shows the relative depression during the storm. The value shown in Figure 14 is the ratio of the cutoff rigidity difference during the event to its original value. The relative depression of −1 means that the corresponding location becomes magnetically connected to the interplanetary magnetic field during the simulated geomagnetic storm. One can see that the general rigidity cutoff patterns have changed mainly in the midlatitude region.
Fig. 14 Example of applying AMPS for rigidity cutoff calculation. The map is calculated for an altitude of 500 km. Left: Rigidity cutoff map calculated for quiet geomagnetic conditions. Right: Depression of the rigidity cutoff during a geomagnetic storm. The calculation was performed for conditions of the geomagnetic storm on 17 March 2015. One can see that the general rigidity cutoff patterns have changed mostly in the midlatitude region (from Tenishev et al., 2021). 
5.6 Mesoscale resolving magnetosphere simulations
While the MHD plasma description has inherent restrictions in describing the microscale processes (see Appendix E for treatment of kinetic processes), BATSRUS, when run with high spatial resolution in key portions of the geospace, can easily resolve the KelvinHelmholtz instability and fluxtransfer events (FTEs) (Kuznetsova et al., 2009), found e.g., at the magnetospheric boundary. Highresolution MHD simulations in the magnetotail can reproduce intricate details of the interchange instability, bursty bulk flows, and other processes (Yu et al., 2017). The adaptive mesh refinement (AMR) guarantees that the run times, while higher for high resolution, remain manageable, as the increase in number of computational cells only increases by about a factor of a few.
An example of a very highresolution simulation is shown in Figure 15. The SWMF/BATSRUS simulation was run with 1/16 R_{E} grid resolution in the tail and magnetopause region in order to resolve small and mesoscale structures in the magnetosphere. The results demonstrate the formation of Kelvin–Helmholtz vortices at the flanks of the magnetopause in response to the solar wind flow past the magnetic boundary. Furthermore, it was shown that reducing the resistivity in the model led to structuring of the reconnection in the magnetotail and the formation of narrow, elongated flow channels (or bursty bulk flows, Angelopoulos et al., 1994) throughout the width of the tail (Haiducek et al., 2020).
Fig. 15 Results of highresolution SWMF/BATSRUS simulation with 1/16 R_{E} grid resolution in the tail and magnetopause region in order to resolve small and mesoscale structures. It shows the filamentary current structures associated with the presence of both Kelvin–Helmholtz vortices at the flanks of the magnetopause and of bursty structures in the tail. Note that the associated flow velocities are not shown (Dorelli & Buzulukova, personal communications, 2020). 
Fig. 16 An illustration of cusp accelerated ion outflow due to soft electron precipitation and waveparticle interactions as modeled by PWOM; adapted from (Glocer et al., 2018). The top panel shows the PWOM computed ion distribution function along a single field line and demonstrates the ion conic evolution with altitude. Here v_{∥} and v_{⊥} are the parallel and perpendicular velocity and the color contour shows the log of phase space density in normalized, dimensionless, units. The bottom panel shows the global simulation, including the kinetic processes in the cusp using 896 field lines. Here u_{O} and u_{H} are the bulk field aligned velocities for O^{+} and H^{+} in units of km/s, and T_{e} is the electron temperature in Kelvin. 
Figure 15 shows the current density in the equatorial plane during a geomagnetically active period. The filamentary current structures on the magnetopause and in the magnetotail are indicative of Kelvin–Helmholtz instability and mesoscale bursty bulk flows, respectively. The associated flow velocities for these structures are not shown. However, in this simulation it was found that while the main flow direction in the more distant magnetotail continues to be Earthward, the reconnection onset at the boundary of the quasidipolar and taillike magnetosphere creates tailward flows that strengthen at substorm onset (Dorelli & Buzulukova, personal communications, 2020). Such simulations are sufficiently accurate that they can be used to reassess the substorm theories (e.g., Baker et al., 1996; Angelopoulos et al., 2008).
5.7 Ionospheric outflow simulations
Observations show that during geomagnetic storms O^{+} can comprise as much as 40%–80% of the ion density in the nearequatorial magnetosphere inside of 15 R_{⊙} (Lennartsson et al., 1981). As O^{+} can only originate in the ionosphere, its observed presence in the magnetosphere during storms is a clear indicator of the importance of the ionosphere in supplying magnetospheric plasma during space weather events. This is true not only of O^{+}, but it has been estimated that 65% of the H^{+} population during geomagnetic storms may also be of ionospheric origin (Gloeckler & Hamilton, 1987).
The Polar Wind Outflow Model (PWOM) supplies the PW component in SWMF that calculates the transport of plasma from the ionosphere and sets the supply for the magnetosphere. This model solves the gyrotropic transport equations (Gombosi & Nagy, 1989) for multiple ion species from below the F2 peak to much higher altitudes. This model was expanded from solving a single field line to solving multiple field lines in order to reconstruct the global 3D outflow distribution (Glocer et al., 2009b). The ability of the model to represent different critical drivers of ion outflow has also grown in recent years. The inclusion of various treatments of superthermal electron populations (photo, auroral, and secondary electrons) to PWOM has improved the model in comparison with observations (Glocer et al., 2012, 2017). Most recently, PWOM has been expanded to move to a hybrid PIC description above 1000 km while maintaining a fluid description at lower altitudes (see Fig. 16, adapted from Glocer et al., 2018). The latter expansion allows PWOM to treat waveparticle interactions due to processes like Ion Cyclotron Resonant Heating, which is thought to be a major mechanism in creating ion conic distributions and energized O^{+} escape.
Simulations with SWMF are able to track the plasma calculated to escape the ionosphere throughout the magnetosphere using BATSRUS when configured with multifluid or multispecies MHD. Using separate fluids or species for each constituent plasma population enables us to track the impact of ion outflow on the magnetosphere.
SWMF in the configurations described above has been used to study the many impacts on ion outflow. The effect of ion outflow on the magnetic field at geosynchronous orbit and on the cross polar cap potential was studied and found to improve the prediction of the magnetic field as well as lower the cross polar cap potential (Glocer et al., 2009a; Welling & Zaharia, 2012). The contribution of ion outflow to the ring current was examined and found to be a major contributor to the total ring current energy content during storms (Ilie et al., 2015; Welling et al., 2015; Glocer et al., 2018, 2020). These codes, coupled together via the SWMF, have also been used to study ring current ENA observations from the TWINS mission (Ilie et al., 2013). These studies are only a subset of the total number of SWMF studies in this area and should not be taken as an exhaustive list.
5.8 Mesoscale ionosphere simulations
The SWMF/Geospace is an important tool in the analysis of the polar region electrodynamics. The model’s advantage is that it can produce superior spatial coverage for the magnetic disturbances whose observations are limited by the oceans and access to remote locations, and better spatial coverage for the fieldaligned currents than those derived from spaceborne magnetic measurements by the AMPERE project (Anderson et al., 2000).
Figure 17 shows simulated fieldaligned currents (FACs) and ground magnetic perturbations due to a solar wind pressure enhancement. The twoway coupled BATSRUS, CRCM, and IE modules are utilized in this run with 1/8 R_{E} resolution used from the dayside magnetopause to the nearEarth magnetotail to capture the magnetosphere reconfiguration due to compression and the subsequent relaxation. Several hundred virtual magnetometers have been included in the simulation at the locations of real magnetometers, and a uniformly distributed array covering the globe at a resolution of 4° in latitude and 12° in longitude. The left panel in Figure 17 shows the transient FACs during the Preliminary Impulse (PI) phase with the ionospheric convection contours superimposed on top. The location of the Poker Flat magnetometer is denoted by the cyan dot. The H component magnetic perturbation contours calculated from the uniform magnetometer array and a comparison of the time series of the simulated and observed H component perturbation at Poker Flat are shown in the middle and right panels of Figure 17, respectively.
Fig. 17 FACs and ground magnetic perturbations due to shock compression on 17 March 2015 from a highresolution BATSRUS run coupled with CRCM and IE. The simulated magnetic perturbations at Poker Flat were compared with the real magnetometer observations. The results shown here illustrate the capacity of the SWMF to resolve mesoscale features of the magnetospheric dynamics in highresolution MHD (Zou et al., 2017). 
The simulated magnetic perturbations matched both the polarity and the magnitude of the H component perturbation very well, suggesting the coupled models were able to capture the source, propagation, and closure of the compressioninduced mesoscale fieldaligned currents. The results of this coupled geospace run were then used to drive the Global Ionosphere and Thermosphere Model (GITM), and revealed a shortlived mesoscale fast flow channel in the ionosphere with intense Joule heating and sudden ion temperature enhancements (Ozturk et al., 2018). These simulation results provided a valuable explanation for a transient ion upflow event observed by the Poker Flat Incoherent Scatter Radar (PFISR) (Zou et al., 2017).
Other important mesoscale features in the coupled geospace system have also been simulated by using the SWMF, including subauroral polarization streams (SAPS) (Yu et al., 2015), and boundary flows between the Region 1 and Region 2 FACs (Wang et al., 2019).
5.9 Geomagnetic index simulations
For operational space weather forecasters, geomagnetic indices are a standard tool when distributing forecasts and warnings to users representing both spacecraft and ground system (e.g. power network) operators. These include K and Kp, Disturbance Storm Time (Dst) and the related SYMH index, and the Auroral Electrojet (AE, AU, AL, and AO) indexes (Rostoker, 1972; Mayaud, 1980; Menvielle et al., 2010). The validity of the SWMF/Geospace suite can be checked by reconstructing these indices from the simulation output data and comparing them with observations. Dstequivalent outputs have long been a staple of SWMF simulations, serving as a quicklook diagnostic of inner magnetosphere performance (e.g., Zhang et al., 2007; Welling et al., 2011). K, Kp, and AE indexes were added more recently and build off of internal virtual magnetometer capabilities.
Several studies have focused on the performance of the virtual geomagnetic indexes produced by the SWMF. Glocer et al. (2016) explored the local Kindex predictive skill of the SWMF and other models, demonstrating the SWMF’s strong capabilities and reproducing this value. Haiducek et al. (2017), simulated the month of January 2005 using the observed solar wind data as input. The simulation was run with two different grid resolutions. The model was found to predict the ring current index SYMH to good accuracy, with a root mean square error of less than 20 nT (see left panel of Fig. 18). The geomagnetic index Kp performed well during storm time, but predicted larger than observed activity during quiet times. On the other hand, the auroral electrojet index AL was predicted reasonably well on average, but was systematically less negative than the observed values during high geomagnetic activity. While the grid resolution caused only small variations to the results, runs without the inner magnetosphere component were not able to produce the storm dynamics (Haiducek et al., 2017). Haiducek et al. (2020) further explored virtual AL performance during substorm activity, finding substormrelated perturbations to be weaker than observations. Liemohn et al. (2018) continued along the same lines and used nearly three years of simulated geomagnetic index data from the experimental realtime SWMF runs at CCMC. The right panel of Figure 18 shows scatter plots of observed and simulated Dst (hourly averaged SYMH) values for the different simulation setups. They examined different metrics of success, and conclude that the correlation coefficient between the observed and model values was 0.69, the prediction efficiency was 0.41, and the Heidke skill score was 0.57 for an event threshold of −50 nT.
Fig. 18 (Left) Comparison of the storm SYMH index for two storms from January 2005. The thick trace shows the observations, the red and orange traces show the normal resolution and high resolution SWMF/Geospace simulations, and the blue trace shows the SWMF/Geospace simulation run without the inner magnetosphere RCM component (from Haiducek et al., 2017). (Right) Scatter plot of the observed real time Dst time series (vertical axis values) against a SWMF prediction of the Dst (horizontal axis values). The green, orange, and cyan contour lines show the regions of 40, 80, and 120 values within a 2.6 × 2.6nT bin. The diagonal longdashed line shows the ideal datamodel relationship. The two dashed lines show the −50 nT threshold values (from Liemohn et al., 2018). 
Overall, all these studies confirm that the Geospace model represents a reasonably accurate approximation to the real magnetosphere during a large variety of circumstances. This provides confidence in the more detailed predictions, such as local magnetic disturbance levels around the globe or the plasma parameters near the geosynchronous orbits. In CCMCled modeling challenges that focus on geomagnetic index comparisons, the SWMF is consistently among the best of the global models (e.g., Pulkkinen et al., 2013; Rastätter et al., 2013; Glocer et al., 2016). While more accurate models exist for predicting and forecasting geomagnetic indices, in particular those based on machine learning algorithms, the SWMF is, todate, the most accurate reproduction of these indices from a solarwindtoionosphere firstprinciples physicsbased model of the full geospace system.
6 Resolving kinetic scales in global simulations
6.1 MHDEPIC and MHDAEPIC
Kinetic models have been used for a long time to model the inner magnetosphere (cf. Wolf et al., 1982; Buzulukova et al., 2010), the radiation belts (cf. Fok et al., 2008) or the transport and diffusion of energetic particles along field lines (cf. Sokolov et al., 2004). These kinetic models use some simplifying assumptions, such as restricting the motion of particles along field lines and ignoring feedback to the magnetic field, to drastically reduce the computational cost. Solving the full kinetic equations in 7 dimensions (one temporal, 3 spatial and 3 velocity) in a global simulation is simply not feasible on the current or even near future supercomputers.
The idea of coupling MHD and kinetic PIC models has been around for a long time (cf. Sugiyama & Kusano, 2007) but making this work in 3D has been an elusive goal. In fact, many in the community argued that coupling fluid and kinetic models is impossible as they are simply not compatible with each other. The MHD and algorithm experts at Michigan initiated a collaboration with several PIC modelers, including Giovanni Lapenta, Stefano Markidis and Jeremiah Brackbill. It took years of working together to overcome all the obstacles. Some were seemingly simple, like converting units, and still took a long time. Others were much more complicated, such as keeping the PIC model stable and suppressing various instabilities, or avoiding discontinuities developing at the interface of the MHD and PIC regions. We found, for example, that using Hall MHD instead of ideal MHD improves stability, or using hyperbolicparabolic cleaning in addition to the 8wave scheme is necessary to eliminate accumulation of errors near the boundaries of the PIC region.
Eventually our work yielded results: the MHD with embedded PIC (MHDEPIC) model became reality (Daldorff et al., 2014). It took a few more years to efficiently couple the models through the SWMF using a newly developed efficient parallel coupler, allow for different grids and different time steps, allowing for multiple PIC domains and generalizing the fluid model from single fluid (Hall) MHD to multispecies and multifluid MHD, as well as the five and sixmoment fluid equations (cf. Chen & Toth, 2019; Zhou et al., 2020).
Running MHDEPIC simulations for long simulation times also revealed hidden issues with the PIC algorithms that could be avoided in standalone PIC simulations by careful tuning of various parameters, but were plaguing the more complicated MHDEPIC simulations. We overcame these issues by developing the Gauss Law satisfying ECSIM (GLECSIM) algorithm (Chen & Toth, 2019) that conserves energy and charge at the same time. This new PIC algorithm, coupled with the extended MHD code, has finally delivered an accurate and reliable MHDEPIC model.
Toth et al. (2017) showed that the kinetic scales can be artificially changed by changing the mass per charge ratio of the ions and electrons and still obtain correct global solutions as well as correct, but scaled, kinetic solutions. The only limitation is that the modified kinetic scales should still be well separated from the global scales. For example, one can increase the kinetic scales by a factor of f = 16 and thus reduce the computational cost of the PIC model by a factor of f^{4} ~ 65,000. With such scaling it became possible to simulate Earth’s magnetosphere with the MHDEPIC model. Chen et al. (2017) modeled the kinetic reconnection process at the dayside magnetopause of Earth in a global simulation. The model correctly reproduced the properties of flux transfer events (FTEs) and revealed several new insights into the birth, development and final fate of FTEs starting from the kinetic scales and growing to the global scales.
While MHDEPIC with kinetic scaling opened the possibility of combining kinetic modeling with global simulations of Earth’s magnetosphere dynamics, the simulations were still very expensive. This is especially true for the magnetotail, where the reconnection sites can move in a large volume due to the intrinsic dynamics of the reconnection Xlines, as well as to the flopping of the magnetotail caused by the changing solar wind.
To further improve the efficiency of the model, we have developed the MHD with an adaptively embedded PIC (MHDAEPIC) algorithm. The main idea comes from the blockadaptive mesh and the hybrid schemes used in BATSRUS: the PIC grid is decomposed into small blocks that can be activated and deactivated dynamically. While the idea is straightforward, the implementation is not. We had to abandon iPIC3D that uses a single grid, and implement the GLECSIM algorithm into the Adaptive Mesh Particle Simulator (AMPS) code (Tenishev et al., 2021). The resulting MHDAEPIC model can achieve an order of magnitude or even more speedup compared to the MHDEPIC model that uses static PIC domains. We also developed a new PIC code, the Flexible Exascale Kinetic Simulator (FLEKS), to be used in MHDAEPIC. FLEKS is based on the AMReX library (Zhang et al., 2019b, 2020) and it was designed for flexibility and high performance with a stateoftheart semiimplicit PIC algorithm. Novel particle splitting and merging algorithms have been designed for FLEKS to control the number of macroparticles per cell during long MHDAEPIC simulations.
In general, the MHD(A)EPIC model offers a powerful tool to study magnetospheric physics. The latest application is covering the tail reconnection site with an adaptive PIC region so that one can study geomagnetic storms and substorms in a more realistic way. Figure 19 shows an example of an adaptive PIC region, which tracks the motion of the magnetotail reconnection site during a storm simulation.
Fig. 19 MHDAEPIC simulation of a geomagnetic storm. The 2D cuts display a part of the much larger 3D domain. Colors show the X component of the velocity and the white lines are traces of the magnetic field in the meridional plane. The black lines in the two snapshots, separated by 159 minutes of simulation time, indicate the edges of the active PIC regions. As the tail evolves, the active PIC region is continuously adapted to cover the reconnection sites of interest. 
6.2 MHDEPIC results
MHDEPIC has been used to simulate the terrestrial magnetosphere (Chen et al., 2017, 2020; Jordanova et al., 2018), the interaction of Mercury (Chen et al., 2019c) and Mars (Ma et al., 2018b) with the solar wind and the minimagnetosphere of Ganymede (Toth et al., 2016; Zhou et al., 2019, 2020). To demonstrate the capabilities of MHDEPIC, here we show some of the Earth magnetosphere simulation results by Chen et al. (2017).
An overview of the evolution of the dayside magnetopause is shown in Figure 20, which contains the Hall magnetic field B_{y} and the field lines at the meridional plane inside the PIC box. At t = 70 the Hall field extends far away from the X line with roughly the same field strength for each branch. Fifteen seconds later, south of the existing reconnection point, another X line starts to form. At t = 145 s, both X lines can be seen clearly, and a flux ropelike structure forms between the two X lines. At t = 325 s, the flux rope moves away from the top X line and the current sheet between them becomes unstable and a secondary flux rope is generated. During the onehour simulation, flux ropes form near the subsolar point and move toward the poles quasiperiodically.
Fig. 20 Snapshots showing By strength (color) and the projected magnetic field lines in the meridional plane inside the PIC region. The color bar is different in each plot (from Chen et al., 2017). 
Crescent shape electron phase space distribution has been observed near the electron diffusion region at the dayside magnetopause by MMS Burch et al. (2016). The same distribution is also found in the 3D MHDEPIC simulation (see Fig. 21). The phase space distribution of electrons inside a cube region on the dayside magnetopause is shown in Figure 21b. The crescent distribution is found in the V_{y}–V_{x} plane, corresponding to the two velocity components perpendicular to the magnetic field. The crescent hot electrons are drifting along the negative y direction with a speed close to 3000 km/s. The direction of the flow is consistent with the E × B direction, and the velocity of the crescent particles is very close to the MMS observed by Burch et al. (2016). Slightly farther away from the reconnection site, where the Larmor field appears, the ion phase space distribution also presents a crescentlike shape, as is shown in Figure 21c. The crescent ions drift in a positive y direction because E_{x} is negative.
Fig. 21 Crescent electron and ion phase space distributions (a) E_{x} (mV/m) in the meridional plane at t = 3600 s; (b) Normalized electron distribution in V_{y}–V_{x} phase space; and (c) Ion phase space distribution. The phasespace density is normalized (from Chen et al., 2017). 
7 Planetary environments and solar analogs
Space weather phenomena at other solar system objects and at astropheres are the subject of increasing interest (cf., Lilensten et al., 2014; Plainaki et al., 2016; André et al., 2018). More recently, space weather phenomena in astropheres harboring extrasolar planets also became the focus of investigations (cf., Pillet et al., 2019).
The SWMF simulation suite has been used to simulate the space environment of most solar system planets, including Mercury (e.g., Kabin et al., 2000a, 2008; Jia et al., 2015, 2019), Venus (e.g., Bauske et al., 1998; Ma et al., 2013), Mars (e.g., Liu et al., 1999; Bauske et al., 2000; Ma et al., 2002, 2018a; Liemohn, 2006; Regoli et al., 2018), Jupiter (e.g., Cravens et al., 2003; Sarkango et al., 2019), Saturn (e.g., Hansen et al., 2000, 2005; Gombosi & Hansen, 2005; Glocer et al., 2007; Gombosi & Ingersoll, 2010; Jia et al., 2010a, b; Zieger et al., 2010), and Uranus (Toth et al., 2004).
Moreover, the SWMF has been applied to comets (e.g., Gombosi et al., 1996, 1999; Häberli et al., 1997; Hansen et al., 2007; Rubin et al., 2012; Gombosi, 2015; Huang et al., 2016a, b), and planetary moons including Io (e.g., Combi et al., 1998; Kabin et al., 2001), Europa (e.g., Kabin et al., 1999a; Liu et al., 2000; Rubin et al., 2015; Jia et al., 2018; Harris et al., 2021), Ganymede (e.g., Toth et al., 2016; Zhou et al., 2019, 2020), Titan (e.g., Kabin et al., 1999b, 2000b; Nagy et al., 2001; Ma et al., 2007), and Enceladus (e.g., Jia et al., 2010a, b, c).
Finally, the SWMF suite of models has also been applied to the outer heliosphere (e.g., Opher et al., 2007, 2009, 2016, 2017) and astrospheres (e.g., Cohen et al., 2010, 2015, 2020; AlvaradoGómez et al., 2020)
As discussed above, the core MHD code within the SWMF, BATSRUS, can be configured to solve the governing equations of ideal MHD, resistive MHD, semirelativistic MHD, multifluid MHD, MHD with anisotropic pressure, or highorder moment MHD. In addition to the basic equations, there are various source and loss terms included in BATSRUS that change from application to application (see details in Appendix C.1.3). Most relevant to our applications to the giant planet magnetospheres (e.g., Jupiter and Saturn) is the capability of including various massloading processes (ionization, chargeexchange, dissociative recombination, etc.) arising from the internal plasma sources associated with planetary moons (Jupiter: Sarkango et al., 2019; Saturn: Zieger et al., 2010; Jia et al., 2012a, b). In modeling planetary magnetospheres with an ionosphere, BATSRUS is normally coupled to the Ionospheric Electrodynamics (IE) module to simulate magnetosphereionosphere coupling. For planetary objects that do not possess a significant atmosphere/ionosphere, such as Mercury and Ganymede, we have extended the BATSRUS MHD model to include the planetary interior as part of the simulation domain such that the influence of the electrical conductivity of the planetary interior on the space environment can be directly modeled. We have applied such a model successfully to the magnetospheres of (Jia et al., 2015, 2019; Chen et al., 2019c) and Ganymede (Zhou et al., 2019, 2020) where the induction effect of the subsurface conducting region (conducting core in the case of Mercury, and subsurface ocean in the case of Ganymede) plays an important role in the global magnetospheric interaction.
Here we show the results from two Mercury simulations that demonstrate the flexibility and capabilities of SWMF. The first is an extended MHD simulation that takes into account the finite conductivity of Mercury’s interior (Jia et al., 2015, 2019), while the second is an MHDEPIC simulation that takes into account kinetic effects (Chen et al., 2019c).
A unique aspect of Mercury’s interaction system arises from the large ratio of the scale of the planet to the scale of the magnetosphere and the presence of a large size core composed of highly conducting material. Consequently, there is strong feedback between the planetary interior and the magnetosphere, especially under conditions of strong external forcing. In applying the SWMF simulation suite to Mercury, Jia et al. (2015), used the resistive MHD version of BATSRUS to develop a global magnetosphere model in which Mercury’s interior is modeled as layers of different electrical conductivities that electromagnetically couple to the surrounding plasma environment. One particular advantage of this model is its ability to characterize the dynamical response of Mercury to time varying external conditions in a selfconsistent manner, such as the induction effect at the planetary core. To demonstrate this capability, we have performed a series of idealized simulations as well as simulations of MESSENGER events for a wide range of upstream solar wind conditions (Jia et al., 2015, 2019). Our results show that, due to the induction effect, Mercury’s core exerts strong global influences on the way Mercury responds to changes in the external environment, including modifying the global magnetospheric structure and current systems as well as affecting the extent to which the solar wind directly impacts the surface. There results have important implications for understanding the role of space weathering in generating Mercury’s tenuous exosphere (e.g., Jia et al., 2019).
As an example, we present results from the MESSENGER M2 flyby simulation by Jia et al. (2015), to illustrate the global configuration of the modeled magnetosphere. Figure 22a shows the model results in the noonmidnight meridional plane, where color contours of plasma pressure, along with projections of sampled magnetic field lines, are plotted to delineate the magnetospheric configuration. One notable feature of the modeled magnetosphere is the pronounced asymmetry with respect to the planet’s equatorial plane, which arises due to the northward offset of the internal dipole as well as the presence of a strong IMF B_{x}, which is typical at Mercury’s orbit. Some general features of the modeled magnetosphere can be compared directly with MESSENGER observations, such as the location and shape of various important boundaries. For the upstream conditions used in this simulation, the magnetopause and bow shock standoff distances are about 1.5 R_{M} (R_{M} = 2440 km is Mercury’s mean radius) and 1.9 R_{M}, respectively, in accordance with MESSENGER observations of the magnetosphere under similar upstream conditions (Winslow et al., 2013). Also plotted in Figure 22a is the modeled magnetopause boundary (magenta line) identified based on the total current density. For comparison, the empirical magnetopause model of Winslow et al. (2013), constructed based on MESSENGER data, is also plotted. As can be seen, the overall shape of the magnetopause boundary in our model is in general agreement with the databased model.
Fig. 22 Cuts through the simulation of the MESSENGER M2 flyby. (a) XZ cut at Y = 0 (noonmidnight meridian) with color contours of plasma thermal pressure. The thick magenta line shows the modeled magnetopause boundary identified based on the current density, while the white dotted line represents the databased empirical magnetopause model of Winslow et al. (2013). (b) YZ cut at X = 0 (terminator plane) with color contours of the x component of the flow velocity (V_{x}). The horizontal line at Z = 1.3 R_{M} is color coded by the y component of the convectional electric field (E_{y}), from which the cross polar cap potential is calculated (from Jia et al., 2015). 
Figure 22b shows a YZ cut at X = 0 through the simulation in which the color contours represent the xcomponent of the modeled plasma flow velocity (V_{x}) and the lines with arrows are sampled field lines. Key regions of the interaction can be readily identified based on the V_{x} contours. The transition from the ambient solar wind speed of ~ 400 km/s to ~ 200 km/s, which is characteristic of the magnetosheath flow at the terminator, marks the boundary of the bow shock, whereas further inward the transition from the sheath flows to convection flows with much smaller speeds (<~100 km/s) marks the magnetopause boundary. Inside the magnetosphere, flows with negative at high latitudes are the crosspolar cap flows moving in the antisunward direction, while the flows with positive at low latitudes are those return flows convecting from the night side to the dayside. From the results shown in Figure 22b we can obtain the crosspolar cap potential in the model, which provides a global measure of the strength of the coupling between the magnetosphere and the solar wind. The total potential drop in our model is about 25 kV, in reasonable agreement with the 30 kV estimated by Slavin et al. (2009), for this flyby based on MESSENGER observations.
To demonstrate the model’s ability to simulate the induction effect, we examine the global response of Mercury’s magnetosphere to solar wind compression. Figure 23 shows the results extracted from a timedependent simulation for the Highly Compressed Magnetosphere (HCM) event observed by MESSENGER on 23 November 2011, which was produced by the passage of a CME (Jia et al., 2019). Figure 23a shows the zcomponent of the magnetic field perturbations, B_{1z}, which result from various current systems, including the ChapmanFerraro currents, the tail current sheet, and the induction currents in the core, all of which are discernible in Figure 23b. As expected for this HCM event, both the magnetopause and the tail current sheet are displaced very close to the surface. The subsolar magnetopause standoff distance in the simulation is ~ 1.12 R_{M}, which is in excellent agreement with the distance of 1.13 R_{M} determined by Slavin et al. (2014), for this event based on MESSENGER observations. For this CME event, the current sheet on the night side almost reaches the planetary surface with its inner edge at only ~ 1.1 R_{M}. As illustrated by the cyan and yellow colors at the core boundary in Figure 23b, strong currents flowing in the direction as indicated by the magenta arrows are induced to prevent the external magnetic perturbations from penetrating into the conducting core. On the day side, the induced currents, together with the intensified ChapmanFerraro currents, produce the strong positive B_{1z} perturbations, which are present throughout Mercury’s resistive mantle (Fig. 23a). The positive B_{1z} perturbations in the dayside magnetosphere and inside the mantle reach amplitudes between ~150 and 200 nT. On the night side, the crosstail currents generate negative B_{1z} perturbations planetward of the current sheet with an average amplitude of −100 nT. The intensification and displacement of the tail current sheet toward the planet also induce currents flowing at the core boundary that act to negate the effect of external variations. By driving the simulation with different upstream solar wind conditions, Jia et al. (2019) conducted a systematic numerical experiment to establish global context for interpreting the HCM events observed by MESSENGER during its entire mission. Their results also provide a quantitative assessment of the relative importance of the shielding effect from induction and the erosion effect from magnetopause reconnection.
Fig. 23 A closeup view of the simulated magnetospheric configuration for the Highly Compressed Magnetosphere (HCM) event observed by MESSEGER on 23 November 2011. (a) Color contours of the zcomponent of the magnetic field perturbation (B_{1z} in nT) shown in the XZ plane and also on a 3D sphere of radius ~0.8 R_{M} corresponding to the coremantle boundary. (b) Same as (a) but for the current density in the ydirection (J_{y} in nA/m^{2}), indicating enhanced magnetopause and tail currents as well as the induction currents at the core in response to solar wind compression. In both panels, the black lines with arrows show projections of sampled magnetic field lines onto the XZ plane and the magenta lines with arrows show sampled streamlines of the induction currents generated at the core. The green circle of radius 1 R_{M} represents the planetary surface (from Jia et al., 2019). 
As established by numerous MESSENGER observations, the dynamics of Mercury’s magnetosphere is driven predominantly by magnetic reconnection, largely due to the IMF and solar wind conditions typically present in the inner heliosphere. Recognizing the importance of reconnection at Mercury, we have adapted the MHDEPIC model to simulate Mercury’s magnetosphere so that reconnection can be treated using a physicsbased model (Chen et al., 2019c), rather than through numerical or ad hoc resistivity as done in MHD models. Figure 24 shows the results from the MHDEPIC Mercury simulation by Chen et al. (2019c), where a PIC box is embedded in the magnetotail to study tail dynamics and various asymmetries as observed by MESSENGER. While the upstream solar wind is held steady during this simulation, the magnetotail as simulated by the PIC code exhibits very dynamic behavior. A series of plasmoids of varying size form in the tail, and both tailwardmoving and planetwardmoving plasmoids are found in the simulation (Fig. 24b). The modeled plasmoids in the nearMercury tail (X = 2–3 R_{M}) have an average size of 0.3 R_{M} in diameter, which is in accordance with the estimate based on MESSENGER observations in this region (Sun et al., 2020).
Fig. 24 Results from a 3D MHDEPIC simulation of Mercurys magnetosphere. (a) 3D view of the simulated magnetosphere. Colors in the equatorial plane represent plasma density (in amu/cm^{3}), whereas solid lines indicate field lines. The red box shows the embedded PIC region in the magnetotail. (b) A closeup view of the PIC solution showing contours of electron pressure (in Log_{10} nPa) and field lines in the noonmidnight meridian (from Chen et al., 2019c). 
While the imposed solar wind flow is symmetric about the Sun–Mercury line, due to kinetic effects significant dawndusk asymmetries develop in the magnetotail in the MHDEPIC simulation, such as a thicker tail current sheet and higher plasma density and pressure on the dawn side. The occurrence and onset location of tail reconnection in the MHDEPIC simulation also exhibit a strong preference for the dawn side, such that almost all dipolarization fronts and highspeed plasma flows arising from tail reconnection concentrate in the dawn sector. These simulation predictions are remarkably consistent with MESSENGER observations (e.g., Sun et al., 2016; Poh et al., 2017). The Mercury simulation presented here along with other published applications (e.g., Toth et al., 2016; Chen et al., 2017; Zhou et al., 2019, 2020) have demonstrated that MHDEPIC is capable of capturing both local and global physics of the magnetosphere, and, therefore, provides an excellent tool for studying solar windmagnetosphere coupling in planetary magnetospheres and for interpreting observations from spacecraft missions.
8 Current and future directions
8.1 Machine learning
The emergence of computational space physics at the turn of the 21st century was made possible by a collaboration between space physicists, applied mathematicians, computer and computational scientists. At the beginning of the 2020s, we are again witnessing a scientific revolution, similarly to the one a quarter century ago. A new scientific discipline is emerging to offer unprecedented opportunities for advancing many research fields, including space weather: Artificial Intelligence (AI), including Machine Learning (ML) and Computer Vision (CV), can help computers “learn” how to find a needle in the haystack, and help us identify new connections between seemingly unrelated phenomena. Moreover, AI can greatly improve the assimilation of observations into computational models and to quantify the uncertainty of model results. The new challenge is how to integrate this new methodology with our existing space weather modeling capabilities.
A number of userfriendly and free libraries make it possible to apply ML tools to a large variety of problems. Software such as TensorFlow from Google (Abadi et al., 2015), AWS from Amazon (Herrero et al., 2011), Azure from Microsoft (Dudley & Duchene, 2010), or PyTorch from Facebook’s AI Research Lab (Paszke et al., 2019), have allowed a proliferation of ML applications to space physics problems. A couple of years ago the SWMF team reached out to data science, machine learning and computer vision experts at the University of Michigan and we forged a new collaboration to bring advanced data science methods to space weather modeling. In this section we show two specific examples that demonstrate the potential of these new approaches.
8.1.1 Total Electron Content (TEC) maps
The ionospheric total electron content (TEC) is arguably the most utilized physical parameter in ionospheric research in the GNSS era. The TEC maps provide us information about the ionospheric density structures and their evolution, and are also of practical importance since they can be used to estimate the GNSS signal delay due to the ionospheric plasma content between a receiver and a GNSS satellite. Recently, we applied advanced machine learning methods to the forecast of global ionospheric total electron content (TEC) maps (GIM) using maps from one of the International GNSS Services (IGS) centers, i.e., the Center for Orbit Determination in Europe (CODE). Spherical harmonic (SH) fitting is often used in constructing the GIM map. We applied an LSTM neural network method Hochreiter & Schmidhuber, 1997) to forecast the 256 SH coefficients, which are then used to construct the GIM maps (Liu et al., 2020). The model results show that the first/second hour TEC root mean square error (RMSE) is 1.27/2.20 TECU during storm time and 0.86/1.51 TECU during quiet time. Comparing with the CODE GIMS, the RMSE of the LSTM prediction is 1.06/1.84 TECU for the 1st /2nd hour, while the RMSE errors from the IRI2016 and NeQuick2 models are around 9.21/5.5 TECU, respectively (Liu et al., 2020). Moreover, typical largescale ionospheric structures, such as equatorial ionization anomaly and storm enhanced density are well reproduced in the predicted TEC maps during storm time. The ML model performs well in predicting global TEC when compared to two empirical models (IRI2016 and NeQuick2, see Fig. 25).
Fig. 25 Global TEC maps with 6 h interval under storm (13 October) conditions in 2016. The predicted TEC maps are from our LSTM/NN and corresponding ones from the CODE GIM, along with their differences, are, respectively, given in the left, middle, and right panels of each figure. The predictions shown are 1hour ahead. The unit of the color contour is TECU (10^{1} 6electrons/m^{2}) (from Liu et al., 2020). 
The IGS GIM maps are constructed based on a few hundred IGS stations with limited spatial resolution. Critical mesoscale ionospheric structures that cause the most severe GNSS scintillations, i.e., the equatorial plasma bubbles, are smoothed out during the SH fitting procedure. Therefore, we also developed an innovative GIM construction method called VISTA (Video Imputation with SoftImpute, Temporal smoothing and Auxiliary data) (Sun et al., 2021). Several extensions of existing matrix completion algorithms have been utilized to achieve TEC map reconstruction, accounting for spatial smoothness and temporal consistency while preserving important multiscale structures of the TEC maps. This method utilizes the Madrigal TEC data based on over 5000 GNSS receivers and is able to overcome large missing data over the oceans. This newly proposed algorithm is targeted but not restricted to the temporal TEC map reconstruction. In the future we will use VISTA to process all historical Madrigal TEC data and then use the fully reconstructed maps to perform TEC forecast and data assimilation.
8.1.2 Neural network predictions of solar flares
Forecasting large solar flares with machine learning (ML) is at the heart of space weather prediction: Flare radiation and any information about the occurrence of a solar eruption event is carried at the speed of light, hence true forecasting is required. The first significant increase of solar energetic particle fluxes can take place within an hour after the flare. The US national space weather forecasting goal is to provide physicsbased hourly space weather forecasts for validity periods of up to 120 h. This goal can only be achieved if we can forecast solar flares and CMEs.
Our approach is to break down the flare forecasting problem into a series of increasingly challenging ML/computer vision steps. First, we applied cuttingedge classification methodology to obtain a flare probability index that jumps from small values (<0.3) to near unity about one day before a large (M/X class) solar flare takes place (Chen et al., 2019b) (see the yellow curves in Fig. 26). The ML algorithm for the flare probability index utilizes SDO/HMIbased SHARP parameters, physicallyinsightful summaries of active region photospheric magnetic fields. Thus, while this MLbased approach adds no new physics, it represents SOLSTICE’s first step toward improving event forecasts with crossdisciplinary efforts.
Fig. 26 The ML strong flare (M/X class) probability index exhibits a dramatic increase about a day before the strong flare occurs on the Sun (Chen et al., 2019b; Wang et al., 2020). The ML methodology also forecasts the GOES Xray peak intensity of the first strong flare (Jiao et al., 2020). 
Second, we developed an innovative mixed LSTM regression model, which combines binary classification of flaring and LSTM regression for flare intensity, that was used to predict the flare onset jointly with the maximum solar flare intensity observed by the GOES satellites within a 24h time window (Jiao et al., 2020). The predicted intensity peaks are shown by light blue curves in Figure 26. While the results are quite encouraging, the method needs to be trained on larger data sets and improved in terms of computational efficiency and model interpretability.
8.2 Data digestion and assimilation, ensemble modeling and uncertainty quantification
Physicsbased modeling, while very successful in many applications, has several inherent limitations. Often the initial and boundary conditions are not fully known, so the equations to solve are not fully specified. The computational cost of the model may become prohibitive for high fidelity and/or wellresolved models. Here we describe a few of the challenges and future directions.
Solar simulations are driven by the boundary conditions applied near the surface of the Sun. While we have reasonable estimates of density, temperature and radial magnetic field at a large part of the surface, the other plasma quantities are not well understood. The lineofsight velocity, for example, can be measured with Doppler shifts, but the other components are not known. The radial magnetic field on the back side of the Sun and near the poles is not well observed in general. The transverse components of the magnetic field can be obtained from vector magnetograms, however, there is an ambiguity for the sign of these components, and the observational errors in the transverse components are very significant. Using the “observed” magnetic field vector as a direct boundary condition for the AWSoM results in substantial unphysical flows due to the errors in the magnetic field. We are currently working on improving the data ingestion algorithm for the vector magnetogram information so that we will, hopefully, be able to initiate CMEs directly from the evolving magnetogram instead of inserting flux ropes, which is our current method.
Another way of using data in physicsbased models is data assimilation. Data assimilation has the potential to significantly improve model performance, as it has been successfully done in terrestrial weather forecasting. To allow for the sparsity of observations of the Sun–Earth system, however, a different data assimilation method needs to be employed than the typical ensemble Kalman filter used in terrestrial forecasting. Presently we are developing a model that combines physics based modeling with data assimilation and uncertainty quantification (UQ). The model will start from the Sun with an ensemble of simulations that span the uncertain observational and model parameters based on a comprehensive UQ analysis. At the end the model will provide a probabilistic forecast of the space weather impacts. While the concept is simple, finding the optimal algorithm that produces the best prediction with minimal uncertainty is a complex and very challenging task that requires developing, implementing and perfecting novel data assimilation and uncertainty quantification methods.
8.3 Opensource Development
Earlier in this paper we described how a large interdisciplinary group of researchers have developed, with sustained effort, the firstprinciples SWMF that is capable of modeling and forecasting space weather and other space physics phenomena. To make the SWMF impactful, it needs to be used by the space physics community.
From the beginning, We have made the SWMF (Toth et al., 2005, 2012) available to the whole community via a user license. Users can obtain the full source code at http://csem.engin.umich.edu/tools/swmf with all scripts and documentation and use it for their research with minimal restrictions. In addition, BATSRUS and the SWMF have been available for runsonrequest through the Community Coordinated Modeling Center (CCMC) at NASA Goddard through https://ccmc.gsfc.nasa.gov/index.php. CCMC made the SWMF, and many other models, accessible to a wide user community who may not have access to large computer resources and/or are not expert users who can configure and run a complex model.
Some parts of the SWMF and related software have been truly opensource for a while. The Global Ionosphere Thermosphere Model (GITM) has been opensource at https://github.com/~aaronjridley/ since 2012. The Space Science library for Python (spacepy), available from https://github.com/spacepy/spacepy since 2012, has become one of the best free visualization and analysis tools for the SWMF output. More recently, VisAnaMatlab at https://github.com/henry2004y/VisAnaMatlab and VisAnaJulia at https://github.com/henry2004y/VisAnaJulia, visualization and analysis tools written by Hongyang Zhou for the Matlab and Julia languages, respectively, have been made opensource too.
Finally, the core SWMF model was also released in 2020 under a noncommercial opensource license at https://github.com/mstemquda. MSTEMQUDA contains the full core of the SWMF and the BATSRUS, RIM, RCM, and RBE models. In addition, it contains a new Python library, swmfpy. We expect to add the CIMI inner magnetosphere model in the near future. The MSTEMQUDA repository is an uptodate mirror of the repositories developed at Michigan. Both the master and stable branches are available.
Making a major part of the SWMF truly opensource opens a new era in the use and development of space physics and space weather model development. We are hopeful that it will lead to more use, faster and more reliable model development and productive collaboration in the community.
9 Concluding remarks
Over the last decades most scientific disciplines have undergone a major revolution, and the science behind space weather is no exception. A few decades ago, observations and theory were the two pillars of scientific discovery. Since then, the explosive advancements in computer hardware, software, numerical algorithms and data assimilation methods have made computational space physics a third pillar of space weather science.
The emergence of computational space physics was made possible by a close collaboration between space physicists, applied mathematicians, computer and computational scientists. But the formation of tightly integrated research efforts did not happen overnight: It takes years to educate researchers from diverse disciplines to understand each other’s terminology, basic concepts, and methodology well enough to create a breakthrough product.
The SWMF development was started in the 1990s. Over the last three decades funding agencies and the University of Michigan invested over $50M and about 200 personyear effort in this project. Maintaining and developing a worldclass modeling framework requires collaboration of space scientists, mathematicians, numerical and computer scientists, and the space weather user community. Such large research environments take time to build and require continued investments both in intellectual and computational capabilities.
In this paper we described the present state of the SWMF and its main component, BATSRUS. We also outlined its history and future directions. Today, SWMF, BATSRUS and all their simulation and analysis tools constitute a cutting edge capability that is available to the space weather community.
Acknowledgments
This work was supported by NASA DRIVE Science Center grant 80NSSC20K0600 and NASA MMS grant 80NSSC19K0564, NASA LWS grants 80NSSC20K0185, 80NSSC20K0190, 80NSSC20K1778 and 80NSSC17K0681, NASA SSW grant 80NSSC20K0854, NASA HSR 80NSSC20K1313, NASA 80NSSC21K0047, the NSF PREEVENTS grant 1663800 and NSF SWQU grant PHY2027555. The authors thank Drs. John Dorelli and Natalia Buzulukova for providing their unpublished result to be shown in the paper (Fig. 15). The authors also thank Ms. Deborah Eddy for making the paper more readable.
Appendix A: Fundamentals of BATSRUS
The BATSRUS code was originally developed in the mid 1990s when there was a major national initiative to utilize the new transition from vector machines to massively parallel architectures. There were three principles guiding this development: (i) apply the latest advances in computational fluid dynamics to MHD, (ii) utilize the emerging adaptive mesh refinement (AMR) technology and (iii) create a data structure that is truly scalable to a very large number of CPU cores.
The emergence of computational space physics at the turn of the 21st century was made possible by a close collaboration between space physicists, applied mathematicians, computer and computational scientists. But the formation of tightly integrated research efforts did not happen overnight: It takes years to educate researchers from diverse disciplines to understand each other’s terminology, basic concepts, and methodology well enough to create a breakthrough product.
A.1 8Wave Riemann Solver
The first step of the BATSRUS development was to attack a fundamental roadblock to highresolution magnetohydrodynamics (MHD) simulations. Highresolution schemes are based on the conservative form of the governing equations:
where is the vector of conserved quantities defined by
where ρ is mass density, u_{x}, u_{y} and u_{z} are the three components of the plasma bulk flow velocity vector, , while is the magnetic field vector and ε is the total energy density
Here ε_{hd} denotes the hydrodynamic energy density, while γ is the specific heat ratio and μ_{0} is the permeability of vacuum. The flux tensor, , can be written as
Finally, is a “source” vector, containing the terms that cannot be expressed in divergence form:
The “source term” given by equation (A.5) can be handled two different ways. One can directly apply Maxwell’s equation that expresses the absence of monopoles resulting in a identity. However, setting to zero results in a degenerate eigensystem for equation (A.1) (cf. Roe & Balsara, 1996). Due to this degeneration even advanced MHD codes solve only the hydrodynamic part of the MHD equations with highresolution methods and advance the magnetic field separately (cf. Lyon et al., 2004; Clarke, 2010; Zhang et al., 2019a). In a groundbreaking paper, Powell (1997) proposed a Riemann solver that formally keeps the term in equation (A.5) and makes it a passively convected quantity. This method resolves the degeneracy of the eigensystem of equation (A.1) and results in an 8th wave that carries information about the discontinuity in the normal component of the magnetic field. This socalled 8wave scheme ensures the solenoidity condition to truncation accuracy (Powell, 1997) and it makes it possible to formulate the MHD problem in a way that makes it suitable for highresolution schemes. Toth (2000) carefully evaluated the various methods that constrain and concluded that the 8wave approach performs as well as the alternative methods generally applied by the computational MHD community.
The complete numerical algorithm solving the MHD equations with the 8wave scheme was published by Powell et al. (1999). It gives a detailed description of the full eigensystem of the MHD equations together with a space physics example. The publication of the Powell et al. (1999) paper created an avalanche of negative reactions from the space physics MHD modeling community, once again proving George Barnard Shaw’s sarcastic comment: “All great truths begin as blasphemies.” The criticism culminated in a paper by Raeder (2000) that tried to discredit the 8wave scheme. The subsequent comment and reply exchange (Gombosi et al., 2000; Raeder, 2000) stopped the open criticism, but the underlying skepticism from some competitors still lingers even today.
A.2 Adaptive mesh refinement
BATSRUS uses a simple and effective blockbased adaptive mesh refinement (AMR) technique (Stout et al., 1997). The approach closely follows that first developed for twodimensional gas dynamics calculations by Berger & Jameson (1985); Berger & Colella (1989). This blockbased tree data structure is advantageous for several reasons. One of the primary advantages is the ease with which the grid can be adapted. If, at some point in the calculation, a particular region of the flow is deemed to be sufficiently interesting, better resolution of that region can be attained by refining a block, and inserting the eight finer blocks that result from this refinement into the data structure. Removing refinement in a region is equally easy. Decisions as to where to refine and coarsen are made based on either geometric considerations or on comparison of local flow quantities to threshold values.
The governing equations are integrated to obtain volumeaveraged solution quantities within computational cells. The computational cells are embedded in regular structured blocks of equalsized cells. The blocks are geometrically selfsimilar. Solution data associated with each block are stored in standard indexed array data structures, making it straightforward to obtain solution information from neighboring cells within a block. Note that the data on each block can be associated with any one of a multitude of coordinate systems including Cartesian, curvilinear, and more.
Computational grids are composed of many selfsimilar blocks. Although each block within a grid has the same datastorage requirements, blocks can be of different sizes in terms of the volume of physical space they occupy. Starting with an initial mesh consisting of blocks of equal size (that is, uniform resolution), spatial adaptation is performed by dividing and coarsening appropriate solution blocks. In regions requiring increased cell resolution, a parent block is refined by dividing itself into eight children, or offspring. Each of the eight octants of a parent block becomes a new block with the same number of cells as the parent, which doubles cell resolution in the region of interest. Conversely, in overresolved regions, the refinement process reverses; eight children coarsen and coalesce into a single parent block. Thus, cell resolution reduces by a factor of 2. Multigridtype restriction and prolongation operators are used to evaluate the solution on all blocks created by the coarsening and division processes, respectively.
When a 3D block is refined, it is split into eight octants (see Fig. A.1). Each octant forms a block with the same number of cells as the original block, but the resolution is increased by a factor of two. The resulting grid structure is an octree of blocks, and the equations are solved at the finest level only, i.e. on the leaves of the tree (see Fig. A.2).
Fig. A.1 Selfsimilar blocks illustrating the double layer of ghost cells for both coarse and fine blocks (from Gombosi et al., 2003). 
Fig. A.2 Solution blocks of the BATSRUS computational mesh with three refinement levels (from Gombosi et al., 2003). 
The hierarchical data structure and selfsimilar blocks simplify domain decomposition and enable good loadbalancing, a crucial element for truly scalable computing. For explicit time stepping (all blocks use the same timestep) natural loadbalancing occurs by distributing the blocks equally among the processors. For more complicated timestepping schemes the loadbalancing is more challenging, but it still can be done on a block by block basis. We achieve additional optimization by ordering the blocks using the PeanoHilbert or Morton spacefilling curves to minimize interprocessor communication. The selfsimilar nature of the solution blocks also means that serial performance enhancements apply to all blocks and that finegrained algorithm parallelization is possible. The algorithm’s parallel implementation is so pervasive that even the grid adaptation performs in parallel.
A.3 BATSRUS Performance
For most computational models that involve the solution of partial differential equations (PDEs), domain decomposition (i.e., partitioning the problem by dividing the computational domain into subdomains, and farming the subdomains off onto separate cores) is a natural and, in many cases, the most practical approach to parallelization. The blockbased AMR solver was designed from the ground up with a view to achieving very high performance on massively parallel architectures (Stout et al., 1997). The hierarchical data structure and selfsimilar blocks make domain decomposition of the problem almost trivial and readily enable good loadbalancing, a crucial element for truly scalable computing. A natural load balancing is accomplished by simply distributing the blocks equally amongst the processors. The parallel implementation of the algorithm has been carried out to such an extent that even the grid adaptation is performed in parallel.
Figure A.3 shows the weak scaling (how the solution time varies with the number of processors for a fixed problem size per processor) of BATSRUS on several supercomputers from 8 up to more than 200,000 cores. Recently, BATSRUS has been further developed to use hybrid MPI and OpenMP parallelism that allows scaling to beyond 500,000 cores (Zhou & Tóth, 2020).
Fig. A.3 The cell update rates as a function of number of cores for the BATSRUS model. The problem size scales in proportion to the number of parallel processes. The dotted lines represent linear scaling. 
Appendix B: SWMF
The Space Weather Modeling Framework (SWMF) was developed in the early 2000s (Toth et al., 2005) with a combination of support from the DoD MURI (Multidisciplinary University Research Initiatives) and NASA Earth and Space Sciences HPCC (High Performance Computing and Communications) programs. This development closely followed the development path of BATSRUS: its intellectual leadership came from a tightly integrated team of senior university faculty in computer science, software engineering, applied mathematics and space plasma physics, while the actual development was carried out by a group of earlytomid career scientists with help from postdocs and graduate students. The availability of significant stable funding for over a decade was a critical element of the success of the SWMF development (the funding came just as the BATSRUS development resources were winding down).
B.1 Structure
Figure B.1 shows the structure of the SWMF. There are over a dozen components or physics domains. In an actual simulation one can use any meaningful subset of the components.
Fig. B.1 The original physics modules of the SWMF (Toth et al., 2005). This was the first successful coupling (De Zeeuw et al., 2004) of a gyrokinetic ring current model (Harel et al., 1981; Wolf et al., 1991; Sazykin et al., 2002; Toffoletto et al., 2003) to a global MHD model describing the magnetosphere. 
If the simulation starts from the Sun, it is typically driven by solar magnetogram data and flare/CME observations. Simulations restricted to magnetospheric components are usually driven by the solar wind data obtained by satellites upstream of the Earth, for example ACE, Wind or Geotail. We also use the F10.7 solar flux for some of the empirical relationships in the ionosphere and thermosphere models.
The SWMF has a layered architecture (see Fig. B.2). The top layer is the user interface. The second layer contains the control module, which is responsible for distributing the active components over the parallel machine, executing the models, and coupling them at the specified frequencies. The third layer contains the physics domain components. Each component can have multiple physics models. Each component version consists of a physics model with a wrapper and one or more couplers. The wrapper is an interface with the control module, while each coupler is an interface with another component. The physics models can also be compiled into standalone executables. The fourth and lowest layer contains the shared library and the utilities that can be used by the physics models as well as by the SWMF core.
Fig. B.2 The layered architecture of the SWMF (from Toth et al., 2012). 
B.2 Couplers
The SWMF couples together the various models at regular intervals, based on either simulation time or iteration number. The relevant physical quantities are passed with efficient MPI communication. In addition to transferring the data, SWMF has to transform between coordinate systems, take care of unit conversions, and interpolate between different grids. Often the models are moving or rotating relative to each other so that the mapping has to be recalculated every coupling time. A further complication arises for adaptive grids that may change between two couplings. SWMF includes utilities to take care of coordinate transformations and interpolation between various grids.
Since the models use widely different grids and time steps, coupling through a simple interface may be very challenging, especially when the flow is slower than the fast magnetosonic speed. A possible solution is to overlap the models. For example the inner boundary of the inner heliosphere model is provided by the solar corona model at 20 R_{⊙}, while the solar corona obtains its outer boundary conditions from inner heliosphere module at 24 R_{⊙}. The overlap serves as a buffer to suppress numerical artifacts due to the differences between the spatial and temporal resolutions.
In some cases the coupling between the physics models requires some complicated and expensive calculations. For example the inner magnetosphere and the radiation belt models require passing the magnetic field geometry and the plasma state along the closed magnetic field lines of the global magnetosphere model. Tracing magnetic field lines is challenging because the global magnetosphere grid is large and distributed over many processors. SWMF uses highly parallel and efficient schemes for tracing multiple field lines (De Zeeuw et al., 2004; Glocer et al., 2009a) that provide mapping information, integrate quantities along the lines, or extract state variables and positions along the lines.
B.3 Original SWMF modules
Figure B.1 shows the components of the original SWMF and the models that can represent these components. Several components were represented by the BATSRUS code. Since the SWMF is compiled into a single library, the components cannot contain modules, external subroutines or functions with identical names. An automated script ensured that BATSRUS codes representing various components could be compiled together and they could be configured and run with different parameters. The original models of SWMF were the following:
B.3.1 Solar Corona (SC)
The Solar Corona (SC) (represented by BATSRUS) domain started at the photosphere and extended to a few solar R_{⊙}. The MHD equations were solved with empirical heating functions, heat conduction, and radiative cooling on a corotating spherical grid with highly stretched radial coordinates to capture the transition region (Cohen et al., 2007; Downs et al., 2010).
B.3.2 Eruptive Events (EE)
The Eruptive Event generator component is responsible for creating a CME. This was achieved with empirical models that insert an unstable flux rope into the steady solar corona solution, or insert an arcade and apply shearing motion at the lower boundary of the corona model (Roussev et al., 2003; Manchester et al., 2004b).
B.3.3 Inner Heliosphere (IH)
The Inner Heliosphere model originally extended from about 20 R_{⊙} to the orbit of the Earth and has been later extended to include the planets. BATSRUS solves the ideal or twotemperature MHD equations on a Cartesian grid in either corotating or inertial frame, and it can model the propagation of CMEs from the Sun to the Earth (Lugaz et al., 2005b; Toth et al., 2005; Manchester et al., 2006).
B.3.4 Global Magnetosphere (GM)
The Global Magnetosphere domain surrounds the Earth and it extends about 30 R_{E} toward the Sun, a few hundred R_{E} toward the magnetotail, and about 60 R_{E} in the orthogonal directions. BATSRUS solves the MHD equations on a Cartesian or spherical grid. As an alternative, the Tsyganenko (1989) empirical model can provide the magnetic field as a function of observed solar wind parameters and planetary indexes.
B.3.5 Inner Magnetosphere (IM)
The Inner Magnetosphere model consists of the closed magnetic field line region around the Earth. The Rice Convection Model (RCM) (Harel et al., 1981; Wolf et al., 1991; Sazykin et al., 2002; Toffoletto et al., 2003) solves for the bounce averaged and isotropic but energy resolved particle distribution of electrons and various ions. This was the first successful coupling (De Zeeuw et al., 2004), of a gyrokinetic ring current model to a global MHD model describing the magnetosphere.
B.3.6 Radiation Belts (RB)
The Radiation Belt domain overlaps with IM but it models the relativistic electrons. The RBE (Fok et al., 2008) model solves the bounceaveraged Boltzmann equation.
B.3.7 Ionospheric Electrodynamics (IE)
The Ionospheric Electrodynamics model is a two dimensional heightintegrated spherical surface at a nominal ionospheric altitude (at around 110 km for the Earth). The Ridley Ionosphere Model (RIM) (Ridley et al., 2004) code uses the fieldaligned currents to calculate particle precipitation and conductances based on empirical relationships, and then it solves for the electric potential on a 2D spherical grid.
B.3.8 Upper Atmosphere (UA)
The Upper Atmosphere contains the thermosphere and the ionosphere extending from around 90 km to about 600 km altitude for the Earth. The GITM (Ridley et al., 2006) code solves the equations of multispecies hydrodynamics with viscosity, thermal conduction, chemical reactions, ionneutral friction, source terms due to solar radiation, etc. on a spherical grid in a corotating frame. The MSIS (Hedin, 1991) and IRI Bilitza (2001) empirical models provide statistical average states for the upper atmosphere and ionosphere, respectively.
B.4 Additional Simulation Tools
In addition to BATSRUS, the SWMF includes several worldclass models that provide one of the most advanced space weather simulation capabilities ranging from kinetic and mesoscales to global description of the space environment. Here we briefly summarize the most important simulation/postprocessing tools included in the SWMF suite.
B.4.1 EEGGL
The Eruptive Event Generator using the GibsonLow configuration tool (EEGGL, Jin et al., 2017a) is the first community model (available at the CCMC) to simulate magnetically driven CMEs. It is an automated tool for finding the flux rope parameters Gibson & Low (1998) to reproduce observed CME events (e.g., Jin et al., 2017a; Borovikov et al., 2017b). The solar magnetogram is first used to specify the inner boundary condition of the magnetic field for AWSoM(R), which is then employed to generate an ambient solar wind solution. Simultaneously, the input magnetogram and the observed CME speed are used by EEGGL to determine the Gibson & Low (1998) flux rope parameters. With the derived parameters, a Gibson & Low (1998) flux rope is inserted into the ambient solar wind to initiate the CME event. The various parameters of the Gibson & Low (1998) model are carefully selected to reproduce the observed CME source region and speed. The user can change parameters, such as helicity and initial orientation, to experiment with the properties of the resulting eruption. Presently, we are working on more eruptive event generator tools, which employ different physical processes (such as the Titov & Démoulin, 1999 mechanism) to initiate solar eruptions.
B.4.2 MFLAMPA
The Multiple Field Line Advection Model for Particle Acceleration solves the kinetic equation for solar energetic particles along a multitude of interplanetary magnetic field lines originating from the Sun (Sokolov et al., 2004; Borovikov et al., 2018). It is seamlessly coupled to AWSoMR and EEGGL and can therefore account for the temporal evolution of field lines as the CME moves outward from the Sun. The diffusion coefficient used in MFLAMPA is selfconsistently calculated from the energy densities of the Alfvénic turbulence simulated by AWSoMR. Together, MFLAMPA, AWSoMR and EEGGL provide a highperformance, selfconsistent description of the solar corona, inner heliosphere, and solar energetic particle distribution to study solar storms and their impact on the inner heliosphere.
B.4.3 AMPS
The Adaptive Mesh Particle Simulator (AMPS) is a highperformance kinetic Monte Carlo code originally developed for modeling neutral planetary environments, where it was used to solve the Boltzmann equation accounting for particle collisions, internal degrees of freedom, and chemical reactions (Tenishev et al., 2008, 2021). AMPS employs AMR mesh with cutcells for discretizing the simulated domain. The implemented cutcells methodology allows the code to simulate gas flows around objects with arbitrarily complex surface geometry like, e.g., nuclei of comets (Tenishev et al., 2016).
The distinct feature of AMPS is the ability to model twophase environments, where a dust phase is simulated concurrently with the ambient gas or plasma. In such a simulation, the electric charge of the dust grains varies according to the ambient plasma conditions, and the size and/or chemical composition of a dust grain can change affected by, e.g., sublimation of volatiles carried the grain (Tenishev et al., 2011).
AMPS was extended to simulating energetic charged particle transport in the inner heliosphere and the Earths magnetosphere (Tenishev et al., 2005, 2018). This model can be used to describe a broad range of suprathermal particle populations including magnetospheric particles with energies exceeding ~1 keV/nucleon, solar energetic particles in the MeV to GeV range, or galactic cosmic rays with energies above ~ 100 MeV/nucleon. Recently, AMPS was extended by adding an implicit PIC capability. Now it can be used for simulating various plasma phenomena either as a standalone modeling tool or coupled to other components of the SWMF.
AMPS is coupled to several components of the SWMF, allowing multiscale and multiphysics simulations. Specifically, twoway coupling has been developed with the Global Magnetosphere (GM) and Outer Heliosphere (OH) modules. A oneway coupling procedure is implemented to couple AMPS to the Solar Corona (SC) and Inner Heliosphere (IH) components of the SWMF for SEP simulations.
B.4.4 iPiC3D
iPiC3D is a parallel highperformance implicit ParticleinCell (PIC) code (cf., Markidis et al., 2010). It solves the full set of Maxwell’s equations for the electromagnetic fields coupled with the equations of motion for electrons and ions on 3D Cartesian grids. The discretization is based on the implicit moment PiC (IMPiC) method that employs an implicit time integration for the electric field, then the magnetic field is updated from the induction equation, finally the particles are moved with a simple iterative scheme (Brackbill & Forslund, 1982, 1986; Lapenta et al., 2006; Brackbill & Lapenta, 2008). The main advantage of iPiC3D is that it is capable of taking larger grid cell sizes and timesteps and thus making the coupled simulation affordable on today’s supercomputers. There are still open questions about the use of implicit PiC codes. The bottom line is that if one wants to resolve Debye scale phenomena the use of expensive explicit PiC codes are necessary. However, if one is mainly interested in reconnection and other space plasma phenomena the use of implicit PiC codes is not only justified, but also necessary (cf., Ricci et al., 2002).
B.4.5 FLEKS
The FLexible Exascale Kinetic Simulator (FLEKS) is a new particleincell (PIC) code that is designed for the MHD with adaptively embedded PIC (MHDAEPIC) simulations. FLEKS uses the Gauss’s Law satisfying EnergyConserving SemiImplicit Method (GLECSIM) (Chen & Toth, 2019) as the base PIC solver. Novel particle splitting and merging algorithms have been designed to control the number of macroparticles per cell during a long MHDAEPIC simulation. The particle splitting algorithm improves statistical representation and reduces noise in the cells with low macroparticle number, while the particle merging algorithm alleviates the load imbalance and speeds up simulations. The FLEKS grid is Cartesian, but the active PIC region is not limited to be a box anymore since any Cartesian cells can be switched on or off at any point of the simulation. FLEKS uses the highperformance parallel data structures provided by the AMReX library (Zhang et al., 2019b, 2020) to store the fields and also the particles.
Appendix C: Physics
BATSRUS has a layered modular software architecture to handle several applications with a single base code (see Fig. D.2). The state variables of the equation system are defined by the equation modules, while the rest of the application dependent details are implemented into user modules. A configuration script is used to select the equation and user modules that are compiled together with the code. There are currently dozens of equation and user modules (obviously not all combinations are possible) which means that BATSRUS can be configured for quite a few different applications. In addition to the basic equations, there are various source terms that change from application to application: collisions, charge exchange, chemistry, photoionization, recombination, radiative losses, etc. The boundary and initial conditions vary greatly as well.
C.1 Conservation Laws in BATSRUS
BATSRUS can be configured to solve the governing equations of ideal and resistive MHD (Powell et al., 1999), semirelativistic (Gombosi et al., 2002), anisotropic (Meng et al., 2012), Hall (Toth et al., 2008), multispecies (Ma et al., 2002) and multifluid (Glocer et al., 2009c) extended magnetofluid equations (XMHD) and more recently nonneutral multifluid plasmas (Huang et al., 2019).
Table C.1 summarizes the various extended MHD conservation laws that can be solved by BATSRUS.
C.1.1 Extended MHD Equations
BATSRUS can solve many approximations to the loworder velocity moments of the Boltzmann equations (we refer the interested readers to the literature (cf., Burgers, 1969; Schunk & Nagy, 1980; Gombosi & Rasmussen, 1991; Gombosi, 1998; Shumlak & Loverich, 2003; Huang et al., 2019). The governing equations for species “s” can be written as
where ρ and denote the mass density and the velocity vector, respectively, and q and m are the charges and masses of the particles. For the pressure tensor we used the CGL approximation (Chew et al., 1956): , where is the identity matrix, is the unit vector along the magnetic field direction, is the pressure along the parallel direction of the magnetic field and p_{⊥} is the pressure in the perpendicular direction. The scalar pressure can be written as p = (p_{∥} +2p_{⊥})/3. BATSRUS has the capability to solve the full equation system (C.1) or reduce it and only solve for the scalar pressure, p.
Equation (C.1) can be obtained from the Boltzmann equation by considering the infinite series of velocity moments, called Maxwell’s equation of change (cf., Gombosi, 1994):
Here is the velocity distribution function of species “s” (expressed in terms of the random velocity, ), is a physical quantity of a single particle of species “s” dependent on the random velocity. 〈 〉 denotes averaging over the entire random velocity space. The order of in the random velocity defines the order of the velocity moment equation. For instance, (zerothorder moment equation) results in the continuity equation, describing the conservation of mass. The firstorder velocity moment equations are obtained by using and they express the conservation of momentum. The secondorder velocity moment equations are obtained by using . There is one zeroorder, three firstorder and six secondorder moment equations (due to the symmetric nature of the diad).
It is important to note that equation (C.2) leads to an infinite number of velocity moment equations. The “villain” is the third term on the left hand side of equation (C.2), . If is nth order in velocity, the term is the (n + 1)th velocity moment. In other words, the transport equation for the nth velocity moment contains the divergence of the (n + 1)th moment, resulting in an infinite series of partial differential equations.
The infinite series of velocity moment equations must be closed some way to obtain a closed set of differential equations. There are a number of closures in the literature (cf., Chapman, 1916; Enskog, 1917; Grad, 1949; Levermore, 1996). The simplest (and most popular) closures either neglect the thirdorder velocity moments (the heat flow), or express a highorder velocity moment in terms of lower moments (cf., Grad, 1949). Equation (C.1) was obtained by neglecting the heat flow tensor and using the CGL approximation (Chew et al., 1956) for the pressure tensor. In this approximation there are six velocity moments we solve for: ρ_{s}, the three components of and the two pressure components, and . For this reason this is called the six moment approximation.
The electric () and magnetic fields () are obtained from Maxwell’s equations:
where ϵ_{0} is the vacuum permittivity, μ_{0} is the vacuum permeability, is the speed of light, is the total charge density and is the current density. Equations (C.3c) and (C.3d) are constraints on the initial conditions and analytically these conditions are preserved. Numerically, however, this is not guaranteed to hold. BATSRUS uses a variety of methods to enforce the solenoidal magnetic field condition (for more details see Toth et al., 2012).
It is important to point out that in the multifluid formulation the electric current density depends on the charge averaged, and not the mass averaged, ion velocity. This can be seen by looking at the definition of :
where Z_{s} is the ionization state of a given ion species and is the charge averaged ion velocity. Note that in general and the two vectors can be quite different. BATSRUS takes into account the full definition of and thus selfconsistently accounts for the different velocities of the various ion species (Glocer et al., 2009c). This is different from the approximate solution applied in the LFM code (cf. Wiltberger et al., 2010; Merkin, 2011) that assumes that the macroscopic plasma velocity in the direction perpendicular to the magnetic field coincides with the electrical drift velocity and therefore is the same for all ion species.
Extended magnetohydrodynamics (XMHD) makes two fundamentally important assumptions: (i) electrons are assumed to be massless and (ii) charge neutrality is assumed at all scales. These two assumptions lead to the generalized Ohm’s law:
In a singleion plasma the electron velocity is resulting in the motional electric field plus the Hall term. The second term in equation (C.5) is the ambipolar electric field. It is interesting to note that the parallel (field aligned) component of the electric field is
where is the parallel gradient operator. In equation (C.6) the first term describes the parallel ambipolar electric field while the second term represents adiabatic focusing. BATSRUS has the capability to solve various XMHD approximations, from ideal MHD to resistive, Hall, anisotropic pressure, multispecies and multifluid limits. A more detailed description of these capabilities can be found in Toth et al. (2012).
C.1.2 Six Moment Equations
A recent addition to the BATSRUS equation set is the sixmoment approximation (Huang et al., 2019). This approximation solves the full set of equations (C.1) and (C.3) without neglecting the electron mass and assuming charge neutrality. Consequently there is no Ohm’s law to express the electric field and we need to solve the full set of electron fluid equations and the full set of Maxwell’s equations. In this approximation the fastest wave mode is the light wave. Since the speed of light typically well exceeds the typical MHD wave speeds, one can artificially reduce it to allow larger time steps and more efficient computation. An additional benefit is that the whistler wave speed is also limited by this reduced speed of light. The sixmoment equations describe several phenomena that are not captured by simpler MHD equations, Hall physics, relativistic limit of fast and whistler waves, net charge, anisotropy of both electron and ion pressures, etc. An additional benefit is that one can have multiple species with positive and negative charges, including multiple electron fluids or negatively charged dust. See Huang et al. (2019) for details.
C.1.3 Source terms
The collision terms in the transport equations describe the various physical processes that transfer mass, momentum and energy between various ionized or neutral species. These terms represent the underlying physics that enable us to model the interaction of space plasma flows with planets (cf. Ma et al., 2004, 2013, 2018b; Sarkango et al., 2019), planetary moons (Rubin et al., 2015; Jia et al., 2018; Harris et al., 2021, cf.), comets (Gombosi et al., 1996; Huang et al., 2016b, cf.), and other objects of interest. The collision term describes the rate of change of the distribution function due to interaction between various species. In BATSRUS we consider the following processes:

elastic collisions,
photoionization and impact ionization (using the BeerLambert Law),
charge transfer, and
recombination.
Next, we discuss the contributions of these processes to the collision term. We make the following simplifying asssumptions:
All particles are assumed to lack any internal degrees of freedom;
Energy thresholds of various processes (such as chemical reactions, ionization thresholds, etc.) will be neglected;
All neutral species are considered cold (T_{n} = 0) and are assumed to move with the same bulk velocity, .
These simplifications limit the scope of our approximations, but our methodology still provides useful insights into collisional effects in space plasmas.
In the present approximation all particles are assumed to possess no intrinsic degrees of freedom, therefore all inelastic collisions change the identity of a particle. These reactions result in ionization, charge transfer, or recombination.
Elastic collisions. Elastic collisions do not change the identity of particles, but do change the momentum and energy of individual particles. The effects of these collisions is described in the general framework of the relaxationtime approximation (cf. Bhatnagar et al., 1954; Burgers, 1969; Gombosi, 1994). The main idea behind this approximation is the recognition that collisions drive all gas components toward equilibrium. Since equilibrium phasespace distributions are Maxwellians, the cumulative effect of elastic collisions can be formally described by gradually replacing the present distribution function (F_{s}) with the appropriate Maxwellian, F_{s(st)} (cf. Gombosi, 1994):
In expression (C.7) the subscript “t” refers to all species other than “s”, and τ_{st} is a “relaxation time” characterizing how fast the distribution function F_{s} approaches equilibrium due to collisions between particles of types “s” and “t”. Equation (C.7) means that “s” and “s′” collisions may drive particles “s” toward two different equilibria: however, in steadystate equilibrium all species will reach the same bulk velocity and temperature.
The relaxation timescale, τ_{st}, can be different for different species. For instance, electrons relax toward equilibrium faster than ions in ionelectron collisions. In practice, the momentum transfer collision frequency, is used instead of the relaxation time. The momentum transfer collision frequency includes a massdependent factor that accounts for the efficiency of momentum transfer in an elastic collision:
The parameters of the Maxwellian, , are chosen in a way that mass, momentum and energy are conserved while the gas is driven toward equilibrium (Burgers, 1969; Gombosi, 1994):
where
In these expressions, k_{B} is the Boltzmann constant and the kinetic temperature is defined as
Equations (C.10) through (C.13) describe the number density of species “s”, the drift velocity of species “s” with respect to the center mass of fluids “s” and “t”, and the “stagnation temperature” of species “s”. It should be noted that and in general .
Ionization. There are four primary ionization processes to be considered: photoionization, impact ionization by superthermal electrons, impact ionization by energetic ions, and finally impact ionization by energetic neutrals. These ionization processes create new charge, therefore we consider them separately from the charge transfer reactions.
The ionization process converts a particle from the thermal neutral population to one of the charged particle species. Since the neutral gas is assumed to be cold (T_{n} = 0) the net ionization source can be approximated by the following expression:
where is the sum of the photoionization and impact ionization frequencies of species “s′”, is the density of particles producing charged particles of type “s”. Throughout this paper the charge state of particles “s′” is one less than the charge state of particles “s”.
Charge Exchange. Charge exchange transfers an electron from one particle to an other (an example is the accidentally resonant reaction). Although there is a transfer of electrons between two heavy particles, in most cases each particle tends to retain its original kinetic energy. Here we limit our consideration to singly charged ions and we consider the following general charge exchange reaction: S + M^{+} → S^{++}M. The ion, S^{+}, is referred to as species “s”, while particles S are species “s′”. In our approximation the neutral particles form a cold gas, therefore one can write the net rate of change of the phasespace distribution function of particles “s” is the following:
Here and are charge exchange rates. The first term describes the loss of particles “s” due to charge exchange with all neutral species, while the second term describes the creation of new “s” particles by charge exchange with “s′”type particles.
Recombination. Recombination removes a positive and a negative charge from the system. It represents a sink for electrons and for particles “s” and a source for particles “s′”. This leads to the following loss rate for ions “s”:
where α_{s} is the recombination coefficient and n_{e} is the electron density. Equation (C.16) also gives the source term for species “s′” (naturally with positive sign).
Combined Collision Term. Next, we combine the collision terms for all processes discussed above and combine equations (C.7), (C.14), (C.15) and (C.16) to obtain:
In BATSRUS we take the appropriate velocity moments of equation (C.17) (corresponding to the actual approximation used for the governing equations).
C.2 Coupled MHD turbulence
The ad hoc elements can be eliminated from the solar corona model by assuming that the coronal plasma is heated by the dissipation of Alfvén wave turbulence (cf. Sokolov et al., 2013). The dissipation itself is caused by the nonlinear interaction between oppositely propagating waves (e.g., Hollweg, 1986).
Within coronal holes, there are no closed magnetic field lines, hence, there are no oppositely propagating waves. Instead, a weak reflection of the outward propagating waves locally generates sunward propagating waves as quantified by van der Holst et al. (2014). The small power in these locally generated (and almost immediately dissipated) inward propagating waves leads to a reduced turbulence dissipation rate in coronal holes, naturally resulting in the bimodal solar wind structure. Another consequence is that coronal holes look like cold black spots in the EUV and Xray images, while closed field regions are hot and bright. Active regions, where the wave reflection is particularly strong, are the brightest in this model (see Oran et al., 2013; Sokolov et al., 2013; van der Holst et al., 2014).
As has been shown by Jacques (1977), the Alfvén waves exert an isotropic pressure on the plasma. The relation between the wave pressure and wave energy density is p_{A} = (w_{+} + w_{−})/2, where w_{±} are the energy densities for the turbulent waves propagating along the magnetic field vector (w_{+}) or in the opposite direction (w_{−}). The Wentzel (1926), Kramers (1926), Brillouin (1926) approximation (WKB) is used to derive the equations that govern the transport of Alfvén waves, which may be reformulated in terms of the wave energy densities. Dissipation of Alfvén waves is the physical process that drives the solar wind and heats the coronal plasma.
Alfvén wave dissipation occurs when two counterpropagating waves interact. Alfvén wave reflection from steep density gradients is the physical process that results in local wave reflection, thus maintaining a source of both types of waves. In order to describe this wave reflection we go beyond the WKB approximation that assumes that the wavelength is much smaller than spatial scales of the background variations.
The equation describing the propagation, dissipation, and reflection of Alfvén turbulence has been derived by van der Holst et al. (2014):
where V_{A} is the Alfvén velocity, while Γ_{±} and are the reflection coefficient and the dissipation rate, respectively. Finally, with the help of the dissipation rate of Alfvén turbulence one can express the ion and electron heating rates (van der Holst et al., 2014; Gombosi et al., 2018).
In this model there are only two free parameters: (i) the Poynting flux of Alfvén waves leaving the photosphere (), and (ii) the transverse correlation length of Alfvén turbulence (L_{⊥}). Our solar corona model assumes that and (cf. Gombosi et al., 2018).
C.3 Gyrokinetic models
C.3.1 Kinetic PWOM
The original polar wind model that PWOM is based on solved the field aligned gyroptropic transport equations for each ion species as described by Gombosi & Nagy (1989). After PWOM was incorporated at the PW component of the SWMF it expanded to a global polar wind model, but retained it’s fluid nature (Glocer et al., 2009b). Given the importance of kinetic processes to many ionospheric outflow mechanisms beyond the polar wind, multiple steps were taken to include these processes in PWOM.
The initial expansion to kinetic processes came with the inclusion of superthermal electrons whose energy is much greater than the thermal energy of ionospheric electrons (0.3 eV). These superthermal electrons are either photoelectrons generated by the photoionization of the neutral atmosphere, precipitating electrons of magnetospheric origin (auroral electrons/polar rain), or secondary electrons generated by other energetic electrons impacting neutral particles. To encorporate these superthermal electrons into PWOM, we split the electron population into thermal and superthermal components that must statisify charge neutrality and current conservation (Glocer et al., 2012):
Here subscripts e and α indicate the thermal and superthermal electrons, respectively. Once the superthermal electron population is known, the thermal population is determined from the above equations as well as the thermal electron energy equation (not shown), which also includes the energy deposition due to collisions between the thermal and superthermal populations. At that point the ambipolar field is determined and the ion solution can be updated, including the effect of additional sources due to impact ionization. In PWOM we have used three approaches to specifying the superthermal electron population. First, in Glocer et al. (2012) we used the output of a twostream calculation of the photoelectron source together with a collisionless kinetic mapping. Second, in Glocer et al. (2017) we coupled PWOM to the kinetic STET code and thereby obtained the superthermal electron solution by solving the Boltzmann equation presented by Khazanov et al. (1994) as:
where the constant β = 1.7 × 10^{−8} eV^{1/2} cm^{−1} s, the superthermal differential flux is given by ϕ = ϕ(t, E, μ, s), the kinetic energy is E, and cosine of the local pitch angle is provide by μ. We also have s defined as the distance along the field line, the magnetic field is given by B and F is the force associated with the parallel electric field and Q is production rate of superthermal electrons, and 〈S〉 represents the collision operators.
The approach to including superthermal electrons as a true kinetic population as described by Glocer et al. (2017) is the most complete and physical approach, but it can be computationally expensive. Therefore, Glocer et al. (2018) included the option to use the twostream approach (cf., Nagy & Banks, 1970; Solomon et al., 1988; Lummerzheim & Lilensten, 1994; Barakat & Schunk, 2001; Schunk & Nagy, 2009). The twostream approach has energy dependence and transport/collisional effects, but does not include the effects of pitch angle diffusion or trapping. It is however dramatically faster than the fully kinetic approach and therefore represents an acceptable compromise between physical completeness and computational efficiency for many problems. However, specific problems that rely on detailed kinetic electron effects such as trapping will still require the more comprehensive treatment.
The inclusion of superthermal electrons allows PWOM to treat only some of the outflow mechanisms, but many other processes require the inclusion of kinetic ions. Most prominant of these processes are waveparticle interactions, which drive ion acceleration in the cusp and auroral zones. Motivated by this, Glocer et al. (2018) expanded the PWOM code to include kinetic ions based on a gyroaveraged PIC approach at high altitudes while keeping the fluid approach at low altitudes for computational efficiency. In the high altitude PIC region each macro particle in PWOM for a species “i” is advanced by solving the gyroaveraged particle equation of motion given by:
where m_{i} specifies the mass, v_{i} specifies the velocity, t specifies the time, and q_{i} is the ion charge. The external forces are given by the parallel electric field E_{∥}, and gravity. In this equation, is the particle’s first adiabatic invariant specified by
At the interface between the low altitude fluid region and the high altitude PIC region information is exchanged. PWOM uses the fluid solution in the last fluid computational cell to sample particles for the PIC region, and the first computational cell of the PIC region is used to compute moments and set boundary conditions for the fluid domain. Collisions are included in both the fluid and PIC regions with the fluid collisional terms provided by Burger’s fully linear approximation (Burgers, 1969) while in the PIC region collisions are included using the Monte Carlo approach described by Takizuka & Abe (1977) and modified by Nanbu & Yonemura (1998) modified by and to allow for particles with variable statistical weights. Wave particle interactions in the PIC region are implemented using the approach described by Retterer et al. (1987). This method includes the heating by randomly perturbing the perpendicular velocity of the macroparticles with the variance determined by a diffusion coefficient, which depends on the wave power. Formulations of these diffusion coefficients are given by Crew et al. (1990), Barakat & Barghouthi (1994).
C.3.2 Dynamic Global Core Plasma Model
Being the coldest (~1 eV) magnetospheric population within the magnetosphere, the plasmasphere’s evolution is dominated by advection via drift and refilling via ionospheric outflow at mid and lowlatitudes. The Dynamic Global Core Plasma Model (DGCPM, Rasmussen et al., 1993; Ober et al., 1997; Liemohn et al., 2004; Borovsky et al., 2014) captures these dynamics by solving a continuity equation for the total flux tube content, N:
is the horizontal bulk velocity of the cold plasmasphere fluid (set by the local drift). A dipole magnetic field is assumed, electric potential is a required input and typically obtained via an empirical model. and represent the net source and loss of plasma from/into a given flux tube, respectively. On the day side, ionospheric plasma is assumed to fill flux tubes until saturation density, N_{S}, is reached:
Saturation values are a function of radial distance and determined empirically Carpenter & Anderson (1992). The filling time constant, τ_{fill}, has a configurable value but defaults to 6.7 days. The loss term includes simple loss into the ionosphere at either end of the flux tube:
The loss time constant, τ_{loss}, is set to 3 days.
Within the SWMF, couplings to other models provide more realistic electric fields and allow for the exploration of the impact of the plasmasphere on the global magnetosphere. DGCPM can obtain electric field from the IE module, opening up a greater range of empirical and firstprinciplesbased electric fields (Borovsky et al., 2014; Ridley et al., 2014). The density of the plasmasphere can be passed to the GM component following the same algorithm used to couple ring current density and pressure (Glocer et al., 2020). It is coupled individually to the HEIDI drift physics model, discussed below (Liemohn et al., 2004), and used extensively to explore ring currentplasmasphere interactions (Liemohn et al., 2006; Liemohn & Jazowski, 2008; Ridley et al., 2014). When at least three modules are used (GM, IE, and PS), the effect of plasmasphere drainage plumes on dayside reconnection can be explored in a selfconsistent manner.
C.3.3 Rice Convection Model
The Rice Convection Model (RCM, Wolf, 1974; Toffoletto et al., 2003; Sazykin et al., 2002) is a guiding center drift model of a set of assumedisotropic populations, each with a given energy invariant, λ, such that,
where W is the kinetic energy of the particles within the population and is the magnetic flux tube volume per unit magnetic flux. The RCM then solves for the evolution of the flux tube density content, η, for each energy invariant population, by solving the continuity equation,
where η is a function of time and space and is defined by,
Here, is the full electromagnetic drift velocity of the population and is given by,
where q is the charge of the individual particles that constitute the energy invariant population. The electric and magnetic fields must be prescribed via external models or empirical relations. Equation (C.28) is solved via an ionospheric grid where each grid point represents the foot point of a magnetospheric flux tube.
The RCM version integrated into the SWMF is unique in its configuration. Flux tube volume, , is obtained from the GM component, as are initial and boundary conditions for η for each energy invariant population. Electric potential is obtained via the IE solution. The total density and pressure at each RCM grid point is handed to GM where it is treated as a source term to the MHD values along each flux tube, nudging the MHD solution towards the RCM solution. This is done via,
where p_{GM} and p_{RCM} refer to the total fluid pressure as calculated by the respective modules and τ_{couple} is a time constant, typically set to 10 s. If other inner magnetosphere (IM) models are used the p_{RCM} term is replaced by the appropriate pressure obtained from the IM model. Equation (C.31) may also be used to return density values alongside pressure. The coupled version of RCM has no explicit source terms outside of advection through the outer boundary; a simple decay is added with a 10 h efolding rate. It should be noted that standalone RCM versions have evolved different capability sets than the one described here (e.g., Yang et al., 2014; Chen et al., 2019a).
C.3.4 The Ring Current Atmospheres Interaction Model (RAM) family
The original Ring current Atmospheres interaction Model (RAM) is a fully kinetic bounceaveraged drift model of the ring current (Fok et al., 1993; Jordanova et al., 1994; Liemohn et al., 1998). It solves the kinetic equation to yield the bounceaveraged distribution function as a function of azimuth, radial distance, energy and pitch angle for several species,
where angle brackets denote average values over the bounce period, Q_{s} is the distribution function for species s, R_{0} is the radial direction in the equatorial plane, “phi” is azimuth, E is particle’s kinetic energy, and μ is the cosine of the particle’s equatorial pitch angle. Finally, h is defined via,
where a and b are mirror points for a particle of a given μ along magnetic field line B(s) with field strength of B_{m} at the mirror point. The relationship between RAMlike models and the RCM (see Sect. C.3.3) is described by Heinemann & Wolf (2001).
RAM has spawned many “child” codes, three of which are integrated into the SWMF and each with its own unique capabilities. These include The Hot Electron Ion Drift Integrator (HEIDI), the Comprehensive Inner Magnetosphere Model (CIMI), and the RAM with a SelfConsistent Magnetic field (RAMSCB) codes. These models all fall into the same class of bounceaveraged models, but have unique implementations as well as differences in variable choice, grid formulation, and source terms that provide strengths for particular problems. Each is integrated into the SWMF such that fields throughout the ring current domain and plasma conditions about the outer boundary are obtained from GM and IE components; pressure and density values are returned to GM following the approach given by equation (C.31).
HEIDI expands upon the original RAM model in several ways. Its usage has mostly focused on largescale dynamics of the inner magnetospheric pressure and current systems, including tracking all source and loss terms (Liemohn et al., 1999, 2002). A key development was the inclusion of selfconsistent electric field calculations (Liemohn et al., 2004), allowing for the analysis of conductance influences on the ring current (Liemohn et al., 2005), and the eventual inclusion of selfconsistent conductance calculation from electron precipitation (Perlongo et al., 2017). A broad set of geocoronal models is available within HEIDI, allowing for deep investigations of ring current decay (Ilie et al., 2013). Also, HEIDI includes a robust definition of nondipolar drift suitable for an arbitrary magnetic field description (Ilie et al., 2012) and can now account for the effects of the inductive electric field as well, giving a more dynamic picture of ring current development (Liu & Ilie, 2021). It runs within the SWMF with RIM and work is in progress toward full coupling with BATSRUS.
The CIMI model represents the combination of two earlier models, the Comprehensive Ring Current Model (CRCM) and the Radiation Belt Environment (RBE) code, to create a complete model of the plasmasphere, ring current, and radiation belt populations of the inner magnetosphere (Fok et al., 2014). The earlier RBE and CRCM codes were coupled into the SWMF (Glocer et al., 2009c, 2013), and CIMI uses an improved version of these couplers that allows it to couple to BATSRUS configured with single fluid MHD, anisotropic MHD, and multifluid MHD (Glocer et al., 2018, 2020). The CIMI grid is located at the magnetic field ionospheric foot points, which simplifies the calculation of the E × B drift as the ionospheric potential does not need to be mapped to the equator. In addition to using different spatial coordinates, CIMI solves a version of equation (C.32) recast in a different set of velocity space coordinates. Namely, the version of CIMI in SWMF uses μ_{a} and K coordinates that correspond to the first and second invariants of motion. This approach lets CIMI represent the advection portion of the transport in a conservation form which can be treated with standard finite volume methods. Specifically, the advection portion of the code can be written as (Fok et al., 2021):
where λ_{i} and ϕ_{i} are latitude and azimuthal angle of the field line foot point, F_{s} is the bounce averaged distribution function times a Jacobian for a particular μ_{a} and K (see Fok et al., 2021, for details). In addition to advection on the left hand side of equation (C.34), wave diffusion terms are often included on the right hand side, shown here only as S_{i} to represent waveparticle interactions with Chorus, hiss, and EMIC waves, which are critical to modeling local acceleration of radiation belt electrons as well as scattering into the loss cone and subsequent precipitation. Additional components of S_{i} include charge exchange loss, loss cone loss, and the effect of Coulomb collisions. Complete descriptions of these terms are given in Fok et al. (2021) along with new forms of equation (C.34) in different coordinates that will eventually be included in the SWMF version. Finally, it is interesting to note that when CIMI is coupled to BATSRUS, the magnetic and electric fields are naturally selfconsistent. Meng et al. (2013) demonstrated that the magnetic field calculated in BATSRUS using the coupled CIMI pressure is consistent with a force balance. Similarly, the pressure feedback from CIMI drives currents in BATSRUS that contribute the ionospheric potential and hence the convection in CIMI.
RAM with SelfConsistent Magnetic field (RAMSCB) uses an Euler potential representation of the magnetic field to achieve a self consistent magnetic field configuration inside the model’s domain (Zaharia et al., 2006; Jordanova et al., 2010). It has an energy range of approximately 100 eV to 500 keV. Loss terms include charge exchange, Coulomb collisions and atmospheric loss at low altitudes. The RAM model was updated to use nondipolar field geometries (Jordanova et al., 2006, 2010). This improvement allows for integration of the 3D force balance magnetic field model (SCB, (Zaharia et al., 2004; Zaharia, 2008). This model balances the force with the divergence of the general pressure tensor to calculate the magnetic field configuration within its domain. The domain ranges from near the Earth’s surface, where the field is assumed dipolar, to the shell created by field lines passing through the equatorial plane at a radial distance of 6.5 R_{E}. Anisotropic pressure both at the outer boundary and inside the code’s domain is required and is provided by RAM. By relying on anisotropic pressure calculated by RAM, the force balance model creates a more stretched, more realistic field than isotropic MHD models that do not capture the ring current pressure build up and are typically very dipolar within 6.6 R_{E}. Initial coupling of these two codes is detailed in Zaharia et al. (2005), Jordanova et al. (2006), and Zaharia et al. (2006) details about the full coupling can be found in Zaharia et al. (2010). RAM provides anisotropic pressure to the 3D equilibrium code, which in return calculates the field aligned integrals required by RAM to calculate particle drift paths. The addition of self consistency creates significant differences in the ring current drift paths (Jordanova et al., 2006) and a depression in the nightside magnetic field (Zaharia et al., 2006). RAMSCB continues to see improvements in its algorithms and implementation related to robustness, efficiency, and performance during extreme driving (Engel et al., 2019). A comprehensive discussion of its coupling within the SWMF is given by Welling et al. (2018).
C.4 Energetic particle models
C.4.1 SEP Models
The acceleration and transport processes of energetic particles in interplanetary space is described by the focused transport equation in which the particle’s gyrophase is averaged out and the particle’s motion is reduced to the guiding center’s motion along the magnetic field and diffusion due to magnetic turbulence (Skilling, 1971; Kóta, 2000; Kóta et al., 2005; Qin et al., 2006; Zhang et al., 2009; Zhao et al., 2016, 2017)
where f is the particle’s distribution function, μ is the cosine of the particle’s pitch angle, D_{μμ} is the pitch angle diffusion coefficient, is solar wind velocity, p is the particle’s momentum in the solar wind frame, s is the direction along the magnetic field, v is the particle’s speed, and Q is the source term.
In the diffusive limit, where the distribution function is assumed to be isotropic, the focused transport equation reduces to the original Parker (1965) equation
where is the diffusion tensor along the magnetic field, with D being the scalar diffusion coefficient. In the Lagrangian coordinates advecting with the background plasma, the above governing equation is reduced to
where ρ(s,t) and B(s,t) are plasma density and total magnetic field magnitude along the magnetic field lines. The 3D problem is then reduced to a set of independent 1D problems on those timeevolving Lagrangian grids (Sokolov et al., 2004). In MFLAMPA the plasma and turbulence parameters along the magnetic field lines are extracted dynamically from the the BATSRUS models (the SC, IH and OH components). The diffusion process is treated as pitch angle scattering of the particles by the magnetic Alfvén waves calculated self consistently within the BATSRUS simulation.
C.4.2 GCR Models
Transport of GCRs in the heliosphere is affected by the solar modulation, which results in a dynamical change of the GCR energy spectrum and anisotropy as they propagate in the heliosphere (e.g., Vainio et al., 2008). The theory of modulation is based on solving the Parker (1965) equation with cosmic ray drift
where f(r,p,t) is the omnidirectional distribution function of GCRs, p the particle’s momentum, r the heliocentric distance, the solar wind velocity, the symmetric part of the diffusion tensor, the pitchangle averaged guiding center drift velocity (cf. Gombosi, 1998) and Q defines the source of the GCRs.
Equation (C.38) describes the main transport processes: (1) diffusion of particles due to their scattering off magnetic inhomogeneities, (2) convection in the outstreaming solar wind, (3) two types of drifts: the gradientcurvature drift in the regular heliospheric magnetic field, and drift along the heliospheric current sheet, and (4) adiabatic energy losses in the expanding solar wind. These processes are defined by the geometrical structure, polarity, strength, and the level of turbulence in the IMF and solar wind, which are ultimately driven by variable solar activity, leading to the temporal variability of the modulation on different timescales.
C.5 Simulating virtual magnetic observatories
There is an increasing number of ground magnetometer stations that provide magnetic field measurements. These observations can be directly compared to simulated observations obtained from SWMF simulations. The large number and high cadence of ground observations necessitates development of a fast and accurate way to generate synthetic magnetometer observations from our simulation results. Here we describe the algorithm that is used in the SWMF suite.
C.5.1 BiotSavart integral for currents in the magnetosphere
The contribution from the magnetospheric current system to the ground magnetic field is given by the BiotSavart integral (neglecting the displacement current):
where is the point where we calculate the synthetic magnetic field, is the contribution to this field from the magnetosphere currents, R_{M} is radius of the magnetosphereionosphere boundary (in our case the inner boundary of the magnetosphere model, R_{M} ≈ 2.5 R_{E} in most simulations) and is the volume element. The current density at a point in the simulation domain is expressed from Ampère’s law and the gradient operators and differentiate over coordinates and , respectively. The gradient of the inverse distance function is
Now the BiotSavart integral (Eq. (C.39)) can be written as
Next, we expand the double vector product using the , , and identities:
If the observation point, , is outside the simulation domain the last integral is zero. Since we are considering magnetic perturbations on the ground – that is outside the simulation region – this last integral can be neglected. Finally, we introduce the unit vector (note, that positive points into the domain of integration) and apply Gauss’ theorem to equation (C.42):
where dS is an area element on the spherical surface, . Note that equation (C.43) replaces the effect of all currents in the simulated magnetosphere with a surface current, and magnetic surface charge, . For the special case when equation (54) reduces to
This result agrees with the wellknown property of a potential field at the center of a sphere that equals the average of the field over a spherical surface (see Jackson, 1975) as long as the field given by equation (C.43) is created by currents located outside the sphere.
C.5.2 Magnetic field perturbations caused by fieldaligned currents
Another source of geomagnetic variations, which is particularly significant at high geomagnetic latitudes, is the magnetic field produced by the currents connecting the magnetosphereionosphere boundary, to the ionosphere and closing there. The currents in the gap region, between and the ionosphere, are fieldaligned, which means that the current density is parallel or antiparallel to the terrestrial magnetic field: , where represents a point in the gap region. This assumption allows us to derive the magnetic field from the field aligned currents.
Through the boundary of each surface element at the MI boundary, dS, a flux tube may be traced to the ionosphere boundary. The total current enclosed by this flux tube is
The field line (and current line) is described by , where is the path length element, and . Expressing in terms of the element of radial coordinate in the gap region, , one can express the equation for the current line in terms of dR′:
The magnetic field line given by equation (C.46) should be integrated from R_{M} down to R_{I}, starting from each point at the MI interface. This way a multitude of field lines, are obtained in the domain. Equation (C.43) can now be generalized to account for the BiotSavart integral from the multitude of field aligned currents, dI, in the gap region,
Note, that the integral over dR′ in equation (C.47) is a complicated vector function. However, its value only depends on and , and consequently, for any given computational grid and set of surface points, this function can be calculated only once (at the beginning of the simulation). After this, the contribution from the currents in the gap region, similarly to that from the magnetosphere currents, is given by the surface integral over the MI interface. We also note that only the derivatives of tangential components of the magnetospheric field at the MI interface contribute to the radial component of the current density. Therefore, only the magnetic field at the MI interface contributes to the surface magnetic variation, but not its radial gradient.
Equation (C.47) is a very important result for computational efficiency. It says that the BiotSavart integral along magnetic field lines going through the gap region can be written as the fieldaligned current multiplied by a constant that only depends on the field line that is approximated and the location of the point where the magnetic perturbation is calculated. These constants can be precalculated and stored (properly distributed among the processors) so that the integrals become a simple weighted sum, which is much faster to calculate than the integrals.
The electric current is a sum of contributions from separate current bundles (flux tubes), therefore one can use the BiotSavart integral to obtain the perturbation magnetic field at a surface point, , produced by the flux tube current, I, for a flux tube described by :
where the magnetic field line passing through the point at the MI interface is parameterized in the gap region with the radial coordinate, (, and ).
The expression for the magnetic field line in the gap region greatly simplifies if the Earth’s magnetic field is described in the dipole approximation. In this case the differential equation for the magnetic field line can be easily solved by projecting this vector equation on the direction of and on the two perpendicular directions. The solution is
where we used the notation, . This expression can be further simplified if we use one of the standard geocentric coordinate systems with the zaxis aligned with the direction of the magnetic dipole moment (the magnetic axis), such as MAG or SM (see Fränz & Harper, 2002).
C.6 Geomagnetic indexes
Geomagnetic indexes, including Dst, Kp, and the AE family of values, are a regular product of the SWMF Geospace. Dst is approximated via a single BiotSavart integral of all currents flowing in the GM component (typically, BATSRUS at Earth, see Sect. C.5.1). Kp and AE indexes leverage virtual magnetometer stations to more closely reflect the calculation of their realworld counterparts.
Realworld Kp is the average of localK index values calculated from 13 midlatitude ground observatories, rounded to the nearest third. LocalK is a range index of the maximum minus the minimum disturbance at a single observatory over set threehour windows (0–3 UT, 3–6 UT, etc), scaled to a 9value integer via a semilogarithmic transformation. The scale factors are station specific; additional adjustments are made for season. In the SWMF, the calculation is optimized for simplicity and performance. For the Kp calculation, 24 stations are used, spread equally about local time and all placed at a constant geomagnetic latitude of 60°, and are scaled such that K = 9 corresponds to ≥600 nT. Virtual Kp is then the average of these 24 stations, rounded to the nearest third. Because the value is written to file during the simulation, set threehour windows are not possible. Instead, a rolling threehour window that ends at the current simulation time is used. Therefore, SWMF virtual Kp range windows line up exactly with the realworld index when the universal time hour is a multiple of 3. No seasonal adjustments are made for virtual Kp. The use of the 24station approach provides a minor improvement in predictive performance compared to using the realworld stations and scalings.
The Auroral Electrojet (AE) family of indexes are highlatitude data products. The realworld indexes are the product of the geomagnetic northsouth perturbations obtained from 13 observatories. AU is the maximum perturbation of the 13 stations as a function of time (reported minutely); AL is the minimum; AE is the difference of AU and AL, and AO is the average of AU and AL. The calculation of virtual AE follows this closely except for the location of the contributing stations: 24 stations at a constant magnetic latitude of 70° evenly spaced in local time are used. Further work is required to optimize the location and number of stations used for virtual AE indexes.
C.7 New Diagnostics
The optically thin solar corona emits in the XUV, visible and IR wavelengths. Currently, there is only one mission (Parker Solar Probe) that takes insitu, local measurements of the solar corona, all information of the global corona measured by groundbased observatories (visible and IR) or via spacecrafts carrying remotesensing instrumentation (XUV). With synthetic observations of the solar corona one can evaluate the solar corona model performance, decompose and analyse physical process and the formation the radiation output. Within the SWMF the solar corona can be visualized via synthetic narrowband imaging (lineofsight, LOS) or synthetic spectra (SPECTRUM).
C.7.1 LineofSight Images
SWMF has the capability of generating various synthetic lineofsight (LOS) plots, such as EUV images. The response of each pixel of the image is treated as a LOS integral of a function f through the plasma:
where l follows along the LOS.
The LOS algorithm is implemented in the following parallel way: For each LOS ray and each grid block we determine first the segment of the ray that intersects the block. Then for each ray and block the function f is trilinear interpolated along the LOS and integrated according to equation (C.50) using a trapezoidal rule. The step size of the integration is proportional to the cell size of the block. Once the integration is finished for all blocks, we add for each LOS ray all integrals over the block segments via MPI reduce.
C.7.2 SPECTRUM
The Spectral Calculations for Global Space Plasma Modeling (SPECTRUM) code (Szente et al., 2019) calculates emissions from the optically thin solar corona by combining AWSoM(R) simulation results with the CHIANTI database (Dere et al., 1997, 2019). Dopplershifted, nonthermal line broadening due to lowfrequency Alfvén waves and anisotropic proton and isotropic electron temperatures can be individually taken into account during the calculations. The synthetic spectral calculations can then be used for model validation, for interpretation of solar observations, and for forward modeling purposes. SPECTRUM is implemented within the Space Weather Modeling Framework (SWMF) and is publicly available.
SPECTRUM is a postprocessing tool within the SWMF: it processes output after a simulation is completed. It is a standalone Fortran code that can process output files originating from any global coronal model, assuming that the data set is formatted appropriately. Currently, SPECTRUM can handle either a Cartesiangrid or the BATSRUS unstructured grid. SPECTRUM uses the same LOS integration technique as described in Section C.7.1.
The spectral calculation is performed the following way. It is assumed that the ion emissions coming from a given volume element dV follow a Gaussian profile, centered at wavelength λ_{0} with line width Δλ:
The total flux at an instrument’s detector at distance d is the sum of all the emission along the LOS from each volume element:
where is the density of the emitting ions, A_{ji} is the Einstein coefficient and ν_{ij} is the transition frequency from j to i. Rewriting the expression into separate density and temperature dependent terms (details in Szente et al., 2019):
G(T_{e}, N_{e}) is the contribution function, slowly varying with density and strongly dependent on temperature. The contribution function for each ion is calculated and saved into a lookup table using tables and procedures from CHIANTI v9 (Dere et al., 1997, 2019) from SolarSoft.
SPECTRUM takes two input files, one is the tabulated contributionfunction values and the BATSRUS output saved from an AWSoM simulation. The following calculations are performed on a linebyline basis, for each cell along the lineofsight. First we apply Dopplershift to the line center:
where u_{LOS} is the lineofsight bulk plasma velocity (positive toward the observer) and c_{light} is the speed of light.
The line width is the sum of instrumental broadening, and a thermal and a nonthermal component:
Thermal broadening is calculated considering the contribution along the lineofsight direction of the anisotropic temperature, the nonthermal component is due to the lowfrequency Alfvén wave contribution along the lineofsight.
The SPECTRUM code output is synthetic spectra or synthetic spectral images. An example is shown in Figures C.1 and C.2.
Fig. C.1 Synthetic spectra (green) compared to observation (black) taken by Hinode/EIS at 12 November 2007 12:32:02 UT of the Northern coronal hole during the Carrington Rotation 2063 (from Szente et al., 2019). Line profiles of Fe XI 201.734 Å, Fe XIII 201.121 Å and 202.044 Å are closely predicted in intensity, width and Dopplershift; while there is no line in the CHIANTI database between 201.5 and 201.6 Å, which explains the missing peak in the model’s result. 
Fig. C.2 Doppler map of line Fe XIII 202.044 Å is obtained from synthetic spectral image using the magnetic boundary of Carrington Rotation 2082, as could be perceived from Earth after the lineofsight integration of individual emission lines (from Szente et al., 2019). Blueshifts show regions where the solar wind is moving toward the observer, redshifts are present where plasma is moving away from the observer. The denser the plasma, the more dominant its effect is on the overall integrated Doppler shift of the observed line. 
Appendix D: Algorithms
Continuous development of the numerical algorithms is a necessity in order to maintain stateoftheart numerical models. Large interdisciplinary teams provide an ideal environment to learn about and adopt the best algorithms in a wide range of applications. Over the last two decades the algorithms of the SWMF and the models in it have improved tremendously. This section highlights some of the most important developments.
D.1 Advanced spatial discretization methods
The numerical error of the solution depends on several factors. In general the numerical error ε generated by the spatial discretization (which usually dominates) in a single time step at a given grid cell can be approximately written as
where k is some coefficient depending on the numerical method, x is the grid resolution and n is the spatial order of the scheme. There are at least three ways to reduce numerical error: reduce the coefficient k, reduce the grid resolution x or increase the order of the scheme n.
The most straightforward approach is to increase grid resolution. Doing this uniformly over the whole computational domain is very expensive. In fact, the computational cost scales roughly with for a threedimensional simulation, because the number of grid cells is and the time step Δt has to be kept proportional to Δx. A much better approach is to increase the grid resolution only where it is necessary. Thanks to its CFD heritage, BATSRUS was born with blockadaptive mesh refinement. This algorithm allows refining the grid where necessary, and coarsen it where possible. Using blockbased adaptation instead of cellbased adaptation (or fully unstructured grids) has distinct advantages for high performance massively parallel codes. In the past 20 years, the original blockadaptive grid implementation has been improved, extended, and in fact, completely rewritten into the Block Adaptive Tree Library (BATL) (Toth et al., 2012). BATL can use an arbitrary number of ghost cells, works in 1, 2 or 3 dimensions, allows for nonCartesian grids, allows grid adaptation based on geometric and physicsbased criteria, and it has very efficient algorithms that scale well to a large number of CPU cores.
Another way to change grid resolution is to use nonCartesian grids. For example, a spherical grid naturally has smaller cells in the longitude and latitude directions near the surface of the central body (Sun, planet, moon) than further away, which is advantageous for typical applications. A further refinement is to use a nonuniform grid spacing by applying some nonlinear stretching. A typical example is to make the grid points linear in the logarithm of the radius instead of the radius itself. Using ln(r) as a generalized coordinate will increase the radial resolution near the central body, which is usually beneficial. One can in fact use custom designed coordinate mapping to resolve specific regions, for example AWSoM uses a special stretching to concentrate cells around the transition region of the Sun. Combining generalized coordinates and adaptive mesh refinement provides great flexibility in using the optimal grid for a given problem.
The original version of BATSRUS used a second order total variation diminishing (TVD) scheme, which was stateoftheart in the 1990s (Powell et al., 1999). But computational fluid dynamics has evolved since then. Inspired by other codes, such as LFM (Lyon et al., 2004), we decided to extend BATSRUS to use higher order schemes. The space physics applications require the solution of complicated partial differential equations typically in three spatial dimensions. The solutions often contain discontinuities, such as shocks or current sheets. TVD schemes excel in maintaining monotonic profiles across shock waves, but at discontinuities the TVD scheme falls back to the first order upwind scheme, which means that the accuracy is only linearly improving with the reduction of the grid cell size.
In search of a suitable highorder scheme, we had the following requirements:
Minimal oscillations near discontinuities.
Conservative scheme that gives correct jump conditions.
High order at grid resolution changes and high order for nonCartesian grids.
Small stencil to allow for small grid blocks.
Only moderately more expensive than the second order TVD scheme.
General method that is high order for various system of equations including nonlinear terms.
To meet these requirements is very challenging. We looked at existing codes and explored the options published in the literature and presented at meetings. It is important to note that higher than second order accurate finite volume schemes require a high order accurate integral (quadrature) of the fluxes over the cell faces and a high order accurate quadrature of the source terms in the cell volume, which makes them rather complicated and expensive in multidimensional simulations. The LFM (Lyon et al., 2004) and GAMERA (Zhang et al., 2019a) codes, for example, are only higher than second order accurate in the finite difference sense for linear systems of equations. In addition, the use of a second order accurate update of the induction equation renders the overall scheme to be second order accurate only when the magnetic field plays an important role. Nevertheless, for the linearly high order donor cell algorithm the coefficient k is small in equation (D.1), which makes the LFM/GAMERA scheme exceptionally accurate, although still second order only (n = 2). After considerable experimentation, we have opted for a conservative finite difference scheme based on the fifth order accurate monotonicity preserving (MP5) limiter.
We have developed a new 5th order scheme (Chen et al., 2016) that satisfies all the requirements listed above. It is 5th order accurate for all terms in the MHD equations, it works for Cartesian and nonCartesian grids alike, and it remains globally 5th order accurate with adaptive mesh refinement included. The stencil is quite compact, so only three ghost cells are needed, which means that the grid blocks can be as small as 6 × 6 × 6 cells, which allows flexible adaptation. Using a third order RungeKutta scheme, the 5th order scheme is only about three times more expensive than the twostage 2nd order TVD method. We can further reduce the computational cost by restricting the 5th order scheme to a part of the computational domain. In Figure D.1, we computed the STEREO images for the three Fe emission lines 171 Å, 195 Å, and 284 Å. The top row is for the AWSoM solar wind model using a 2nd order scheme, while the middle row is for using the 5th order scheme, which gives more detail and compares favorably with the observations (bottom panels).
Fig. D.1 Comparison of synthesized EUV images of the model with observational STEREO A/EUVI images. The columns are from left to right for 171 Å, 195 Å, and 284 Å. Top panels: synthesized EUV images for 2nd order scheme. Middle panels: synthesized EUV images for 5th order scheme. Bottom panels: observational STEREO A/EUVI images. The observation time is 7 March 2011 20:00 UT. 
D.2 Advanced Time Integration Methods
Most numerical models employ an explicit time stepping scheme, where the values at the next time step in a given grid cell are calculated from the current values in the vicinity of this cell, the stencil. Explicit schemes are simple and fast, but the time step Δt is limited by the CourantFriedrichLewy (CFL) condition:
i.e., it cannot exceed the distance of the neighboring cells Δx divided by the fastest characteristic wave speed c_{max} of the system of equations. The proportionality factor C depends on the numerical scheme, but it is typically less than unity. The CFL condition is a simple but fundamental consequence of causality. When the solution changes due to the fastest waves, there is no magical bullet, the explicit method is optimal. In many cases, however, the solution changes at a much lower rate, because the fastest modes are not present in the solution. For example, the magnetosphere of the Earth typically changes due to changes in the solar wind and not due to propagating fast magnetosonic waves.
The BATSRUS code has immediately benefited from the aerospace CFD expertise: a simple but incredibly efficient way to accelerate convergence to a steady state solution is local time stepping. Each grid cell takes the largest time step allowed by the Courant condition. While propagating time at different rates in different grid cells is not physical, the final steady state will be still correct, as it finds a balance of the divergence of the fluxes and the source terms, where the time step is a simple multiplier that makes no difference. Local time stepping has been used routinely in the aerospace community, but was virtually unknown in the space physics community. In combination with adaptive mesh refinement, BATSRUS routinely obtains exact or approximate steady state solutions 10 to 1000 times faster than the simple explicit method on a static grid (Toth et al., 2012).
The Courant condition due to a fast wave mode is a specific example of stiffness. Stiffness means that the partial differential equations contain potentially large terms that happen to cancel each other out. Simple explicit time integration will blow up if the time step exceeds some restrictive limit. Implicit time integration offers a way to speed up the calculation: the fluxes and source terms are calculated from the values based on the next time step. Obviously, the values at the next time step are not yet known, hence the name, implicit. Typically, an implicit time integration scheme requires to solve a system of equations. The simplest example is a stiff source term, for example collisional terms among multiple species. Since these source terms do not involve spatial derivatives, one can solve the equation for state variables (for example densities) at each grid point independently, which is why this is called the pointimplicit method. When the stiff terms involve spatial derivatives, for example heat conduction, the system of equations involve all the grid cells together. Typically, we employ an iterative scheme to solve a linearized system. Since the rest of the equations are solved explicitly, this method is called semiimplicit. Finally, one may solve the full system of equations implicitly with an iterative scheme, which alleviates all the stability restrictions. Solving a large system of equations is, of course, computationally expensive. It only makes sense if the time step can be increased sufficiently to beat the efficiency of the explicit method. The time step of an implicit scheme is always limited by accuracy considerations. The various implicit schemes in BATSRUS originate from an interdisciplinary project of applied mathematicians, computer scientists and plasma physicists in the Netherlands in the 1990s (Toth et al., 1998).
One does not have to choose a certain time integration scheme for the whole computational domain. In fact, in some applications the pointimplicit method is used only where the stiff source terms are present (for example the collisional terms are only important in the ionosphere of Mars), the semiimplicit scheme may be limited to the region where the Hall term is important (for example the magnetotail), and the implicit scheme for the full set of equations may also be combined with the explicit method in an adaptive manner based on the stability constraint for a given time step (Toth et al., 2006). In practice, we are using all of these schemes in various combinations. The optimal choice depends on the application and it can be orders of magnitude faster than the simple explicit time stepping. Figure D.2 shows how BATSRUS implements the various time integration schemes in a hierarchical manner. This allows using the different schemes independently or combined for the various equation sets and applications.
Fig. D.2 The layered software structure of BATSRUS. The arrows point from the module that is using data or methods from the other module. There are multiple versions of the equation and user modules. The various time stepping schemes are independent of the details of the equations being solved (from Toth et al., 2012). 
A particularly interesting application of advanced time stepping algorithms is for particleincell (PIC) codes. Explicit PIC models are limited by the Courant condition for light waves. In addition, the grid resolution is also limited: the model has to resolve the Debye length, which can be exceedingly small compared to the scales of the full system. The usual remedy is to artificially increase the electron mass and reduce the speed of light, but explicit PIC remains extremely expensive even with these tricks. Using an implicit algorithm removes both the spatial and temporal limitations: the semiimplicit particleincell algorithm (Brackbill & Forslund, 1982) can use arbitrarily large grid cells and the time step is limited by the Courant condition for the thermal velocities of the particles rather than the speed of light.
This advance was a game changer for modeling large systems with a kinetic model. A further major improvement was the energy conserving semiimplicit method (ECSIM) developed by Lapenta (2017). Energy conservation is crucial to maintain the longterm stability of PIC models. Before ECSIM, the original remedy was applying some smoothing on the electric field that required experimentation with the amount of smoothing and resulted in somewhat diffused solutions. ECSIM, however, did not take care of enforcing Gauss’ law that the divergence of the electric field equals the net charge density, which limited its applicability. This problem was solved by Chen et al. (2018), who developed the Gauss’ Law satisfying ECSIM algorithm. GLECSIM is now the workhorse PIC algorithm used in the various kinetic models (iPIC3D, AMPS and FLEKS) in the SWMF.
D.3 Hybrid schemes
The SWMF allows applying different models in different physics domains. This flexibility is crucial to model a complex multiscale and multiphysics system. A similar approach is used in BATSRUS to apply different numerical methods (for example high order vs. second order scheme) or even different physics (for example Hall MHD vs. ideal MHD) in different regions to achieve optimal performance. The regions can be selected using geometry and/or physicsbased criteria. We have developed a general library that can define regions in the computational domain using addition and subtraction of simple geometrical objects (boxes, cylinders, spheres, cones, paraboloids, shells, rings, etc.) This allows the user to define regions of complicated shape from the input parameter file. When the region is dynamic, it can be defined by some local physical quantities based on the numerical solution. For example, one can use some threshold for the current density to apply adaptive mesh refinement. The two approaches can also be used in combination, for example the mesh refinement based on current density can be restricted to a certain part of the magnetotail. These capabilities are now available for a variety of options: adaptive mesh refinement, high vs. low order scheme, explicit vs. implicit schemes, Hall term, resistivity, heat conduction, viscosity, conservative vs. nonconservative energy equation, semirelativistic correction, etc. To minimize numerical artifacts at the interfaces, we allow for a linear tampering in the critical parameter when applicable, for example the coefficient of the Hall term is 0 outside the Hall region, 1 inside, and varies linearly from 0 to 1 at the interface of a finite width.
Using different schemes with different computational costs in the computational domain poses new requirements for the load balancing algorithm. Our approach is to assign a type for each grid block based on the combination of numerical schemes used. Blocks of the same type use the same combination of schemes, so their computational cost is similar. Then we load balance the various types of grid blocks independently. As long as there are enough grid blocks to fill the CPU cores, this approach works well.
D.4 Achieving and maintaining high performance
Both BATSRUS and the SWMF were designed to achieve high performance on massively parallel super computers. BATSRUS uses a blockadaptive grid for multiple reasons: it makes load balancing simple, it provides fixed loop sizes over the grid cells that can be easily optimized by the compiler, and the amount of data associated with each grid block can fit into the cache memory. The original design of BATSRUS was based on the Message Passing Interface (MPI) library, which is still the most used parallelization tool on current supercomputers. All these design features improve performance and parallel scaling. In fact, in 1997 BATSRUS achieved 13 Gflops on 512 cores of a Cray T3D computer. At that time this was among the largest supercomputers and >10 Gflop performance with excellent parallel scaling was a heroic achievement. The current version demonstrated nearly perfect weak scaling up to 250,000 cores of a Cray supercomputer, as shown in Figure A.3. While this is very respectable and more than sufficient for current supercomputers, we have to prepare the code for future architectures with even more cores and less memory per core. One reason we cannot run the code on even more cores is that the data structure describing the blockadaptive grid keeps growing with the problem size. To mitigate this issue, we have implemented a hybrid MPI+OpenMP parallelization (Zhou & Tóth, 2020): MPI is used to parallelize over the CPU nodes, while OpenMP is used to use multithreading over the cores on a single node. This means that large data structures can be shared by multiple threads, which reduces the memory use substantially. Using this hybrid parallelizaton, BATSRUS could run up to 500,000 cores of the Blue Waters super computer while still maintaining excellent performance.
The next frontier is porting large simulation codes to GPUs. Using support from NSF, we have started to work on porting BATSRUS to GPUs. Our current approach is using the OpenACC library to parallelize loops over grid cells and run them on separate GPU threads. This work is in a preliminary phase now, but we already have some simple tests running on one GPU.
The SWMF was designed to be as light weight as possible. The models can run serially or concurrently and synchronization is only performed when necessary (Toth, 2006). The current SWMF also supports models running with OpenMP: each model can use a different number of threads. Typically one thread per core is used, but hyperthreading is also supported. Coupling between the models also needs to be efficient, especially when a large amount of information is exchanged frequently. We have developed efficient and flexible coupling libraries that allow direct parallel coupling between two massively parallel models. As Figure D.3 shows, the coupled BATSRUS and iPIC3D models scale to 32,000 cores with minimal loss of efficiency. The two models exchange the MHD quantities calculated from the PIC distribution function for every MHD grid cell inside the PIC region every time step.
Fig. D.3 The cell and particle update rates as a function of number of cores for the SWMF running the twoway coupled BATSRUS and iPIC3D models. The problem size scales in proportion to the number of parallel processes. The dotted lines represent linear scaling. 
Appendix E: MHDEPIC and MHDAEPIC
The magnetohydrodynamics with embedded particleincell (MHDEPIC) algorithm allows global MHD simulations performed with the kinetic physics properly handled by a PIC model in regions of interest (Daldorff et al., 2014). The PIC domain can cover the regions where kinetic effects are most important, such as reconnection sites. In the newly developed MHD with adaptively embedded PIC (MHDAEPIC) algorithm, the PIC domain consists of small blocks that can be adaptively activated and deactivated to cover the dynamically changing regions of interest. While keeping the expensive PIC model restricted to a small region, the BATSRUS code can efficiently handle the rest of the computational domain where the MHD or Hall MHD description is sufficient. Since the PIC model is able to describe the electron behavior selfconsistently, our coupled MHDEPIC and MHDAEPIC models are well suited for investigating the nature and role of magnetic reconnection in space weather phenomena.
Figure E.1 illustrates the overall flow of the coupling algorithm between BATSRUS and iPIC3D (Daldorff et al., 2014), while Figure E.2 shows the spatial discretization of the coupling. It is important to point out that the BATSRUS – iPIC3D coupling via SWMF is genuinely 2way: all physical quantities are selfconsistently advanced by both codes and the relevant information is fully exchanged in every timestep.
Fig. E.1 The overall flow of MHDPIC coupling. At t = 0 the MHD code sends the MHD state inside and around the PIC region to the PIC code. Both the MHD and PIC codes then advance by one or more time steps until both models reach the next coupling time. Information is exchanged both ways, but this time the PIC code only uses the MHD solution as a boundary condition, while the MHD code overwrites its solution in the PIC region with the PIC solution. This process continues until they reach the final simulation time or until the PIC region is removed (after Daldorff et al., 2014). 
Fig. E.2 Spatial discretization of the MHDEPIC coupling. The Cartesian grid of the PIC region is indicated with the darker gray area. The lighter gray area shows the ghost cell/node region of the PIC grid. The large red dots are node values obtained from the MHD variables. The small red dots illustrate particles created in the ghost cells of the PIC grid. The small red squares are the ghost cell centers of the PIC grid where the magnetic field is set from the MHD solution. The black dots indicate the MHD cell centers where the solution is obtained from the PIC code. The MHD grid can be either Cartesian or spherical (after Chen et al., 2017). 
MHDEPIC (Daldorff et al., 2014; Chen et al., 2017; Toth et al., 2017; Zhou & Tóth, 2020) uses the node centered number densities, velocities, pressures and magnetic field (large red dots in Fig. E.2) to create the macroparticles inside the ghost cells of the PIC grid, as illustrated by the small red dots in the light gray area in Figure E.2. The particles that leave the PIC region (dark gray area in the figure) are discarded at the end of the PIC time step. New macroparticles are generated for each species in each ghost cell of the PIC domain with the appropriate (bi)Maxwellian distribution functions using the MHD solution. The locations of the new particles are random with a uniform distribution over the ghost cell. For each ghost cell the corresponding number density, velocity and pressure are linearly interpolated from the surrounding MHD values (large red dots in the figure) to the given location. In this twoway coupled method the MHD values in the cell centers covered by the nodes of the PIC grid (black dots in Fig. E.2) are fully overwritten by the PIC solution. The magnetic field can simply be interpolated from the PIC field. For the other MHD variables MHDEPIC (Daldorff et al., 2014; Chen et al., 2017; Toth et al., 2017; Zhou & Tóth, 2020) and takes various moments of the distribution function represented by the macroparticles.
The SWMF and the BATSRUS codes require a Fortran compiler and the MPI library. The iPiC3D code requires a C++ compiler and the MPI library and optionally the parallel HDF5 library. Very good scaling up to 32K MPI processes on 1024 XE nodes of the Blue Waters supercomputer were achieved, as shown in Figures A.3 and D.3.
Currently, SWMF, BATSRUS and iPiC3D mainly use pure MPI parallelism that works fine up to 1024 XE nodes and 32 thousand MPI processes (see Figs. A.3 and D.3). While this is more than sufficient for most applications, the codes encounter some limitations when running with 65K and more MPI processes. Using OpenMP parallelism on the nodes can reduce memory usage and the number of messages sent between the MPI processes. The OpenMP+MPI approach does not improve the performance relative to the pure MPI parallelization, but the hybrid approach allows running larger problems on a larger number of nodes (Zhou & Tóth, 2020).
References
 Abadi M, Agarwal A, Barham P, Brevdo E, Chen Z, et al. 2015. TensorFlow: Largescale machine learning on heterogeneous systems. Software available from tensorflow.org, URL http://tensorflow.org/. [Google Scholar]
 Abbett WP, Fisher GH. 2003. A coupled model for the emergence of active region magnetic flux into the solar corona. Astrophys J 582: 475. https://doi.org/10.1086/344613. [Google Scholar]
 Abbett WP, Mikić Z, Linker JA, McTiernan JM, Magara T, Fisher G. 2004. The photospheric boundary of SuntoEarth coupled models. J Atmos SolarTerr Phys 66: 1257–1270. https://doi.org/10.1016/j.jastp.2004.03.016. [Google Scholar]
 AlvaradoGómez JD, Drake JJ, Garraffo C, Cohen O, Poppenhaeger K, Yadav RK, Moschou SP. 2020. An Earthlike stellar wind environment for proxima centauri c. Astrophys J 902(1): L9. https://doi.org/10.3847/20418213/abb885. [Google Scholar]
 Anderson BJ, Takahashi K, Toth BA. 2000. Sensing global Birkeland currents with Iridium engineering magnetometer data. Geophys Res Lett 27(24): 4045–4048. https://doi.org/10.1029/2000GL000094. [Google Scholar]
 André N, Grande M, Achilleos N, Barthélémy M, Bouchemit M, et al. 2018. Virtual Planetary Space Weather Services offered by the Europlanet H2020 Research Infrastructure. Planet Space Sci 150: 50–59. https://doi.org/10.1016/j.pss.2017.04.020. [NASA ADS] [CrossRef] [Google Scholar]
 Angelopoulos V, Kennel CF, Coroniti FV, Pellat R, Kivelson MG, Walker RJ, Russell CT, Baumjohann W, Feldman WC, Gosling JT. 1994. Statistical characteristics of bursty bulk flow events. J Geophys Res 99(21): 257. https://doi.org/10.1029/94ja01263. [Google Scholar]
 Angelopoulos V, McFadden JP, Larson D, Carlson CW, Mende SB, et al. 2008. Tail reconnection triggering substorm onset. Science 321(5891): 931–935. https://doi.org/10.1126/science.1160495. [NASA ADS] [CrossRef] [Google Scholar]
 Antiochos SK, DeVore CR, Klimchuk JA. 1999. A model for solar coronal mass ejections. Astrophys J 510: 485–493. https://doi.org/10.1086/306563. [Google Scholar]
 Badavi FF, Nealy JE, Wilson JW. 2011. The Low Earth Orbit validation of a dynamic and anisotropic trapped radiation model through ISS measurements. Adv Space Res 48(8): 1441–1458. https://doi.org/10.1016/j.asr.2011.06.009. [Google Scholar]
 Baker D, Kanekal S. 2008. Solar cycle changes, geomagnetic variations, and energetic particle properties in the inner magnetosphere. J Atmos SolarTerr Phys 70(2–4): 195–206. https://doi.org/10.1016/j.jastp.2007.08.031. [Google Scholar]
 Baker DN, Balstad R, Bodeau JM, Cameron E, Fennell JF, et al. 2009. Severe Space Weather EventsUnderstanding Societal and Economic Impacts Workshop Report. National Academies Press, Washington, DC. ISBN 9780309138116. https://doi.org/10.17226/12643. [Google Scholar]
 Baker DN, Pulkkinen TI, Angelopoulos V, Baumjohann W, McPherron RL. 1996. Neutral line model of substorms: Past results and present view. J Geophys Res 101(A6): 12975–13010. https://doi.org/10.1029/95JA03753. [Google Scholar]
 Barakat A, Schunk R. 2001. Effects of wave “particle interactions on the dynamic behavior of the generalized polar wind. J Atmos SolarTerr Phys 63(1): 75–83. https://doi.org/10.1016/S13646826(00)001061. [Google Scholar]
 Barakat AR, Barghouthi IA. 1994. The effect of waveparticle interactions on the polar wind O^{+}. Geophys Res Lett 21(21): 2279–2282. https://doi.org/10.1029/94gl01701. [Google Scholar]
 Bauske R, Nagy AF, Dezeeuw DL, Gombosi TI, Powell KG. 2000. 3D multiscale mass loaded MHD simulations of the solar wind interaction with Mars. Adv Space Res 26: 1571–1575. https://doi.org/10.1016/s02731177(00)001058. [Google Scholar]
 Bauske R, Nagy AF, Gombosi TI, De Zeeuw DL, Powell KG, Luhmann JG. 1998. A threedimensional MHD study of solar wind mass loading processes at Venus: Effects of photoionization, electron impact ionization, and charge exchange. J Geophys Res 103(A10): 23625–23638. https://doi.org/10.1029/98JA01791. [Google Scholar]
 Berger MJ, Colella P. 1989. Local adaptive mesh refinement for shock hydrodynamics. J Comput Phys 82(1): 67–84. https://doi.org/10.1016/00219991(89)900351. [Google Scholar]
 Berger MJ, Jameson A. 1985. Automatic adaptive grid refinement for the Euler equations. AIAA J 23(4): 561–568. https://doi.org/10.2514/3.8951. [Google Scholar]
 Bhatnagar PL, Gross EP, Krook M. 1954. A model for collision processes in gases I. Phys Rev 94: 511–525. https://doi.org/10.1103/physrev.94.511. [Google Scholar]
 Bilitza D. 2001. International reference ionosphere 2000. Radio Sci 36(2): 261–275. https://doi.org/10.1029/2000RS002432. [CrossRef] [Google Scholar]
 Bolduc L. 2002. GIC observations and studies in the HydroQubec power system. J Atmos SolarTerr Phys 64(16): 1793–1802. https://doi.org/10.1016/S13646826(02)001281. [Google Scholar]
 Borovikov D, Sokolov IV, Manchester WB, Jin M, Gombosi TI. 2017a. Eruptive event generator based on the GibsonLow magnetic configuration. J Geophys Res 122(8): 7979–7984. https://doi.org/10.1002/2017JA024304. [Google Scholar]
 Borovikov D, Sokolov IV, Roussev II, Taktakishvili A, Gombosi TI. 2018. Toward a quantitative model for simulation and forecast of solar energetic particle production during gradual events. I. Magnetohydrodynamic background coupled to the SEP Model. Astrophys J 864(1): 88. https://doi.org/10.3847/15384357/aad68d. [Google Scholar]
 Borovikov D, Tenishev V, Gombosi TI, Guidoni SE, DeVore CR, Karpen JT, Antiochos SK. 2017b. Electron Acceleration in Contracting Magnetic Islands during Solar Flares. Astrophys. J. 835(1): 48. https://doi.org/10.3847/15384357/835/1/48. [Google Scholar]
 Borovsky JE, Welling DT, Thomsen MF, Denton MH. 2014. Longlived plasmaspheric drainage plumes: Where does the plasma come from? J Geophys Res 119(8): 6496–6520. https://doi.org/10.1002/2014JA020228. [Google Scholar]
 Brackbill JU, Forslund D. 1982. An implicit method for electromagnetic plasma simulation in two dimensions. J Comput Phys 46: 271–308. https://doi.org/10.1016/00219991(82)90016X. [Google Scholar]
 Brackbill JU, Forslund DW. 1986. Simulation of lowfrequency electromagnetic phenomena in plasmas. In: Multiple Time Scales. Brackbill JU, Cohen BI, (Eds.) Academic Press, New York, NY. pp. 271–310. [Google Scholar]
 Brackbill JU, Lapenta G. 2008. Magnetohydrodynamcis with implicit plasma simulation. Comm Comput Phys 4: 433–456. [Google Scholar]
 Brecht S, Lyon J, Fedder J, Hain K. 1981. A simulation study of EastWest IMF effects on the magnetosphere. Geophys Res Lett 8: 397–400. https://doi.org/10.1029/GL008i004p00397. [Google Scholar]
 Brecht S, Lyon J, Fedder J, Hain K. 1982. A timedependent threedimensional simulation of the Earth’s magnetosphere: Reconnection events. J Geophys Res 87(A8): 6098–6108. https://doi.org/10.1029/ja087ia08p06098. [Google Scholar]
 Brillouin L. 1926. La mécanique ondulatoire de Schrödinger: une méthode générale de resolution par approximations successives. Comptes Rendus de l’Acadmie des Sciences 183: 24–26. [Google Scholar]
 Burch JL, Torbert RB, Phan TD, Chen LJ, Moore TE, et al. 2016. Electronscale measurements of magnetic reconnection in space. Science 352: 6290. https://doi.org/10.1126/science.aaf2939. [Google Scholar]
 Burgers JM. 1969. Flow equations for composite gases. Academic Press, New York. [Google Scholar]
 Buzulukova N, Fok MC, Pulkkinen A, Kuznetsova M, Moore TE, Glocer A, Brandt PC, Toth G, Rastätter L. 2010. Dynamics of ring current and electric fields in the inner magnetosphere during disturbed periods: CRCMBATSRUS coupled model. J Geophys Res 115(A05): 210. https://doi.org/10.1029/2009JA014621. [Google Scholar]
 Camporeale E. 2019. The challenge of machine learning in Space Weather: Nowcasting and forecasting. Space Weather 17: 1166–1207https://doi.org/10.1029/2018sw002061. [CrossRef] [Google Scholar]
 Carpenter D, Anderson RR. 1992. An ISEE/whistler model of equatorial electron density in the magnetosphere. J Geophys Res 97(A2): 1097. https://doi.org/10.1029/91ja01548. [Google Scholar]
 Chapman S. 1916. On the Law of Distribution of Molecular Velocities, and on the Theory of Viscosity and Thermal Conduction, in a NonUniform Simple Monatomic Gas. Philos Trans Royal Soc A 216: 279–348. https://doi.org/10.1098/rsta.1916.0006. [Google Scholar]
 Chen MW, Lemon CL, Hecht J, Sazykin S, Wolf RA, Boyd A, Valek P. 2019a. Diffuse auroral electron and ion precipitation effects on RCM E comparisons with satellite data during the 17 March 2013 storm. J. Geophys. Res. 124(6): 4194–4216. https://doi.org/10.1029/2019JA026545. [Google Scholar]
 Chen Y, Manchester WB, Hero AO, Toth G, DuFumier B, Zhou T, Wang X, Zhu H, Sun Z, Gombosi TI. 2019b. Identifying Solar Flare Precursors Using Time Series of SDO/HMI Images and SHARP Parameters. Space Weather 17(10): 1404–1426. https://doi.org/10.1029/2019SW002214. [Google Scholar]
 Chen Y, Meng XL, Wang X, van Dyk DA, Marshall HL, Kashyap VL. 2018. Calibration concordance for astronomical instruments via multiplicative shrinkage. J Am Stat Assoc 114: 1–34. https://doi.org/10.1080/01621459.2018.1528978. [Google Scholar]
 Chen Y, Toth G. 2019. Gauss’s law satisfying energyconserving semiimplicit particleincell method. J Comput Phys 386: 632–652. https://doi.org/10.1016/j.jcp.2019.02.032. [Google Scholar]
 Chen Y, Toth G, Cassak P, Jia X, Gombosi TI, Slavin JA, Markidis S, Peng IB, Jordanova VK, Henderson MG. 2017. Global threedimensional simulation of Earth’s dayside reconnection using a twoway coupled magnetohydrodynamics with embedded particleincell model: initial results. J Geophys Res 122(10): 10318–10335. https://doi.org/10.1002/2017JA024186. [Google Scholar]
 Chen Y, Toth G, Gombosi TI. 2016. A fifthorder finite difference scheme for hyperbolic equations on blockadaptive curvilinear grids. J Comput Phys 305: 604–621. https://doi.org/10.1016/j.jcp.2015.11.003. [Google Scholar]
 Chen Y, Tóth G, Hietala H, Vines SK, Zou Y, Nishimura Y, Silveira MVD, Guo Z, Lin Y, Markidis S. 2020. Magnetohydrodynamic with embedded particleincell simulation of the geospace environment modeling dayside kinetic processes challenge event. Earth Space Sci 7(11): e2020EA001331. https://doi.org/10.1029/2020EA001331. [Google Scholar]
 Chen Y, Tóth G, Jia X, Slavin JA, Sun W, Markidis S, Gombosi TI, Raines JM. 2019c. Studying dawndusk asymmetries of Mercury’s magnetotail using MHDEPIC simulations. J Geophys Res 124(11): 8954–8973. https://doi.org/10.1029/2019JA026840. [Google Scholar]
 Chew GF, Goldberger ML, Low FE. 1956. The Boltzmann equation and the onefluid hydromagnetic equations in the absence of particle collisions. Proc R Soc Lond A 236(1204): 112–118. https://doi.org/10.1098/rspa.1956.0116. [Google Scholar]
 Clarke DA. 2010. On the Reliability of ZEUS3D. Astrophys J Suppl 187(1): 119–134. https://doi.org/10.1088/00670049/187/1/119. [Google Scholar]
 Cohen O, Drake JJ, Kashyap VL, Sokolov IV, Gombosi TI. 2010. The impact of hot Jupiters on the spindown of their host stars. Astrophys J Lett 723(1): L64–L67. https://doi.org/10.1088/20418205/723/1/L64. [Google Scholar]
 Cohen O, Garraffo C, Moschou SP, Drake JJ, AlvaradoGómez JD, Glocer A, Fraschetti F. 2020. The space environment and atmospheric Joule heating of the habitable zone exoplanet TOI700d. Astrophys J 897(1): 101. https://doi.org/10.3847/15384357/ab9637. [Google Scholar]
 Cohen O, Ma Y, Drake JJ, Glocer A, Garraffo C, Bell JM, Gombosi TI. 2015. The interaction of Venuslike, Mdwarf planets with the stellar wind of their host star. Astrophys J 806(1): 41. https://doi.org/10.1088/0004637X/806/1/41. [Google Scholar]
 Cohen O, Sokolov IV, Roussev II, Arge CN, Manchester WB, et al. 2007. A semiempirical magnetohydrodynamical model of the solar wind. Astrophys J Lett 654: L163–L166. https://doi.org/10.1086/511154. [Google Scholar]
 Combi MR, Kabin K, Gombosi T, De Zeeuw D, Powell K. 1998. Io’s plasma environment during the Galileo flyby: Global threedimensional MHD modeling with adaptive mesh refinement. J Geophys Res 103(A5): 9071–9081. https://doi.org/10.1029/98JA00073. [Google Scholar]
 Cranmer SR, Van Ballegooijen AA. 2010. Can the solar wind be driven by magnetic reconnection in the Sun’s magnetic carpet? Astrophys J 720(1): 824–847. https://doi.org/10.1088/0004637X/720/1/824. [Google Scholar]
 Cravens T, Waite J, Gombosi T, Lugaz N, Gladstone G, Mauk B, MacDowall R. 2003. Implications of Jovian Xray emission for magnetosphereionosphere coupling. J Geophys Res 108: 1465. https://doi.org/10.1029/2003JA010050. [Google Scholar]
 Crew GB, Chang T, Retterer JM, Peterson WK, Gurnett DA, Huff RL. 1990. Ion cyclotron resonance heated conics: Theory and observations. J Geophys Res 95(A4): 3959–3985. https://doi.org/10.1029/JA095iA04p03959. [Google Scholar]
 Daldorff LK, Toth G, Gombosi TI, Lapenta G, Amaya J, Markidis S, Brackbill JU. 2014. Twoway coupling of a global Hall magnetohydrodynamics model with a local implicit particleincell model. J Comput Phys 268: 236–254. https://doi.org/10.1016/j.jcp.2014.03.009. [NASA ADS] [CrossRef] [Google Scholar]
 De Zeeuw DL, Sazykin S, Wolf RA, Gombosi TI, Ridley AJ, Toth G. 2004. Coupling of a global MHD code and an inner magnetospheric model: Initial results. J Geophys Res 109(A12): A12219. https://doi.org/10.1029/2003ja010366. [Google Scholar]
 Dere K, Landi E, Mason H, MonsignoriFossi B, Young P. 1997. CHIANTI  an atomic database for emission lines. Astron Astrophys Suppl Ser 125(1): 149–173. https://doi.org/10.1051/aas:1997368. [Google Scholar]
 Dere KP, Del Zanna G, Young PR, Landi E, Sutherland RS. 2019. CHIANTI – An atomic database for emission lines. XV. Version 9, Improvements for the XRay Satellite Lines. Astrophys J Suppl 241(2): 22. https://doi.org/10.3847/15384365/ab05cf. [Google Scholar]
 Downs C, Roussev II, van der Holst B, Lugaz N, Sokolov IV, Gombosi TI. 2010. Toward a realistic thermodynamic magnetohydrodynamic model of the global solar corona. Astrophys J 712(2): 1219–1231. https://doi.org/10.1088/0004637X/712/2/1219. [Google Scholar]
 Dudley RJ, Duchene N. 2010. Microsoft Azure: Enterprise Application Development. Packt Publishing, Birmingham, UK. ISBN 1849680981. [Google Scholar]
 Engel MA, Morley SK, Henderson MG, Jordanova VK, Woodroffe JR, Mahfuz R. 2019. Improved Simulations of The Inner Magnetosphere During High Geomagnetic Activity With the RAM SCB Model. J. Geophys. Res. 124(6): 4233–4248. https://doi.org/10.1029/2018JA026260. [Google Scholar]
 Enskog D. 1917. Kinetische Theorie der Vorgänge in mässig verdünnten Gasen. Ph.D. thesis, University of Uppsala, Sweden. Also published in Kungliga Suenska vetenskapsakademiens handlingar, 63(4), 1922. [Google Scholar]
 Fang F, Manchester WB IV, Abbett WP, van der Holst B. 2012a. Buildup of magnetic shear and free energy during flux emergence and cancellation. Astrophys J 754: 15. https://doi.org/10.1088/0004637X/754/1/15. [Google Scholar]
 Fang F, Manchester WB IV, Abbett WP, van der Holst B. 2012b. Dynamic coupling of convective flows and magnetic field during flux emergence. Astrophys J 745: 37. https://doi.org/10.1088/0004637X/745/1/37. [Google Scholar]
 FEMA. 2019. 2019 National Threat and Hazard Identification and Risk Assessment (THIRA). Tech. Rep. FEMA2019508c. Federal Emergency Management Administration, Washington, DC. [Google Scholar]
 Feng X, Xiang C, Zhong D, Zhou Y, Yang L, Ma X. 2014a. SIPCESE MHD model of solar wind with adaptive mesh refinement of hexahedral meshes. Comput Phys Commun 185(7): 1965–1980. https://doi.org/10.1016/j.cpc.2014.03.027. [Google Scholar]
 Feng X, Zhang M, Zhou Y. 2014b. A new threedimensional solar wind model in spherical coordinates with a sixcomponent grid. Astrophys J Suppl 214(1): 6. https://doi.org/10.1088/00670049/214/1/6. [Google Scholar]
 Fok MC, Buzulukova NY, Chen SH, Glocer A, Nagai T, Valek P, Perez JD. 2014. The Comprehensive Inner MagnetosphereIonosphere Model. J Geophys Res 119(9): 7522–7540. https://doi.org/10.1002/2014JA020239. [Google Scholar]
 Fok MC, Kang SB, Ferradas CP, Buzulukova NB, Glocer A, Komar CM. 2021. New Developments in the Comprehensive Inner MagnetosphereIonosphere (CIMI) Model. J Geophys Res 126: e2020JA028987. https://doi.org/10.1029/2020JA028987. [Google Scholar]
 Fok MC, Kozyra JU, Nagy AF, Rasmussen CE, Khazanov GV. 1993. A decay model of equatorial ring current and the associated aeronomical consequences. J Geophys Res 98(19): 381. [Google Scholar]
 Fok MH, Horne RE, Meredith NP, Glauert SA. 2008. Radiation Belt Environment model: Application to space weather nowcasting. J Geophys Res 113: A03S08. https://doi.org/10.1029/2007JA012558. [Google Scholar]
 Fränz M, Harper D. 2002. Heliospheric coordinate systems. Planet Space Sci 50: 217–233. https://doi.org/10.1016/S00320633(01)001192. [Google Scholar]
 Gibson S, Low BC. 1998. A timedependent threedimensional magnetohydrodynamic model of the coronal mass ejection. Astrophys J 493: 460–473. https://doi.org/10.1086/305107. [Google Scholar]
 Glocer A, Fok M, Meng X, Toth G, Buzulukova N, Chen S, Lin K. 2013. CRCM + BATSRUS twoway coupling. J Geophys Res 118(4): 1635–1650. https://doi.org/10.1002/jgra.50221. [Google Scholar]
 Glocer A, Gombosi TI, Toth G, Hansen KC, Ridley AJ, Nagy A. 2007. The Polar Wind Outflow Model: Saturn Results. J Geophys Res 112(A01): 304. https://doi.org/10.1029/2006JA011755. [Google Scholar]
 Glocer A, Khazanov G, Liemohn M. 2017. Photoelectrons in the quiet polar wind. J Geophys Res 122(6): 6708–6726. https://doi.org/10.1002/2017JA024177. [Google Scholar]
 Glocer A, Kitamura N, Toth G, Gombosi T. 2012. Modeling solar zenith angle effects on the polar wind. J Geophys Res 117: https://doi.org/10.1029/2011JA017136. [Google Scholar]
 Glocer A, Rastätter L, Kuznetsova M, Pulkkinen A, Singer HJ, et al. 2016. Communitywide validation of geospace model local Kindex predictions to support model transition to operations. Space Weather 14(7): 469–480. https://doi.org/10.1002/2016SW001387. [CrossRef] [Google Scholar]
 Glocer A, Toth G, Fok M, Gombosi T, Liemohn M. 2009a. Integration of the radiation belt environment model into the space weather modeling framework. J Atmos SolarTerr Phys 71(16): 1653–1663. https://doi.org/10.1016/j.jastp.2009.01.003. [Google Scholar]
 Glocer A, Toth G, Fok MC. 2018. Including Kinetic Ion Effects in the Coupled Global Ionospheric Outflow Solution. J Geophys Res 123(4): 2851–2871. https://doi.org/10.1002/2018JA025241. [Google Scholar]
 Glocer A, Toth G, Gombosi T, Welling D. 2009b. Modeling ionospheric outflows and their impact on the magnetosphere, initial results. J Geophys Res 114: A05216. https://doi.org/10.1029/2009JA014053. [Google Scholar]
 Glocer A, Tóth G, Ma Y, Gombosi T, Zhang JC, Kistler LM. 2009c. Multifluid blockadaptivetree solar wind Roetype upwind scheme: Magnetospheric composition and dynamics during geomagnetic stormsInitial results. J Geophys Res 114(A12): A12203. https://doi.org/10.1029/2009ja014418. [Google Scholar]
 Glocer A, Welling D, Chappell CR, Toth G, Fok MC, et al. 2020. A case study on the origin of nearEarth plasma. J Geophys Res 125(11): e2020JA028205. https://doi.org/10.1029/2020JA028205. [Google Scholar]
 Gloeckler G, Hamilton DC. 1987. AMPTE ion composition results. Phys Scripta T18: 73–84. https://doi.org/10.1088/00318949/1987/T18/009. [Google Scholar]
 Godunov SK. 1959. A difference scheme for numerical computation of discontinuous solutions of hydrodynamic equations. Mat Sb 47(3): 271–306 (in Russian). [Google Scholar]
 Gombosi T, De Zeeuw DL, Powell K, Ridley A, Sokolov I, Stout Q, Toth G. 2003. Adaptive mesh refinement MHD for global space weather simulations. In: Space Plasma Simulation, Vol. 615 in Lecture Notes in Physics. Büchner J, Dum CT, Scholer M, (Eds.) Springer, BerlinHeidelbergNew York. pp. 251. https://doi.org/10.1007/3540270396_36. [Google Scholar]
 Gombosi T, Powell K, van Leer B. 2000. Comment on “Modeling the magnetosphere for northward IMF: Effects of electrical resistivity” by. J Raeder J Geophys Res 105(A6): 13141–13147. https://doi.org/10.1029/1999ja000342. [Google Scholar]
 Gombosi T, Rasmussen C. 1991. Transport of gyrationdominated space plasmas of thermal origin I – Generalized transport equations. J Geophys Res 96: 7759–7778. https://doi.org/10.1029/91JA00012. [Google Scholar]
 Gombosi TI. 1994. Gaskinetic Theory. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/CBO9780511524943. [Google Scholar]
 Gombosi TI. 1998. Physics of the space environment. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/CBO9780511529474. [Google Scholar]
 Gombosi TI. 2015. Physics of cometary magnetospheres. Geophys Monograph 207: 169–188. https://doi.org/10.1002/9781118842324.ch10. [Google Scholar]
 Gombosi TI, Baker DN, Balogh A, Erickson PJ, Huba JD, Lanzerotti LJ. 2017. Anthropogenic Space Weather. Space Sci Rev 1–55: https://doi.org/10.1007/s1121401703575. [Google Scholar]
 Gombosi TI, De Zeeuw DL, Häberli RM, Powell KG. 1996. Threedimensional multiscale MHD model of cometary plasma environments. J Geophys Res 101: 15233–15252. https://doi.org/10.1029/96JA01075. [Google Scholar]
 Gombosi TI, Hansen KC. 2005. Saturn’s variable magnetosphere. Science 307: 1224–1226. https://doi.org/10.1126/science.1108226. [Google Scholar]
 Gombosi TI, Hansen KC, De Zeeuw DL, Combi MR, Powell KG. 1999. MHD simulation of comets: the plasma environment of comet HaleBopp. Earth Moon Planet. 79(1/3): 179–207. https://doi.org/10.1023/A:1006289418660. [Google Scholar]
 Gombosi TI, Ingersoll AP. 2010. Saturn: Atmosphere, ionosphere, and magnetosphere. Science 327(5972): 1476–1479. https://doi.org/10.1126/science.1179119. [Google Scholar]
 Gombosi TI, Nagy A. 1989. Timedependent modeling of field aligned currentgenerated ion transients in the polar wind. J Geophys Res 94: 359–369. https://doi.org/10.1029/ja094ia01p00359. [Google Scholar]
 Gombosi TI, Toth G, De Zeeuw DL, Hansen KC, Kabin K, Powell KG. 2002. Semirelativistic magnetohydrodynamics and physicsbased convergence acceleration. J Comput Phys 177: 176–205. https://doi.org/10.1006/jcph.2002.7009. [Google Scholar]
 Gombosi TI, van der Holst B, Manchester WB, Sokolov IV. 2018. Extended MHD modeling of the steady solar corona and the solar wind. Liv Rev Sol Phys 15(1): 4. https://doi.org/10.1007/s4111601800144. [Google Scholar]
 Grad H. 1949. On the kinetic theory of rarefied gases. Commun Pure Appl Math 2: 331–407. https://doi.org/10.1002/cpa.3160020403. [Google Scholar]
 Häberli RM, Gombosi TI, De Zeuuw DL, Combi MR, Powell KG. 1997. Modeling of cometary Xrays caused by solar wind minor ions. Science 276: 939–942. https://doi.org/10.1126/science.276.5314.939. [NASA ADS] [CrossRef] [Google Scholar]
 Haiducek JD, Welling DT, Ganushkina NY, Morley SK, Ozturk DS. 2017. SWMF global magnetosphere simulations of January 2005: Geomagnetic indices and crosspolar cap potential. Space Weather 15(12): 1567–1587. https://doi.org/10.1002/2017SW001695. [CrossRef] [Google Scholar]
 Haiducek JD, Welling DT, Morley SK, Ganushkina NY, Chu X. 2020. Using multiple signatures to improve accuracy of substorm identification. J Geophys Res 125(4): e2019JA027559. https://doi.org/10.1029/2019JA027559. [Google Scholar]
 Hansen KC, Bagdonat T, Motschmann U, Alexander C, Combi MR, Cravens TE, Gombosi TI, Jia Y, Robertson IP. 2007. The plasma environment of comet 67P/ChuryumovGerasimenko throughout the Rosetta main mission. Space Sci Rev 128: 133–166. https://doi.org/10.1007/s1121400691426. [Google Scholar]
 Hansen KC, Gombosi TI, DeZeeuw DL, Groth CPT, Powell KG. 2000. A 3D Global MHD Simulation of Saturn’s Magnetosphere. Adv Space Res 26(10): 1681–1690. https://doi.org/10.1016/s02731177(00)000788. [Google Scholar]
 Hansen KC, Ridley AJ, Gombosi TI, Hospodarsky G. 2005. Global simulations of Saturn’s magnetosphere at the time of Cassini approach. Geophys Res Lett 32(20): L20S06. https://doi.org/10.1029/2005GL022835. [Google Scholar]
 Harel M, Wolf RA, Reiff PH, Spiro RW, Burke WJ, Rich FJ, Smiddy M. 1981. Quantitative simulation of a magnetospheric substorm 1, Model logic and overview. J Geophys Res 86: 2217–2241. https://doi.org/10.1029/ja086ia04p02217. [Google Scholar]
 Harris C, Jia X, Slavin J, Tóth G, Huang Z, Rubin M. 2021. Multifluid MHD simulations of Europas plasma interaction under different magnetospheric conditions. J Geophys Res e2020JA028888. https://doi.org/10.1029/2020JA028888. [Google Scholar]
 Hayashi K. 2013. An MHD simulation model of timedependent global solar corona with temporally varying solarsurface magnetic field maps. J Geophys Res 118: 6889–6906. https://doi.org/10.1002/2013JA018991. [Google Scholar]
 Hedin AE. 1987. MSIS86 thermospheric model. J Geophys Res 92: 4649. https://doi.org/10.1029/ja092ia05p04649. [NASA ADS] [CrossRef] [Google Scholar]
 Hedin AE. 1991. Extension of the MSIS Thermosphere Model into the middle and lower atmosphere. J Geophys Res 96(A2): 1159–1172. https://doi.org/10.1029/90JA02125. [NASA ADS] [CrossRef] [Google Scholar]
 Heinemann M, Wolf RA. 2001. Relationships of models of the inner magnetosphere to the Rice Convection Model. J Geophys Res 106(A8): 15545–15554. https://doi.org/10.1029/2000ja000389. [Google Scholar]
 Herrero JL, Lucio F, Carmona P. 2011. Web services and web components. In: 2011 7th International Conference on Next Generation Web Services Practices. Salamanca, Spain. pp. 164–169. https://doi.org/10.1109/NWeSP.2011.6088171. [Google Scholar]
 Hill C, DeLuca C, Balaji V, Suarez M, da Silva A, the ESMF Joint Specification Team. 2004. The architecture of the Earth system modeling framework. Comput Sci Eng 6(1): 18–28. https://doi.org/10.1109/MCISE.2004.1255817. [Google Scholar]
 Hochreiter S, Schmidhuber J. 1997. Long shortterm memory. Neural Comput 9(8): 1735–1780. https://doi.org/10.1162/neco.1997.9.8.1735. [Google Scholar]
 Hollweg JV. 1986. Transition region, corona, and solar wind in coronal holes. J Geophys Res 91: 4111–4125. https://doi.org/10.1029/JA091iA04p04111. [NASA ADS] [CrossRef] [Google Scholar]
 Howard RA, Moses JD, Vourlidas A, Newmark JS, Socker DG, et al. 2008. Sun Earth Connection Coronal and Heliospheric Investigation (SECCHI). Space Sci Rev 136: 67–115. https://doi.org/10.1007/s1121400893414. [Google Scholar]
 Huang Z, Toth G, Gombosi T, Bieler A, Combi M, et al. 2016a. A possible mechanism for the formation of magnetic field dropouts in the coma of 67P/Churyumov–Gerasimenko. MNRAS 462(1): S468–S475. https://doi.org/10.1093/mnras/stw3118 [Google Scholar]
 Huang Z, Toth G, Gombosi TI, Jia X, Rubin M, et al. 2016b. Fourfluid MHD simulations of the plasma and neutral gas environment of comet 67P/ChuryumovGerasimenko near perihelion. J Geophys Res 121(5): 4247–4268. https://doi.org/10.1002/2015JA022333. [Google Scholar]
 Huang Z, Toth G, van der Holst B, Chen Y, Gombosi TI. 2019. A sixmoment multifluid plasma model. J Comput Phys 387: 134–153. https://doi.org/10.1016/j.jcp.2019.02.023. [Google Scholar]
 Ilie R, Liemohn MW, Toth G, Ganushkina NY, Daldorff LKS. 2015. Assessing the role of oxygen on ring current formation and evolution through numerical experiments. J Geophys Res 120(6): 4656–4668. https://doi.org/10.1002/2015ja021157. [Google Scholar]
 Ilie R, Liemohn MW, Toth G, Skoug RM. 2012. Kinetic model of the inner magnetosphere with arbitrary magnetic field. J Geophys Res 117(A4): A04208. https://doi.org/10.1029/2011JA017189. [Google Scholar]
 Ilie R, Skoug R, Funsten H, Liemohn M, Bailey J, Gruntman M. 2013. The impact of geocoronal density on ring current development. J Atmos SolarTerr Phys 99: 92–103. https://doi.org/10.1016/j.jastp.2012.03.010. [Google Scholar]
 Illing RME, Hundhausen AJ. 1985. Observation of a coronal transient from 1.2 to 6 solar radii. J Geophys Res 90: 275–282. https://doi.org/10.1029/JA090iA01p00275. [Google Scholar]
 Jackson JD. 1975. Classical Eelectrodynamics. John Wiley and Sons, New York. [Google Scholar]
 Jacques SA. 1977. Momentum and energy transport by waves in the solar atmosphere and solar wind. Astrophys J 215: 942–951. https://doi.org/10.1086/155430. [Google Scholar]
 Janhunen P. 1996. GUMICS3: A global ionospheremagnetosphere coupling simulation with high ionospheric resolution. In: Guyenne TD, Hilgers A, (Eds.). Proceedings of the ESA 1996 Symposium on Environment Modelling for SpaceBased Applications, vol. 392 of ESA Special Publication, pp. 233–239. ESA SP392. [Google Scholar]
 Jia X, Hansen KC, Gombosi TI, Kivelson MG, Toth G, DeZeeuw DL, Ridley AJ. 2012a. Magnetospheric configuration and dynamics of Saturn’s magnetosphere: A global MHD simulation. J Geophys Res 117: A05225. https://doi.org/10.1029/2012JA017575. [Google Scholar]
 Jia X, Kivelson MG, Gombosi TI. 2012b. Driving Saturn’s magnetospheric periodicities from the upper atmosphere/ionosphere. J Geophys Res 117: A04215. https://doi.org/10.1029/2011JA017367. [Google Scholar]
 Jia X, Kivelson MG, Khurana KK, Kurth WS. 2018. Evidence of a plume on Europa from Galileo magnetic and plasma wave signatures. Nature Astron 2: 459–464. https://doi.org/10.1038/s415500180450z. [NASA ADS] [CrossRef] [Google Scholar]
 Jia X, Slavin JA, Gombosi TI, Daldorff LKS, Toth G, van der Holst B. 2015. Global MHD simulations of Mercury’s magnetosphere with coupled planetary interior: Induction effect of the planetary conducting core on the global interaction. J Geophys Res 120(6): 4763–4775. https://doi.org/10.1002/2015JA021143. [NASA ADS] [CrossRef] [Google Scholar]
 Jia X, Slavin JA, Poh G, DiBraccio GA, Toth G, Chen Y, Raines JM, Gombosi TI. 2019. MESSENGER observations and global simulations of highly compressed magnetosphere events at Mercury. J Geophys Res 124(1): 229–247. https://doi.org/10.1029/2018JA026166. [CrossRef] [Google Scholar]
 Jia YD, Russell CT, Khurana KK, Ma YJ, Kurth W, Gombosi TI. 2010a. Interaction of Saturn’s magnetosphere and its moons: 3. Time variation of the Enceladus plume. J Geophys Res 115: A12243. https://doi.org/10.1029/2010JA015534. [Google Scholar]
 Jia YD, Russell CT, Khurana KK, Ma YJ, Najib D, Gombosi TI. 2010b. Interaction of Saturn’s magnetosphere and its moons: 2. Shape of the Enceladus plume. J Geophys Res 115: A04,215. https://doi.org/10.1029/2009JA014873. [Google Scholar]
 Jia YD, Russell CT, Khurana KK, Toth G, Leisner JS, Gombosi TI. 2010c. Interaction of Saturn’s magnetosphere and its moons: 1. Interaction between corotating plasma and standard obstacles. J Geophys Res 115: A04214. https://doi.org/10.1029/2009JA014630. [Google Scholar]
 Jiao Z, Sun H, Wang X, Manchester W, Gombosi T, Hero A, Chen Y. 2020. Solar flare intensity prediction with machine learning models. Space Weather 18(7): e2020SW002440. https://doi.org/10.1029/2020SW002440. [CrossRef] [Google Scholar]
 Jin M, Manchester WB, van der Holst B, Sokolov I, Toth G, Mullinix RE, Taktakishvili A, Chulaki A, Gombosi TI. 2017a. Dataconstrained Coronal Mass Ejections in a Global Magnetohydrodynamics Model. Astrophys J 834(2): 173. https://doi.org/10.3847/15384357/834/2/173. [NASA ADS] [CrossRef] [Google Scholar]
 Jin M, Manchester WB, van der Holst B, Sokolov I, Toth G, Vourlidas A, de Koning CA, Gombosi TI. 2017b. Chromosphere to 1 AU Simulation of the 2011 March 7th Event: A Comprehensive Study of Coronal Mass Ejection Propagation. Astrophys J 834(2): 172. https://doi.org/10.3847/15384357/834/2/172. [NASA ADS] [CrossRef] [Google Scholar]
 Jin M, Petrosian V, Liu W, Nitta NV, Omodei N, et al. 2018. Probing the puzzle of behindthelimb γray flares: DATAdriven simulations of magnetic connectivity and CMEdriven shock evolution. Astrophys J 867(2): 122. https://doi.org/10.3847/15384357/aae1fd. [CrossRef] [Google Scholar]
 Jordanova V, Delzanno G, Henderson M, Godinez H, Jeffery C, et al. 2018. Specification of the nearEarth space environment with SHIELDS. J Atmos SolarTerr Phys 177: 148–159. https://doi.org/10.1016/j.jastp.2017.11.006. [CrossRef] [Google Scholar]
 Jordanova VK, Kozyra JU, Khazanov GV, Nagy AF, Rasmussen CE, Fok MC. 1994. A bounceaveraged kinetic model of the ring current ion population. Geophys Res Lett 21: 2785–2788. https://doi.org/10.1029/94gl02695. [CrossRef] [Google Scholar]
 Jordanova VK, Miyoshi YS, Zaharia S, Thomsen MF, Reeves GD, Evans DS, Mouikis CG, Fennell JF. 2006. Kinetic simulations of ring current evolution during the Geospace Environment Modeling challenge events. J Geophys Res 111(A11): A11S10. https://doi.org/10.1029/2006JA011644. [CrossRef] [Google Scholar]
 Jordanova VK, Zaharia S, Welling DT. 2010. Comparative study of ring current development using empirical, dipolar, and selfconsistent magnetic field simulations. J Geophys Res 115(A12): A00J11. https://doi.org/10.1029/2010JA015671. [Google Scholar]
 Kabin K, Combi MR, Gombosi TI, De Zeeuw DL, Hansen KC, Powell KG. 2001. Io’s magnetospheric interaction: an MHD model with daynight asymmetry. Planet Space Sci 49: 337–344. https://doi.org/10.1016/S00320633(00)001550. [CrossRef] [Google Scholar]
 Kabin K, Combi MR, Gombosi TI, Nagy AF, De Zeeuw DL, Powell KG. 1999a. On Europa’s magnetospheric interaction: An MHD simulation of the E4 flyby. J Geophys Res 104(A9): 19983–19992. https://doi.org/10.1029/1999JA900263. [CrossRef] [Google Scholar]
 Kabin K, Gombosi TI, De Zeeuw DL, Powell KG. 2000a. Interaction of Mercury with the solar wind. Icarus 143: 397–406. https://doi.org/10.1006/icar.1999.6252. [NASA ADS] [CrossRef] [Google Scholar]
 Kabin K, Gombosi TI, De Zeeuw DL, Powell KG, Israelevich PL. 1999b. Interaction of the Saturnian magnetosphere with Titan: Results of a 3D MHD simulation. J Geophys Res 104(A2): 2451–2458. https://doi.org/10.1029/1998JA900080. [CrossRef] [Google Scholar]
 Kabin K, Heimpel MH, Rankin R, Aurnou JM, GomezPerez N, Paral J, Gombosi TI, Zurbuchen TH, Koehn PL, DeZeeuw DL. 2008. Global MHD modeling of Mercury’s magnetosphere with applications to the MESSENGER mission and dynamo theory. Icarus 195(1): 1–15. https://doi.org/10.1016/j.icarus.2007.11.028. [NASA ADS] [CrossRef] [Google Scholar]
 Kabin K, Israelevich PL, Ershkovich AI, Neubauer FM, Gombosi TI, De Zeeuw DL, Powell KG. 2000b. Titan’s magnetic wake: Atmospheric or magnetospheric interaction. J Geophys Res 105: 10761–10770. https://doi.org/10.1029/2000JA900012. [CrossRef] [Google Scholar]
 Keppens R, Teunissen J, Xia C, Porth O. 2021. MPIAMRVAC: A parallel, gridadaptive PDE toolkit. Comput Math Appl 81: 316–333. https://doi.org/10.1016/j.camwa.2020.03.023. [CrossRef] [Google Scholar]
 Khazanov G, Neubert T, Gefan G. 1994. A unified theory of ionosphereplasmasphere transport of suprathermal electrons. IEEE Trans Plasma Sci 22(2): 187–198. https://doi.org/10.1109/27.279022. [CrossRef] [Google Scholar]
 Kóta J. 2000. Diffusion of energetic particles in focusing fields. J Geophys Res 105: 2403–2412. https://doi.org/10.1029/1999ja900469. [NASA ADS] [CrossRef] [Google Scholar]
 Kóta J, Manchester WB, Jokipii JR, de Zeeuw DL, Gombosi TI. 2005. Simulation of SEP Acceleration and Transport at CMEdriven Shocks. In: AIP Conference Proceedings, Vol. 781, AIP, pp. 201–206. https://doi.org/10.1063/1.2032697. [CrossRef] [Google Scholar]
 Kramers HA. 1926. Wellenmechanik und halbzahlige Quantisierung. Zeitschrift für Physik 39(10): 828–840. https://doi.org/10.1007/BF01451751. [CrossRef] [Google Scholar]
 Kuznetsova MM, Hesse M, Rastaetter L, Taktakishvili A, Toth G, De Zeeuw DL, Ridley A, Gombosi TI. 2007. Multiscale modeling of magnetospheric reconnection. J Geophys Res 112(A10): A10210. https://doi.org/10.1029/2007JA012316. [CrossRef] [Google Scholar]
 Kuznetsova MM, Sibeck DG, Hesse M, Wang Y, Rastaetter L, Toth G, Ridley A. 2009. Cavities of weak magnetic field strength in the wake of FTEs: Results from global magnetospheric MHD simulations. Geophys Res Lett 36(10): L10104. https://doi.org/10.1029/2009GL037489. [CrossRef] [Google Scholar]
 Lapenta G. 2017. Exactly energy conserving semiimplicit particle in cell formulation. J Comput Phys 334: 349–366. https://doi.org/10.1016/j.jcp.2017.01.002. [CrossRef] [Google Scholar]
 Lapenta G, Brackbill JU, Ricci P. 2006. Kinetic approach to microscopicmacroscopic coupling in space and laboratory plasmas. Phys. Plasmas 13(055): 904. https://doi.org/10.1063/1.2173623. [CrossRef] [Google Scholar]
 LeBoeuf JN, Tajima T, Kennel CF, Dawson JM. 1981. Global simulations of the threedimensional magnetosphere. Geophys Res Lett 8: 257–260. https://doi.org/10.1029/gl008i003p00257. [CrossRef] [Google Scholar]
 Leka K, Barnes G. 2018. Solar flare forecasting: Present methods and challenges. In: Extreme Events in Geospace. Buzulukova N, (Ed.), Elsevier, Amsterdam, The Netherlands. pp. 65–98. https://doi.org/10.1016/B9780128127001.000030. [CrossRef] [Google Scholar]
 Lemen JR, Title AM, Akin DJ, Boerner PF, Chou C, et al. 2012. The Atmospheric Imaging Assembly (AIA) on the Solar Dynamics Observatory (SDO). Sol Phys 275(1–2): 17–40. https://doi.org/10.1007/s1120701197768. [Google Scholar]
 Lennartsson W, Sharp RD, Shelley EG, Johnson RG, Balsiger H. 1981. Ion composition and energy distribution during 10 magnetic storms. J Geophys Res 86: 4628–4638. https://doi.org/10.1029/JA086iA06p04628. [CrossRef] [Google Scholar]
 Levermore C. 1996. Moment closure hierarchies for kinetic theories. J Stat Phys 83: 1021–1065. https://doi.org/10.1007/BF02179552. [CrossRef] [MathSciNet] [Google Scholar]
 Liemohn M, Ganushkina NY, De Zeeuw DL, Rastaetter L, Kuznetsova M, Welling DT, Toth G, Ilie R, Gombosi TI, van der Holst B. 2018. Realtime SWMF at CCMC: assessing the Dst output from continuous operational simulations. Space Weather 16(10): 1583–1603. https://doi.org/10.1029/2018SW001953. [CrossRef] [Google Scholar]
 Liemohn M, Ridley A, Gallagher D, Ober D, Kozyra J. 2004. Dependence of plasmaspheric morphology on the electric field description during the April 17, 2002 magnetic storm. J Geophys Res 109: A03209. https://doi.org/10.1029/2003JA010304. [Google Scholar]
 Liemohn MW. 2006. Introduction to special section on Results of the National Science Foundation Geospace Environment Modeling Inner Magnetosphere/Storms Assessment Challenge. J Geophys Res 111: A11S01. https://doi.org/10.1029/2006JA011970. [Google Scholar]
 Liemohn MW, Jazowski M. 2008. Ring current simulations of the 90 intense storms during solar cycle 23. J Geophys Res 113: A00A17. https://doi.org/10.1029/2008JA013466. [CrossRef] [Google Scholar]
 Liemohn MW, Khazanov GV, Kozyra JU. 1998. Banded electron structure formation in the inner magnetosphere. Geophys Res Lett 25: 877–880. https://doi.org/10.1029/98gl00411. [CrossRef] [Google Scholar]
 Liemohn MW, Kozyra JU, Clauer CR, Khazanov GV, Thomsen MF. 2002. Adiabatic energization in the ring current and its relation to other source a nd loss terms. J Geophys Res 107(A4): https://doi.org/10.1029/2001JA000243. [Google Scholar]
 Liemohn MW, Kozyra JU, Jordanova VK, Khazanov GV, Thomsen MF, Cayton TE. 1999. Analysis of early phase ring current recovery mechanisms during geomagnetic storms. Geophys Res Lett 26(18): 2845–2848. https://doi.org/10.1029/1999gl900611. [CrossRef] [Google Scholar]
 Liemohn MW, Ridley AJ, Brandt PC, Gallagher DL, Kozyra JU, Ober DM, Mitchell DG, Roelof EC, DeMajistre R. 2005. Parametric analysis of nightside conductance effects on inner magnetospheric dynamics for the 17 April 2002 storm. J Geophys Res 110: A12S22. https://doi.org/10.1029/2005JA011109. [CrossRef] [Google Scholar]
 Liemohn MW, Ridley AJ, Kozyra JU, Gallagher DL, Thomsen MF, Henderson MG, Denton MH, Brandt PC, Goldstein J. 2006. Analyzing electric field morphology through datamodel comparisons of the GEM IM/S Assessment Challenge events. J Geophys Res 111: A11S11. https://doi.org/10.1029/2006JA011700. [Google Scholar]
 Lilensten J, Coates AJ, Dehant V, de Wit TD, Horne RB, Leblanc F, Luhmann J, Woodfield E, Barthélemy M. 2014. What characterizes planetary space weather? Astron Astrophys Rev 22(1): https://doi.org/10.1007/s0015901400796. [CrossRef] [Google Scholar]
 Linker JA, Mikić 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 104: 9809–9830. https://doi.org/10.1029/1998ja900159. [NASA ADS] [CrossRef] [Google Scholar]
 Linker JA, Mikić Z, Schnack DD. 1994. Modeling coronal evolution. European Space Agency, Estes Park, CO. pp. 249–252. [Google Scholar]
 Lionello R, Linker JA, Mikič 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]
 Liu J, Ilie R. 2021. The effects of inductive electric field on the spatial and temporal evolution of inner magnetospheric ring current. J Geophys Res Space Phys 126: e2020JA028554. https://doi.org/10.1029/2020JA028554. [Google Scholar]
 Liu L, Zou S, Yao Y, Wang Z. 2020. Forecasting global ionospheric tec using deep learning approach. Space Weather 18(11): e2020SW002501. https://doi.org/10.1029/2020SW002501. [Google Scholar]
 Liu Y, Nagy A, Kabin K, Combi M, DeZeeuw D, Gombosi T, Powell K. 2000. Twospecies, 3D, MHD simulation of Europa’s interaction with Jupiter’s magnetosphere. Geophys Res Lett 27(12): 1791–1794. https://doi.org/10.1029/1999GL003734. [CrossRef] [Google Scholar]
 Liu Y, Nagy AF, Groth CPT, De Zeeuw DL, Gombosi TI, Powell KG. 1999. 3D Multifluid MHD studies of the solar wind interaction with Mars. Geophys Res Lett 26(17): 2689–2692. https://doi.org/10.1029/1999GL900584. [CrossRef] [Google Scholar]
 Lugaz N, Manchester WB, Gombosi TI. 2005a. The evolution of CME density structures. Astrophys J 627: 1019–1030. https://doi.org/10.1086/430465. [NASA ADS] [CrossRef] [Google Scholar]
 Lugaz N, Manchester WB, Gombosi TI. 2005b. Numerial Simulation of the Interaction of Two Coronal Mass Ejections from Sun to Earth. Astrophys J 634: 651–662. https://doi.org/10.1086/491782. [NASA ADS] [CrossRef] [Google Scholar]
 Lugaz N, Manchester WB, Roussev II, Gombosi TI. 2008. Observational evidence of CMEs Interacting in the Inner Heliosphere as Inferred from MHD Simulations. J Atmos SolarTerr Phys 70(2–4): 598–604. https://doi.org/10.1016/j.jastp.2007.08.033. [NASA ADS] [CrossRef] [Google Scholar]
 Lugaz N, Vourlidas A, Roussev II, Morgan H. 2009. SolarTerrestrial Simulation in the STEREO Era: The 24–25 January 2007 Eruptions. Sol Phys 256: 269–284. https://doi.org/10.1007/s1120700993394. [CrossRef] [Google Scholar]
 Luhmann JG, Solomon SC, Linker JA, Lyon JG, Mikic Z, Odstrcil D, Wang W, Wiltberger M. 2004. Coupled model simulation of a SuntoEarth space weather event. J Atmos SolarTerr Phys 66: 1243–1256. https://doi.org/10.1016/j.jastp.2004.04.005. [CrossRef] [Google Scholar]
 Lummerzheim D, Lilensten J. 1994. Electron transport and energy degradation in the ionosphere: evaluation of the numerical solution, comparison with laboratory experiments and auroral observations. Ann Geophys 12(10/11): 1039–1051. https://doi.org/10.1007/s0058599410397. [NASA ADS] [CrossRef] [Google Scholar]
 Lundstedt H, Wintoft P. 1994. Prediction of geomagnetic storms from solar wind data with the use of a neural network. Ann Geophys 12(1): 19–24. https://doi.org/10.1007/s0058599400192. [CrossRef] [Google Scholar]
 Lyon JG, Fedder J, Huba J. 1986. The effect of different resistivity models on magnetotail dynamics. J Geophys Res 91(A7): 8057–8064. https://doi.org/10.1029/ja091ia07p08057. [CrossRef] [Google Scholar]
 Lyon JG, Fedder JA, Mobarry CM. 2004. The Lyon–Fedder–Mobarry (LFM) global MHD magnetospheric simulation code. J Atmos SolarTerr Phys 66(15–16): 1333–1350. https://doi.org/10.1016/j.jastp.2004.03.020. [CrossRef] [Google Scholar]
 Ma Y, Nagy AF, Sokolov IV, Hansen KC. 2004. Threedimensional, multispecies, high spatial resolution MHD studies of the solar wind interaction with Mars. J Geophys Res 109: A07211. https://doi.org/10.1029/2003JA010367. [Google Scholar]
 Ma Y, Fang X, Halekas JS, Xu S, Russell CT, et al. 2018a. The Impact and Solar Wind Proxy of the 2017 September ICME Event at Mars. Geophys Res Lett 45(15): 7248–7256. https://doi.org/10.1029/2018GL077707. [CrossRef] [Google Scholar]
 Ma Y, Nagy AF, Hansen KC, De Zeeuw DL, Gombosi TI, Powell K. 2002. Threedimensional multispecies MHD studies of the solar wind interaction with Mars in the presence of crustal fields. J Geophys Res 107(A10): 1282. https://doi.org/10.1029/2002JA009293. [CrossRef] [Google Scholar]
 Ma Y, Nagy AF, Russell C, Strangeway RJ, Wei HY, Toth G. 2013. A Global MultiSpecies SingleFluid MHD Study of the Plasma Interaction around Venus. J Geophys Res 118: 321. https://doi.org/10.1029/2012JA018265. [CrossRef] [Google Scholar]
 Ma Y, Russell CT, Toth G, Chen Y, Nagy AF, et al. 2018b. Reconnection in the Martian Magnetotail: HallMHD with embedded particleincell simulations. J Geophys Res 123(5): 3742–3763. https://doi.org/10.1029/2017JA024729. [CrossRef] [Google Scholar]
 Ma YJ, Nagy AF, Toth G, Cravens TE, Russell CT, et al. 2007. 3D global multispecies HallMHD simulation of the Cassini T9 flyby. Geophys Res Lett 34(24): L24S10. https://doi.org/10.1029/2007GL031627. [Google Scholar]
 Manchester W, Gombosi TI, De Zeeuw DL, Sokolov IV, Roussev II, Powell KG, Kóta J, Toth G, Zurbuchen TH. 2005. Coronal Mass Ejection shock and sheath structures relevant to particle acceleration. Astrophys J 622: 1225–1239. https://doi.org/10.1086/427768. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Gombosi TI, De Zeeuw DL, Fan Y. 2004a. Eruption of a Buoyantly Emerging Magnetic Flux Rope. Astrophys J 610: 588. https://doi.org/10.1086/421516. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Gombosi TI, Roussev II, Zeeuw DLD, Sokolov IV, Powell KG, Toth G, Opher M. 2004b. Threedimensional MHD simulation of a fluxrope driven CME. J Geophys Res 109(A1): 1102–1119. https://doi.org/10.1029/2002JA009672. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Kozyra JU, Lepri ST, Lavraud B. 2014a. Simulation of magnetic cloud erosion during propagation. J Geophys Res 119: 1–16. https://doi.org/10.1002/2014JA019882. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Ridley AJ, Gombosi TI, De Zeeuw D. 2006. Modeling the SunEarth Propagation of a Very Fast CME. Adv Space Res 38(2): 253–262. https://doi.org/10.1016/j.asr.2005.09.044. [CrossRef] [Google Scholar]
 Manchester WB, van der Holst B. 2017. The interaction of coronal mass ejections with alfvénic Turbulence. In: Journal of Physics Conference Series. Vol. 900of Journal of Physics Conference Series, 012015. https://doi.org/10.1088/17426596/900/1/012015. [CrossRef] [Google Scholar]
 Manchester WB, van der Holst B, Lavraud B. 2014b. Flux rope evolution in interplanetary coronal mass ejections: the 13 May 2005 event. Plasma Phys Controlled Fusion 56(6): 1–11. https://doi.org/10.1088/07413335/56/6/064006. [CrossRef] [Google Scholar]
 Manchester WB, van der Holst B, Toth G, Gombosi TI. 2012. The coupled evolution of electrons and ions in coronal mass ejectiondriven shocks. Astrophys J 756(1):81. https://doi.org/10.1088/0004637X/756/1/81. [NASA ADS] [CrossRef] [Google Scholar]
 Manchester WB, Vourlidas A, Toth G, Lugaz N, Roussev II, Sokolov IV, Gombosi TI, Zeeuw DLD, 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]
 Markidis S, Lapenta G, RizwanUddin. 2010. Multiscale simulations of plasma with iPIC3D. Math Comput Simul 80(7): 1509–1519. https://doi.org/10.1016/j.matcom.2009.08.038. [CrossRef] [Google Scholar]
 Mayaud PN. 1980. Derivation, Meaning, and Use of Geomagnetic Indices. In: Vol. 22 of Geophysical Monograph Series, American Geophysical Union, Washington, DC. https://doi.org/10.1029/GM022. [Google Scholar]
 Meng X, Toth G, Glocer A, Fok MC, Gombosi TI. 2013. Pressure anisotropy in global magnetospheric simulations: Coupling with ring current models. J Geophys Res 118: 5639. https://doi.org/10.1002/jgra.50539. [CrossRef] [Google Scholar]
 Meng X, Toth G, Gombosi TI. 2012. Classical and semirelativistic magnetohydrodynamics with anisotropic ion pressure. J. Comput. Phys. 231: 3610–3622. https://doi.org/10.1016/j.jcp.2011.12.042. [CrossRef] [Google Scholar]
 Menvielle M, Iyemori T, Marchaudon A, Nosé M. 2010. Geomagnetic indices. Geomagnetic observations and models. Springer, Netherlands. pp. 183–228. https://doi.org/10.1007/9789048198580_8. [Google Scholar]
 Merkin VG. 2011. Effects of ionospheric O + on the magnetopause boundary wave activity. AIP Conf Proc 1320: 208–212. https://doi.org/10.1063/1.3544326. [CrossRef] [Google Scholar]
 Mukhopadhyay A, Jia X, Welling D, Liemohn M. 2021. Global magnetohydrodynamic simulations: Performance quantification of magnetopause distances and convection potential prediction. Earth Space Sci Open Archive, pp. 1–22, in press, https://doi.org/10.3389/fspas.2021.637197. [Google Scholar]
 Mukhopadhyay A, Welling DT, Liemohn MW, Ridley AJ, Chakraborty S, Anderson BJ. 2020. Conductance model for extreme events: impact of auroral conductance on space weather forecasts. Space Weather 18(11): e2020SW002551. https://doi.org/10.1029/2020SW002551. [CrossRef] [Google Scholar]
 Nagy AF, Banks PM. 1970. Photoelectron fluxes in the ionosphere. J Geophys Res 75: 6260. [CrossRef] [Google Scholar]
 Nagy AF, Liu Y, Hansen KC, Kabin K, Gombosi T, Combi M, De Zeeuw DL, Powell KG, Kliore A. 2001. The interaction between the magnetosphere of Saturn and Titan’s ionosphere. J Geophys Res 106: 6151–6160. https://doi.org/10.1029/2000JA000183. [CrossRef] [Google Scholar]
 Nanbu K, Yonemura S. 1998. Weighted Particles in Coulomb Collision Simulations Based on the Theory of a Cumulative Scattering Angle. J Comput Phys 145(2): 639–654. https://doi.org/10.1006/jcph.1998.6049. [CrossRef] [Google Scholar]
 Ober DM, Horwitz JL, Gallagher DL. 1997. Formation of density troughs embedded in the outer plasmasphere by subauroral ion drift events. J Geophys Res 102(A7): 14595–14602. https://doi.org/10.1029/97JA01046. [CrossRef] [Google Scholar]
 Odstrčil D. 2003. Modeling 3D solar wind structures. Adv Space Res 32(4): 497–506. [NASA ADS] [CrossRef] [Google Scholar]
 Odstrčil D, Pizzo VJ. 2009. Numerical Heliospheric Simulations as Assisting Tool for Interpretation of Observations by STEREO Heliospheric Imagers. Sol Phys 259: 297–309. https://doi.org/10.1007/s112070099449z. [NASA ADS] [CrossRef] [Google Scholar]
 Ogino T. 1986. A threedimensional MHD simulation of the interaction of the solar wind with the Earth’s magnetosphere: The generation of fieldaligned currents. J Geophys Res 91: 6791–6806. https://doi.org/10.1029/ja091ia06p06791. [CrossRef] [Google Scholar]
 Opher M, Bibi FA, Toth G, Richardson JD, Izmodenov VV, Gombosi TI. 2009. A strong, highlytilted interstellar magnetic field near the Solar System. Nature 462: 1036–1038. https://doi.org/10.1038/nature08567. [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Opher M, Drake J, Swisdak M, Zieger B, Toth G. 2017. The twist of the draped interstellar magnetic field ahead of the heliopause: Aa magnetic reconnection driven rotational discontinuity. Astrophys J Lett 839(1): L12. https://doi.org/10.3847/20418213/aa692f. [NASA ADS] [CrossRef] [Google Scholar]
 Opher M, Drake JF, Zieger B, Swisdak M, Toth G. 2016. Magnetized jets driven by the Sun: The structure of the heliosphere revisited – Updates. Phys Plasmas 23(5): 056501. https://doi.org/10.1063/1.4943526. [NASA ADS] [CrossRef] [Google Scholar]
 Opher M, Stone EC, Gombosi TI. 2007. The orientation of the local interstellar magnetic field. Science 316(5826): 875–878. https://doi.org/10.1126/science.1139480. [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Oran R, van der Holst B, Landi E, Jin M, Sokolov IV, Gombosi TI. 2013. A global wavedriven magnetohydrodynamic solar model with a unified treatment of open and closed magnetic field topologies. Astrophys J 778: 176–195. https://doi.org/10.1088/0004637X/778/2/176. [NASA ADS] [CrossRef] [Google Scholar]
 Ozturk DS, Zou S, Ridley AJ, Slavin JA. 2018. Modeling study of the geospace system response to the solar wind dynamic pressure enhancement on 17 March 2015. J Geophys Res 123(4): 2974–2989. https://doi.org/10.1002/2017JA025099. [CrossRef] [Google Scholar]
 Parker EN. 1965. The passage of energetic charged particles through interplanetary space. Planet Space Sci 13: 9–49. https://doi.org/10.1016/00320633(65)901315. [NASA ADS] [CrossRef] [Google Scholar]
 Paszke A, Gross S, Massa F, Lerer A, Bradbury J, et al. 2019. PyTorch: An imperative style, highperformance deep learning library. Adv Neural Inform Process Syst 8026–8037. [Google Scholar]
 Perlongo NJ, Ridley AJ, Liemohn MW, Katus RM. 2017. The effect of ring current electron scattering rates on magnetosphereionosphere coupling. J Geophys Res 122(4): 4168–4189. https://doi.org/10.1002/2016ja023679. [CrossRef] [Google Scholar]
 Pillet VM, Hill F, Hammel H, de Wijn AG, Gosain S, et al. 2019. Synoptic studies of the sun as a key to understanding stellar astrospheres. Bull AAS 51(3): Retrieved from https://baas.aas.org/pub/2020n3i110. [Google Scholar]
 Plainaki C, Lilensten J, Radioti A, Andriopoulou M, Milillo A, et al. 2016. Planetary space weather: Scientific aspects and future perspectives. J Space Weather Space Clim 6: A31. https://doi.org/10.1051/swsc/2016024. [CrossRef] [Google Scholar]
 Poedts S, Kochanov A, Lani A, Scolini C, Verbeke C, et al. 2020. The virtual space weather modelling centre. J Space Weather Space Clim 10: 14. https://doi.org/10.1051/swsc/2020012. [CrossRef] [Google Scholar]
 Poh G, Slavin JA, Jia X, Raines JM, Imber SM, Sun WJ, Gershman DJ, DiBraccio GA, Genestreti KJ, Smith AW. 2017. Mercury’s crosstail current sheet: Structure, Xline location and stress balance. Geophys Res Lett 44(2): 678–686. https://doi.org/10.1002/2016GL071612. [NASA ADS] [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]
 Powell K, Roe P, Linde T, Gombosi T, 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]
 Powell KG. 1997. An approximate Riemann solver for magnetohydrodynamics. Upwind and HighResolution Schemes. Springer, Berlin Heidelberg. pp. 570–583. https://doi.org/10.1007/9783642605437_23. [CrossRef] [Google Scholar]
 Pulkkinen A, Rastätter L, Kuznetsova M, Singer H, Balch C, et al. 2013. Communitywide validation of geospace model ground magnetic field perturbation predictions to support model transition to operations. Space Weather 11(6): 369–385. https://doi.org/10.1002/swe.20056. [CrossRef] [Google Scholar]
 Pulkkinen TI, Baker DN, Pellinen RJ, Büchner J, Koskinen HEJ, Lopez RE, Dyson RL, Frank LA. 1992. Particle scattering and current sheet stability in the geomagnetic tail during the substorm growth phase. J. Geophys. Res. 97(A12): 19283–19297. https://doi.org/10.1029/92JA01189. [CrossRef] [Google Scholar]
 Pulkkinen TI, Ganushkina NY, Tanskanen EI, Kubyshkina M, Reeves GD, Thomsen MF, Russell CT, Singer HJ, Slavin JA, Gjerloev J. 2006. Magnetospheric current systems during stormtime sawtooth events. J Geophys Res 111(A11): xxx. https://doi.org/10.1029/2006JA011627. [CrossRef] [Google Scholar]
 Qin G, Zhang M, Dwyer JR. 2006. Effect of adiabatic cooling on the fitted parallel mean free path of solar energetic particles. J Geophys Res 111(A8): A08101. https://doi.org/10.1029/2005JA011512. [Google Scholar]
 Raeder J. 2000. Reply. J Geophys Res 105: 13149–13153. [CrossRef] [Google Scholar]
 Raeder J, Berchem J, AshourAbdalla M. 1996. The importance of small scale processes in global MHD simulations: Some numerical experiments. In: The Physics of Space Plasmas. Vol. 14, Chang T, Jasperse JR, (Eds.) MIT Cent. for Theoret. Geo/Cosmo Plasma Phys, Cambridge, Mass. pp. 403. [Google Scholar]
 Raeder J, Walker RJ, AshourAbdalla M. 1995. The structure of the distant geomagnetic tail during long periods of northward IMF. Geophys Res Lett 22(4): 349–352. https://doi.org/10.1029/94gl03380. [CrossRef] [Google Scholar]
 Rasmussen C, Guiter SM, Thomas SG. 1993. Twodimensional model of the plasmasphere: Refilling time constants. Planet Space Sci 41: 35. [CrossRef] [Google Scholar]
 Rastätter L, Kuznetsova MM, Glocer A, Welling D, Meng X, et al. 2013. Geospace environment modeling 2008–2009 challenge: Dst index. Space Weather 11(4): 187–205. https://doi.org/10.1002/swe.20036. [CrossRef] [Google Scholar]
 Rastätter L, Tóth G, Kuznetsova MM, Pulkkinen AA. 2014. CalcDeltaB: An efficient postprocessing tool to calculate groundlevel magnetic perturbations from global magnetosphere simulations. Space Weather 12(9): 553–565. https://doi.org/10.1002/2014SW001083. [CrossRef] [Google Scholar]
 Reames DV. 1999. Particle acceleration at the Sun and in the heliosphere. Space Sci Rev 90(3–4): 413–491. https://doi.org/10.1023/A:1005105831781. [NASA ADS] [CrossRef] [Google Scholar]
 Regoli LH, Dong C, Ma Y, Dubinin E, Manchester WB, Bougher SW, Welling DT. 2018. Multispecies and multifluid MHD approaches for the study of ionospheric escape at Mars. J Geophys Res 123(9): 7370–7383. https://doi.org/10.1029/2017JA025117. [CrossRef] [Google Scholar]
 Retterer JM, Chang T, Crew GB, Jasperse JR, Winningham JD. 1987. Monte Carlo modeling of ionospheric oxygen acceleration by cyclotron resonance with broadband electromagnetic turbulence. Phys. Rev. Lett 59(1): 148–151. https://doi.org/10.1103/PhysRevLett.59.148. [CrossRef] [Google Scholar]
 Ricci P, Lapenta G, Brackbill JU. 2002. GEM reconnection challenge: Implicit kinetic simulations with the physical mass ratio. Geophys Res Lett 29(23): 3–1–3–4. https://doi.org/10.1029/2002GL015314. [CrossRef] [Google Scholar]
 Ridley AJ, Deng Y, Toth G. 2006. The global ionospherethermosphere model. J Atmos SolarTerr Phys 68(8): 839–864. https://doi.org/10.1016/j.jastp.2006.01.008. [NASA ADS] [CrossRef] [Google Scholar]
 Ridley AJ, Dodger AM, Liemohn MW. 2014. Exploring the efficacy of different electric field models in driving a model of the plasmasphere. J Geophys Res 119(6): 4621–4638. https://doi.org/10.1002/2014JA019836. [CrossRef] [Google Scholar]
 Ridley AJ, Gombosi TI, DeZeeuw DL. 2004. Ionospheric control of the magnetosphere: Conductance. Ann Geophys 22: 567–584. https://doi.org/10.5194/angeo225672004. [CrossRef] [Google Scholar]
 Roe PL. 1981. Approximate Riemann solvers, parameter vectors, and difference schemes. J Comput Phys 43: 357–372. https://doi.org/10.1006/jcph.1997.5705. [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Roe PL, Balsara DS. 1996. Notes on the eigensystem of magnetohydrodynamics. SIAM J Appl Math 56(1): 57–67. https://doi.org/10.1137/S003613999427084X. [CrossRef] [Google Scholar]
 Rostoker G. 1972. Geomagnetic indices. Rev Geophys 10(4): 935. https://doi.org/10.1029/RG010i004p00935. [CrossRef] [Google Scholar]
 Roussev II. 2008. Eruptive events in the solar atmosphere: new insights from theory and 3D numerical modelling. J Contemp Phys 49(4): 237–254. https://doi.org/10.1080/00107510802366658. [CrossRef] [Google Scholar]
 Roussev II, Forbes TG, Gombosi TI, Sokolov IV, DeZeeuw DL, Birn J. 2003. A threedimensional flux rope model for coronal mass ejections based on a loss of equilibrium. Astrophys J Lett 588: L45–L48. https://doi.org/10.1086/375442. [NASA ADS] [CrossRef] [Google Scholar]
 Roussev II, Lugaz N, Sokolov IV. 2007. New physical insight on the changes in magnetic topology during coronal mass ejections: Case studies for the 2002 April 21 and August 24 Events. Astrophys J Lett 668: L87–L90. https://doi.org/10.1086/522588. [NASA ADS] [CrossRef] [Google Scholar]
 Roussev II, Sokolov IV. 2006. Models of solar eruptions: Recent advances from theory and simulations. Geophys Monograph 165: 89–102. https://doi.org/10.1029/165gm10. [Google Scholar]
 Roussev II, Sokolov IV, Forbes TG, Gombosi TI, Lee MA, Sakai JI. 2004. A numerical model of a coronal mass ejection: Shock development with implications for the acceleration of GeV protons. Astrophys J Lett 605: L73–L76. https://doi.org/10.1086/392504. [CrossRef] [Google Scholar]
 Rubin M, Hansen KC, Combi MR, Daldorff LKS, Gombosi TI, Tenishev VM. 2012. KelvinHelmholtz instabilities at the magnetic cavity boundary of comet 67P/ChuryumovGerasimenko. J Geophys Res 117(A06): 227. https://doi.org/10.1029/2011JA017300. [NASA ADS] [CrossRef] [Google Scholar]
 Rubin M, Jia X, Altwegg K, Combi MR, Daldorff LKS, et al. 2015. Selfconsistent multifluid MHD simulations of Europa’s exospheric interaction with Jupiter’s magnetosphere. J Geophys Res 120. https://doi.org/10.1002/2015JA021149. [Google Scholar]
 Sachdeva N, van der Holst B, Manchester WB, Tóth G, Chen Y, et al. 2019. Validation of the Alfvén Wave Solar Atmosphere Model (AWSoM) with Observations from the Low Corona to 1 au. Astrophys J 887(1): 83. https://doi.org/10.3847/15384357/ab4f5e. [NASA ADS] [CrossRef] [Google Scholar]
 Sarkango Y, Jia X, Toth G. 2019. Global MHD simulations of the response of Jupiter’s magnetosphere and ionosphere to changes in the solar wind and IMF. J Geophys Res 124: https://doi.org/10.1029/2019JA026787. [Google Scholar]
 Sazykin S, Wolf RA, Spiro RW, Gombosi TI, De Zeeuw DL, Thomsen MF. 2002. Interchange instability in the inner magnetosphere associated with geosynchronous particle flux decreases. Geophys Res Lett 29(10): 881–884. https://doi.org/10.1029/2001GL014416. [CrossRef] [Google Scholar]
 Schunk RW, Nagy AF. 1980. Ionospheres of the terrestrial planets. Rev Geophys Space Phys 18: 813–852. [NASA ADS] [CrossRef] [Google Scholar]
 Schunk RW, Nagy AF. 2009. Ionospheres: Physics, Plasma Physics, and Chemistry. Cambridge University Press, Cambridge, UK. https://doi.org/10.1017/cbo9780511635342. [CrossRef] [Google Scholar]
 Shumlak U, Loverich J. 2003. Approximate Riemann solver for the twofluid plasma model. J Comput Phys 187(2): 620–638. https://doi.org/10.1016/S00219991(03)001517. [CrossRef] [Google Scholar]
 Siscoe GL, Crooker NU, Erickson GM, Sonnerup BUÖ, Siebert KD, Weimer DR, White WW, Maynard NC. 2000. Global geometry of magnetospheric currents inferred from MHD Simulations. Geophys Monograph 118: 41–52. https://doi.org/10.1029/GM118p0041. [Google Scholar]
 Skilling J. 1971. Cosmic rays in the Galaxy: Convection or diffusion? Astrophys J 170: 265. [NASA ADS] [CrossRef] [Google Scholar]
 Slavin JA, Acuna MH, Anderson BJ, Baker DN, Benna M, et al. 2009. MESSENGER observations of magnetic reconnection in Mercury’s magnetosphere. Science 324(5927): 606–610. https://doi.org/10.1126/science.1172011. [NASA ADS] [CrossRef] [Google Scholar]
 Slavin JA, DiBraccio GA, Gershman DJ, Imber SM, Poh GK, et al. 2014. MESSENGER observations of Mercury’s dayside magnetosphere under extreme solar wind conditions. J Geophys Res 119(10): 8087–8116. https://doi.org/10.1002/2014JA020319. [NASA ADS] [CrossRef] [Google Scholar]
 Sokolov IV, Roussev II, Gombosi TI, Lee MA, Kóta J, Forbes TG, Manchester WB IV, Sakai JI. 2004. A new field line advection model for solar particle acceleration. Astrophys J Lett 616: L171–L174. https://doi.org/10.1086/426812. [NASA ADS] [CrossRef] [Google Scholar]
 Sokolov IV, van der Holst B, Manchester WB, Ozturk DCS, Szente J, Taktakishvili A, Tóth G, Jin M, Gombosi TI. 2021. Threadedfieldlines model for the low solar corona powered by the Alfvén wave turbulence. Astrophys J 908(1): 172–183. https://doi.org/10.3847/15384357/abc000. [CrossRef] [Google Scholar]
 Sokolov IV, van der Holst B, Oran R, Downs C, Roussev II, Jin M, Manchester WB, Evans RM, Gombosi TI. 2013. Magnetohydrodynamic waves and coronal heating: Unifying empirical and MHD Turbulence Models. Astrophys J 764: 23. https://doi.org/10.1088/0004637X/764/1/23. [NASA ADS] [CrossRef] [Google Scholar]
 Solomon SC, Hays PB, Abreu VJ. 1988. The auroral 6300Å emission: Observations and modeling. J Geophys Res 93(A9): 9867. https://doi.org/10.1029/ja093ia09p09867. [CrossRef] [Google Scholar]
 Stout QF, Zeeuw DLD, Gombosi TI, Groth CPT, Marshall HG, Powell KG. 1997. Proceedings of the 1997 ACM/IEEE Conference on Supercomputing (CDROM) – Supercomputing ‘97. ACM Press. https://doi.org/10.1145/509593.509650. [Google Scholar]
 Sugiyama T, Kusano K. 2007. Multiscale plasma simulation by the interlocking of magnetohydrodynamic model and particleincell kinetic model. J Comput Phys 227(2): 1340–1352. https://doi.org/10.1016/j.jcp.2007.09.011. [CrossRef] [Google Scholar]
 Sun H, Hua Z, Ren J, Zou S, Sun Y, Chen Y. 2021. Matrix completion methods for the total electron content video reconstruction. Arxiv: http://arxiv.org/abs/2012.01618. [Google Scholar]
 Sun WJ, Fu SY, Slavin JA, Raines JM, Zong QG, Poh GK, Zurbuchen TH. 2016. Spatial distribution of Mercury’s flux ropes and reconnection fronts: MESSENGER observations. J Geophys Res 121(8): 7590–7607. https://doi.org/10.1002/2016JA022787. [CrossRef] [Google Scholar]
 Sun WJ, Slavin JA, Dewey RM, Chen Y, DiBraccio GA, Raines JM, Jasinski JM, Jia X, AkhavanTafti M. 2020. MESSENGER Observations of Mercury’s nightside magnetosphere under extreme solar wind conditions: Reconnectiongenerated structures and steady convection. J Geophys Res 125(3): e2019JA027490. https://doi.org/10.1029/2019JA027490. [Google Scholar]
 Szente J, Landi E, Manchester WB, Toth G, van der Holst B, Gombosi TI. 2019. SPECTRUM: Synthetic Spectral Calculations for Global Space Plasma Modeling. Astrophys J Suppl 242(1): 1. https://doi.org/10.3847/15384365/ab16d0. [CrossRef] [Google Scholar]
 Takizuka T, Abe H. 1977. A binary collision model for plasma simulation with a particle code. J Comput Phys 25(3): 205–219. https://doi.org/10.1016/00219991(77)900997. [CrossRef] [Google Scholar]
 Tenishev V, Borovikov D, Combi MR, Sokolov I, Gombosi T. 2018. Toward development of the energetic particle radiation nowcast model for assessing the radiation environment in the altitude range from that used by the commercial aviation in the troposphere to LEO, MEO, and GEO. In: 2018 Atmospheric and Space Environments Conference. AIAA. https://doi.org/10.2514/6.20183650 [Google Scholar]
 Tenishev V, Combi M, Davidsson B. 2008. A global kinetic model for cometary comae: The evolution of the coma of the Rosetta target comet ChuryumovGerasimenko throughout the mission. Astrophys J 685(1): 659–677. https://doi.org/10.1086/590376. [NASA ADS] [CrossRef] [Google Scholar]
 Tenishev V, Combi M, Sokolov IV, Roussev II, Gombosi TI. 2005. Numerical studies of the solar energetic particle transport and acceleration, AIAA, Toronto, Ontario, Canada. https://doi.org/10.2514/6.20054928. [Google Scholar]
 Tenishev V, Combi MR, Rubin M. 2011. Numerical simulation of dust in a cometary coma: Application to comet 67P/ChuryumovGerasimenko. Astrophys J 732(2): 104. https://doi.org/10.1088/0004637x/732/2/104. [NASA ADS] [CrossRef] [Google Scholar]
 Tenishev V, Fougere N, Borovikov D, Combi MR, Bieler A, et al. 2016. Analysis of the dust jet imaged by Rosetta VIRTISM in the coma of comet 67P/ChuryumovGerasimenko on April 12, 2015. Mon Not R Astron Soc 462(Suppl 1): S370–S375. https://doi.org/10.1093/mnras/stw2793. [CrossRef] [Google Scholar]
 Tenishev V, Shou Y, Borovikov D, Lee Y, Fougere N, Michael A, Combi MR. 2021. Application of the Monte Carlo method in modeling dusty gas, dust in plasma, and energetic ions in planetary, magnetospheric, and heliospheric environments. J Geophys Res 126(2): e2020JA028242. https://doi.org/10.1029/2020JA028242. [CrossRef] [Google Scholar]
 Titov VS, Démoulin P. 1999. Basic topology of twisted magnetic configurations in solar flares. Astron Astrophys 351: 707–720. [Google Scholar]
 Toffoletto F, Sazykin S, Spiro R, Wolf R. 2003. Inner magnetospheric modeling with the Rice Convection Model. Space Sci Rev 107: 175–196. https://doi.org/10.1023/A:1025532008047. [CrossRef] [Google Scholar]
 Toth G. 1996. A general code for modeling MHD flows on parallel computers: Versatile advection code. Astroph Lett Comm 34: 245–250. [Google Scholar]
 Toth G. 2000. The ∇·B = 0 constraint in shock capturing magnetohydrodynamic codes. J Comput Phys 161(2): 605–652. https://doi.org/10.1006/jcph.2000.6519. [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Toth G. 2006. Flexible, efficient and robust algorithm for parallel execution and coupling of components in a framework. Comput Phys Commun 174(10): 793–802. https://doi.org/10.1016/j.cpc.2005.12.017. [CrossRef] [Google Scholar]
 Toth G, Chen Y, Gombosi TI, Cassak P, Markidis S, Peng IB. 2017. Scaling the ion inertial length and its implications for modeling reconnection in global simulations. J Geophys Res 122(10): 10336–10355. https://doi.org/10.1002/2017JA024189. [CrossRef] [Google Scholar]
 Toth 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. [CrossRef] [Google Scholar]
 Toth G, De Zeeuw DL, Gombosi TI, Powell KG. 2006. A parallel explicit/implicit time stepping scheme on blockadaptive grids. J Comput Phys 217: 722–758. https://doi.org/10.1016/j.jcp.2006.01.029. [NASA ADS] [CrossRef] [Google Scholar]
 Toth G, Jia X, Markidis S, Peng B, Chen Y, et al. 2016. Extended Magnetohydrodynamics with Embedded ParticleinCell Simulation of Ganymede’s Magnetosphere. J Geophys Res 121(2): 1273–1293. https://doi.org/10.1002/2015JA021997. [CrossRef] [Google Scholar]
 Toth G, Keppens R, Botchev MA. 1998. Implicit and semiimplicit schemes in the Versatile Advection Code: Numerical tests. Astron Astrophys 332: 1159–1170. [Google Scholar]
 Toth G, Kovács D, Hansen KC, Gombosi TI. 2004. 3D MHD Simulations of the Magnetosphere of Uranus. J Geophys Res 109(A11): A11210. https://doi.org/10.1029/2004JA010406. [NASA ADS] [CrossRef] [Google Scholar]
 Toth G, Ma YJ, Gombosi TI. 2008. Hall magnetohydrodynamics on block adaptive grids. J Comput Phys 227: 6967–6984. https://doi.org/10.1016/j.jcp.2008.04.010. [NASA ADS] [CrossRef] [Google Scholar]
 Toth G, Meng X, Gombosi TI, Rastätter L. 2014. Predicting the time derivative of local magnetic perturbations. J Geophys Res 119. https://doi.org/10.1002/2013JA019456. [Google Scholar]
 Toth G, Sokolov IV, Gombosi TI, Chesney DR, Clauer CR, et al. 2005. Space weather modeling framework: A new tool for the space science community. J Geophys Res 110(A12): 12226–12237. https://doi.org/10.1029/2005JA011126. [NASA ADS] [CrossRef] [Google Scholar]
 Toth 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]
 Tsyganenko NA. 1989. Global quantitative models of the geomagnetic field in the cislunar magnetosphere for different disturbance levels. Planet Space Sci 35: 1347. [CrossRef] [Google Scholar]
 Tsyganenko NA. 1995. Modeling the Earth’s magnetospheric magnetic field confined within a realistic magnetopause. J Geophys Res 100: 5599–5612. [CrossRef] [Google Scholar]
 Tsyganenko NA. 2002a. A model of the near magnetosphere with a dawndusk asymmetry – 1. Mathematical Structure. J Geophys Res 107: SMP 121–SMP 1215. https://doi.org/10.1029/2001JA000219. [CrossRef] [Google Scholar]
 Tsyganenko NA. 2002b. A model of the near magnetosphere with a dawndusk asymmetry – 2. Parameterization and fitting to observations. J Geophys Res 107: SMP 101–SMP 1017. https://doi.org/10.1029/2001JA000220. [CrossRef] [Google Scholar]
 Tsyganenko NA, Sitnov MI. 2005. Modeling the dynamics of the inner magnetosphere during strong geomagnetic storms. J Geophys Res 110(A9): 3208. https://doi.org/10.1029/2004JA010798. [CrossRef] [Google Scholar]
 Tsyganenko NA, Stern D. 1996. Modeling the global magnetic field of the largescale birkeland current systems. J Geophys Res 101(27): 187. [NASA ADS] [CrossRef] [Google Scholar]
 Usadi A, Kageyama A, Watanabe K, Sato T. 1993. A global simulation of the magnetosphere with a long tail: Southward and northward interplanetary magnetic field. J Geophys Res 98: 7503–7517. [NASA ADS] [CrossRef] [Google Scholar]
 Usmanov AV. 1993. A global numerical 3D MHD model of the solar wind. Sol Phys 146: 377–396. https://doi.org/10.1007/BF00662021. [CrossRef] [Google Scholar]
 Usmanov AV, Goldstein ML, Besser BP, Fritzer JM. 2000. A global MHD solar wind model with WKB Alfvén waves: Comparison with Ulysses data. J Geophys Res 105: 12675–12695. https://doi.org/10.1029/1999JA000233. [NASA ADS] [CrossRef] [Google Scholar]
 Vainio R, Desorgher L, Heynderickx D, Storini M, Flückiger E, et al. 2008. Dynamics of the Earth’s particle radiation environment. Space Sci Rev 147: 187–231. https://doi.org/10.1007/s1121400994967. [NASA ADS] [CrossRef] [Google Scholar]
 van der Holst B, Jacobs C, Poedts S. 2007. Simulation of a Breakout Coronal Mass Ejection in the Solar Wind. Astrophys. J. 671(1): L77–L80. https://doi.org/10.1086/524732. [NASA ADS] [CrossRef] [Google Scholar]
 van der Holst B, Manchester W IV, Sokolov IV, Toth G, Gombosi TI, DeZeeuw D, Cohen O. 2009. Breakout coronal mass ejection or streamer blowout: The bugle effect. Astrophys J 693(2): 1178–1187. https://doi.org/10.1088/0004637X/693/2/1178. [NASA ADS] [CrossRef] [Google Scholar]
 van der Holst B, Manchester WB, Frazin RA, Vásquez AM, Toth G, Gombosi TI. 2010. A datadriven, twotemperature solar wind model with Alfvén waves. Astrophys J 725(1): 1373–1383. https://doi.org/10.1088/0004637X/725/1/1373. [NASA ADS] [CrossRef] [Google Scholar]
 van der Holst B, Sokolov IV, Meng X, Jin M, Manchester WB, Toth G, Gombosi TI. 2014. Alfvén Wave Solar Model (AWSoM): Coronal heating. Astrophys J 782(2): 81. https://doi.org/10.1088/0004637X/782/2/81. [NASA ADS] [CrossRef] [Google Scholar]
 van Leer B. 1973. Towards the ultimate conservative difference scheme. I. The quest of monoticity, Vol. 18of Lecture Notes in Physics, Berlin Springer Verlag. pp. 163–168. [CrossRef] [Google Scholar]
 van Leer B. 1974. Towards the ultimate conservative difference scheme. II. Monoticity and conservation combined in a secondorder scheme. J Comput Phys 14(4): 361–370. https://doi.org/10.1016/00219991(74)900199. [NASA ADS] [CrossRef] [Google Scholar]
 van Leer B. 1977a. Towards the ultimate conservative difference scheme. III. Upstreamcentered finitedifference schemes for ideal compressible flow. J Comput Phys 23(3): 263–275. https://doi.org/10.1016/00219991(77)900948. [NASA ADS] [CrossRef] [Google Scholar]
 van Leer B. 1977b. Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection. J Comput Phys 23(3): 276–299. https://doi.org/10.1016/00219991(77)90095x. [NASA ADS] [CrossRef] [Google Scholar]
 van Leer B. 1979. Towards the ultimate conservative difference scheme. V. A secondorder Sequel to Godunov’s method. J Comput Phys 32(1): 101–136. https://doi.org/10.1016/00219991(79)901451. [NASA ADS] [CrossRef] [Google Scholar]
 Wang X, Chen Y, Toth G, Manchester WB, Gombosi TI, Hero AO, Jiao Z, Sun H, Jin M, Liu Y. 2020. Predicting solar flares with machine learning: Investigating solar cycle dependence. Astrophys J 895(1): 3. https://doi.org/10.3847/15384357/ab89ac. [CrossRef] [Google Scholar]
 Wang Z, Zou S, Coppeans T, Ren J, Ridley A, Gombosi T. 2019. Segmentation of SED by boundary flows associated with westward drifting partial ring current. Geophys Res Lett 46(14): 7920–7928. [CrossRef] [Google Scholar]
 Washimi H, Tanaka T. 1996. 3D Magnetic field and current system in the heliosphere. Space Sci Rev 78: 85–94. [NASA ADS] [CrossRef] [Google Scholar]
 Watanabe K, Sato T. 1990. Global simulation of the solar windmagnetosphere interaction: The importance of its numerical validity. J Geophys Res 95: 75–88. [CrossRef] [Google Scholar]
 Weigel RS, Klimas AJ, Vassiliadis D. 2003. Solar wind coupling to and predictability of ground magnetic fields and their time derivatives. J Geophys Res 108(A7): 1298. https://doi.org/10.1029/2002JA009627. [CrossRef] [Google Scholar]
 Welling D. 2019. Magnetohydrodynamic models of B and their use in GIC estimates, American Geophysical Union (AGU). pp. 43–65. https://doi.org/10.1002/9781119434412.ch3. [Google Scholar]
 Welling D, Love J, Rigler EJ, Oliveira D, Komar C. 2020. Numerical simulations of the geospace response to the arrival of a perfect interplanetary coronal mass ejection. Earth Space Sci Open Archive 16. https://doi.org/10.1002/essoar.10502106.1. [Google Scholar]
 Welling DT, Jordanova VK, Glocer A, Toth G, Liemohn MW, Weimer DR. 2015. The twoway relationship between ionospheric outflow and the ring current. J Geophys Res 120(6): 4338–4353. https://doi.org/10.1002/2015JA021231. [Google Scholar]
 Welling DT, Jordanova VK, Zaharia SG, Glocer A, Toth G. 2011. The effects of dynamic ionospheric outflow on the ring current. J Geophys Res 116(A2): A00J19. https://doi.org/10.1029/2010JA015642. [Google Scholar]
 Welling DT, Toth G, Jordanova VK, Yu Y. 2018. Integration of RAMSCB into the Space Weather Modeling Framework. J Atmos SolarTerr Phys 177: 160–168. https://doi.org/10.1016/j.jastp.2018.01.007. [Google Scholar]
 Welling DT, Zaharia SG. 2012. Ionospheric outflow and cross polar cap potential: What is the role of magnetospheric inflation? Geophys Res Lett 39(L23): 101. https://doi.org/10.1029/2012GL054228. [Google Scholar]
 Wentzel G. 1926. Eine Verallgemeinerung der Quantenbedingungen für die Zwecke der Wellenmechanik. Zeitschrift für Physik 38(6): 518–529. https://doi.org/10.1007/BF01397171. [Google Scholar]
 White WW, Siscoe GL, Erickson GM, Kaymaz Z, Maynard NC, Siebert KD, Sonnerup BUÖ, Weimer DR. 1998. The magnetospheric sash and the crosstail S. Geophys Res Lett 25(10): 1605–1608. [Google Scholar]
 Wiltberger M, Lotko W, Lyon JG, Damiano P, Merkin V. 2010. Influence of cusp O+ outflow on magnetotail dynamics in a multifluid MHD model of the magnetosphere. J Geophys Res 115: A00J05. https://doi.org/10.1029/2010JA015579. [Google Scholar]
 Winglee RM. 1998. Multifluid simulations of the magnetosphere: The identification of the geopause and its variation with IMF. Geophys Res Lett 25: 4441–4444. [Google Scholar]
 Winglee RM, Lewis W, Lu G. 2005. Mapping of the heavy ion outflows as seen by IMAGE and multifluid global modeling for the 17 April 2002 storm. J Geophys Res 110: A12S24. https://doi.org/10.1029/2004JA010909. [Google Scholar]
 Winslow RM, Anderson BJ, Johnson CL, Slavin JA, Korth H, Purucker ME, Baker DN, Solomon SC. 2013. Mercury’s magnetopause and bow shock from MESSENGER Magnetometer observations. J Geophys Res 118(5): 2213–2227. https://doi.org/10.1002/jgra.50237. [Google Scholar]
 Wolf R. 1974. Calculations of magnetospheric electric fields. In: Magnetospheric Physics. McCormac BM (Ed.), D. Reidel Publishing, Hingham, MA. pp. 167–177. [Google Scholar]
 Wolf RA, Harel M, Spiro RW, Voigt G, Reiff PH, Chen CK. 1982. Computer simulation of inner magnetospheric dynamics for the magnetic storm of July 29, 1977. J Geophys Res 87: 5949–5962. https://doi.org/10.1029/JA087iA08p05949. [Google Scholar]
 Wolf RA, Spiro RW, Rich FJ. 1991. Extension of the Rice Convection Model into the highlatitude ionosphere. J Atmos SolarTerr Phys 50: 817–829. [Google Scholar]
 Wu CC, Walker R, Dawson JM. 1981. A threedimensional MHD model of the Earth’s magnetosphere. Geophys Res Lett 8: 523–526. https://doi.org/10.1029/GL008i005p00523. [Google Scholar]
 Yang J, Toffoletto FR, Wolf RA. 2014. RCME simulation of a thin arc preceded by a northsouthaligned auroral streamer. Geophys Res Lett 41(8): 2695–2701. https://doi.org/10.1002/2014GL059840. [Google Scholar]
 Yu Y, Cao J, Fu H, Lu H, Yao Z. 2017. The effects of bursty bulk flows on globalscale current systems. J Geophys Res 122(6): 6139–6149. https://doi.org/10.1002/2017JA024168. [Google Scholar]
 Yu Y, Jordanova V, Zou S, Heelis R, Ruohoniemi M, Wygant J. 2015. Modeling subauroral polarization streams during the 17 March 2013 storm. J Geophys Res 120(3): 1738–1750. [Google Scholar]
 Yu Y, Ridley AJ. 2008. Validation of the space weather modeling framework using groundbased magnetometers. Space Weather 6(S05): 002. https://doi.org/10.1029/2007SW000345. [Google Scholar]
 Yu Y, Ridley AJ, Welling DT, Toth G. 2010. Including gap region fieldaligned currents and magnetospheric currents in the MHD calculation of groundbased magnetic field perturbations. J Geophys Res 115(A08): 207. https://doi.org/10.1029/2009JA014869. [Google Scholar]
 Zaharia S. 2008. Improved Euler potential method for threedimensional magnetospheric equilibrium. J Geophys Res 113(A8). https://doi.org/10.1029/2008JA013325. [Google Scholar]
 Zaharia S, Cheng CZ, Maezawa K. 2004. 3D Forcebalanced magnetospheric configurations. Ann Geophys 22: 251–266. [Google Scholar]
 Zaharia S, Jordanova VK, Thomsen MF, Reeves GD. 2006. Selfconsistent modeling of magnetic fields and plasmas in the inner magnetosphere: Application to a geomagnetic storm. J Geophys Res 111(A11): A11S14. https://doi.org/10.1029/2006JA011619. [Google Scholar]
 Zaharia S, Jordanova VK, Welling DT, Toth G. 2010. Selfconsistent inner magnetosphere simulation driven by a global MHD model. J Geophys Res 115(A12): 228. https://doi.org/10.1029/2010JA015915. [Google Scholar]
 Zaharia S, Thomsen MF, Birn J, Denton MH, Jordanova VK, Cheng CZ. 2005. Effect of stormtime plasma pressure on the magnetic field in the inner magnetosphere. Geophys Res Lett 32(L03): 102. [Google Scholar]
 Zhang B, Sorathia KA, Lyon JG, Merkin VG, Garretson JS, Wiltberger M. 2019a. GAMERA: A Threedimensional Finitevolume MHD Solver for Nonorthogonal Curvilinear Geometries. Astrophys J Suppl 244(1): 20. https://doi.org/10.3847/15384365/ab3a4c. [Google Scholar]
 Zhang J, Liemohn MW, De Zeeuw DL, Borovsky JE, Ridley AJ, et al. 2007. Understanding stormtime ring current development through datamodel comparisons of a moderate storm. J Geophys Res 112(A04): 208. https://doi.org/10.1029/2006JA011846. [Google Scholar]
 Zhang M, Qin G, Rassoul H. 2009. Propagation of solar energetic particles in threedimensional interplanetary magnetic fields. Astrophys J 692(1): 109–132. https://doi.org/10.1088/0004637x/692/1/109. [NASA ADS] [CrossRef] [Google Scholar]
 Zhang W, Almgren A, Beckner V, Bell J, Blaschke J, et al. 2019b. AMReX: A framework for blockstructured adaptive mesh refinement. J Open Source Softw 4(37): 1370. https://doi.org/10.21105/joss.01370. [Google Scholar]
 Zhang W, Myers A, Gott K, Almgren A, Bell J. 2020. AMReX: BlockStructured Adaptive Mesh Refinement for Multiphysics Applications. http://arxiv.org/abs/2009.12009. [Google Scholar]
 Zhao L, Zhang M, Rassoul HK. 2016. Double power laws in the eventintegrated solar energetic particle spectrum. Astrophys J 821(1): 62. https://doi.org/10.3847/0004637x/821/1/62. [CrossRef] [Google Scholar]
 Zhao L, Zhang M, Rassoul HK. 2017. The effects of interplanetary transport in the eventintergrated solar energetic particle spectra. Astrophys J 836(1): 31. https://doi.org/10.3847/15384357/836/1/31. [Google Scholar]
 Zhou H, Tóth G. 2020. Efficient OpenMP parallelization to a complex MPI parallel magnetohydrodynamics code. J Parallel Distributed Comput 139: 65–74. https://doi.org/10.1016/j.jpdc.2020.02.004. [Google Scholar]
 Zhou H, Tóth G, Jia X, Chen Y. 2020. Reconnectiondriven dynamics at Ganymede’s Upstream Magnetosphere: 3D Global Hall MHD and MHDEPIC Simulations. J Geophys Res 125(8): e28162. https://doi.org/10.1029/2020JA028162. [Google Scholar]
 Zhou H, Tóth G, Jia X, Chen Y, Markidis S. 2019. Embedded kinetic simulation of Ganymede’s magnetosphere: Improvements and inferences. J Geophys Res 124: 5441–5460. https://doi.org/10.1029/2019JA026643. [Google Scholar]
 Zieger B, Hansen KC, Gombosi TI, De Zeeuw DL. 2010. Periodic plasma escape from the massloaded Kronian magnetosphere. J Geophys Res 115: A8: https://doi.org/10.1029/2009JA014951. [Google Scholar]
 Zou S, Ozturk D, Varney R, Reimer A. 2017. Effects of sudden commencement on the ionosphere: PFISR observations and global MHD simulation. Geophys Res Lett 44(7): 3047–3058. https://doi.org/10.1002/2017GL072678. [Google Scholar]
Cite this article as: Gombosi TI, Chen Y, Glocer A, Huang Z, Jia X, et al. 2021. What sustained multidisciplinary research can achieve: The space weather modeling framework. J. Space Weather Space Clim. 11, 42. https://doi.org/10.1051/swsc/2021020.
All Tables
All Figures
Fig. 1 Overview of the BATSRUS multiphysics code. 

In the text 
Fig. 2 Schematic diagram of the Space Weather Modeling Framework. The SWMF and its core models are open source (https://github.com/MSTEMQUDA), while the full SWMF is available via registration under a user license (http://csem.engin.umich.edu/tools/smmf). 

In the text 
Fig. 3 Overview of the AWSoM and AWSoMR physics. They solve XMHD equations with separate ion and electron temperatures. The energy densities of parallel and antiparallel propagating turbulence that are selfconsistently coupled to each other and to the plasma are solved together with the XMHD equations. Heat conduction and radiative cooling are also taken into account. The turbulence is powered by the Poynting flux leaving the solar photosphere. 

In the text 
Fig. 4 Illustration of the models (components within SWMF) and couplings in the SWMF/Geospace configuration. Arrows denote the information that is passed between the components (adapted from Haiducek et al., 2017). 

In the text 
Fig. 5 Grid configurations for BATSRUS within the SWMF Geospace. The left and right hand panels illustrate the grid configuration of the operational Geospace model versions 1 and 2, respectively (from Haiducek et al., 2017). 

In the text 
Fig. 6 Background corona and solar wind solutions with the AWSoMR model for solar minimum and maximum conditions. The background solar wind is driven by an outward going photospheric turbulent energy flux per unit magnetic field of 1 MJ m^{−2} s^{−1} Tesla^{−1} (CR2209) and by 0.45 MJ m^{−2} s^{−1} Tesla^{−1} (CR2123). 

In the text 
Fig. 7 Three examples of Gibson & Low (1998) flux ropes with different size and magnetic strength parameters. Panels (a)–(f) and (g)–(i) show, respectively, flux ropes specified with radii of 0.8 and 0.6 R_{s}. Strength parameters are set to 0.6 for model run (a)–(c) and 2.25 for (d)–(i). The left column shows the initial configuration of the flux ropes with blue and red isosurfaces showing, respectively, the ratios of 0.3 and 2.5 of the mass density of the CME model divided by that of the preevent corona. The middle column shows the resulting CME evolution at t = 20 min. Here, magnetic field lines are colored red, grayshaded and green to illustrate the flux rope, largescale helmet streamers, and magnetic fields surrounding active regions and open flux. Color contour images show the ratio of the mass density of the CME divided by that of the preevent corona. The right column shows modelproduced SOHO/LASCO white light images, where the total brightness is normalized by dividing by that of the preevent background solar wind (from Jin et al., 2017a). 

In the text 
Fig. 8 CMEdriven EUV waves in the simulation (left) and in the corresponding SDO/AIA observation (right). Both the simulation and observation images are produced by a triratio running difference method. The tricolor channels are AIA 211 Å (red), AIA 193 Å (green), and AIA 171 Å (blue). The ratio in each channel is identically scaled to 1 ± 0.2 for both observation and simulation (from Jin et al., 2017b). 

In the text 
Fig. 9 Comparison showing a general agreement between the whitelight observations from SOHO LASCO C2 (top left) and STEREOB COR1 (top right) and the respective synthesized whitelight images from the simulation (bottom). The color contours show the relative total brightness changes compared to the preevent background level (from Jin et al., 2018). 

In the text 
Fig. 10 1 AU results of the EEGGL simulation of the 12 July 2012 CME event simulated at CCMC. Shown are the simulated and observed plasma quantities plotted with dashed and solid lines, respectively. From top to bottom are the magnetic field component B_{x}, B_{y} and B_{z}, the mass density, and the Earthdirected velocity V_{x}. Simulation results are shifted 10 h to match the shock arrival. We find good agreement, with the exception of the B_{x}–B_{y} rotation and the excessive trailing velocity. 

In the text 
Fig. 11 Distribution of the energetic particles (> 10 MeV) along the extracted magnetic field lines at 10 min (left panel) and 20 min (right panel) after the eruption of CME. The flux is in the unit of particle flux unit (pfu, particles/cm^{2}/s/sr). The green isosurface represents the leading edge of the CME. 

In the text 
Fig. 12 Example of gyroradii of particles with 1 MeV < E < 16 MeV during quiet geomagnetic conditions. Left: Gyroradius map in the equatorial plane. Right: Gyroradius map in the meridional plane (X = 0). The gyroradii of these particles can be as large as tens of R_{E}. Here, Xaxis is directed toward the Sun, and Yaxis is in the equatorial plane, and Zaxis is such that the frame of reference is righthanded. The freespace energy spectrum of the simulated energetic particles is taken from Badavi et al. (2011). 

In the text 
Fig. 13 Example of the calculated density of energetic protons with energies 1 MeV < E < 100 MeV) in geospace. Both the SEP’s energy spectrum and geomagnetic parameters are taken for quiet conditions. The figure demonstrates that the topology of the SEPs population in the geospace is affected by the particles’ gyromotion. Here, Xaxis is directed toward the Sun, Yaxis in the equatorial plane, and Zaxis is such that the coordinate frame is righthanded. 

In the text 
Fig. 14 Example of applying AMPS for rigidity cutoff calculation. The map is calculated for an altitude of 500 km. Left: Rigidity cutoff map calculated for quiet geomagnetic conditions. Right: Depression of the rigidity cutoff during a geomagnetic storm. The calculation was performed for conditions of the geomagnetic storm on 17 March 2015. One can see that the general rigidity cutoff patterns have changed mostly in the midlatitude region (from Tenishev et al., 2021). 

In the text 
Fig. 15 Results of highresolution SWMF/BATSRUS simulation with 1/16 R_{E} grid resolution in the tail and magnetopause region in order to resolve small and mesoscale structures. It shows the filamentary current structures associated with the presence of both Kelvin–Helmholtz vortices at the flanks of the magnetopause and of bursty structures in the tail. Note that the associated flow velocities are not shown (Dorelli & Buzulukova, personal communications, 2020). 

In the text 
Fig. 16 An illustration of cusp accelerated ion outflow due to soft electron precipitation and waveparticle interactions as modeled by PWOM; adapted from (Glocer et al., 2018). The top panel shows the PWOM computed ion distribution function along a single field line and demonstrates the ion conic evolution with altitude. Here v_{∥} and v_{⊥} are the parallel and perpendicular velocity and the color contour shows the log of phase space density in normalized, dimensionless, units. The bottom panel shows the global simulation, including the kinetic processes in the cusp using 896 field lines. Here u_{O} and u_{H} are the bulk field aligned velocities for O^{+} and H^{+} in units of km/s, and T_{e} is the electron temperature in Kelvin. 

In the text 
Fig. 17 FACs and ground magnetic perturbations due to shock compression on 17 March 2015 from a highresolution BATSRUS run coupled with CRCM and IE. The simulated magnetic perturbations at Poker Flat were compared with the real magnetometer observations. The results shown here illustrate the capacity of the SWMF to resolve mesoscale features of the magnetospheric dynamics in highresolution MHD (Zou et al., 2017). 

In the text 
Fig. 18 (Left) Comparison of the storm SYMH index for two storms from January 2005. The thick trace shows the observations, the red and orange traces show the normal resolution and high resolution SWMF/Geospace simulations, and the blue trace shows the SWMF/Geospace simulation run without the inner magnetosphere RCM component (from Haiducek et al., 2017). (Right) Scatter plot of the observed real time Dst time series (vertical axis values) against a SWMF prediction of the Dst (horizontal axis values). The green, orange, and cyan contour lines show the regions of 40, 80, and 120 values within a 2.6 × 2.6nT bin. The diagonal longdashed line shows the ideal datamodel relationship. The two dashed lines show the −50 nT threshold values (from Liemohn et al., 2018). 

In the text 
Fig. 19 MHDAEPIC simulation of a geomagnetic storm. The 2D cuts display a part of the much larger 3D domain. Colors show the X component of the velocity and the white lines are traces of the magnetic field in the meridional plane. The black lines in the two snapshots, separated by 159 minutes of simulation time, indicate the edges of the active PIC regions. As the tail evolves, the active PIC region is continuously adapted to cover the reconnection sites of interest. 

In the text 
Fig. 20 Snapshots showing By strength (color) and the projected magnetic field lines in the meridional plane inside the PIC region. The color bar is different in each plot (from Chen et al., 2017). 

In the text 
Fig. 21 Crescent electron and ion phase space distributions (a) E_{x} (mV/m) in the meridional plane at t = 3600 s; (b) Normalized electron distribution in V_{y}–V_{x} phase space; and (c) Ion phase space distribution. The phasespace density is normalized (from Chen et al., 2017). 

In the text 
Fig. 22 Cuts through the simulation of the MESSENGER M2 flyby. (a) XZ cut at Y = 0 (noonmidnight meridian) with color contours of plasma thermal pressure. The thick magenta line shows the modeled magnetopause boundary identified based on the current density, while the white dotted line represents the databased empirical magnetopause model of Winslow et al. (2013). (b) YZ cut at X = 0 (terminator plane) with color contours of the x component of the flow velocity (V_{x}). The horizontal line at Z = 1.3 R_{M} is color coded by the y component of the convectional electric field (E_{y}), from which the cross polar cap potential is calculated (from Jia et al., 2015). 

In the text 
Fig. 23 A closeup view of the simulated magnetospheric configuration for the Highly Compressed Magnetosphere (HCM) event observed by MESSEGER on 23 November 2011. (a) Color contours of the zcomponent of the magnetic field perturbation (B_{1z} in nT) shown in the XZ plane and also on a 3D sphere of radius ~0.8 R_{M} corresponding to the coremantle boundary. (b) Same as (a) but for the current density in the ydirection (J_{y} in nA/m^{2}), indicating enhanced magnetopause and tail currents as well as the induction currents at the core in response to solar wind compression. In both panels, the black lines with arrows show projections of sampled magnetic field lines onto the XZ plane and the magenta lines with arrows show sampled streamlines of the induction currents generated at the core. The green circle of radius 1 R_{M} represents the planetary surface (from Jia et al., 2019). 

In the text 
Fig. 24 Results from a 3D MHDEPIC simulation of Mercurys magnetosphere. (a) 3D view of the simulated magnetosphere. Colors in the equatorial plane represent plasma density (in amu/cm^{3}), whereas solid lines indicate field lines. The red box shows the embedded PIC region in the magnetotail. (b) A closeup view of the PIC solution showing contours of electron pressure (in Log_{10} nPa) and field lines in the noonmidnight meridian (from Chen et al., 2019c). 

In the text 
Fig. 25 Global TEC maps with 6 h interval under storm (13 October) conditions in 2016. The predicted TEC maps are from our LSTM/NN and corresponding ones from the CODE GIM, along with their differences, are, respectively, given in the left, middle, and right panels of each figure. The predictions shown are 1hour ahead. The unit of the color contour is TECU (10^{1} 6electrons/m^{2}) (from Liu et al., 2020). 

In the text 
Fig. 26 The ML strong flare (M/X class) probability index exhibits a dramatic increase about a day before the strong flare occurs on the Sun (Chen et al., 2019b; Wang et al., 2020). The ML methodology also forecasts the GOES Xray peak intensity of the first strong flare (Jiao et al., 2020). 

In the text 
Fig. A.1 Selfsimilar blocks illustrating the double layer of ghost cells for both coarse and fine blocks (from Gombosi et al., 2003). 

In the text 
Fig. A.2 Solution blocks of the BATSRUS computational mesh with three refinement levels (from Gombosi et al., 2003). 

In the text 
Fig. A.3 The cell update rates as a function of number of cores for the BATSRUS model. The problem size scales in proportion to the number of parallel processes. The dotted lines represent linear scaling. 

In the text 
Fig. B.1 The original physics modules of the SWMF (Toth et al., 2005). This was the first successful coupling (De Zeeuw et al., 2004) of a gyrokinetic ring current model (Harel et al., 1981; Wolf et al., 1991; Sazykin et al., 2002; Toffoletto et al., 2003) to a global MHD model describing the magnetosphere. 

In the text 
Fig. B.2 The layered architecture of the SWMF (from Toth et al., 2012). 

In the text 
Fig. C.1 Synthetic spectra (green) compared to observation (black) taken by Hinode/EIS at 12 November 2007 12:32:02 UT of the Northern coronal hole during the Carrington Rotation 2063 (from Szente et al., 2019). Line profiles of Fe XI 201.734 Å, Fe XIII 201.121 Å and 202.044 Å are closely predicted in intensity, width and Dopplershift; while there is no line in the CHIANTI database between 201.5 and 201.6 Å, which explains the missing peak in the model’s result. 

In the text 
Fig. C.2 Doppler map of line Fe XIII 202.044 Å is obtained from synthetic spectral image using the magnetic boundary of Carrington Rotation 2082, as could be perceived from Earth after the lineofsight integration of individual emission lines (from Szente et al., 2019). Blueshifts show regions where the solar wind is moving toward the observer, redshifts are present where plasma is moving away from the observer. The denser the plasma, the more dominant its effect is on the overall integrated Doppler shift of the observed line. 

In the text 
Fig. D.1 Comparison of synthesized EUV images of the model with observational STEREO A/EUVI images. The columns are from left to right for 171 Å, 195 Å, and 284 Å. Top panels: synthesized EUV images for 2nd order scheme. Middle panels: synthesized EUV images for 5th order scheme. Bottom panels: observational STEREO A/EUVI images. The observation time is 7 March 2011 20:00 UT. 

In the text 
Fig. D.2 The layered software structure of BATSRUS. The arrows point from the module that is using data or methods from the other module. There are multiple versions of the equation and user modules. The various time stepping schemes are independent of the details of the equations being solved (from Toth et al., 2012). 

In the text 
Fig. D.3 The cell and particle update rates as a function of number of cores for the SWMF running the twoway coupled BATSRUS and iPIC3D models. The problem size scales in proportion to the number of parallel processes. The dotted lines represent linear scaling. 

In the text 
Fig. E.1 The overall flow of MHDPIC coupling. At t = 0 the MHD code sends the MHD state inside and around the PIC region to the PIC code. Both the MHD and PIC codes then advance by one or more time steps until both models reach the next coupling time. Information is exchanged both ways, but this time the PIC code only uses the MHD solution as a boundary condition, while the MHD code overwrites its solution in the PIC region with the PIC solution. This process continues until they reach the final simulation time or until the PIC region is removed (after Daldorff et al., 2014). 

In the text 
Fig. E.2 Spatial discretization of the MHDEPIC coupling. The Cartesian grid of the PIC region is indicated with the darker gray area. The lighter gray area shows the ghost cell/node region of the PIC grid. The large red dots are node values obtained from the MHD variables. The small red dots illustrate particles created in the ghost cells of the PIC grid. The small red squares are the ghost cell centers of the PIC grid where the magnetic field is set from the MHD solution. The black dots indicate the MHD cell centers where the solution is obtained from the PIC code. The MHD grid can be either Cartesian or spherical (after Chen et al., 2017). 

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.