© F. dos Santos Prol et al., Published by EDP Sciences 2025

1 Introduction

Improving the representation of the topside ionosphere and plasmasphere has been a recent focus of several well-established models. Notable examples include recent updates to the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) (Richmond et al., 1992), the International Reference Ionosphere (IRI) (Bilitza et al., 2022), and NeQuick 2 (Nava et al., 2008). For example, Cai et al. (2022) extended the TIE-GCM topside ionosphere to cover altitudes from 700 to 1,200 km. The IRI model has also incorporated new options for representing the topside ionosphere and plasmasphere (Bilitza et al., 2024; Servan-Schreiber et al., 2025) by integrating the Global Core Plasmaspheric Model (GCPM) (Gallagher et al., 2000) and the empirical model from Ozhogin et al. (2012). Similarly, NeQuick 2 underwent revision, with recent studies exploring new ways to represent empirical parameters that govern ionospheric scale height (Themens et al., 2018; Pezzopane et al., 2023) along with efforts to incorporate the Gallagher plasmasphere model into NeQuick 2 (Cueto et al., 2007). These advances are driven by the need to address the mismodeling of the topside ionosphere and plasmasphere, which can contribute up to 50–60% of the total ionospheric delay in Global Navigation Satellite Systems (GNSS) (Yizengaw et al., 2008; Gulyaeva & Gallagher, 2007).

To better account for the ionospheric delay corrections in GNSS measurements, significant progress has been made using radio-occultation (RO) techniques, e.g., Pignalberi et al. (2022); Smirnov et al. (2023). RO observations have shown that ionospheric electron density can be effectively characterized from the peak height (hmF2) up to approximately 800 km using a linear scale height with Chapman or Semi-Epstein functions (Prol et al., 2018; Pignalberi et al., 2020). However, accurately representing electron density beyond this altitude remains challenging. Recent studies have highlighted the potential of using the Van Allen Probes data for modeling the topside ionosphere up to the plasmaspause altitudes (Prol et al., 2022), with observations extending above 8,000 km. A key limitation of this approach, however, is the data gap between RO and Van Allen Probes measurements, particularly in the transition region. Although the model provided accurate representations, it did not establish a unified analytical formulation that integrated the ionosphere and plasmasphere. Instead, the study emphasized the need for a dedicated model and additional data sources to achieve a more comprehensive representation of the topside ionosphere and plasmasphere.

This study introduces a new formulation for representing the topside ionosphere and plasmasphere. Designed based on the NeQuick 2 model, the proposed approach ensures compatibility of the ionospheric observations with Van Allen Probes measurements in the plasmasphere while keeping the representation of the lower topside ionosphere. The primary objective is to improve the known underestimate of electron density of NeQuick 2 (Cherniak & Zakharenkova, 2016; Kashcheyev & Nava, 2019) while preserving a realistic plasmaspheric structure and its morphology along the geomagnetic field lines. The following sections outline the method developed, present an evaluation using reference data from 2008 to 2024, and present conclusions with a summary of findings.

2 Material and method

The model developed in this study is named IonoMAP (Ionospheric Mapping). Figure 1 presents a block diagram that outlines the key parameters of the model. The primary inputs are position and time, while the outputs include electron density profiles, global 3D electron density grids, and total electron content (TEC) along the line of sight between GNSS satellites and ground- or satellite-based receivers. The implementation is based on the NeQuick 2 source code with several modifications. In the following, we detail these differences and the enhancements introduced in the new formulation.

Figure 1 Block-level diagram of the main parameters (gray boxes) estimated with IonoMAP. Black boxes provide the input and output parameters.

2.1 IonoMAP model formulation

The solar flux index F10.7 is obtained from the continuously updated apf107.dat table, available on the IRI website.¹ An empirical formula is then applied to convert F10.7 into the sunspot number, as specified by the Radiocommunication Sector of the International Telecommunication Union (International Telecommunication Union, Radiocommunication Sector, 2023):

$$R12=\sqrt{167273.0+(F107-63.7)1123.6}-408.99,$$(1)

where the yearly-smoothed sunspot number is adopted.

The MOdified DIP (MODIP) latitude is obtained as (Rawer, 1963):

