© T.I. Gombosi et al., Published by EDP Sciences 2021

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.

Because of society’s reliance on the electrical grid, the internet, high-frequency 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 low-frequency events, they have the potential to cause crippling long-term damage. In fact, the risk of high impact damages due to space weather fits the profile of a market-changing 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 first-principles models requires space physics expertise for the various sub-domains and advanced numerical algorithms. Since the sub-domain 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 high-performance computing is required. Clearly, developing a first-principles space weather model requires sustained multi-disciplinary collaboration of space physicists, applied mathematicians, computer scientists and software engineers.

Presently there are only a couple of physics-based 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 BATS-R-US (Block Adaptive-Tree Solar-wind Roe-type Upwind Scheme) and SWMF. The fundamentals of the BATS-R-US 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 physics-based 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₂, O, O₂, He, Ar, and H. Low-order 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 space-borne 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 large-scale current systems with parametrized empirical functions, and the parameter values are found through least-squares 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 Black-box 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 storm-time) 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 open-access, robust, and effective software tools. Typically, the models are custom-made 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 Physics-based 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. Physics-based 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 low-frequency 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 Hydro-Quebec 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 wide-spread 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 (SIP-AMR-CESE MHD Model), Feng et al., 2014a, b). More general-use 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, MPI-AMRVAC (Keppens et al., 2021), KU Leuven’s European heliospheric forecasting information asset (EUHFORIA, Pomoell & Poedts, 2018) and the University of Michigan’s BATS-R-US (Powell et al., 1999; Toth et al., 2012), model.

3 The origins of BATS-R-US 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 high-order 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, BATS-R-US (Powell et al., 1999). Figure 1 summarizes the present capabilities of BATS-R-US; the algorithms are discussed in detail in Appendix D.

The BATS-R-US (Powell et al., 1999; Toth et al., 2012) is a versatile, high-performance, 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), semi-relativistic (Gombosi et al., 2002), anisotropic (Meng et al., 2012), Hall (Toth et al., 2008), multispecies (Ma et al., 2002), and multi-fluid (Glocer et al., 2009c), extended magnetofluid equations (XMHD) and, most recently, non-neutral multifluid plasmas (Huang et al., 2019). BATS-R-US is used to model several physics domains (see Fig. 2). The efficiency of BATS-R-US is crucial to reach faster than real-time performance with the SWMF while maintaining high resolution in the domains of interest.

In a number of fields in which computer-based modeling of complex, multi-scale, multi-physics 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 high-performance computational capability to simulate the space-weather 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 high-performance computational capability to simulate the space-weather 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 open-source license (https://github.com/MSTEM-QUDA). 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 user-friendly web interface. The user specifies the domains and the driving input parameters, and the CCMC runs-on-request 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 Solar-atmosphere Model (AWSoM/AWSoM-R) and the SWMF/Geospace Model.

AWSoM/AWSoM-R (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⁴ 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 magnetosphere-ionosphere 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/AWSoM-R 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 counter-propagating waves as the underlying cause of the turbulence energy cascade (e.g., Cranmer & Van Ballegooijen, 2010), which transports the energy of turbulence from the large-scale motions across the inertial range of the turbulence spatial scale to short-wavelength perturbations. The latter can be efficiently damped due to wave-particle 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 self-consistently incorporates low-frequency 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 field-aligned Alfvén speed gradients and the vorticity of the background. In addition, the two populations counter-stream 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.

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²” 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²/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 BATS-R-US and it solves the same equations as the solar corona model, but on a Cartesian grid in either co-rotating 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 Threaded-Field-Line Model and AWSoM-R

In the transition region the plasma temperature increases some two orders of magnitude over ~10² km, resulting in a temperature gradient of ~10⁴ K/km. To resolve this gradient, 3D numerical simulations require sub-kilometer 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 Threaded-Field-Line 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 time-dependent simulations. The AWSoM model with TFLM inner boundary conditions is called AWSoM-R, 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 BATS-R-US 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 quasi-dipolar 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 BATS-R-US 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). BATS-R-US supplies near-body field-aligned 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 BATS-R-US to set the plasma tangential velocity at the inner boundary. The RCM receives its initial and boundary field and plasma conditions from BATS-R-US as well as electric field from RIM. It returns total plasma pressure and density to BATS-R-US 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 BATS-R-US and RIM and solves for the energetic electron population in the radiation belts.

The couplings default to 5-second (GM-IE) and 10-second (all other) frequency; faster coupling frequencies are required under extreme driving or when high-frequency output is produced (Welling et al., 2020). While the explicit couplings are shown, the self-consistent nature of multi-model SWMF simulations produces implicit couplings. For example, while region-2 Birkeland currents are not explicitly passed from IM to GM physics modules, the improved pressure gradients in BATS-R-US due to the pressure coupling from RCM drives region-2 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 time-dependent 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. BATS-R-US 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 near-body resolution of 1/4 R_E and ~2 million grid cells with 1/8 R_E maximum resolution, respectively.

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 large-scale 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 coupled-model 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 Biot-Savart 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 post-processing (Rastätter et al., 2014), the SWMF/Geospace combines information from the IE and GM models on-the-fly 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 self-consistent 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 in-situ 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 physics-based and empirical models to predict ground magnetic perturbations (Pulkkinen et al., 2013; Glocer et al., 2016). The physics-based 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 global-scale, meso-scale and micro-scale processes in a self-consistent manner, and integrate these processes to form a truly multi-scale space weather simulation capability. In addition, significant validation efforts have been made by a variety of comparisons with both in-situ and remote-sensing 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 quiet-time 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 STEREO-A (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/AWSoM-R 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 AWSoM-R 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⁻³) 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).

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 pre-event conditions. We also have a SWMF component (EE), which is a physics-based extended MHD model (BATS-R-US) 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 physics-based EE model only works in a stand-alone 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).

Magnetically-driven CMEs were first modeled with the SWMF suite in the early 2000s. First, the distorted spheromac-type Gibson & Low (1998) (GL) unstable flux-rope 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 physics-based Sun-to-Earth 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 in-situ observations at 1 AU and coronagraph images from LASCO C2 and C3. An automated tool, the Eruptive Event Generator using Gibson-Low (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 runs-on-request service to provide CME simulations.

Representative results from EEGGL-driven 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 (model-derived) SOHO/LASCO white light images. The color scale shows the white light total brightness divided by that of the pre-event 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.

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⁴ K and density n = 2 × 10¹⁷ m⁻³. 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⁻³ 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 CME-driven shocks and compressional waves.

In addition to EUV images, our model also allows us to make synthetic Thomson-scattered 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 running-difference images are able to reproduce the observed typical three-part 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 pile-up 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 CME-driven 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.

EEGGL was designed to provide data-drive 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 L1-observed plasma conditions (shown with solid lines) resulting from the Earth-directed CME that occurred on 12 July 2012. Here, time-series data are shown (top to bottom) for the Cartesian components of the magnetic field, mass density and Earth-directed 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 miss-matched 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 storm-driving B_z, which allows the model to successfully drive a magnetospheric simulation, while issues with flux rope rotation and stream-interaction remain to be addressed. This EEGGL-driven simulation was performed on demand at the CCMC where the model outputs are available to the public.

5.4 Solar energetic particle simulations

The acceleration of energetic particles in a CME-driven shock and the subsequent transport processes are modeled using the M-FLAMPA 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). M-FLAMPA 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 BATS-R-US simulations.

Figure 11 shows the application of M-FLAMPA 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 steady-state 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 flux-rope 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 AWSoM-R, EEGGL, and M-FLAMPA 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²/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.

Figure 11 demonstrates the capability of using the self-consistent physics-based 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 gyro-motion of the energetic particles must be employed. The effect of the gyro-motion 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 iso-surface in geospace.

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 time-backward 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 mid-latitude region.

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), BATS-R-US, when run with high spatial resolution in key portions of the geospace, can easily resolve the Kelvin-Helmholtz instability and flux-transfer events (FTEs) (Kuznetsova et al., 2009), found e.g., at the magnetospheric boundary. High-resolution 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 high-resolution simulation is shown in Figure 15. The SWMF/BATS-R-US 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).

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 re-assess 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 near-equatorial 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 wave-particle 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 BATS-R-US when configured with multi-fluid or multi-species 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 field-aligned currents than those derived from spaceborne magnetic measurements by the AMPERE project (Anderson et al., 2000).

