We present a community effort to assess how open science can advance heliophysics and space weather modeling. Open science has the potential to enhance the quality and pace of scientific discovery, but its application to scientific modeling requires more careful consideration with respect to open data and open software guidelines, as complex scientific models are not ordinary software. We gathered feedback from modeling teams worldwide through a living survey and discussion sessions at the Open Science Workshop in College Park, USA, in 2024, and the COSPAR ISWAT Initiative Working Meeting in Cape Canaveral, USA, in 2025. We complement these findings with lessons learned from almost 25 years of experience at the Community Coordinated Modeling Center in enabling open use of models. We identify key roadblocks in current open science practices and guidelines and offer recommendations for future progress. Our findings are organized into four overlapping themes: open use of models and simulation results, open validation, open development, and open collaboration. An essential outcome of the discussion is the need for model developers and users to speak with a united voice and promote the role of models in future open science efforts. We introduce a new cross-domain community initiative called Heliophysics Open Modeling Environment (HOME), which will be integrated as an overarching activity within the COSPAR ISWAT Initiative. HOME will serve as a platform for modelers and model users to work together, facilitate community modeling, improve the scientific return on modeling investment, and advance innovation in heliophysics and space weather.

Abstract Space weather poses an important but under‐quantified threat to society. While severe geomagnetic storms are recognized as potential global catastrophes, their socio‐economic impacts remain poorly quantified. We present a novel physics‐engineering‐economic framework that links geophysical drivers to power grid geoelectric fields, transformer vulnerability, and macroeconomic consequences. Using the United States as an example, we estimate daily economic losses for a 250‐year geomagnetic storm from transformer thermal heating of 1.81 billion USD (95 percent confidence interval: 1.65 to 1.96 billion USD), disrupting power for approximately 5.1 million people and 135,000 businesses. These estimates are conservative lower bounds, reflecting only transformer thermal heating effects and excluding voltage collapse, cascading failures, and restoration costs. The true societal risk could be substantially higher. Nonetheless, this contribution provides the first nationwide end‐to‐end coupling from space physics to potential macroeconomic loss, with quantified uncertainties. Our results demonstrate that coupled socio‐economic modeling of space weather is feasible and essential, and the framework is scalable and transferable, offering a template for assessing space weather risk to critical infrastructure in other countries.

The US National Aeronautics and Space Administration (NASA) DRIVE Center for Geospace Storms (CGS) hosted a workshop for graduate students conducting research in geospace science in November 2024. The research relevant to geospace has traditionally been separated into the Magnetosphere, including Magnetosphere-Ionosphere coupling, and the Ionosphere-Thermosphere-Mesosphere disciplines. This is apparent in sections of the American Geophysical Union and the funding structure of the United States National Science Foundation (NSF) and NASA. Following that structure, the geospace science community is divided into sub-communities separated along the discipline boundaries. However, it is now widely accepted that geospace is a tightly coupled system whose parts are strongly interconnected. The future Heliophysics workforce, particularly those focused on geospace, must have an understanding of how the regions of geospace influence and are influenced by each other. The goal of the workshop hosted by CGS was to bring together graduate students to improve their understanding across all regions of geospace with an emphasis on how to use geospace modeling to address compelling research questions throughout the domain. We report on the content of the workshop and the evaluation of achievement of these goals.

Cloud computing has gained substantial momentum across diverse applications in recent years, notably in scientific computing, collaborative research, and large-scale machine learning operations. Its integration of data and code within a unified system facilitates swift data transfer and sharing among various research groups. However, despite its prominence in research, cloud computing usage in education is still limited beyond computer science courses. Embracing this technological shift presents an opportunity for graduate students and early-career researchers to familiarize themselves with these tools, contributing to open research and facilitating global collaboration. In this work, we explore from a user perspective the use of cloud computing in two NASA projects, particularly the Center for Geospace Storms (CGS) and Heliocloud, shedding light on how these initiatives can benefit the scientific community. By bridging higher education with academic and research environments through workshops and tutorials, these efforts can play a pivotal role in educating the next generation of researchers.