$$\tan \mu =\frac{I}{\sqrt{\cos \varphi }},$$(2)

where the magnetic inclination I is computed using the International Geomagnetic Reference Field (IGRF) model. The symbol ϕ refers to the input geographic latitude.

The parameters derived by Jones & Gallet (1962), formally known as the coefficients of the Consultative Committee on International Radio (CCIR), are used to calculate the critical frequency of the F2 layer (foF2) and the propagation factor for 3000 km (M(3000)F2). Details on utilizing the CCIR coefficients are provided in Forsythe et al. (2024). The F2-layer peak altitude (hmF2) is determined based on M(3000)F2, as described in Nava et al. (2008).

The critical parameters of the E- and F1-layers are computed as described in NeQuick 2, including hmE, hmF1, foE, and foF1. Similarly, the bottomside thickness parameters, including the BE_bot, BE_top, B1_bot, B1_top, and B2_bot, are derived using the same approach as NeQuick 2. The topside thickness parameter (B2_top) is also computed using the NeQuick 2 formulation, as follows:

$$B{2}_{\mathrm{top}}=\mathrm{kB}{2}_{\mathrm{bot}}\left[1+\frac{\mathrm{rg}(h-\mathrm{hmF}2)}{{\mathrm{rkB}2}_{\mathrm{bot}}+g(h-\mathrm{hmF}2)}\right],$$(3)

where h stands for altitude and

$$B{2}_{\mathrm{bot}}=\frac{0.385\mathrm{NmF}2}{(\mathrm{d}N/\mathrm{d}h{)}_{\max }},$$(4)

$$\ln \left[\right(\mathrm{d}N/\mathrm{d}h{)}_{\max }]=-3.467+1.714 \ln (\mathrm{foF}2)+2.02 \ln (M\left(3000\right)F2),$$(5)

$$k=3.22-0.0538 \mathrm{foF}2-0.00664 \mathrm{h}\mathrm{mF}2+0.113\frac{\mathrm{hmF}2}{B{2}_{\mathrm{bot}}}+0.00257 R12,$$(6)

with r = 100 and g = 0.125.

To ensure a smooth “fadeout” of the E and F1 layers near their peak heights and to prevent artificial secondary maxima, NeQuick 2 employs semi-Epstein layers with a ξ function for electron density computation. This function ensures an exponential decay between the different layers. A key distinction of the current model is its direct use of the Epstein formulation, eliminating the need for the ξ function and its associated iterative process.

The Epstein layer is described as follows:

$$N_{\mathrm{Epstein}}=\frac{4N_{\max }}{{\left(1+\exp \left(\frac{h-h_{\max }}{B}\right)\right)}^2}\exp \left(\frac{h-h_{\max }}{B}\right),$$(7)

where N_max is the electron density at the peak of the corresponding layer, h_max is the altitude at the peak of the layer, and B is the thickness parameter of the corresponding layer.

To control the transition between Epstein layers, we utilized the sigmoid function:

$$\sigma =\frac{1}{1+e^y},$$(8)

where y is defined to vary linearly from −3π to 3π within the altitudes. To achieve this, we formulate the following equation

$$y=3\pi +\frac{(h_{\min }-h)}{h_{\max }-h_{\min }}6\pi ,$$(9)

where the main idea is that when y = −3π, which corresponds to h = h_max, the value of σ approaches unity, thereby assigning greater weight to the upper atmospheric layers. Otherwise, when y = 3π, which corresponds to h = h_min, σ tends toward zero, reducing the weight assigned to the lower atmospheric layers. Additionally, we impose that σ = 1 for altitudes above h_max, and σ = 0 for altitudes below h_min. Table 1 summarizes the values adopted for each sigmoid function.

The electron density is then given as:

$$N_e=N_{E,\mathrm{bot}}+N_{E,\mathrm{top}}\left(1-{\sigma }_{E,F1}\right)+N_{F1,\mathrm{bot}}\left({\sigma }_{E,F1}\right)+N_{F1,\mathrm{top}}\left(1-{\sigma }_{F1,F2}\right)+N_{F2,\mathrm{bot}}\left({\sigma }_{F1,F2}\right)+N_{\mathrm{top}},$$(10)

