The core of the code is the FORTRAN ionospheric dynamics engine. At this stage of development, this is numerically unchanged from the original release of the SAMI2 model, though the handling of some variables has been updated to accommodate compilation using GNU compilers (e.g., gfortran team, 2022). SAMI2 solves for the chemical and dynamical evolution of seven ion species in this environment (H⁺, He⁺, N⁺, O⁺, N 2 + , NO⁺, and O 2 + ). The temperature equation is solved for three ion species (H⁺, He⁺ and O⁺) and for the electrons. Ion inertia is included in the ion momentum equation for motion along the geomagnetic field. This is important in modeling the topside ionosphere where the plasma transitions from collisional to collisionless. SAMI2 uses a nonorthogonal, nonuniform, fixed grid. The grid is designed to optimize the numerical mesh so that the spatial resolution decreases with increasing altitude. The plasma is transported along the magnetic field using a semi-implicit transport algorithm, and transverse to the geomagnetic field using a finite volume method in conjunction with the donor cell method (Huba, 2003). The numerical solutions are well documented in Huba et al. (2000). A brief summary follows.

The SAMI2 model simulates the production, motion, and loss of ions along a two-dimensional slice of Earth’s ionosphere, as shown in Figure 1. This slice is aligned with magnetic field lines as calculated for an offset tilted dipole field. The continuity, momentum, and temperature equations for ions and electrons are solved. The model is initialized and driven by empirical models, as discussed in Section 2.2. A series of scaling factors can be used to alter the magnitude of these empirical values through the namelist file. In general, the model is run for 24 h before modelled values are output to files. This is done to clear transients from the system.

The sami2py software builds on the modular nature of the SAMI2 model. In the original release, SAMI2 used four key empirical models to prime the ionospheric solutions: NRLMSISe-00 (Picone, 2002) to provide the neutral atmosphere, EUVAC (Richards et al., 1994) to provide the EUV spectrum, HWM-93 to provide neutral winds (Hedin et al., 1993b,a), and the Fejer-Scherliess model of low-latitude E×B drifts (Scherliess and Fejer, 1999). sami2py updates these component models, whose acronyms are defined below, to the latest versions and includes the older versions as optional inputs. Additionally, the number of scalable parameters has been expanded. A full list of the available models and scalable parameters is included in Table 1.

The Naval Research Laboratory Mass Spectrometer and Incoherent Scatter radar (NRLMSIS) model is a semi-empirical model representing multiple decades of neutral atmospheric measurements, including mass spectrometer, radar, and satellite drag data (Picone, 2002). The version implemented in sami2py is a modification of the extended version of the model released in 2000 (NRLMSISe-00). During the solar minimum between cycles 23 and 24, record low densities in the thermosphere were observed through satellite drag measurements Emmert et al. (2010) and direct measurement of neutral pressure density (Haaser et al., 2010). These measurements were outside of the underlying database used to construct the model. Solomon et al. (2010) suggested that anomalously low Extreme Ultraviolet (EUV) radiation during this period resulted in a much cooler thermosphere than expected from the radio flux proxy for solar activity (F_10.7). Since F_10.7 rather than EUV is used to drive the thermospheric model, Klenzing et al. (2013) implemented a scalar factor for the exospheric temperature in their empirical study of altered electrodynamics during extreme solar minima. The SAMI2 model already allows users to scale the resultant density profiles independently for each species after NRLMSISe-00 has run. The modification implemented here adds the capability to scale the exospheric temperature directly in NRLMSISe-00 in addition to constantly scaling each species. An example of the effect of this reduced temperature run is shown in Figure 2.

The Extreme Ultraviolet for Aeronomic Calculations (EUVAC) model provides a calculation of the EUV flux as a function of the solar radio flux proxy F_10.7 (Richards et al., 1994). For SAMI2, the model is used to calculate the photo-ionization rate of the ionosphere. While the implementation is unchanged from the SAMI2 1.00 release, a scalar parameter has been added to the code to allow sensitivity studies for directly changing the total photo-ionization rate.

The Horizontal Wind Model (HWM) provides a statistical view of neutral winds gathered from world-wide Fabry-Perot Interferometers, Incoherent Scatter Radars, satellites, and rockets (Drob et al., 2015). The latest version (HWM14) is incorporated as the default, thought users can run numerical experiments with HWM07 (Drob et al., 2008) and HWM93 as options.

The Fejer-Scherliess model of E×B drift climatology (e.g., Scherliess and Fejer, 1999) provides the two-dimensional drifts perpendicular to the magnetic field lines as a function of local time, solar activity, day of year, and longitude. This is done through cubic spline fits to data from the Jicamarca Incoherent Scatter Radar and the Atmospheric Explorer E satellite. The model is unchanged in the sami2py implementation. As in SAMI2, scalar parameters allow users to directly change the magnitude and offset of the drifts.