Figure 17 shows simulated field-aligned currents (FACs) and ground magnetic perturbations due to a solar wind pressure enhancement. The two-way coupled BATSRUS, CRCM, and IE modules are utilized in this run with 1/8 R_E resolution used from the dayside magnetopause to the near-Earth 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.

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 compression-induced meso-scale field-aligned currents. The results of this coupled geospace run were then used to drive the Global Ionosphere and Thermosphere Model (GITM), and revealed a short-lived meso-scale 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 meso-scale 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 SYM-H 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. Dst-equivalent outputs have long been a staple of SWMF simulations, serving as a quick-look 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 K-index 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 SYM-H 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 substorm-related 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 real-time SWMF runs at CCMC. The right panel of Figure 18 shows scatter plots of observed and simulated Dst (hourly averaged SYM-H) 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.

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 CCMC-led 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, to-date, the most accurate reproduction of these indices from a solar-wind-to-ionosphere first-principles physics-based model of the full geospace system.

6 Resolving kinetic scales in global simulations

6.1 MHD-EPIC and MHD-AEPIC

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 hyperbolic-parabolic cleaning in addition to the 8-wave scheme is necessary to eliminate accumulation of $\bar{\nabla }\cdot \bar{\mathbf{B}}$ errors near the boundaries of the PIC region.

Eventually our work yielded results: the MHD with embedded PIC (MHD-EPIC) 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 multi-species and multi-fluid MHD, as well as the five- and six-moment fluid equations (cf. Chen & Toth, 2019; Zhou et al., 2020).

Running MHD-EPIC simulations for long simulation times also revealed hidden issues with the PIC algorithms that could be avoided in stand-alone PIC simulations by careful tuning of various parameters, but were plaguing the more complicated MHD-EPIC simulations. We overcame these issues by developing the Gauss Law satisfying ECSIM (GL-ECSIM) 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 MHD-EPIC 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⁴ ~ 65,000. With such scaling it became possible to simulate Earth’s magnetosphere with the MHD-EPIC 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 MHD-EPIC 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 X-lines, 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 (MHD-AEPIC) algorithm. The main idea comes from the block-adaptive mesh and the hybrid schemes used in BATS-R-US: 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 GL-ECSIM algorithm into the Adaptive Mesh Particle Simulator (AMPS) code (Tenishev et al., 2021). The resulting MHD-AEPIC model can achieve an order of magnitude or even more speed-up compared to the MHD-EPIC model that uses static PIC domains. We also developed a new PIC code, the Flexible Exascale Kinetic Simulator (FLEKS), to be used in MHD-AEPIC. FLEKS is based on the AMReX library (Zhang et al., 2019b, 2020) and it was designed for flexibility and high performance with a state-of-the-art semi-implicit PIC algorithm. Novel particle splitting and merging algorithms have been designed for FLEKS to control the number of macro-particles per cell during long MHD-AEPIC 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.

6.2 MHD-EPIC results

MHD-EPIC 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 mini-magnetosphere of Ganymede (Toth et al., 2016; Zhou et al., 2019, 2020). To demonstrate the capabilities of MHD-EPIC, 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 rope-like 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 one-hour simulation, flux ropes form near the subsolar point and move toward the poles quasi-periodically.

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 MHD-EPIC 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 crescent-like shape, as is shown in Figure 21c. The crescent ions drift in a positive y direction because E_x is negative.

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; Alvarado-Gó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, semi-relativistic MHD, multi-fluid MHD, MHD with anisotropic pressure, or high-order 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 mass-loading processes (ionization, charge-exchange, 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 magnetosphere-ionosphere 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 MHD-EPIC 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 noon-midnight 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.

Figure 22b shows a YZ cut at X = 0 through the simulation in which the color contours represent the x-component 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 $V_x$ at high latitudes are the crosspolar cap flows moving in the anti-sunward direction, while the flows with positive $V_x$ 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 time-dependent 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 z-component of the magnetic field perturbations, B_1z, which result from various current systems, including the Chapman-Ferraro 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 stand-off 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 Chapman-Ferraro 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 cross-tail 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.

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 MHD-EPIC model to simulate Mercury’s magnetosphere so that reconnection can be treated using a physics-based 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 MHD-EPIC 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 tailward-moving and planetward-moving plasmoids are found in the simulation (Fig. 24b). The modeled plasmoids in the near-Mercury 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).

While the imposed solar wind flow is symmetric about the Sun–Mercury line, due to kinetic effects significant dawn-dusk asymmetries develop in the magnetotail in the MHD-EPIC 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 MHD-EPIC simulation also exhibit a strong preference for the dawn side, such that almost all dipolarization fronts and high-speed 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 MHD-EPIC is capable of capturing both local and global physics of the magnetosphere, and, therefore, provides an excellent tool for studying solar wind-magnetosphere 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 user-friendly 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 IRI-2016 and NeQuick-2 models are around 9.21/5.5 TECU, respectively (Liu et al., 2020). Moreover, typical large-scale 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).

The IGS GIM maps are constructed based on a few hundred IGS stations with limited spatial resolution. Critical meso-scale 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 multi-scale 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 physics-based 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 cutting-edge 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/HMI-based SHARP parameters, physically-insightful summaries of active region photospheric magnetic fields. Thus, while this ML-based approach adds no new physics, it represents SOLSTICE’s first step toward improving event forecasts with cross-disciplinary efforts.

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 24-h 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

Physics-based 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 well-resolved 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 line-of-sight 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 physics-based 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 Open-source Development

Earlier in this paper we described how a large interdisciplinary group of researchers have developed, with sustained effort, the first-principles 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, BATS-R-US and the SWMF have been available for runs-on-request 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 open-source for a while. The Global Ionosphere Thermosphere Model (GITM) has been open-source 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 open-source too.

Finally, the core SWMF model was also released in 2020 under a non-commercial open-source license at https://github.com/mstem-quda. MSTEM-QUDA contains the full core of the SWMF and the BATS-R-US, 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 MSTEM-QUDA repository is an up-to-date mirror of the repositories developed at Michigan. Both the master and stable branches are available.

Making a major part of the SWMF truly open-source 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 person-year effort in this project. Maintaining and developing a world-class 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, BATS-R-US. We also outlined its history and future directions. Today, SWMF, BATS-R-US 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 PRE-EVENTS grant 1663800 and NSF SWQU grant PHY-2027555. 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 BATS-R-US

The BATS-R-US 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 8-Wave Riemann Solver

The first step of the BATS-R-US development was to attack a fundamental roadblock to high-resolution magnetohydrodynamics (MHD) simulations. High-resolution schemes are based on the conservative form of the governing equations:

$$\frac{\mathrm{ }\mathrm{\partial }\bar{\mathbf{U}}}{\mathrm{\partial }t}+{\left(\bar{\nabla }\cdot \overline{\overline{\mathbf{F}}}\right)}^T=\overline{\mathbf{S}}$$(A.1)

where $\bar{\mathbf{U}}$ is the vector of conserved quantities defined by

$$\bar{\mathbf{U}}={\left(\rho ,\rho u_x,\rho u_y,\rho u_z,B_x,B_y,B_z,\epsilon \right)}^T$$(A.2)

where ρ is mass density, u_x, u_y and u_z are the three components of the plasma bulk flow velocity vector, $\bar{\mathbf{u}}$, while $\bar{\mathbf{B}}=\{B_x,B_y,B_z\}$ is the magnetic field vector and ε is the total energy density

$$\epsilon ={\epsilon }_{\mathrm{hd}}+\frac{\bar{\mathbf{B}}\cdot \bar{\mathbf{B}}}{2{\mu }_0}=\frac{p}{\gamma -1}+\rho \frac{\bar{\mathbf{u}}\cdot \bar{\mathbf{u}}}{2}+\frac{\bar{\mathbf{B}}\cdot \bar{\mathbf{B}}}{2{\mu }_0}.$$(A.3)

Here ε_hd denotes the hydrodynamic energy density, while γ is the specific heat ratio and μ₀ is the permeability of vacuum. The flux tensor, $\overline{\overline{\mathbf{F}}}$, can be written as

$$\overline{\overline{\mathbf{F}}}={\left(\begin{array}{c}\rho \bar{\mathbf{u}}\\ \rho \bar{\mathbf{u}}\bar{\mathbf{u}}+\left(p+\frac{\bar{\mathbf{B}}\cdot \bar{\mathbf{B}}}{2{\mu }_0}\right)\overline{\overline{\mathbf{I}}}-\frac{\bar{\mathbf{B}}\bar{\mathbf{B}}}{{\mu }_0}\\ \begin{array}{c}\bar{\mathbf{u}}\bar{\mathbf{B}}-\bar{\mathbf{B}}\bar{\mathbf{u}}\\ \bar{\mathbf{u}}\left(\epsilon +p+\frac{\bar{\mathbf{B}}\cdot \bar{\mathbf{B}}}{2{\mu }_0}\right)-\frac{\left(\bar{\mathbf{u}}\cdot \bar{\mathbf{B}}\right)\bar{\mathbf{B}}}{{\mu }_0}\end{array}\end{array}\right)}^T.$$(A.4)

Finally, $\bar{\mathbf{S}}$ is a “source” vector, containing the terms that cannot be expressed in divergence form:

$$\bar{\mathbf{S}}=-\bar{\nabla }\cdot \bar{\mathbf{B}}\left(\begin{array}{l}0\\ \bar{\mathbf{B}}\\ \bar{\mathbf{u}}\\ \bar{\mathbf{u}}\cdot \bar{\mathbf{B}}\end{array}\right).$$(A.5)

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 $\bar{\mathbf{S}}\equiv 0$ identity. However, setting $\bar{\mathbf{S}}$ 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 high-resolution 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 $\bar{\nabla }\cdot \bar{\mathbf{B}}$ 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 so-called 8-wave 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 high-resolution schemes. Toth (2000) carefully evaluated the various methods that constrain $\bar{\nabla }\cdot \bar{\mathbf{B}}$ and concluded that the 8-wave 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 8-wave 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 8-wave 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

BATS-R-US uses a simple and effective block-based adaptive mesh refinement (AMR) technique (Stout et al., 1997). The approach closely follows that first developed for two-dimensional gas dynamics calculations by Berger & Jameson (1985); Berger & Colella (1989). This block-based 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 volume-averaged solution quantities within computational cells. The computational cells are embedded in regular structured blocks of equal-sized cells. The blocks are geometrically self-similar. 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 self-similar blocks. Although each block within a grid has the same data-storage 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 over-resolved regions, the refinement process reverses; eight children coarsen and coalesce into a single parent block. Thus, cell resolution reduces by a factor of 2. Multigrid-type 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).

The hierarchical data structure and self-similar blocks simplify domain decomposition and enable good load-balancing, a crucial element for truly scalable computing. For explicit time stepping (all blocks use the same time-step) natural load-balancing occurs by distributing the blocks equally among the processors. For more complicated time-stepping schemes the load-balancing 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 Peano-Hilbert or Morton space-filling curves to minimize inter-processor communication. The self-similar nature of the solution blocks also means that serial performance enhancements apply to all blocks and that fine-grained algorithm parallelization is possible. The algorithm’s parallel implementation is so pervasive that even the grid adaptation performs in parallel.

A.3 BATS-R-US 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 block-based 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 self-similar blocks make domain decomposition of the problem almost trivial and readily enable good load-balancing, 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 BATS-R-US on several supercomputers from 8 up to more than 200,000 cores. Recently, BATS-R-US has been further developed to use hybrid MPI and OpenMP parallelism that allows scaling to beyond 500,000 cores (Zhou & Tóth, 2020).

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 BATS-R-US: 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 early-to-mid 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 BATS-R-US 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.

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 stand-alone 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.

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 BATS-R-US 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 BATS-R-US 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 BATS-R-US) 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 co-rotating 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. BATS-R-US solves the ideal or two-temperature MHD equations on a Cartesian grid in either co-rotating 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. BATS-R-US 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 bounce-averaged Boltzmann equation.

B.3.7 Ionospheric Electrodynamics (IE)

The Ionospheric Electrodynamics model is a two dimensional height-integrated 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 field-aligned 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 multi-species hydrodynamics with viscosity, thermal conduction, chemical reactions, ion-neutral 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 BATS-R-US, the SWMF includes several world-class models that provide one of the most advanced space weather simulation capabilities ranging from kinetic and meso-scales 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 Gibson-Low 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 M-FLAMPA

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 AWSoM-R 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 M-FLAMPA is self-consistently calculated from the energy densities of the Alfvénic turbulence simulated by AWSoM-R. Together, M-FLAMPA, AWSoM-R and EEGGL provide a high-performance, self-consistent 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 high-performance 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 cut-cells for discretizing the simulated domain. The implemented cut-cells 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 two-phase 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 stand-alone modeling tool or coupled to other components of the SWMF.

AMPS is coupled to several components of the SWMF, allowing multi-scale and multi-physics simulations. Specifically, two-way coupling has been developed with the Global Magnetosphere (GM) and Outer Heliosphere (OH) modules. A one-way 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 high-performance implicit Particle-in-Cell (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 time-steps 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 particle-in-cell (PIC) code that is designed for the MHD with adaptively embedded PIC (MHD-AEPIC) simulations. FLEKS uses the Gauss’s Law satisfying Energy-Conserving Semi-Implicit Method (GL-ECSIM) (Chen & Toth, 2019) as the base PIC solver. Novel particle splitting and merging algorithms have been designed to control the number of macro-particles per cell during a long MHD-AEPIC simulation. The particle splitting algorithm improves statistical representation and reduces noise in the cells with low macro-particle 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 high-performance parallel data structures provided by the AMReX library (Zhang et al., 2019b, 2020) to store the fields and also the particles.

Appendix C: Physics

BATS-R-US 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 BATS-R-US 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, photo-ionization, recombination, radiative losses, etc. The boundary and initial conditions vary greatly as well.

C.1 Conservation Laws in BATS-R-US

BATS-R-US can be configured to solve the governing equations of ideal and resistive MHD (Powell et al., 1999), semi-relativistic (Gombosi et al., 2002), anisotropic (Meng et al., 2012), Hall (Toth et al., 2008), multispecies (Ma et al., 2002) and multi-fluid (Glocer et al., 2009c) extended magnetofluid equations (XMHD) and more recently non-neutral multifluid plasmas (Huang et al., 2019).

Table C.1 summarizes the various extended MHD conservation laws that can be solved by BATS-R-US.

C.1.1 Extended MHD Equations

BATS-R-US can solve many approximations to the low-order 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

$$\frac{\mathrm{\partial }{\rho }_s}{\mathrm{\partial }t}+\bar{\nabla }\cdot \left({\rho }_s{\bar{\mathbf{u}}}_s\right)=0$$(C.1a)

$$\frac{\mathrm{\partial }{\rho }_s{\bar{\mathbf{u}}}_s}{\mathrm{\partial }t}+\bar{\nabla }\cdot \left[{\rho }_s{\bar{\mathbf{u}}}_s {\bar{\mathbf{u}}}_s+p_{s_{\perp }}\overline{\overline{\mathbf{I}}}+(p_{s_{\parallel }}-p_{s_{\perp }}) \bar{\mathbf{b}} \bar{\mathbf{b}}\right]=\frac{q_s}{m_s}{\rho }_s(\bar{\mathbf{E}}+{\bar{\mathbf{u}}}_s\times \bar{\mathbf{B}})$$(C.1b)

$$\frac{\mathrm{\partial }{p}_{s_{\parallel }}}{\mathrm{\partial }t}+\bar{\nabla }\cdot \left(p_{s_{\parallel }}{\bar{\mathbf{u}}}_s\right)=-2p_{s_{\parallel }}\bar{\mathbf{b}}\cdot (\bar{\mathbf{b}}\cdot \bar{\nabla }) {\bar{\mathbf{u}}}_s$$(C.1c)

$$\frac{\mathrm{\partial }{p}_{s_{\perp }}}{\mathrm{\partial }t}+\bar{\nabla }\cdot \left(p_{s_{\perp }}{\bar{\mathbf{u}}}_s\right)=-p_{s_{\perp }}(\bar{\nabla }\cdot {\bar{\mathbf{u}}}_s)+p_{s_{\perp }}\bar{\mathbf{b}}\cdot (\bar{\mathbf{b}}\cdot \bar{\nabla }) {\bar{\mathbf{u}}}_s$$(C.1d)

where ρ and $\bar{\mathbf{u}}$ 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): $\overline{\overline{\mathbf{P}}}=p_{\perp }\overline{\overline{\mathbf{I}}}+(p_{\parallel }-p_{\perp }) \bar{\mathbf{b}} \bar{\mathbf{b}}$, where $\overline{\overline{\mathbf{I}}}$ is the identity matrix, $\bar{\mathbf{b}}$ is the unit vector along the magnetic field direction, $p_{\mathrm{P}}$ 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. BATS-R-US 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):

$$\begin{array}{ll}\frac{\mathrm{\partial }\langle {\mathcal{M}}_s{\mathcal{F}}_s\rangle }{\mathrm{\partial }t}& +\bar{\nabla }\cdot \left({\bar{\mathbf{u}}}_s\langle {\mathcal{M}}_s{\mathcal{F}}_s\rangle \right)+\bar{\nabla }\cdot \langle {{\bar{\mathbf{c}}}_s\mathcal{M}}_s{\mathcal{F}}_s\rangle +\langle {\mathcal{F}}_s \left[\right({\bar{\mathbf{c}}}_s\cdot \bar{\nabla }\left){\bar{\mathbf{u}}}_s\right]\cdot {\bar{\nabla }}_{c_s}{\mathcal{M}}_s\rangle \\ & +\left(\frac{\mathrm{\partial }{\bar{\mathbf{u}}}_s}{\mathrm{\partial }t}+\left[\right({\bar{\mathbf{u}}}_s\cdot \bar{\nabla }\left){\bar{\mathbf{u}}}_s\right]\right)\cdot \langle {\mathcal{F}}_s {\bar{\nabla }}_{c_s}{\mathcal{M}}_s\rangle -\langle {\mathcal{F}}_s ({\bar{\mathbf{a}}}_s\cdot {\bar{\nabla }}_{c_s}){\mathcal{M}}_s\rangle =\langle {\mathcal{M}}_s{\left(\frac{\delta {\mathcal{F}}_s}{\delta t}\right)}_{\mathrm{coll}}\rangle \end{array}$$(C.2)