The adiabatic invariants (M, J, Φ) and the related invariants (M, K, L*) have been established as effective coordinate systems for describing radiation belt dynamics at a theoretical level, and through numerical techniques, can be paired with in situ observations to order phase‐space density. To date, methods for numerical techniques to calculate adiabatic invariants have focused on empirical models such the Tsyganenko models TS05, T96, and T89. In this work, we develop methods based on numerical integration and variable step size iteration for the calculation of adiabatic invariants, applying the method to the Lyon‐Fedder‐Mobarry (LFM) global magnetohydrodynamics (MHD) simulation code, with optional coupling to the Rice Convection Model (RCM). By opening the door to adiabatic invariant modeling with MHD magnetic fields, the opportunity for exploratory modeling work of radiation belt dynamics is enabled. Calculations performed using LFM are cross‐referenced with the same code applied to the T96 and TS05 Tsyganenko models evaluated on the LFM grid. Important aspects of the adiabatic invariant calculation are reviewed and discussed, including (a) sensitivity to magnetic field model used, (b) differences in the problem between quiet and disturbed geomagnetic states, and (c) the selection of key parameters, such as the magnetic local time step size for drift shell determination. The rigorous development and documentation of this algorithm additionally acts as preliminary step for future thorough reassessment of in situ phase‐space density results using alternative magnetic field models.

The reprocessing of radiation belt electron flux measurements into phase space density (PSD) as a function of the adiabatic invariants is a widely-used method to address major questions regarding electron energization and loss in the outer radiation belt. In this reprocessing, flux measurements j (α, E) at local pitch angles α, energies E, and optionally magnetometer measurements B, are combined with a global magnetic field model to express the phase space density f (L*) in terms of the third invariant Φ ∝ 1/L* at fixed first and second invariants M and K. While the general framework of the calculation is agreed upon, implementation details vary amongst the literature, and the issue of magnetic field model dependence is rarely addressed. This work reviews the steps of the calculation with lists of commonly used implementation options. For the first time, analysis is presented to display the effect of doing the calculation with different implementation options and with different backing models (including both empirical and MHD-driven models). The results are summarized to inform evaluation of existing results and future efforts calculating and analyzing radiation belt electron phase space density. Three events are analyzed, and while differences are found, the primary structural interpretations of the phase space density analysis exhibit model independence.

Space weather poses an important but under-quantified threat to critical infrastructure, the economy, and society. While extreme geomagnetic storms are recognized as potential global catastrophes, their socio-economic impacts remain poorly quantified. Here we present a novel physics-engineering-economic framework that links geophysical drivers of geomagnetic storms to power grid geoelectric fields, transformer vulnerability, and macroeconomic consequences. Using the United States as an example, we estimate daily economic losses from transformer thermal heating of 1.37 billion USD (95 percent confidence interval: 1.16 to 1.58 billion USD) for a 100-year geomagnetic storm, with power outages affecting 4.1 million people and 101,000 businesses. A 250-year event could raise losses to 2.09 billion USD per day (95 percent confidence interval: 1.84 to 2.34 billion USD), disrupting power for more than 6 million people and 155,000 businesses. Crucially, the framework is scalable and transferable, offering a template for assessing space weather risk to critical infrastructure in other countries. This integrative approach provides the first end-to-end quantification of space weather socio-economic impacts, bridging space physics through to policy-relevant metrics. Our results demonstrate that coupled socio-economic modeling of space weather is both feasible and essential, enabling evidence-based decision making in infrastructure protection and global risk management.

Space weather risk assessment is constrained by the lack of available asset information needed to model geomagnetically induced currents (GICs) in electricity transmission infrastructure. We propose a systematic method that enables risk analysts to collect their own open-source substation data. Using a web browser platform for annotation, we convert OpenStreetMap (OSM) substation locations into high-resolution, component-level mappings of electricity transmission assets. We convert an initial 1,313 high-voltage (>=230 kV) substations to 52,273 components using low-altitude, satellite, and Street View imagery accessed through Google Earth, identifying 7,949 transformers. Compared to the OSM baseline, this approach provides detailed insights on voltage levels and substation configurations. We then construct a geospatial GIC network for the Tennessee Valley Authority (TVA) region, comparing May 2024 results with the University of Illinois Urbana-Champaign 150-bus (UIUC150) synthetic network and with measured ground GICs at 13 monitoring devices. The transformer types at unannotated substations and the grounding resistances are unknown, so we sample both across a Monte Carlo ensemble. This gives a median TVA 95th-percentile peak ground GIC of 29.1 A, with a 90 percent confidence interval of 23.9-36.1 A. The UIUC150 network yields a 95th-percentile peak ground GIC of 35.8 A under the same forcing, falling within this interval, and the modeled time series broadly capture the temporal morphology of the geomagnetic storm at the monitoring sites. This method shows promise for spatially explicit, screening-level GIC assessment without requiring access to operator data.