where the electron density value refers to the corresponding thickness parameter. For example, N_E,bot is calculated based on BE_bot. Notably, when the E or F1 layer is absent, its contribution is set to zero.

2.2 New formulation for the topside ionosphere and plasmasphere

The topside function provided by equations (3)–(7) is suited to describe the topside ionosphere. The model developed in this work is intended to accommodate this topside ionosphere formulation with a general representation of the plasmasphere given by Van Allen Probes data. All empirically defined numerical values are based on the values found by Prol et al. (2022).

As shown in Prol et al. (2022), the electron density measured by the Van Allen Probes can be approximated using a linear scale height, which is equivalent to describing the logarithm of the electron density as linear decay. Accordingly, the electron density of the plasmasphere N_Plas is given by a linear logarithmic relation:

$$\mathrm{lo}{\mathrm{g}}_{10}\left(N_{\mathrm{Plas}}\right)=P_0+\mathrm{d}{P}_0(h-h{P}_0),$$(11)

where P₀ is the starting point of the logarithmic decay, represented as:

$$P_0=\mathrm{lo}{\mathrm{g}}_{10}\left(N_{\mathrm{top}}^{1500}\right).$$(12)

The quantity $N_{\mathrm{top}}^{1500}$ denotes the electron density of the topside model defined in the previous section. The reference altitude is set at 1500 km to mark an average location of the transition region between the topside ionosphere and plasmasphere at the equator.

The gradient of the exponential decay is provided by:

$$\mathrm{d}{P}_0=\frac{\mathrm{lo}{\mathrm{g}}_{10}\left(N_{\mathrm{pp}}\right)-P_0}{h_{\mathrm{pp}}-h{P}_0},$$(13)

being h_pp the altitude of the plasmapause and N_pp the electron density at the plasmapause. Based on the Van Allen Probe data (Prol et al., 2022), we have observed that the plasmapause at the equator occurs roughly at 30000 km, with an average electron density of N_pp = 1.01 × 10⁸ el/m³. The plasmapause location at the equator is propagated by the geomagnetic field lines, considering the cosine of the MODIP latitude:

$$h_{\mathrm{pp}}=25 000\cos \left(\mathrm{modip}\right)+5000,$$(14)

in which h_pp refers to the altitude of the plasmapause in km. At the geomagnetic equator, equation (14) defines the plasmapause with around the L-shell of 5.7 (L_pp). It should be noted that the L_pp value of 5.7 is quite common for geomagnetically quiet conditions, but is not valid during geomagnetic storm conditions, and incorporating the activity-dependent plasmapause models will be performed in future studies. Furthermore, around the geomagnetic poles, where the magnetic field lines are open and the plasmasphere is not present, the physical meaning of the h_pp parameter is different and does not represent the plasmapause location; rather, it signifies the altitude at which the ionospheric content becomes almost non-existent due to the vertical outflow.

The $N_{\mathrm{top}}^{1500}$ parameter is obtained by propagating a virtual point of the electron density at the equator through the geomagnetic field lines. This virtual point of the electron density at the equator is identified as $N_{\mathrm{top},\mathrm{eq}}^{1500}$. The $N_{\mathrm{top},\mathrm{eq}}^{1500}$ parameter is obtained for a specific position (lat,lon,alt) with equations (3) to (7) as if the point is located at the equator and altitude of 1500 km. The underlying relation between $N_{\mathrm{top}}^{1500}$ and $N_{\mathrm{top},\mathrm{eq}}^{1500}$ is described as:

$$N_{\mathrm{top}}^{1500}=N_{\mathrm{top},\mathrm{eq}}^{1500}\cos (\mathrm{modip}{)}^2+10^8,$$(15)

where the virtual point of the electron density in the equator is identified as $N_{\mathrm{top},\mathrm{eq}}^{1500}$. In equation (15), 10⁸ el/m³ is summed to keep the background ionosphere with a reasonable existence.

The electron density at the topside ionosphere and plasmasphere is computed as:

$$N_{\mathrm{top}}^{\mathrm{new}}=N_{\mathrm{top}}(1-{\sigma }_{\mathrm{top},\mathrm{plas}})+N_{\mathrm{plas}}\left({\sigma }_{\mathrm{top},\mathrm{plas}}\right),$$(16)