Here ${\mathcal{F}}_s$ is the velocity distribution function of species “s” (expressed in terms of the random velocity, ${\bar{\mathbf{c}}}_s$), ${\mathcal{M}}_s$ 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 ${\mathcal{M}}_s$ in the random velocity defines the order of the velocity moment equation. For instance, ${\mathcal{M}}_s=m_s$ (zeroth-order moment equation) results in the continuity equation, describing the conservation of mass. The first-order velocity moment equations are obtained by using ${\mathcal{M}}_s=m_s{\bar{\mathbf{c}}}_s$ and they express the conservation of momentum. The second-order velocity moment equations are obtained by using ${\mathcal{M}}_s=m_s {\bar{\mathbf{c}}}_s {\bar{\mathbf{c}}}_s$. There is one zero-order, three first-order and six second-order moment equations (due to the symmetric nature of the ${\bar{\mathbf{c}}}_s {\bar{\mathbf{c}}}_s$ 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), $\bar{\nabla }\cdot \langle {\bar{\mathbf{c}}}_s {\mathcal{M}}_s{\mathcal{F}}_s\rangle$. If ${\mathcal{M}}_s$ is n-th order in velocity, the term $\langle {\bar{\mathbf{c}}}_s {\mathcal{M}}_s{\mathcal{F}}_s\rangle$ is the (n + 1)-th velocity moment. In other words, the transport equation for the n-th 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 third-order velocity moments (the heat flow), or express a high-order 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 ${\bar{\mathbf{u}}}_s$ and the two pressure components, $p_{\parallel }$ and $p_{\perp }$. For this reason this is called the six moment approximation.

The electric ($\bar{\mathbf{E}}$) and magnetic fields ($\bar{\mathbf{B}}$) are obtained from Maxwell’s equations:

$$\frac{\mathrm{\partial }\bar{\mathbf{B}}}{\mathrm{\partial }t}+\bar{\nabla }\times \bar{\mathbf{E}}=0$$(C.3a)

$$\frac{\mathrm{\partial }\bar{\mathbf{E}}}{\mathrm{\partial }t}-c^2 \bar{\nabla }\times \bar{\mathbf{B}}=-c^2{\mu }_0 \bar{\mathbf{j}}$$(C.3b)

$$\bar{\nabla }\cdot \bar{\mathbf{E}}=\frac{{\rho }_c}{{ϵ}_0}$$(C.3c)

$$\bar{\nabla }\cdot \bar{\mathbf{B}}=0$$(C.3d)

where ϵ₀ is the vacuum permittivity, μ₀ is the vacuum permeability, $c=1/\sqrt{{ϵ}_0{\mu }_0}$ is the speed of light, ${\rho }_c={\sum }_s\mathrm{}(q_s/m_s){\rho }_s$ is the total charge density and $\bar{\mathbf{j}}={\sum }_s\mathrm{}(q_s/m_s) {\rho }_s\overline{{\mathbf{u}}_s}$ 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. BATS-R-US 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 $\bar{\mathbf{j}}$:

$$\bar{\mathbf{j}}=e\left(\sum _{\mathrm{s}=\mathrm{ions}}\mathrm{}{Z}_s n_s{\bar{\mathbf{u}}}_s-n_e{\bar{\mathbf{u}}}_s\right)=e n_e ({\bar{\mathbf{u}}}_{+}-{\bar{\mathbf{u}}}_e)$$(C.4)

where Z_s is the ionization state of a given ion species and ${\bar{\mathbf{u}}}_{+}={\sum }_{\mathrm{s}}\mathrm{}{Z}_s(n_s/n_e) {\bar{\mathbf{u}}}_s$ is the charge averaged ion velocity. Note that in general ${\bar{\mathbf{u}}}_{+}\ne \bar{\mathbf{u}}$ and the two vectors can be quite different. BATS-R-US takes into account the full definition of ${\bar{\mathbf{u}}}_{+}$ and thus self-consistently 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:

$$\bar{\mathbf{E}}=-{\mathbf{u}}_e\times \bar{\mathbf{B}}-\frac{1}{e{n}_e}\bar{\nabla }\cdot \left[p_{e_{\perp }}\overline{\overline{\mathbf{I}}}+(p_{e_{\parallel }}-p_{e_{\perp }})\bar{\mathbf{b}} \bar{\mathbf{b}}\right].$$(C.5)

In a single-ion plasma the electron velocity is ${\bar{\mathbf{u}}}_e={\bar{\mathbf{u}}}_i-\bar{\mathbf{j}}/\left(e{n}_e\right)$ 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

$$E_{\parallel }=\bar{\mathbf{b}}\cdot \bar{\mathbf{E}}=-\frac{{\bar{\nabla }}_{\parallel }{p}_{e_{\parallel }}}{e{n}_e}+\frac{p_{e_{\parallel }}-p_{e_{\perp }}}{e{n}_e}\frac{{\bar{\nabla }}_{\parallel }B}{B}$$(C.6)

where ${\bar{\nabla }}_{\parallel }=\bar{\mathbf{b}}\cdot \bar{\nabla }$ 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. BATS-R-US 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 BATS-R-US equation set is the six-moment 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 six-moment 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 BATS-R-US we consider the following processes:

• elastic collisions,

• photoionization and impact ionization (using the Beer-Lambert 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, ${\bar{\mathbf{u}}}_n$.

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 relaxation-time 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 phase-space 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):

$${\left(\frac{\delta F_s}{\delta t}\right)}_{\mathrm{el}}=\sum _{t=\mathrm{all}}\mathrm{}\frac{F_{s\left(\mathrm{st}\right)}-F_s}{{\tau }_{\mathrm{st}}}.$$(C.7)

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 steady-state 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 ion-electron collisions. In practice, the momentum transfer collision frequency, ${\overline{\nu }}_{\mathrm{st}}$ is used instead of the relaxation time. The momentum transfer collision frequency includes a mass-dependent factor that accounts for the efficiency of momentum transfer in an elastic collision:

$$\frac{1}{{\tau }_{\mathrm{st}}}=\frac{m_s+m_t}{m_t}{\overline{\nu }}_{\mathrm{st}}.$$(C.8)

The parameters of the Maxwellian, $F_{s\left(\mathrm{st}\right)}$, are chosen in a way that mass, momentum and energy are conserved while the gas is driven toward equilibrium (Burgers, 1969; Gombosi, 1994):

$$F_{s\left(\mathrm{st}\right)}=n_s{\left(\frac{m_s}{2\pi k_{\mathrm{B}}{T}_{s\left(\mathrm{st}\right)}}\right)}^{3/2}\times \exp \left[-\frac{m_s}{2k_{\mathrm{B}}{T}_{s\left(\mathrm{st}\right)}}{\left({\bar{\mathbf{v}}}_s+{\bar{\mathbf{u}}}_s-{\bar{\mathbf{u}}}_{\mathrm{st}}\right)}^2\right]$$(C.9)

where

$$n_s=\int \mathrm{} \underset{\mathrm{\infty }}{\int }\mathrm{} \int \mathrm{}{F}_s(t,\bar{\mathbf{r}},{\bar{\mathbf{v}}}_s) d^3{v}_s$$(C.10)

$${\bar{\mathbf{u}}}_{\mathrm{st}}=\frac{m_t{\bar{\mathbf{u}}}_t+m_s{\bar{\mathbf{u}}}_s}{m_s+m_t}$$(C.11)

$$T_{s\left(\mathrm{st}\right)}=T_s+\frac{m_s{m}_t}{(m_s+m_t{)}^2}\left[2(T_t-T_s)+\frac{m_t}{3k_{\mathrm{B}}}({\bar{\mathbf{u}}}_t-{\bar{\mathbf{u}}}_s{)}^2\right].$$(C.12)

In these expressions, k_B is the Boltzmann constant and the kinetic temperature is defined as

$$T_s=\frac{m_s}{3n_s{k}_{\mathrm{B}}}\int \mathrm{} \underset{\mathrm{\infty }}{\int }\mathrm{} \int \mathrm{}{v}_s^2{F}_s(t,\bar{\mathbf{r}},{\bar{\mathbf{v}}}_s) d^3{v}_s$$(C.13)

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 $u_{\mathrm{st}}=u_{\mathrm{ts}}$ and in general $T_{s\left(\mathrm{st}\right)}\ne T_{t\left(\mathrm{ts}\right)}$.

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:

$${\left(\frac{\delta F_s}{\delta t}\right)}_{\mathrm{ion}}={\nu }_{s\mathrm{\prime }}^{\mathrm{io}}{n}_{s\mathrm{\prime }} {\delta }^3\left({\bar{\mathbf{u}}}_s+{\bar{\mathbf{v}}}_s-{\bar{\mathbf{u}}}_n\right)$$(C.14)

where ${\nu }_{s\mathrm{\prime }}^{\mathrm{io}}$ is the sum of the photoionization and impact ionization frequencies of species “s′”, $n_{s\mathrm{\prime }}$ 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 ${\mathrm{O}}^{+}+\mathrm{H}\rightleftharpoons \mathrm{O}+{\mathrm{H}}^{+}$ 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 phase-space distribution function of particles “s” is the following:

$${\left(\frac{\delta F_s}{\delta t}\right)}_{\mathrm{cx}}=-F_s\left(\sum _{t\mathrm{\prime }=\mathrm{neutrals}}\mathrm{}{k}_{\mathrm{st}\mathrm{\prime }}{n}_{t\mathrm{\prime }}\right)+\left(\sum _{t=\mathrm{ions}}\mathrm{}{k}_{\mathrm{ts}\mathrm{\prime }} n_t\right) n_{s\mathrm{\prime }}{\delta }^3\left({\bar{\mathbf{u}}}_s+{\bar{\mathbf{v}}}_s-{\bar{\mathbf{u}}}_n\right).$$(C.15)

Here $k_{\mathrm{ts}\mathrm{\prime }}$ and $k_{\mathrm{ts}}$ 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”:

$${\left(\frac{\delta F_s}{\delta t}\right)}_{\mathrm{rec}}=-{\alpha }_s{n}_e{F}_s$$(C.16)

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:

$$\frac{\delta F_s}{\delta t}=-\left(\sum _{t=\mathrm{all}}\mathrm{}\frac{m_s+m_t}{m_t}{\overline{\nu }}_{\mathrm{st}}+\sum _{t\mathrm{\prime }=\mathrm{neutrals}}\mathrm{}{k}_{\mathrm{st}\mathrm{\prime }}{n}_{t\mathrm{\prime }}+{\alpha }_s{n}_e\right)F_s+\left({\nu }_{s\mathrm{\prime }}^{\mathrm{io}}+\sum _{t=\mathrm{ions}}\mathrm{}{k}_{\mathrm{ts}\mathrm{\prime }} n_t\right)n_{s\mathrm{\prime }}{\delta }^3\left({\bar{\mathbf{u}}}_s+{\bar{\mathbf{v}}}_s-{\bar{\mathbf{u}}}_n\right)+\sum _{t=\mathrm{all}}\mathrm{}\frac{m_s+m_t}{m_t}{\overline{\nu }}_{\mathrm{st}}{F}_{s\left(\mathrm{st}\right)}.$$(C.17)

In BATS-R-US 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 X-ray 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 counter-propagating 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):

$$\frac{\mathrm{\partial }{w}_{\pm }}{\mathrm{\partial }t}+\bar{\nabla }\cdot \left[(\bar{\mathbf{u}}\pm {\bar{\mathbf{V}}}_{\mathrm{A}})w_{\pm }\right]+\frac{w_{\pm }}{2}\left(\bar{\nabla }\cdot \bar{\mathbf{u}}\right)=-{\mathrm{\Gamma }}_{\pm }{w}_{\pm }\mp \mathcal{R}\sqrt{w_{-}{w}_{+}}$$(C.18)

where V_A is the Alfvén velocity, while Γ_± and $\mathcal{R}$ 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 ($P_{{\mathrm{A}}_{\odot }}$), and (ii) the transverse correlation length of Alfvén turbulence (L_⊥). Our solar corona model assumes that $P_{{\mathrm{A}}_{\odot }}\sim B_{\odot }$ and $L_{\perp }\sim B_{\odot }^{-1/2}$ (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):

$$n_e+n_{\alpha }=\sum _i\mathrm{}{n}_i$$(C.19)

$$n_e{u}_e+n_{\alpha }{u}_{\alpha }=\left(\sum _i\mathrm{}{n}_i{u}_i-\frac{j}{e}\right).$$(C.20)

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 two-stream 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:

$$\frac{\beta }{\sqrt{E}}\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }t}+\mu \frac{\mathrm{\partial }\varphi }{\mathrm{\partial }s}-\frac{1-{\mu }^2}{2}\left(-\frac{F}{E}+\frac{1}{B}\frac{\mathrm{\partial }B}{\mathrm{\partial }s}\right)\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }\mu }+\mathrm{EF\mu }\frac{\mathrm{\partial }}{\mathrm{\partial }E}\left(\frac{\varphi }{E}\right)=q+\langle S\rangle$$(C.21)

where the constant β = 1.7 × 10⁻⁸ eV^1/2 cm⁻¹ 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 two-stream approach (cf., Nagy & Banks, 1970; Solomon et al., 1988; Lummerzheim & Lilensten, 1994; Barakat & Schunk, 2001; Schunk & Nagy, 2009). The two-stream 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 wave-particle 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 gyro-averaged particle equation of motion given by:

$$m_i\frac{\mathrm{\partial }{v}_{i_{\parallel }}}{\mathrm{\partial }t}-q_i{E}_{\parallel }+\frac{G m_i{M}_{\mathrm{planet}}}{r^2}+{\mu }_{a_i}\frac{\mathrm{\partial }B}{\mathrm{\partial }s}=0$$(C.22)

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, ${\mu }_{a_i}$ is the particle’s first adiabatic invariant specified by

$${\mu }_{a_i}=\frac{m_i{v}_{\perp }^2}{2B}.$$(C.23)

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 $\bar{\mathbf{E}}\times \bar{\mathbf{B}}$ drift and refilling via ionospheric outflow at mid- and low-latitudes. 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:

$$\frac{\mathrm{\partial }N}{\mathrm{\partial }t}=\mathcal{S}-\mathcal{L}-{\bar{\mathbf{u}}}_{\perp }\cdot \bar{\nabla }N$$(C.24)

${\bar{\mathbf{u}}}_{\perp }$ is the horizontal bulk velocity of the cold plasmasphere fluid (set by the local $\bar{\mathbf{E}}\times \bar{\mathbf{B}}$ drift). A dipole magnetic field is assumed, electric potential is a required input and typically obtained via an empirical model. $\mathcal{S}$ and $\mathcal{L}$ 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:

$$\mathcal{S}=\frac{N_S-N\left(t\right)}{{\tau }_{\mathrm{fill}}}.$$(C.25)

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:

$$\mathcal{L}=\frac{N\left(t\right)}{{\tau }_{\mathrm{loss}}}.$$(C.26)

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 first-principles-based 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 current-plasmasphere 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 self-consistent 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 assumed-isotropic populations, each with a given energy invariant, λ, such that,

$$\lambda =W{\mathcal{V}}^{2/3}$$(C.27)

where W is the kinetic energy of the particles within the population and $\mathcal{V}$ 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,

$$\frac{\mathrm{\partial }\eta }{\mathrm{\partial }t}-{\bar{v}}_{\mathrm{D}}\cdot \bar{\nabla }\eta =0$$(C.28)

where η is a function of time and space and is defined by,

$$\eta =\int \frac{n\mathrm{d}s}{B}=n\mathcal{V}.$$(C.29)

Here, ${\bar{v}}_{\mathrm{D}}$ is the full electromagnetic drift velocity of the population and is given by,

$${\bar{v}}_{\mathrm{D}}=\frac{\bar{\mathbf{E}}\times \bar{\mathbf{B}}}{B^2}+\frac{\lambda \bar{\mathbf{B}}\times \bar{\nabla }{\mathcal{V}}^{-2/3}}{q{B}^2}$$(C.30)

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, $\mathcal{V}$, 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,

$${p\mathrm{\prime }}_{\mathrm{GM}}^{}=p_{\mathrm{GM}}+\min \left(1,\frac{\mathrm{d}t}{{\tau }_{\mathrm{couple}}}\right)(p_{\mathrm{RCM}}-p_{\mathrm{GM}})$$(C.31)

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 e-folding rate. It should be noted that stand-alone 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 bounce-averaged 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 bounce-averaged distribution function as a function of azimuth, radial distance, energy and pitch angle for several species,