An alternative E×B is provided for users wanting to investigate alternate drift climatologies. Since the model is constrained to a local series of flow tubes in a single magnetic meridian, the alternate model is incorporated as a series of Fourier coefficients that are user-specified that describe a function of Solar Local Time (SLT), as shown in Eq. 1. E × B total S L T = Σ i = 1 10 C i 0 ⁡ cos i π S L T 12 + C i 1 ⁡ sin i π S L T 12 (1) This allows users with direct measurements to create a localized drift model. Examples of this type of usage are presented in Klenzing et al. (2013) and Smith and Klenzing (2022). An additional input file to the FORTRAN code names exb.inp was added so that the localized model can be changed without recompiling the FORTRAN engine. Note that this creates a function that averages to zero over all local times, ensuring that there is no net upward or downward drift over the course of a day.

The sami2py Python code wraps the compiled SAMI2 FORTRAN engine (see Figure 3) in a standardized Python package. It provides an interface for users to directly update the namelist and E×B input files via keywords, and returns the results in an xarray.Dataset object (Hoyer and Hamman, 2017).

The core SAMI2 code in sami2py is compatible with FORTRAN 90 and is suitable for compilation under multiple compilers. The variable parameters, such as geographic location, solar activity, and season, are input via a namelist file, and the resulting modelled parameters are sent to binary output files. An additional exb.inp file is included to generate alternate E×B drift models via a Fourier series over solar local time. The sami2py code provides a user interface to both the input namelist files (through the sami2py.run_model method) and the output binaries (through the sami2py.Model class).

The method sami2py.run_model allows the user to directly run the compiled FORTRAN executable. The namelist that specifies the parameters of the model run can be adjusted via keyword arguments, which are fully documented in the code docstrings and in the detailed documentation that is available in the GitHub repository and online at readthedocs. This includes a user-specified “tag” to quickly describe the run for archival purposes (e.g., “solarmin”). The FORTRAN executable saves each variable as a separate file. By default, this method will move all of the output files, as well as the input namelist and exb.inp files, to an archival directory. All files are grouped under subdirectories by the tag name, longitude, and date in case a user runs multiple dates or locations for the same input conditions.

The sami2py.Model class loads the raw output of the model run. It loads each individual file and reshapes them into a single xarray. Dataset object for convenience of use. This class will also load the namelist info as metadata to allow inspection of input parameters, as well as any custom E×B input that was used. When working within sami2py, this information is stored in the model.MetaData object as a dictionary. The parameters are reshaped as 4D arrays with appropriate coordinates. Examples are shown in the sample code in Section 3. Users may run analysis directly from the Model object or save to a single file.

For portability and reproducability, both data and metadata can be exported to a netCDF4 file (Whitaker et al., 2020) using the to_netcdf method on the model. The metadata will be included as top-level attributes in the output file, documenting how the run was initialized and including both the sami2py version number and commit hash (in case a custom branch based on an official version was created). The netCDF4 versions of the file are constructed to be compatible with pysat.

The pysat ecosystem (Stoneback et al., 2018) has evolved to support management and analysis of a number of data sets throughout the space science community. The core pysat engine provides a framework to manage data sets, including acquisition, archival, and management. As a management tool, it has been used operationally in missions and analysis projects, including the ICON and COSMIC2 missions. A series of libraries has been written to translate between the core pysat commands and individual data sets. This standardization allows pysat to manage the metadata as well.

These files can be integrated into the pysat ecosystem by using the custom sami2py instrument module at pysatModels (Burrell et al., 2022). This package includes a number of other tools to compare observational data with models.

The code here has been tested in linux, Mac, and Windows environments through Github Actions. Each environment is tested through a unit test suite with 97.6% code coverage as of version 0.3.0. The unit tests are configured to use the latest python packages under python 3.9 and 3.10 environments, as well as a version limited to numpy 1.20 under python 3.8. The specific versions used for the core requirements as of the publication of this paper are listed in Table 2.

The sami2py.run_model method and sami2py.Model class provide the core functionality of sami2py. The following code snippet prepares the archive directory, and specifies the time and location for the run as well as declaring custom E×B input.

Note that setting the user archive directory only needs to be run when the package is first installed.

Now that the custom input has been declared and the environment is prepared for archival, the model can now be executed. The time, location, F10.7 and E×B are provided to the sami2py.run_model method. Upon completion the model output is loaded as a sami2py.Model object and archived as a netCDF file.

The following code snippet loads the archived model run, adds a new variable to the data set which consists of the total plasma density, and then plots the total plasma density as a function of local time and altitude with the E×B drift superimposed over the density. Note that by default, the ion density variable (deni) is a four-dimensional object, with one of the dimensions (retrievable as ion) specifies the individual ion species. A summation over this third axis is needed to extract total ion density.

The resulting figure is shown in Figure 4.

One of the applications of sami2py is for early-phase research projects. The user is the broader ionospheric research community who are communicated with on a direct basis with the development team and at conferences such as CEDAR. The operational environment is then considered to be an individual’s work computer.