where the sigmoid between the topside ionosphere and plasmasphere (σ_top,plas) is computed using equations (8)–(9) and considering h_min = 800 km and h_max = 2000 km. The lower boundary of 800 km approximately corresponds to both the upper limit of GNSS Radio Occultation (RO) data and the upper boundary of the ionosphere. Although RO data is not utilized in the present study, this choice facilitates its assimilation into the model in future research. The upper boundary of 2000 km is defined as the initial point of the plasmasphere. Beyond this altitude, we assume that only the plasmaspheric electron content remains. Figure 2 shows an example of an electron density profile obtained by IonoMAP with the new model formulation.

Figure 2 Typical electron density profile representation by IonoMAP.

2.3 Model validation

The validation primarily aims to assess the IonoMAP model’s ability to represent the topside ionosphere and plasmasphere. We compare the estimates from IonoMAP based on the NeQuick 2 topside formulation (Eqs. (3)–(7)) with those from the proposed new formulation (Eqs. (11)–(16)). To achieve this, the validation is conducted using multiple years of data from four datasets: Swarm-A, COSMIC-1, Van Allen Probes, and Global Ionospheric Maps (GIMs). Although our model was primarily developed to represent the ionosphere and plasmasphere under quiet conditions, no validation data were excluded based on geomagnetic activity indicators. Given that several years of data were used, we expect that ionospheric and plasmaspheric disturbances were largely smoothed out in the validation. However, a comprehensive assessment of the model’s performance during disturbed periods should be addressed in future work.

The Swarm-A dataset is obtained from the ESA web platform² and includes calibrated Slant Total Electron Content (STEC) data as level 2 products (TECATMS). We gathered data from the following directories: Level2daily, Entire mission data, TEC, TMS, Sat A. We evaluated data from 2014 to 2024, organized in 30-day intervals. The dataset contains absolute STEC values derived from precise orbit determination (POD) observations made by receivers on the Swarm-A satellite (van den IJssel et al., 2016), which orbits at an altitude of approximately 470 km. We filtered out data with negative or excessively high STEC values, retaining only data within the range of 1 to 200 TECU, where 1 TECU equals 10¹⁶ el/m². An elevation mask of 30 degrees is applied to remove noise. Additionally, short arcs with fewer than 10 observations, or where the daily minimum STEC exceeded 8 TECU, are excluded. The developed ionospheric model is then used to compute the STEC values along the line of sight between the Swarm-A and GNSS satellites. While this approach eliminates the unrealistic approximations often provided by mapping functions, it imposes a significant computational burden, as the model must be run for multiple points along each line of sight. To mitigate this, we subsampled the daily files by selecting only one point out of every 100 points. This has lowered the number of STEC values computed in each TECATMS file by about 100 times.

COSMIC-1 dataset (UCAR COSMIC Program, 2022⁴) is obtained from the COSMIC Data Analysis and Archive Center (CDAAC) and includes calibrated STEC data as level 1b products (podTec). We evaluated data from 2008 to 2018, organized in 30-day intervals. The STEC values are collected by POD receivers onboard all COSMIC-1 satellites, which orbited at an altitude of approximately 800 km. We filtered out data with negative STEC values, retaining only those where STEC > 1 TECU. An elevation mask of 30 degrees is applied to remove noise. Additionally, daily arcs where the minimum STEC exceeded 10 TECU are excluded. The validation is carried out at the STEC level. To reduce computational burden, we reduced each podTec file by selecting only one point out of every 10 points.

Electron density measurements from the Van Allen Probes mission are obtained using the Neural Network-based Upper Hybrid Resonance Determination (NURD) algorithm developed by Zhelavskaya et al. (2017). We evaluated data from 2013 to 2017 collected by probe A (VAP-A). The dataset includes electron density measurements at altitudes up to 5.8 Earth radii. We filtered the data to retain electron density values (N_e) within the range 10⁶ < N_e < 10¹⁰ el/m³ and altitudes (h) between 7,000 and 27,000 km.