$$\begin{array}{ll}\langle \frac{\mathrm{d}{Q}_s}{\mathrm{d}t}\rangle =& \frac{\mathrm{\partial }{Q}_s}{\mathrm{\partial }t}+\frac{1}{R_0^2}\frac{\mathrm{\partial }}{\mathrm{\partial }{R}_0}\left(R_0^2\langle \frac{\mathrm{d}{R}_0}{\mathrm{d}t}{Q}_s\rangle \right)+\frac{\mathrm{\partial }}{\mathrm{\partial }\varphi }\left(\langle \frac{\mathrm{d}\varphi }{\mathrm{d}t}\rangle Q_s\right)\\ & +\frac{1}{\sqrt{E}}\frac{\mathrm{\partial }}{\mathrm{\partial }E}\left(\sqrt{(}E)\langle \frac{\mathrm{d}E}{\mathrm{d}t}\rangle Q_s\right)+\frac{1}{h\left(\mu \right) \mu }\frac{\mathrm{\partial }}{\mathrm{\partial }\mu }\left(h\left(\mu \right) \mu \langle \frac{\mathrm{d}\mu }{\mathrm{d}t}\rangle Q_s\right)=\langle \frac{\mathrm{d}{Q}_s}{\mathrm{d}t}{\rangle }_{\mathrm{loss}}\end{array}$$(C.32)

where angle brackets denote average values over the bounce period, Q_s is the distribution function for species s, R₀ 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,

$$h\left(\mu \right)=\frac{1}{2R_0}{\int }_a^b\mathrm{}\frac{\mathrm{d}s}{\sqrt{1-B\left(s\right)/B_m}}$$(C.33)

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 RAM-like 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 Self-Consistent Magnetic field (RAM-SCB) codes. These models all fall into the same class of bounce-averaged 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 large-scale 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 self-consistent 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 self-consistent 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 non-dipolar 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 BATS-R-US.

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 BATS-R-US configured with single fluid MHD, anisotropic MHD, and multi-fluid 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):

$$\frac{\mathrm{\partial }{F}_s}{\mathrm{\partial }t}+\frac{\langle {\dot{\lambda }}_i\rangle F_s}{\mathrm{\partial }{\lambda }_i}+\frac{\langle {\dot{\varphi }}_i\rangle F_s}{\mathrm{\partial }{\varphi }_i}=S_i$$(C.34)

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 wave-particle 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 self-consistent. Meng et al. (2013) demonstrated that the magnetic field calculated in BATS-R-US using the coupled CIMI pressure is consistent with a force balance. Similarly, the pressure feedback from CIMI drives currents in BATS-R-US that contribute the ionospheric potential and hence the convection in CIMI.

RAM with Self-Consistent Magnetic field (RAM-SCB) 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 3-D force balance magnetic field model (SCB, (Zaharia et al., 2004; Zaharia, 2008). This model balances the $\bar{\mathbf{j}}\times \bar{\mathbf{B}}$ 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). RAM-SCB 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)

$$\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mu v\frac{\mathrm{\partial }f}{\mathrm{\partial }s}+(\bar{\mathbf{u}}\cdot \bar{\nabla })f+\frac{\mathrm{d}p}{\mathrm{d}t}\frac{\mathrm{\partial }f}{\mathrm{\partial }p}+\frac{\mathrm{d}\mu }{\mathrm{d}t}\frac{\mathrm{\partial }f}{\mathrm{\partial }\mu }-\frac{\mathrm{\partial }}{\mathrm{\partial }\mu }\left(D_{\mu \mu }\frac{\mathrm{\partial }f}{\mathrm{\partial }\mu }\right)=Q,$$(C.35)

where f is the particle’s distribution function, μ is the cosine of the particle’s pitch angle, D_μμ is the pitch angle diffusion coefficient, $\bar{\mathbf{u}}$ 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

$$\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+(\bar{\mathbf{u}}\cdot \bar{\nabla })f-\frac{1}{3}(\bar{\nabla }\cdot \bar{\mathbf{u}}\frac{\mathrm{\partial }f}{\mathrm{\partial }\ln p}=\bar{\nabla }\cdot (\overline{\overline{\mathbf{D}}}\cdot \bar{\nabla }f),$$(C.36)

where $\overline{\overline{\mathbf{D}}}=D\mathbf{bb}$ 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

$$\frac{\mathrm{df}}{\mathrm{dt}}+\frac{1}{3}\frac{d\ln \rho }{\mathrm{dt}}\frac{\mathrm{\partial }f}{\mathrm{\partial }\ln p}=B\frac{\mathrm{\partial }}{\mathrm{\partial }s}\left(\frac{D}{B}\frac{\mathrm{\partial }f}{\mathrm{\partial }s}\right),$$(C.37)

where ρ(s,t) and B(s,t) are plasma density and total magnetic field magnitude along the magnetic field lines. The 3-D problem is then reduced to a set of independent 1-D problems on those time-evolving Lagrangian grids (Sokolov et al., 2004). In M-FLAMPA the plasma and turbulence parameters along the magnetic field lines are extracted dynamically from the the BATS-R-US 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 BATS-R-US 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

$$\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+(\bar{\mathbf{u}}\cdot \bar{\nabla })f+\left({\bar{\mathbf{v}}}_{\mathrm{D}}\cdot \bar{\nabla }\right)f-\frac{1}{3}(\bar{\nabla }\cdot \bar{\mathbf{u}})\frac{\mathrm{\partial }f}{\mathrm{\partial }\ln p}=\bar{\nabla }\cdot (\overline{\overline{\mathbf{D}}}\cdot \bar{\nabla }f)Q,$$(C.38)

where f(r,p,t) is the omnidirectional distribution function of GCRs, p the particle’s momentum, r the heliocentric distance, $\bar{\mathbf{u}}$ the solar wind velocity, $\overline{\overline{\mathbf{D}}}$ the symmetric part of the diffusion tensor, ${\bar{\mathbf{v}}}_{\mathrm{D}}$ the pitch-angle 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 out-streaming solar wind, (3) two types of drifts: the gradient-curvature 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 Biot-Savart integral for currents in the magnetosphere

The contribution from the magnetospheric current system to the ground magnetic field is given by the Biot-Savart integral (neglecting the displacement current):

$${\bar{\mathbf{B}}}_{\mathrm{M}}\left({\bar{\mathbf{x}}}_0\right)=\frac{{\mu }_0}{4\pi }\underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}\bar{\mathbf{J}}\left(\bar{\mathbf{x}}\right)\times \frac{{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^3} \mathrm{d}V, \bar{\mathbf{J}}\left(\bar{\mathbf{x}}\right)=\frac{1}{{\mu }_0}{\bar{\nabla }}_{\bar{\mathbf{x}}}\times \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right),$$(C.39)

where ${\bar{\mathbf{x}}}_0$ is the point where we calculate the synthetic magnetic field, ${\bar{\mathbf{B}}}_{\mathrm{M}}\left({\bar{\mathbf{x}}}_0\right)$ is the contribution to this field from the magnetosphere currents, R_M is radius of the magnetosphere-ionosphere boundary (in our case the inner boundary of the magnetosphere model, R_M ≈ 2.5 R_E in most simulations) and $\mathrm{d}V={\mathrm{d}}^3\bar{\mathbf{x}}$ 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 ${\bar{\nabla }}_{\bar{\mathbf{x}}}$ and ${\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}$ differentiate over coordinates $\bar{\mathbf{x}}$ and ${\bar{\mathbf{x}}}_0$, respectively. The gradient of the inverse distance function is

$${\bar{\nabla }}_{\bar{\mathbf{x}}}\frac{1}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}|}=-{\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}\frac{1}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}|}=\frac{{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^3}.$$(C.40)

Now the Biot-Savart integral (Eq. (C.39)) can be written as

$${\bar{\mathbf{B}}}_{\mathrm{M}}\left({\bar{\mathbf{x}}}_0\right)=\frac{1}{4\pi }{\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}\times \underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}\frac{{\bar{\nabla }}_{\bar{\mathbf{x}}}\times \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)}{\left|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}\right|}\mathrm{d}V.$$(C.41)

Next, we expand the double vector product using the ${\bar{\nabla }}_{\bar{\mathbf{x}}}{\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^{-1}={\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}{\bar{\nabla }}_{\bar{\mathbf{x}}}|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^{-1}$, ${\bar{\nabla }}_{\bar{\mathbf{x}}}\cdot \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)=0$, and ${\bar{\nabla }}_{\bar{\mathbf{x}}}^2|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^{-1}=-4\pi {\delta }^3({\bar{\mathbf{x}}}_0-\bar{\mathbf{x}})$ identities:

$${\bar{\mathbf{B}}}_{\mathrm{M}}\left({\bar{\mathbf{x}}}_0\right)=\frac{1}{4\pi }\underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}{\bar{\nabla }}_{\bar{\mathbf{x}}}\left[\bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)\cdot {\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}\left(\frac{1}{\left|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}\right|}\right)\right] \mathrm{d}V-\frac{1}{4\pi }\underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}{\bar{\nabla }}_{\bar{\mathbf{x}}}\cdot \left[\bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right){\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}\left(\frac{1}{\left|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}\right|}\right)\right] \mathrm{d}V-\frac{1}{4\pi }\underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}{\bar{\nabla }}_{\bar{\mathbf{x}}}\cdot \left[{\bar{\nabla }}_{{\bar{\mathbf{x}}}_0}\left(\frac{1}{\left|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}\right|}\right)\bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)\right] \mathrm{d}V+\underset{\left|\bar{\mathbf{x}}\right|>R_{\mathrm{M}}}{\int }\mathrm{}\bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right){\delta }^3\left({\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}\right)\mathrm{d}V.$$(C.42)