The formation of the stormtime ring current is a result of the inward transport and energization of plasma sheet ions. Previous studies have demonstrated that a significant fraction of the total inward plasma sheet transport takes place in the form of bursty bulk flows (BBFs), known theoretically as flux tube entropy-depleted “bubbles.’ However, it remains an open question to what extent bubbles contribute to the buildup of the stormtime ring current. Using the Multiscale Atmosphere Geospace Environment (MAGE) Model, we present a case study of the March 17, 2013 storm, including a quantitative analysis of the contribution of plasma transported by bubbles to the ring current. We show that bubbles are responsible for at least 50\% of the plasma energy enhancement within 6 R$_E$ during this strong geomagnetic storm. The bubbles that penetrate within 6 R$_E$ transport energy primarily in the form of enthalpy flux, followed by Poynting flux and relatively little as bulk kinetic flux. Return flows can transport outwards a significant fraction of the plasma energy being transported by inward flows, and therefore must be considered when quantifying the net contribution of bubbles to the energy buildup. Data-model comparison with proton intensities observed by the Van Allen Probes show that the model accurately reproduces both the bulk and spectral properties of the stormtime ring current. The evolution of the ring current energy spectra throughout the modeled storm is driven by both inward transport of an evolving plasma sheet population and by charge exchange with Earth’s geocorona.

Supporting information for "The contribution of plasma sheet bubbles to stormtime ring current buildup and evolution of the energy composition" by Sciola et al. 2023 submitted to the Journal of Geophysical Research: Space Physics

Thermospheric mass density perturbations are commonly observed during geomagnetic storms and fundamental to upper atmosphere dynamics, but the sources of these perturbations are not well understood. Large neutral density perturbations during storms greatly affect the drag experienced by low Earth orbit. We investigated the thermospheric density perturbations at all latitudes observed along the CHAMP and GRACE satellite trajectories during the August 24–25, 2005 geomagnetic storm. Observations show that large neutral density enhancements occurred not only at high latitudes, but also globally. Large density perturbations were seen in the equatorial regions away from the high‐latitude, magnetospheric energy sources. We used the high‐resolution Multiscale Atmosphere Geospace Environment (MAGE) model to simulate consecutive neutral density changes observed by satellites during the storm. The MAGE simulation, which resolved mesoscale high‐latitude convection electric fields and field‐aligned currents, and included physics‐based specification of auroral precipitation, was contrasted with a standalone ionosphere‐thermosphere simulation driven by a high‐latitude electrodynamics empirical model. The comparison demonstrates that first‐principles representations of highly dynamic and localized Joule heating events in a fully coupled whole geospace model is critical to accurately capture both generation and propagation of traveling atmospheric disturbances (TADs) that produce neutral density perturbations globally. The MAGE simulation shows that larger density peaks in the equatorial region observed by CHAMP and GRACE are the result of TADs generated at high‐latitudes in both hemispheres, and intersect at low‐latitudes. This study reveals the importance of investigating thermospheric density variations at all latitudes in a fully coupled geospace model with sufficiently high resolving power.

Explosive magnetotail activity has long been understood in the context of its auroral manifestations. While global models have been used to interpret and understand many magnetospheric processes, the temporal and spatial scales of some auroral forms have been inaccessible to global modeling creating a gulf between observational and theoretical studies of these phenomena. We present here an important step toward bridging this gulf using a newly developed global magnetosphere-ionosphere model with resolution capturing less than or similar to 30 km azimuthal scales in the auroral zone. In a global magnetohydrodynamic (MHD) simulation of the growth phase of a synthetic substorm, we find the self-consistent formation and destabilization of localized magnetic field minima in the near-Earth magnetotail. We demonstrate that this destabilization is due to ballooning-interchange instability which drives earthward entropy bubbles with embedded magnetic fronts. Finally, we show that these bubbles create localized field-aligned current structures that manifest in the ionosphere with properties matching observed auroral beads.

We study electron injection and energization by bursty bulk flows (BBFs), by tracing electron trajectories using magnetohydrodynamic (MHD) field output from the Lyon-Fedder-Mobarry (LFM) code. The LFM MHD simulations were performed using idealized solar wind conditions to produce BBFs. We show that BBFs can inject energetic electrons of few to 100 keV from the magnetotatail beyond -24 R-E to inward of geosynchronous, while accelerating them in the process. We also show the dependence of energization and injection on the initial relative position of the electrons to the magnetic field structure of the BBF, the initial pitch angle, and the initial energy. In addition, we have shown that the process can be nonadiabatic with violation of the first adiabatic invariant (mu). Further, we discuss the mechanism of energization and injection in order to give generalized insight into the process.