The GIM dataset (Hernández-Pajares et al., 2009) is obtained from GNSS products provided by the International GNSS Service (IGS) and processed by the Center for Orbit Determination (CODE). It includes VTEC measurements modeled using spherical harmonics, assuming the ionosphere as a single-shell layer located at an altitude of 450 km. We analyzed 10 years of data, spanning from 2008 to 2018. For each day, we evaluated data at two time instances: 0- and 12 hours UT. The integration of electron density into VTEC using IonoMAP is performed with a 20 km step size.

3 Results and discussion

The results presented in this section are based on the IonoMAP model developed using equations (1)–(16). The model was implemented in Python. While the implementation is based on the NeQuick 2 source code, there are notable differences (discussed in the previous section) between the developed model and NeQuick 2. Therefore, the results presented here reflect the performance of our internally developed adaptation of NeQuick 2, rather than NeQuick 2 itself. It is important to note that the default topside function in our model is the same as in NeQuick 2, which we refer to as NEQ in the results. The new topside model formulation is referred to as NEW in the results.

Figure 3 shows the electron density representation of the ionosphere produced by IonoMAP, plotted in terms of latitude and altitude. The graphs correspond to December 30, 2013, at 00:00 Universal Time (UT). The left panel displays the ionospheric profile using the NeQuick 2 formulation for the topside, where the crests of the equatorial ionization anomaly (EIA) extend across all altitudes in the topside, even as high as 2000 km. This is due to the fact that the topside thickness parameter in the NeQuick 2 is modeled as a function of the bottomside thickness and therefore extends the ionization crests upward even to unphysical heights. In contrast, the right panel shows the new proposed formulation, where the two crests merge into one, as expected in the equatorial fountain effect.

Figure 3 Topside ionosphere representation based on IonoMAP assuming the NeQuick 2 formulation for the topside (left panel) and the new topside formulation (right panel). The color bar is shown in logarithmic scale for better visibility. The results were obtained on December 30, 2013, at 00:00 UT for a longitude of 180°.

Figure 4 presents a meridional cut obtained using the NeQuick 2 formulation (left panel) and the new proposed model for the topside (right panel). The meridional cut illustrates the ionosphere and plasmasphere at 180°W (left side) and the opposite longitude 0°E (right side). The left side of each graph represents daytime, while the right side represents nighttime. The proposed model follows the geomagnetic field lines and is not influenced by the EIA crests, unlike the current NeQuick 2 formulation.

Figure 4 Meridional cuts of the ionosphere and plasmasphere based on the developed models. The result obtained with NeQuick 2 formulation is shown in the left panel, while the new formulation is shown in the right panel. The color bar is shown in logarithmic scale for better visibility. The representation is provided in geographic coordinates on December 30, 2013, at 00:00 UT for a longitude of 0° and 180°.

Figure 5 presents a time series of the absolute mean and standard deviation of the error between the modeled and reference STEC values, with COSMIC-1 data used as the reference. Data above approximately 800 km is utilized for the STEC computation. The absolute mean error using the IonoMAP with the NeQuick 2 formulation for the topside is about 1 TECU worse than when using the new parametrization. This improvement is evident over several years. Indeed, both models show worse performance between 2012 and 2014, which refers to periods of higher solar activity.

Figure 5 Absolute mean error and standard deviation (STD) of the error for over 10 years of STEC data collected by COSMIC-1.

Figure 6 presents a time series of the absolute mean error and standard deviation of the error, using Swarm-A STEC data as the reference, for the two topside formulations. In both cases, the altitude variations of the satellite are considered to compute the STEC values. The new formulation consistently outperforms the traditional one in both bias and dispersion. In both models, the error increases approaching 2014 and 2024, coinciding with two peaks in solar activity. During periods of low solar activity, the mean error remains below 5 TECU, while it can rise to around 15 TECU during higher solar activity.

Figure 6 Absolute mean error and standard deviation (STD) of the error for over 10 years of STEC data collected by the Swarm-A satellite.

Figure 7 presents an example of the difference between VTEC data derived from IonoMAP and CODE-GIM during high solar activity. Figure 8 shows the same, but for low solar activity. Notice that the bottom side is now considered. Overall, we can see a systematic underestimation of VTEC values by IonoMAP, especially in the EIA crests. The difference is more evident during high solar activity, ranging from −35 to 20 TECU. During low solar activity, IonoMAP better describes the VTEC values, with differences ranging from −10 to 10 TECU.