If the observation point, ${\bar{\mathbf{x}}}_0$, 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 ${\bar{\mathbf{n}}}_{\mathrm{R}}=\bar{\mathbf{x}}/\left|\bar{\mathbf{x}}\right|$ (note, that positive ${\bar{\mathbf{n}}}_{\mathrm{R}}$ points into the domain of integration) and apply Gauss’ theorem to equation (C.42):

$${\bar{\mathbf{B}}}_{\mathrm{M}}\left({\bar{\mathbf{x}}}_0\right)=\underset{\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}}{\int }\mathrm{}\left[{\sigma }_m\left(\bar{\mathbf{x}}\right)\frac{{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^3}+\frac{{\mu }_0}{4\pi }{\mathbf{i}}_m\left(\bar{\mathbf{x}}\right)\times \frac{{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}}{|{\bar{\mathbf{x}}}_0-\bar{\mathbf{x}}{|}^3}\right] \mathrm{d}S,$$(C.43)

where dS is an area element on the spherical surface, $\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}$. Note that equation (C.43) replaces the effect of all currents in the simulated magnetosphere with a surface current, ${\mathbf{i}}_m\left(\bar{\mathbf{x}}\right)=-{\mathbf{n}}_{\mathrm{R}}\times \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)/{\mu }_0$ and magnetic surface charge, ${\sigma }_m\left(\bar{\mathbf{x}}\right)=-{\mathbf{n}}_{\mathrm{R}}\cdot \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)/4\pi$. For the special case when ${\bar{\mathbf{x}}}_0=0$ equation (54) reduces to

$${\bar{\mathbf{B}}}_{\mathrm{M}}\left(0\right)=\underset{\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}}{\int }\mathrm{}\frac{\bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)\mathrm{d}S}{4\pi R_{\mathrm{M}}^2}.$$(C.44)

This result agrees with the well-known 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 field-aligned currents

Another source of geomagnetic variations, which is particularly significant at high geomagnetic latitudes, is the magnetic field produced by the currents connecting the magnetosphere-ionosphere boundary, $\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}$ to the ionosphere and closing there. The currents in the gap region, between $\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}$ and the ionosphere, $R_{\mathrm{I}}\le \left|{\bar{\mathbf{x}}}^{\mathrm{\prime }}\right|\le R_{\mathrm{M}}$ are field-aligned, which means that the current density is parallel or anti-parallel to the terrestrial magnetic field: $\bar{\mathbf{J}}\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right)\left|\right| {\bar{\mathbf{B}}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right)$, where ${\bar{\mathbf{x}}}^{\mathrm{\prime }}$ 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 M-I boundary, dS, a flux tube may be traced to the ionosphere boundary. The total current enclosed by this flux tube is

$$\mathrm{d}I=\frac{1}{{\mu }_0}{\mathbf{n}}_{\mathrm{R}}\cdot \left[{\bar{\nabla }}_{\bar{\mathbf{x}}}\times \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)\right] \mathrm{d}S.$$(C.45)

The field line (and current line) is described by $\mathrm{d}{\bar{\mathbf{x}}}^{\mathrm{\prime }}/\mathrm{d}\mathcal{l}=\pm {\bar{\mathbf{b}}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right)$, where $\mathrm{d}\mathcal{l}$ is the path length element, and ${\mathbf{b}}_0={\bar{\mathbf{B}}}_0/\left|{\bar{\mathbf{B}}}_0\right|$. Expressing $\mathrm{d}\mathcal{l}=\pm \mathrm{d}{R}^{\mathrm{\prime }}/[{\mathbf{n}}_{R^{\mathrm{\prime }}}\cdot {\mathbf{b}}_0({\bar{\mathbf{x}}}^{\mathrm{\prime }}\left)\right]$ in terms of the element of radial coordinate in the gap region, $R^{\mathrm{\prime }}=\sqrt{|{\bar{\mathbf{x}}}^{\mathrm{\prime }}{|}^2}$, one can express the equation for the current line in terms of dR′:

$$\frac{\mathrm{d}{\bar{\mathbf{x}}}^{\mathrm{\prime }}}{\mathrm{d}{R}^{\mathrm{\prime }}}=\frac{{\mathbf{b}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right)}{{\mathbf{n}}_{R^{\mathrm{\prime }}}\cdot {\mathbf{b}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right)}.$$(C.46)

The magnetic field line given by equation (C.46) should be integrated from R_M down to R_I, starting from each point $\bar{\mathbf{x}}$ at the M-I interface. This way a multitude of field lines, ${\bar{\mathbf{x}}}^{\mathrm{\prime }}(\bar{\mathbf{x}},R^{\mathrm{\prime }})$ are obtained in the $R_{\mathrm{I}}\le R^{\mathrm{\prime }}\le R_{\mathrm{M}}$ domain. Equation (C.43) can now be generalized to account for the Biot-Savart integral from the multitude of field aligned currents, dI, in the gap region,

$${\bar{\mathbf{B}}}_{\mathrm{Gap}}\left({\bar{\mathbf{x}}}_0\right)=\frac{1}{4\pi }\underset{\left|\bar{\mathbf{x}}\right|=R_{\mathrm{M}}}{\int }\mathrm{}{\mathbf{n}}_{\mathrm{R}}\cdot \left[{\bar{\nabla }}_{\bar{\mathbf{x}}}\times \bar{\mathbf{B}}\left(\bar{\mathbf{x}}\right)\right]\underset{R_{\mathrm{I}}}{\overset{R_{\mathrm{M}I}}{\int }}\mathrm{}\frac{{\mathbf{b}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right(\bar{\mathbf{x}},R^{\mathrm{\prime }}\left)\right)}{{\mathbf{n}}_{R^{\mathrm{\prime }}}\cdot {\mathbf{b}}_0\left({\bar{\mathbf{x}}}^{\mathrm{\prime }}\right(\bar{\mathbf{x}},R^{\mathrm{\prime }}\left)\right)}\times \frac{{\bar{\mathbf{x}}}_0-{\bar{\mathbf{x}}}^{\mathrm{\prime }}(\bar{\mathbf{x}},R^{\mathrm{\prime }})}{|{\bar{\mathbf{x}}}_0-{\bar{\mathbf{x}}}^{\mathrm{\prime }}(\bar{\mathbf{x}},R^{\mathrm{\prime }}){|}^3}\mathrm{d}{R}^{\mathrm{\prime }}\mathrm{d}S.$$(C.47)

Note, that the integral over dR′ in equation (C.47) is a complicated vector function. However, its value only depends on $\bar{\mathbf{x}}$ and ${\bar{\mathbf{x}}}_0$, 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 M-I interface. We also note that only the derivatives of tangential components of the magnetospheric field at the M-I interface contribute to the radial component of the current density. Therefore, only the magnetic field at the M-I 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 Biot-Savart integral along magnetic field lines going through the gap region can be written as the field-aligned 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 Biot-Savart integral to obtain the perturbation magnetic field at a surface point, ${\bar{\mathbf{x}}}_0$, produced by the flux tube current, I, for a flux tube described by ${\bar{\mathbf{x}}}^{\mathrm{\prime }}({\bar{\mathbf{x}}}_{\mathrm{M}},R^{\mathrm{\prime }})$:

$$\delta =\bar{\mathbf{B}}\left({\bar{\mathbf{x}}}_0\right)=\frac{{\mu }_{0}I}{4\pi }\underset{R_{\mathrm{I}}}{\overset{R_{\mathrm{M}}}{\int }}\mathrm{}\left[\frac{\mathrm{d}{\bar{\mathbf{x}}}^{\mathrm{\prime }}({\bar{\mathbf{x}}}_{\mathrm{M}},R^{\mathrm{\prime }})}{\mathrm{d}{R}^{\mathrm{\prime }}}\times \frac{{\bar{\mathbf{x}}}^{\mathrm{\prime }}({\bar{\mathbf{x}}}_{\mathrm{M}},R^{\mathrm{\prime }})-{\bar{\mathbf{x}}}_0}{\left|{\bar{\mathbf{x}}}^{\mathrm{\prime }}\right({\bar{\mathbf{x}}}_{\mathrm{M}},R^{\mathrm{\prime }})-{\bar{\mathbf{x}}}_0{|}^3}\right]\mathrm{d}{R}^{\text{'}}$$(C.48)

where the magnetic field line passing through the point ${\bar{\mathbf{x}}}_{\mathrm{M}}$ at the M-I interface is parameterized in the gap region with the radial coordinate, $R^{\mathrm{\prime }}=\left|{\bar{\mathbf{x}}}^{\mathrm{\prime }}\right|$ ($R_{\mathrm{I}}\le \left|{\bar{\mathbf{x}}}^{\mathrm{\prime }}\right|\le R_{\mathrm{M}}$, and $\mathrm{d}{\bar{\mathbf{x}}}^{\mathrm{\prime }}=(\mathrm{d}{\bar{\mathbf{x}}}^{\mathrm{\prime }}/\mathrm{d}{R}^{\mathrm{\prime }}) \mathrm{d}{R}^{\mathrm{\prime }}$).

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 ${\bar{\mathbf{e}}}_{\mathrm{M}}$ and on the two perpendicular directions. The solution is

$${\bar{\mathbf{x}}}^{\mathrm{\prime }}\left({\bar{\mathbf{x}}}_{\mathrm{M}},R^{\mathrm{\prime }}\right)=\left\{{\bar{\mathbf{x}}}_{\mathrm{M}}+{\bar{\mathbf{e}}}_{\mathrm{M}}\left[\mathrm{sgn}\left(z_{\mathrm{M}}\right)\sqrt{\frac{R_{\mathrm{M}}^3}{R^{\mathrm{\prime }}}+z_{\mathrm{M}}^2-R_{\mathrm{M}}^2}-z_{\mathrm{M}}\right]\right\}{\left(\frac{R^{\mathrm{\prime }}}{R_{\mathrm{M}}}\right)}^{\frac{3}{2}}.$$(C.49)

where we used the notation, $z_{\mathrm{M}}={\bar{\mathbf{x}}}_{\mathrm{M}}\cdot {\bar{\mathbf{e}}}_{\mathrm{M}}$. This expression can be further simplified if we use one of the standard geocentric coordinate systems with the z-axis 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 Biot-Savart integral of all currents flowing in the GM component (typically, BATS-R-US at Earth, see Sect. C.5.1). Kp and AE indexes leverage virtual magnetometer stations to more closely reflect the calculation of their real-world counterparts.

Real-world Kp is the average of local-K index values calculated from 13 mid-latitude ground observatories, rounded to the nearest third. Local-K is a range index of the maximum minus the minimum disturbance at a single observatory over set three-hour windows (0–3 UT, 3–6 UT, etc), scaled to a 9-value integer via a semi-logarithmic 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 three-hour windows are not possible. Instead, a rolling three-hour window that ends at the current simulation time is used. Therefore, SWMF virtual Kp range windows line up exactly with the real-world index when the universal time hour is a multiple of 3. No seasonal adjustments are made for virtual Kp. The use of the 24-station approach provides a minor improvement in predictive performance compared to using the real-world stations and scalings.

The Auroral Electrojet (AE) family of indexes are high-latitude data products. The real-world indexes are the product of the geomagnetic north-south 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 in-situ, local measurements of the solar corona, all information of the global corona measured by ground-based observatories (visible and IR) or via spacecrafts carrying remote-sensing 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 narrow-band imaging (line-of-sight, LOS) or synthetic spectra (SPECTRUM).

C.7.1 Line-of-Sight Images

SWMF has the capability of generating various synthetic line-of-sight (LOS) plots, such as EUV images. The response $\mathcal{R}$ of each pixel of the image is treated as a LOS integral of a function f through the plasma:

$$\mathcal{R}=\int \mathrm{}f\left(\mathcal{l}\right) \mathrm{d}\mathcal{l}$$(C.50)

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 tri-linear 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). Doppler-shifted, nonthermal line broadening due to low-frequency 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 post-processing tool within the SWMF: it processes output after a simulation is completed. It is a stand-alone 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 Cartesian-grid or the BATS-R-US 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 λ₀ with line width Δλ:

$$\varphi \left(\lambda \right)=\frac{1}{\sqrt{2\pi }\Delta \lambda }{\mathrm{e}}^{-{\left(\frac{\lambda -{\lambda }_0}{2\Delta \lambda }\right)}^2}.$$(C.51)

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:

$$F=\frac{1}{4\pi d^2}{\int }_V\mathrm{}N\left(X_j^{+m}\right)A_{\mathrm{ji}}h{\nu }_{\mathrm{ij}}\mathrm{d}V,$$(C.52)

where $N\left(X_j^{+m}\right)$ is the density of the emitting $X_j^{+m}$ 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):

$$F=\frac{1}{4\pi d^2}{\int }_V\mathrm{}G\left(T_e,N_e\right)N_e^2\mathrm{d}V,$$(C.53)

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 contribution-function values and the BATS-R-US output saved from an AWSoM simulation. The following calculations are performed on a line-by-line basis, for each cell along the line-of-sight. First we apply Doppler-shift to the line center:

$${\lambda }_{\mathrm{shifted}}=\left(1-\frac{u_{\mathrm{LOS}}}{c_{\mathrm{light}}}\right){\lambda }_0,$$(C.54)

where u_LOS is the line-of-sight 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 non-thermal component:

$${\Delta \lambda }^2=\Delta {\lambda }_{\mathrm{instrument}}^2+{\lambda }^2\frac{u_{\mathrm{th}}^2+u_{n\mathrm{th}}^2}{c_{\mathrm{light}}^2}.$$(C.55)

Thermal broadening is calculated considering the contribution along the line-of-sight direction of the anisotropic temperature, the non-thermal component is due to the low-frequency Alfvén wave contribution along the line-of-sight.

The SPECTRUM code output is synthetic spectra or synthetic spectral images. An example is shown in Figures C.1 and C.2.

Appendix D: Algorithms

Continuous development of the numerical algorithms is a necessity in order to maintain state-of-the-art 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

$$\epsilon =k(\Delta x{)}^n$$(D.1)

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 $(\Delta x{)}^{-4}$ for a three-dimensional simulation, because the number of grid cells is $\propto (\Delta x{)}^{-3}$ 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, BATS-R-US was born with block-adaptive mesh refinement. This algorithm allows refining the grid where necessary, and coarsen it where possible. Using block-based adaptation instead of cell-based adaptation (or fully unstructured grids) has distinct advantages for high performance massively parallel codes. In the past 20 years, the original block-adaptive 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 non-Cartesian grids, allows grid adaptation based on geometric and physics-based 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 non-Cartesian 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 non-uniform grid spacing by applying some non-linear 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 BATS-R-US used a second order total variation diminishing (TVD) scheme, which was state-of-the-art 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 BATS-R-US 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 high-order scheme, we had the following requirements:

1. Minimal oscillations near discontinuities.

2. Conservative scheme that gives correct jump conditions.

3. High order at grid resolution changes and high order for non-Cartesian grids.

4. Small stencil to allow for small grid blocks.

5. Only moderately more expensive than the second order TVD scheme.

6. General method that is high order for various system of equations including non-linear 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 multi-dimensional 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 non-Cartesian 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 Runge-Kutta scheme, the 5th order scheme is only about three times more expensive than the two-stage 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).

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 Courant-Friedrich-Lewy (CFL) condition:

$$\Delta t<C\frac{\Delta x}{c_{\max }}$$(D.2)

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 BATS-R-US 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, BATS-R-US 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 point-implicit 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 semi-implicit. 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 BATS-R-US 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 point-implicit 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 semi-implicit 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 BATS-R-US 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.

A particularly interesting application of advanced time stepping algorithms is for particle-in-cell (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 semi-implicit particle-in-cell 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 semi-implicit method (ECSIM) developed by Lapenta (2017). Energy conservation is crucial to maintain the long-term 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. GL-ECSIM 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 multi-scale and multi-physics system. A similar approach is used in BATS-R-US 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 physics-based 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. non-conservative energy equation, semi-relativistic 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 BATS-R-US and the SWMF were designed to achieve high performance on massively parallel super computers. BATS-R-US uses a block-adaptive 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 BATS-R-US 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 BATS-R-US 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 block-adaptive 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 multi-threading 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, BATS-R-US 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 BATS-R-US 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 hyper-threading 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 BATS-R-US 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.

Appendix E: MHD-EPIC and MHD-AEPIC

The magnetohydrodynamics with embedded particle-in-cell (MHD-EPIC) 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 (MHD-AEPIC) 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 BATS-R-US 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 self-consistently, our coupled MHD-EPIC and MHD-AEPIC 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 BATS-R-US 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 BATS-R-US – iPIC3D coupling via SWMF is genuinely 2-way: all physical quantities are self-consistently advanced by both codes and the relevant information is fully exchanged in every time-step.

MHD-EPIC (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 macro-particles 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 macro-particles 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 two-way 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 MHD-EPIC (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 macro-particles.

The SWMF and the BATS-R-US 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, BATS-R-US 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