It is now over three decades since the first paper was published using the code that has come to be known as LFM (Lyon-Fedder-Mobarry). The code, used since then extensively in heliophysics research, had a number of novel features: eighth-order centered spatial differencing, the Partial Donor Cell Method limiter for shock capturing, a non-orthogonal staggered spherical mesh with constrained transport, conservative averaging-reconstruction for axis singularities and the capability to handle multiple ion species. However the computational kernel of the LFM code, designed and optimized for architectures long retired, has aged and is difficult to adapt to the modern multicore era of supercomputing. To carry its legacy forward, we re-envisage the LFM as GAMERA, Grid Agnostic MHD for Extended Research Applications, which preserves the core numerical philosophy of LFM while also incorporating numerous algorithmic and computational improvements. The upgrades in the numerical schemes include accurate grid metric calculations using high-order Gaussian quadrature techniques, high-order upwind reconstruction, non-clipping options for interface values. The improvements in the code implementation includes the use of data structures and memory access patterns conducive to aligned, vector operations and the implementation of hybrid parallelism, using MPI and OMP. Thus, while keeping the best elements of LFM, GAMERA is designed to be a portable and easy-to-use code that provides multi-dimensional MHD simulations in non-orthogonal curvilinear geometries on modern supercomputer architectures.

Field line curvature scattering by the magnetic field structure associated with bursty bulk flows (BBFs) has been studied, using simulated output fields from the Lyon-Fedder-Mobarry global magnetohydrodynamic code for specified solar wind input. There are weak magnetic field strength (B) regions adjacent to BBFs observed in the simulations. We show that these regions can cause strong scattering where the first adiabatic invariant changes by several factors within one equatorial crossing of energetic electrons of a few kiloelectron volts when the BBFs are beyond 10 R-E geocentric in the tail. Scattering by BBFs decreases as they move toward the Earth or when the electron energy decreases. For radiation belt electrons near or inside geosynchronous orbit we demonstrate that the fields associated with BBFs can cause weak scattering where the fractional change of the first invariant (mu(0)) within one equatorial crossing is small, but the change due to several crossings can accumulate. For the weak scattering case we developed a method of calculating the pitch angle diffusion coefficient D-alpha alpha.D-alpha alpha for radiation belt electrons for one particular BBF were calculated as a function of initial energy, equatorial pitch angle, and radial location. These D-alpha alpha values were compared to calculated D-alpha alpha for a dipole field with no electric field. We further compared D-alpha alpha values with that of stretched magnetic fields calculated by Artemyev et al. (2013, https://doi.org/10.5194/angeo-31-1485-2013) at r approximate to 7 R-E. Results show that scattering by BBFs can be comparable to the most highly stretched magnetic field they studied.

During geomagnetic storms the intensities of the outer radiation belt electron population can exhibit dramatic variability. In the main phase electron intensities exhibit deep depletion over a broad region of the outer belt. The intensities then increase during the recovery phase, often to levels that significantly exceed their pre-storm values. In this study we analyze the depletion, recovery and eventual enhancement of radiation belt intensities during the 2013 St. Patrick's day geomagnetic storm. We simulate the evolution of the high-energy electron population comprising the outer radiation belt using our newly-developed test particle radiation belt model (CHIMP) based on a hybrid guiding-center/Lorentz integrator and electromagnetic fields derived from a coupled 3D ring current and global MHD simulation (LFM-RCM). Our approach differs from previous work in that we use MHD information to identify regions of strong, bursty, and azimuthally localized Earthward convection in the magnetotail where test particles are then seeded. In other words, magnetospheric electromagnetic fields inform how test-particles evolve and the plasma flow informs where and when new test-particles are created. This is key to reproducing injection of energized plasma sheet electrons into the outer belt.

Much of plasma heating and transport from the magnetotail into the inner magnetosphere occurs in the form of mesoscale discrete injections associated with sharp dipolarizations of magnetic field (dipolarization fronts). In this study we investigate the mechanisms of ion acceleration at dipolarization fronts in a high-resolution global magnetospheric MHD model (LFM). We use large-scale three-dimensional test-particle simulations (CHIMP) to address the following science questions: 1) what are the characteristic scales of dipolarization regions that can stably trap ions? 2) what role does the trapping play in ion transport and acceleration? 3) how does it depend on particle energy and distance from Earth? 4) to what extent ion acceleration is adiabatic? High-resolution LFM was run using idealized solar wind conditions with fixed nominal values of density and velocity and a southward IMF component of −5 nT. To simulate ion interaction with dipolarization fronts, a large ensemble of test particles distributed in energy, pitch-angle, and gyrophase was initialized inside one of the LFM dipolarization channels in the magnetotail. Full Lorentz ion trajectories were then computed over the course of the front inward propagation from the distance of 17 to 6 Earth radii. A large fraction of ions with different initial energies stayed in phase with the front over the entire distance. The effect of magnetic trapping at different energies was elucidated with a correlation of the ion guiding center and the ExB drift velocities. The role of trapping in ion energization was quantified by comparing the partial pressure of ions that exhibit trapping to the pressure of all trapped ions.