Figure 7 Difference between VTEC values derived by global ionospheric maps of IonoMAP and CODE. Top panel is obtained with the NeQuick 2 formulation for the topside ionosphere, and the bottom panel is obtained with the new formulation. Maps are related to January 2nd, 2013, at 00:00 hours UT.

Figure 8 Difference between VTEC values derived by global ionospheric maps of IonoMAP and CODE. Top panel is obtained with the NeQuick 2 formulation for the topside ionosphere and the bottom panel is obtained with the new formulation. Maps are related to January 1st, 2008, at 00:00 hours UT.

Figure 9 shows the mean error and standard deviation of the error when using CODE VTEC data as a reference for two time instances: 0- and 12-hours UT. In this case, the error includes contributions from both the topside and bottomside formulations. Again, we see that IonoMAP generally underestimates VTEC, with this underestimation closely linked to solar activity. During years of low solar activity (e.g., 2008–2011), the bias remains around zero, while during high solar activity years (2012–2015), the bias reaches up to nearly 25 TECU. The standard deviation of the error also increases with solar activity, though the change is more gradual.

Figure 9 Mean error and standard deviation (STD) of the error of IonoMAP, considering the NeQuick 2 formulation and the new topside formulation. VTEC data obtained through GIM are used as a benchmark for two instances: 0- and 12-hours UT.

Figure 10 presents an evaluation of IonoMAP in the plasmasphere, using electron density values from the Van Allen Probes as a reference. All individual electron density profiles collected by VAP-A are superimposed as a function of altitude. The distributions are shown in terms of the normalized number of occurrences on a logarithmic scale, calculated by dividing the number of times that a specific value of electron density was found on a specific pixel in the 2D histogram by the maximum number of times at that specific height. The top panels display the electron density distributions, while the bottom panels show the electron density errors. It is evident that the error in the NeQuick 2 topside formulation is underestimated by around 10 times. This bias is partially corrected by the new formulation. Furthermore, Van Allen Probes data show two peaks at high altitudes, one of them due to the plasmatrough, where the densities are two orders of magnitude lower than in the plasmasphere. In the developed model, this region is not considered.

Figure 10 2D histogram of profiles from the VAP-A instrument (top left), NeQuick 2 model (top middle) and the new formulation of the topside (top right) as a function of the altitude and logarithm of electron density with base 10. The bottom panels show the corresponding errors of the profiles. The color of the points represents the normalized occurrence in percentage.

Tables 2 and 3 present the quantitative values of the mean and standard deviation of the error for the estimated STEC and VTEC values using IonoMAP with both the NeQuick 2 and the new formulation for the topside ionosphere and plasmasphere. The NeQuick 2 formulation shows an overall underestimation of about 2–4 TECU, as demonstrated by previous investigations by Cherniak & Zakharenkova (2016); Kashcheyev & Nava (2019); Servan-Schreiber et al. (2025). In contrast, the new formulation, although still underestimating, exhibits a bias of around 1–3 TECU, representing an improvement in bias of approximately 28–40% in modeling the topside ionosphere and plasmasphere. In case of the VAP-A dataset, the mean error was 0.2 × 10⁹ el/m³ and 0.13 × 10⁹ el/m³ for the NeQuick 2 formulation and new model, respectively. The standard deviation of the error demonstrates similar precision between the two formulations.

4 Conclusions

A new approach for modeling the topside ionosphere and plasmasphere has been presented. The new model has been evaluated using TEC data from the COSMIC-1 mission, the Swarm-A satellite, and GIMs. Compared to the topside formulation implemented in NeQuick 2, the new formulation for the topside ionosphere and plasmasphere shows improvements of approximately 28–40%, with an underestimation ranging from 1 to 3 TECU, whereas the NeQuick 2 formulation underestimates by around 2–4 TECU. Although the developed model still exhibits an overall underestimation of VTEC, this is also linked to errors in mapping the bottomside and peak parameters. Additionally, the new model shows near-zero bias when compared to the Van Allen Probes in-situ measurements, while the formulation used in the NeQuick 2 model demonstrates a negative bias of approximately 10 times greater. Furthermore, the developed model better represents the plasmasphere’s dependence on geomagnetic field lines.