An example of the early-phase research projects is running sensitivity studies on proposed physical forcing mechanisms. For this paper, an example of an identified user for this application is Klenzing et al. (2013) where the early phase research includes a series of sensitivity studies for proposed modifications to ionospheric drivers under extremely low levels of solar activity. This study was originally conducted using a prototype of the sami2py model written in MatLab, but the functionality applies to the Python version as well. Each empirical model that drives the SAMI2 ion dynamics engine can be modified to reflect proposed changes to the forcing of the ionosphere, including reductions in exospheric temperature for the MSIS model and the direct input of user-specified E×B drift profiles as a function of local time.

Examples of how the ionospheric density changes by altering the E×B drift assumptions are shown in Figures 4, 5. Each plot shows the evolution of the vertical ionospheric density profile over time. The white line plotted above the ionospheric density represents the driving E×B timeseries used in sami2py, with Figure 4 driven by the Fejer-Scherliess model (Scherliess and Fejer, 1999) and Figure 5 driven by climatology measured by the Coupled Ion-Neutral Dynamics Investigation (CINDI) mission of opportunity (Smith and Klenzing, 2022).

The work discussed above has shown how a Python version of SAMI2 will provide a path beyond the current state of the art capabilities for individual research projects. The Python interface for the SAMI2 model also provides a new capability making it easier for more researchers to access and use this model, as well as document results. Moving from a MatLab interface to an open source language improves the accessibility of the work. Incorporation of the resulting modeled data into an xarray. Dataset object improves the usability of the output. The primary requirement for this application at this phase is to ensure that this Python package is open access and works across computer operating systems. We have satisfied the milestones for AUL 3 with the release of sami2py version 0.2.0 in December 2019 (Klenzing et al., 2019).

The AUL four to six milestones require improved documentation and testing of the beta prototype of the model. Changes incorporated since version 0.2.0 include docstrings for all functions, improved Continuous Integration (CI) testing, and improved compatibility with external Python packages, including numpy, xarray, and pysat. The model undergoes continuous integration tests in the GitHub Actions environment with > 97 % coverage, fulfilling simulation in an operational environment. The CI tests are run for Linux, mac, and windows systems to satisfy AUL 5 (demonstration in a relevant context). Additionally, tests for older versions of numpy are included to maintain compliance with NEP029 (Caswell et al., 2019). Since sami2py is being developed on GitHub, it is easily transferred from the development environment to the operational environment (end user’s workstation) across the community. Regular updates are given at community workshops. With the documentation of the code, including the online documentation at readthedocs and the examples within this paper, and the release of version 0.2.5 (Klenzing et al., 2021) all milestones through AUL six have been completed.

AUL level 7 is the Application Prototype of the project. This requires demonstration of the prototype and dissemination of results. Both of these goals are achieved with the release of version 0.3.0 (Klenzing et al., 2022) and the publication of this paper. Improvements to the user interface and code style have been implemented in version 0.3.0 to maintain PyHC standards and improve code maintainability.

For AUL 8 and 9, a finalized project for on-demand usage needs to be released. In the context of this application for sami2py, a series of updates focusing on an improved workflow and code maintainability have been identified. These are demarcated as a future 0.4.0 release. Input from the community will be evaluated alongside these updates as the user base grows.

As an additional demonstration of the prototype, the sami2py module is a central dependency for the growin python module which was written to compute the Rayleigh-Taylor instability (RTI) growth rate. The calculation of the RTI growth rate is central to the development and growth of plumes of depleted plasma, or plasma bubbles, in the bottomside of the equatorial ionosphere. The growin module uses the sami2py module to run the SAMI2 model, archive the output, and load the output into Python data structures (Klenzing et al., 2022). Similar to the example code above, drift measurements are used to create a climatological drift profile from in-situ measurements. These drifts are then passed to sami2py and an ionosphere is simulated with the typical ionospheric indices for the corresponding time period. Subsequently the produced ionospheric plasma densities, drifts, and winds are used to compute flux-tube integrated quantities necessary to compute the RTI growth rate. These growth rates have been previously used to discuss bubble occurrence frequencies obtained from the CINDI (Smith and Klenzing, 2022) and Global Observations of the Limb and Disk (GOLD) (Martinis et al., 2021) missions.

Similar to the previous application, the broader ionospheric research community is the user and will benefit from a Python version of growin and the inclusion of sami2py within it. The feasibility, viability, and expected improvements can all be found within Smith and Klenzing (2022). Thus many of the milestones have been completed for this application through the previously discussed application in Section 4.1. As shown in Table 5, the key additional requirement here is the output of neutral atmospheric data, which is required to perform the RTI calculations. This has been added to sami2py as an optional output. As the other components growin were already within the operational/end user environment, the final AUL is now dependent on the progress of sami2py. Similar to the previous application, the usage of sami2py in the growin package is at an AUL of 7.

Beyond the research community, another user community has been identified but not yet contacted. The code here can also be used as an educational tool as part of a Space Weather of Ionospheric Electrodynamics curriculum. The straightforward and modular nature of the code makes it practical to incorporate into homework or class projects as needed. As this application has been identified, but specific requirements have not been defined and incorporated into the code, this is defined as an AUL 1 project. Work is ongoing, and interested parties should contact the authors to help better refine this project and requirements for these purposes.