We believe that our proposed formulation for the topside ionosphere and plasmasphere is also capable of enhancing performance when applied to the original NeQuick 2 source code. However, this cannot be confirmed based on the current results. Future work may involve the integration of the proposed formulation into the original NeQuick 2 source code, followed by a systematic evaluation of its performance and consistency within the model framework.

Acknowledgments

The authors would like to acknowledge the Joint Study Group 1 “Understanding Ionospheric and Plasmaspheric Processes” of the Global Geodetic Observing System (GGOS) Focus Area on Geodetic Space Weather Research (GSWR) for their support and collaboration. The authors are grateful to the Telecommunications/ICT for Development (T/ICT4D) Laboratory team of the Abdus Salam International Centre for Theoretical Physics (Trieste, Italy), for developing, maintaining, and making freely available the NeQuick 2 model source code that served as a foundation for the present work. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Funding

This work was supported by the Research Council of Finland under the project “Exploring the Upper Limits of the Ionosphere” (EULI), under decision number 362096. Open Access funding provided by the National Land Survey of Finland.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

Swarm data used for validation are publicly available in the ESA Swarm Dissemination Server (https://swarm-diss.eo.esa.int/). COSMIC-1 data used for validation are publicly available in the UCAR portal (https://doi.org/10.5065/ZD80-KD74). The Van Allen Probes obtained by the NURD dataset are available at https://dataservices.gfz-potsdam.de/. The GIM dataset can be obtained through the NASA’s Crustal Dynamics Data Information System (CDDIS) via the portal https://www.earthdata.nasa.gov/centers/cddis-daac. The data for this article were generated by IonoMAP and are available upon request.

References

All Tables

All Figures

	Figure 1 Block-level diagram of the main parameters (gray boxes) estimated with IonoMAP. Black boxes provide the input and output parameters.
In the text

	Figure 2 Typical electron density profile representation by IonoMAP.
In the text

	Figure 3 Topside ionosphere representation based on IonoMAP assuming the NeQuick 2 formulation for the topside (left panel) and the new topside formulation (right panel). The color bar is shown in logarithmic scale for better visibility. The results were obtained on December 30, 2013, at 00:00 UT for a longitude of 180°.
In the text

Figure 4 Meridional cuts of the ionosphere and plasmasphere based on the developed models. The result obtained with NeQuick 2 formulation is shown in the left panel, while the new formulation is shown in the right panel. The color bar is shown in logarithmic scale for better visibility. The representation is provided in geographic coordinates on December 30, 2013, at 00:00 UT for a longitude of 0° and 180°.

In the text

	Figure 5 Absolute mean error and standard deviation (STD) of the error for over 10 years of STEC data collected by COSMIC-1.
In the text

	Figure 6 Absolute mean error and standard deviation (STD) of the error for over 10 years of STEC data collected by the Swarm-A satellite.
In the text

	Figure 7 Difference between VTEC values derived by global ionospheric maps of IonoMAP and CODE. Top panel is obtained with the NeQuick 2 formulation for the topside ionosphere, and the bottom panel is obtained with the new formulation. Maps are related to January 2nd, 2013, at 00:00 hours UT.
In the text

	Figure 8 Difference between VTEC values derived by global ionospheric maps of IonoMAP and CODE. Top panel is obtained with the NeQuick 2 formulation for the topside ionosphere and the bottom panel is obtained with the new formulation. Maps are related to January 1st, 2008, at 00:00 hours UT.
In the text

	Figure 9 Mean error and standard deviation (STD) of the error of IonoMAP, considering the NeQuick 2 formulation and the new topside formulation. VTEC data obtained through GIM are used as a benchmark for two instances: 0- and 12-hours UT.
In the text

	Figure 10 2D histogram of profiles from the VAP-A instrument (top left), NeQuick 2 model (top middle) and the new formulation of the topside (top right) as a function of the altitude and logarithm of electron density with base 10. The bottom panels show the corresponding errors of the profiles. The color of the points represents the normalized occurrence in percentage.
In the text