A new climatological electron density model for supporting space weather services

The ionosphere is the ionized part of the Earth’s atmosphere, ranging from about 60 km up to several Earth radii, whereas the upper part above about 1000 km height up to the plasmapause is usually called the plasmasphere. We present a new three-dimensional electron density model to support space weather services and mitigate propagation errors for trans-ionospheric signals. The model is developed by superposing the Neustrelitz Plasmasphere Model (NPSM) to an ionosphere model composed of separate F and E-layer distributions. It uses the Neustrelitz TEC model (NTCM), Neustrelitz Peak Density Model (NPDM), and the Neustrelitz Peak Height Model (NPHM) for the total electron content (TEC), peak ionization, and peak height information. These models describe the spatial and temporal variability of the key parameters as a function of local time, geographic/geomagnetic location, solar irradiation, and activity. The model is developed to calculate the electron concentration at any given location and time in the ionosphere for trans-ionospheric applications and named the Neustrelitz Electron Density Model (NEDM2020). A comprehensive validation study is conducted against electron density in-situ data from DMSP and Swarm, Van Allen Probes and ICON missions, and topside TEC data from COSMIC/FORMOSAT-3 mission, bottom side TEC data from TOPEX/Poseidon mission, and ground-based TEC data from International GNSS Service (IGS) covering both high and low solar activity conditions. Additionally, the model performance is compared with the 3D electron density model NeQuick2. Our investigation shows that the NEDM2020 performs better than the NeQuick2 compared with the in-situ data from Van Allen Probes and ICON satellites and TEC data from COSMIC and TOPEX/Poseidon missions. When compared with DMSP and IGS TEC data, both NEDM2020 and NeQuick2 perform very similarly.


Introduction
Models of electron density have been developed and utilized for many years. Physics-based models describe the ionospheric plasma by basic equations such as the continuity equation, equations of motion, and energy considering the close interaction of electrons and ions with thermosphere and magnetosphere (e.g., Sojka, 1989). On the other hand, empirical models describe the average behavior of measured variables (e.g., Bilitza, 2001). Ionospheric models can act as background models providing stability and filling data gaps in electron density or TEC reconstructions through data assimilation (e.g., Jakowski et al., 2011b). Powerful assimilative approaches have been developed in recent years by combining actual measurements with background information from empirical or physical models or at least from empirical orthonormal functions related to Chapman profiles (e.g., Angling & Cannon, 2004;Bust et al., 2004;Schunk et al., 2004;Sezen et al., 2013;Elvidge & Angling, 2019;Themens et al., 2020).
Although future developments comprise assimilative datadriven physics-based models, data merging into empirical models is still justified (e.g., Angling & Cannon, 2004) because theoretical models are extremely time-consuming and require numerous input parameters characterizing complex coupling processes with the thermosphere, plasmasphere, and magnetosphere, which are not available in most cases.
Key observables of the ionosphere, such as the peak electron density (NmF2), total electron content (TEC), and the peak electron density height (hmF2), have been measured over more than five decades. Measurement techniques comprise, in particular vertical sounding radar and trans-ionospheric radio beacon techniques. Ground and space-based techniques utilizing dualfrequency signals of Global Navigation Satellite Systems (GNSS) have shown their high efficiency in monitoring TEC and vertical electron density distribution on a global scale already over more than two decades. Consequently, huge data pools are currently available for developing empirical ionospheric models considering special user needs. Thus, the International Reference Ionosphere (IRI), which has been developed for more than five decades by the international community, is a rather complex model that describes different physical quantities and their relationships up to plasmaspheric heights (Bilitza, 2001;Bilitza & Reinisch, 2008). Certified as a reference model by the International Organization for Standardization (ISO), IRI can be used both for physical studies (e.g., Angling & Cannon, 2004) as well as for simulations in engineering tools (e.g., Maltseva et al., 2012;Feltens et al., 2018). Because ionospheric impacts on radio wave propagation depend primarily on the electron density along the ray path, correction models for GNSS users focus on describing only electron density or TEC behavior as a function of geophysical and solar irradiation conditions. Whereas GNSS such as Global Positioning System (GPS) and Beidou correct ionospheric propagation delay by TEC models (Klobuchar, 1987;Yuan et al., 2019), the European Galileo system utilizes NeQuickG Orús Perez et al., 2018) that provides 3D electron densities along ray paths for estimating propagation errors. Whereas TEC models (Jakowski et al., 2011a;Hoque & Jakowski, 2015;Hoque et al., 2018) are robust and, thanks to only a few coefficients easy to handle, electron density models like NeQuick are more complex and require integration of the electron density along used ray paths. However, on the other hand, TEC models require a mapping function to transfer the vertical TEC provided by the model to slant ray paths which is commonly an additional source of errors (Hoque & Jakowski, 2013). Consequently, reliable simulation studies of the ionospheric impact on radio wave propagation require global 3D electron density and geomagnetic field models. This is a strong motivation to develop a climatologic global three-dimensional electron density model that combines accuracy, robustness, and fast computation. Here we present the empirical model NEDM2020 of the global electron density distribution in the ionosphere/plasmasphere system. Like other climatology models, NEDM2020 can be used directly for climatological estimations of the 3D electron density development as a function of the level of solar activity. However, if merged with actual data obtained from different sources like integral or in situ measurements of the electron density, the resulting reconstruction will be able to provide a more realistic 3D plasma distribution following the dynamics of space weather drivers such as solar irradiation, particle precipitation, and the convection electric field. Applied in this way, the model is well suited to be operated in space weather services, e.g., for simulations, climatological forecasts and as a background model for near realtime estimation of propagation errors as principally has been shown in earlier publications (e.g., Jakowski et al., 2005;Berdermann et al., 2014). The practicability of this approach has directly been approved by Prol et al. (2021) recently. The model describes the average behavior of the three-dimensional electron density distribution around the Earth up to GNSS orbit heights as a function of the solar activity characterized by the 10.7 cm solar radio flux index F10.7 as a proxy for ionizing Extreme Ultraviolet (EUV) radiation. To get a 3D solution, here we combine former model approaches describing TEC (Jakowski et al., 2011a), NmF2 (Hoque & Jakowski, 2011), hmF2 (Hoque & Jakowski, 2012), and electron density of the plasmasphere  with good accuracy.
The 3D electron density model Neustrelitz Electron Density Model (NEDM2020) presented here is established by superposing electron density models of the plasmasphere and the ionosphere consisting of separate F and E-layer distributions of electron density. In the following, at first, the model approach is described, and then the approach is justified by comparing model results with satellite in-situ electron density and space-and ground-based TEC data covering different altitude ranges. Additionally, the model results are compared with the NeQuick2 model, which is a standard ionosphere model recommended by the International Telecommunication Union (ITU) for trans-ionospheric applications and (a modified version of NeQuick2 known as NeQuickG is) endorsed by Galileo satellite navigation system for single-frequency users, and thus to balance our expectation and model performance.

Modelling of the Earth's ionosphere and plasmasphere
The ionosphere is the ionized part of the Earth's atmosphere, ranging from about 60 km up to several Earth radii, whereas the upper part above about 1000 km height up to the plasmapause is usually called the plasmasphere (Davies, 1990). The electron density distribution of the Earth's ionosphere can be modeled by adding a Chapman layer (Rishbeth & Garriott, 1969) representing the ionospheric F-layer N eF with an E-layer N eE model superposed by exponential decay functions representing the plasmasphere N eP , The Chapman layer function has been found very convenient for describing the vertical structure of the ionospheric electron density (Rishbeth & Garriott, 1969). Although the Chapman layer is based on the ion production by a one-component gas, a more realistic ionosphere can be obtained by varying different layer parameters such as the maximum ionization N m , corresponding height h m , and the scale height H; thus, a wider variety of ionospheric conditions can be created. The ionospheric electron density distribution N eF along height h can be written by the Chapman layer as, where z = (h À hm)/H F . The atmospheric scale height H F , the peak ionization Nm and the corresponding height of peak ionization hm practically define the shape of the Chapman layer, and thus by tuning these parameters, we can generate electron density profiles with varying shapes. The expression of vertical total electron content (VTEC) can be derived by integrating equation (2) with respect to the height and written as (Hoque & Jakowski, 2008), The ionospheric F region around the peak electron density height of about 250-450 km causes the most significant impact M.M. Hoque et al.: J. Space Weather Space Clim. 2022, 12, 1 on trans-ionospheric radio wave propagation. Therefore, the peak electron density of the F-layer Nm and corresponding height hm are known as key parameters for characterizing the ionosphere. For determining global Nm and hm distribution Hoque & Jakowski (2011  satellite missions and about 60 years of data from about 180 worldwide ionosonde stations. The nonlinear approach describes the Nm and hm dependencies on local time, geographic/geomagnetic location, and solar irradiance and activity with a low number of model coefficients and a few empirically fixed parameters. In the present study, the NPDM and NPHM models are used for describing the F-layer distribution. Therefore, a brief description of NPDM and NPHM is given below.

Peak electron density model
The basic modeling approach includes five major dependencies of Nm on the local time, season, geomagnetic latitude, equatorial anomaly, and solar cycle, and they are combined in a multiplicative way as, The factor F i explicitly contains the model functions and coefficients (for details, we refer to Hoque & Jakowski, 2011). The variation with local time contains the diurnal, semi-diurnal, and ter-diurnal harmonic components. The model function (e.g., F 1 ) describing the local time variation includes the dependency on the solar zenith angle. The solar zenith angle can be described as a function of season (i.e., day of the year) and geographic latitude (Davies, 1990). The so-called summer daytime bite out (BO) effect (i.e., decrease of electron density) in the middle latitude is modeled as an additional term in the local time variation. The seasonal variation of Nm is modeled by two components, namely the annual and the semi-annual variation. The geomagnetic field dependency is given by the geomagnetic latitude dependency and modeled by a half-cosine function having the peak at the geomagnetic equator. A dipole representation of the Earth's magnetic field is used for simplicity. The latitudinal distribution of the peak electron density shows a minimum at the geomagnetic equator together with two maxima on opposite sides of the equator near magnetic latitudes of 15-20°N and S. This phenomenon is called the equatorial or Appleton anomaly. The northward and southward ionization crests are modeled by two Gaussian functions whose peaks are approximated at 16°N and 15°S geomagnetic latitudes. In general, the crests are found prominent during daytime and less prominent during nighttime. Considering this, the half-width of the Gaussian function is modeled in such a way that the width becomes maximum during nighttime and minimum at 14:00 LT and thus generating a visible crest during daytime. The solar activity dependence of Nm is formulated by a slowly increasing exponential function. The nonlinear approach needs 13 model coefficients and a few empirically fixed parameters for describing the Nm dependencies on different geophysical conditions.

Peak electron density height model
Following the same basic approach as Nm, a set of nonlinear equations describing the dependencies of hm on local time, season, geomagnetic field, and solar activity are combined in a multiplicative way (for details, see Hoque & Jakowski, 2012). The local time variation contains the diurnal, semi-diurnal, and ter-diurnal harmonic components as well as dependency on the solar zenith angle. Like Nm, the seasonal variation is modeled by the annual and semi-annual components. The latitudinal distribution of hm shows a maximum at the geomagnetic equator but no visible crest on both sides of the equator. Instead, the hm distribution shows a gradual decrease on both sides of the magnetic equator. The solar activity dependence of hm is modeled as a slowly varying exponential function. The nonlinear approach needs 13 coefficients and a few empirically fixed parameters for global modeling of hm.
Total vertical electron content from ground up to GNSS height (~20,000 km) could vary largely depending on geophysical conditions such as geographic/geomagnetic location, diurnal, seasonal, and solar cycle variation. For determining global VTEC, Jakowski et al. (2011a) developed the Neustrelitz TEC Model (NTCM) based on global TEC data from the Center for Orbit Determination in Europe (CODE). In the present study, the NTCM model is used for describing the global VTEC distribution. Therefore, a brief description of NTCM is given below.

Total electron content model
The major dependencies of TEC are combined in a multiplicative way, similar to the Nm model (see Eq. (4)). The factors explicitly contain the model functions and coefficients whose number varies between 2 and 6 for different factors (for details, see Jakowski et al., 2011a). As discussed for Nm and hm models, the variation with local time includes the diurnal, semi-diurnal, and ter-diurnal harmonic components and expressions that consider the dependency on the irradiation angle of the Sun. In general, the daily maximum of TEC occurs at 14:00 LT instead of at 12:00 LT indicating about 2 h delay in the response of the ionospheric/thermospheric system to the daily solar excitation. Therefore, the phase shift of the diurnal component is approximated at LT = 14:00 h. The seasonal variation is composed of annual and semi-annual components. Similar to the Nm model, the geomagnetic latitude dependency is modeled by a halfcosine function with the maximum at the geomagnetic equator. Equatorial crests are modeled by two Gaussian functions whose peaks are approximated at 16°N and 10°S, respectively. The southern crest is 5°closer to the geomagnetic equator compared to the Nm southern crest. This discrepancy is found in the model input database. The TEC model is based on Global Ionospheric Maps (GIM), whereas the Nm model is based on GNSS radio occultation and ionosonde data. The solar activity dependence is modeled as a linear function of F10.7. The empirical approach describes the VTEC dependencies on different geophysical conditions using 12 model coefficients and a few empirically fixed parameters.
It is worthy of mentioning that the global NTCM as well as NPDM and NPHM are climatological and thus describe the average behavior under quiet geomagnetic conditions. Knowing the peak ionization Nm and VTEC contribution of the F-layer, the corresponding scale height H F can be determined by equation (3) as,

E-layer model
For estimating the electron density profile of the E-layer, we follow the approach published by Ching & Chiu (1973): cos where h is the height, hmE = 110 km is the E-layer peak density height, H E = 10 km is the E-Layer scale height, u and k are geographic latitude and longitude, respectively, d is the declination of the sun, v is the solar zenith angle and R is the monthly smoothed Zurich sunspot number defined by equation (12) in which F10 is the solar radio flux index F10.7.
Due to the small width, the E-layer contribution to vertical TEC can be ignored in ground-based GNSS applications. Nevertheless, we have included the E-layer with regard to future three-dimensional ionospheric electron density reconstructions using GNSS radio occultation techniques. The E-layer is of interest because it plays an important role in terrestrial High Frequency (HF) communication (Hervás et al., 2020).

Modeling of the Earth's plasmasphere
At altitudes well above the F-layer, the ionosphere merges to the so-called plasmasphere or protonosphere at a transition height of about 1000 km where the density of H + ion starts to dominate over the density of O + ion (Böhm et al., 2013). The upper boundary is defined by the plasmapause, at which the plasma density drops by one or two orders of magnitude. The plasmapause is located by the McIlwain parameter L at approximately L = 5 in the equatorial plane of the Earth under quiet geomagnetic conditions and can contract to L-shell values of 2 or 3 during severe disturbances.
Recently  developed an empirical Neustrelitz Plasmasphere Model (NPSM) based on electron density data derived from topside GPS measurements onboard the CHAMP satellite. The NPSM approach consists of two exponential decay functions representing an upper part N ePL where plasmaspheric processes related to plasmapause and magnetosphere dominate, and a lower altitude part N ePh where ionospheric coupling takes place, The electron density distribution in the high altitude plasmasphere is closely related to geomagnetic field lines, and therefore the high-altitude part is modeled as a function of the McIlwain parameter L (see Eq. (14)). However, the electron density distribution in the lower part N ePh is closely related to the altitude, and therefore, it is described by an exponential decay function of height h (see Eq. (15)). The quantity N ePL1 represents the basic electron density of the high-altitude part, and H PL is the corresponding plasma scale height. Similarly, the quantity N ePho denotes the basic electron density of the low altitude part, and H Ph is the corresponding plasma scale height. The quantities N ePL1 , N ePho , H PL , and H Ph are separately modeled as functions of geophysical conditions. For details of modeling we refer to . The plasmapause position is fixed at L PP = 5 assuming geomagnetically quiet conditions, and the corresponding radial distance is calculated by, Equation (16) indicates that the plasmapause location varies significantly with the geomagnetic latitude u m showing the highest value at the equator, whereas it vanishes at the poles. The NPSM approach starts computation at the F2-layer peak height and describes the plasmaspheric electron density from about 1000 km to GNSS orbit heights of about 20,000 km. Our investigation shows that NPSM overestimates the plasmaspheric density, especially during low solar activity conditions. The reason may be that NPSM was developed based on CHAMP topside GPS data during the period 2000-2005 while the solar activity conditions were mostly high . To mitigate the overestimation, we reduced the NPSM contribution multiplying the original N eP by a solar flux dependent factor given by, Factor Nep ¼ 0:5 þ 0:5 F 10 À 60 ð Þ 140 : The fixed parameters in equation (18) are optimized by trial and error methods. The optimization is done by comparing model outputs with TEC and electron density in-situ data recorded onboard low Earth orbit (LEO) satellites. The COSMIC/FOR-MOSAT-3 topside TEC data and the Defense Meteorological Satellite Program (DMSP) in-situ data are used for this purpose.
As representative datasets of low solar activity conditions, we used one-day data from Northern winter (1st January 2008) and one-day data from Northern summer (1st July 2008). Similarly, data from 1st January and 1st July of 2014 are used as representative datasets for high solar activity in Northern winter and summer conditions, respectively. For the selected days, the model-derived TEC and Ne values are compared with corresponding TEC data from COSMIC/FORMOSAT-3 and in-situ data from DMSP satellites, respectively. The fixed parameters in equation (18) are chosen so that the model results roughly follow the representative datasets when visually inspected. Note that COSMIC/FORMOSAT-3 and DMSP data from the selected 4 days are not used in the validation study later. The NPSM coefficients are not updated rather, the model output is reduced by a common factor for the whole vertical profile. Thus, only the basic electron densities of high and low altitude parts (see Eqs. (13)-(15)) are modified, and their scale heights are not changed. Our limited investigation shows that during a low solar activity condition (F10.7 = 80 f.u.), the approach reduces the plasmaspheric content up to 4.5 TECU over the equator and 0.02 TECU over the Poles. During a high solar activity condition (F10.7 = 160 f.u.), the approach reduces the plasmaspheric content up to 1.7 TECU over the equator and 0.01 TECU over the Poles. A common driving parameter, namely the solar radio flux index F10.7 is used for operating the ionospheric key parameter models NTCM, NPDM, and NPHM, and the plasmaspheric density model NPSM. In the present study, we superposed the described ionosphere model based on Chapman functions for the F and E-layers with the plasmasphere model. The vertical electron density distribution at certain geolocation is computed by the following steps.  (2)). 3. E-layer profile (N eE values from ground up to GNSS height) is computed (see Eqs. (6)- (12)). It is worth mentioning that the VTEC contribution from the E-layer is practically negligible in applications (<1 TECU). 4. For computing the final electron density distribution, the NPSM, F-layer, and E-layer profiles are superposed (added).
Note that no separate merging algorithm is needed for superposing ionosphere and plasmasphere models. This is possible since the peak density heights and scale heights of E-layer and F-layer differ largely. The E-layer peak density height (110 km) is far below the F-layer peak density height (>250 km), and the E-layer scale height is relatively small (10 km) compared to the F-layer scale height. Again, on the other hand, the plasmaspheric peak density is about 2 orders of magnitudes smaller than the F-layer peak density, and on the other hand, the plasmaspheric scale height is about 1-2 order of magnitude larger than the F-layer scale height. We see that at around 1000 km, the F-layer electron density reduces sharply to the level of electron density given by the plasmasphere model (see Fig. 1). Figure 1 shows the electron density distribution obtained by the ionospheric E-layer (green curve), F-layer (black curve), and plasmasphere model (blue curve), as well as the combined electron density distribution (red broken line curve) during low and high solar activity conditions represented by F10.7 = 80 and 160 f.u., respectively (see left and right panel). We arbitrarily choose the day of the year, DOY = 1 and universal time UT = 14:00 h and geographic latitudes as 0, 30, 60, and 85°( see bottom to top panel) along the 0°meridian for electron density profile generation. Our investigations show that the superposition works well also at other geophysical conditions represented by different F10.7, DOY, UT, geographic latitude, and longitude values. However, we found that the E-layer model gives a higher electron density value than that of the F-layer at and around Poles. Our investigation shows that the NmF2 model gives very low values at and around Poles which makes the E-layer density dominate over the F-layer density. Due to the lack of NmF2 observations in the Polar region, the NmF2 model could not model the peak density properly. This issue is addressed in Sector 5 where model scopes and scope of improvement are discussed.
Subplots in Figure 1 show that the electron density given by the Chapman layer (black curve) decreases rapidly above the F layer at about 1000 km, whereas the electron density given by NPSM (blue line) decreases slowly, maintaining the plasma density. We see that the developed model maintains the plasma density at upper heights after combining E-and F-layer profiles with the plasmasphere model.
Each subplot of Figure 1 lists estimated scale heights H E , H F , H PL , and H Ph for E-layer, F-layer, plasmasphere L-shell and height dependent layers, respectively. The E-layer scale height is fixed at 10 km, whereas the F-layer scale height is varied between 50 and 90 km. The plasmasphere scale heights H PL and H Ph are derived from respective models (see , and their values (H PL =~4000-4300 km and H ph =~700-1200 km) are 1-2 order of magnitudes larger than the F-layer scale height.
In the following, the quality of NEDM profiles are validated in terms of electron density and TEC estimates at different altitude ranges against reference in-situ and TEC datasets, respectively.

Validation database and data sources
We used satellite-based in-situ electron density and TEC data as well as ground-based TEC data for NEDM2020 validation. For this, electron density in-situ data from DMSP, Swarm (A, B, C satellite), Van Allen Probes and ICON (Ionospheric Connection Explorer) missions, topside TEC data from COSMIC/FORMOSAT-3 mission, bottom side TEC data from TOPEX/Poseidon altimeter mission and ground-based TEC data from International GNSS Service (IGS) are used as reference datasets. Additionally, the peak electron density data from COSMIC radio occultation (RO) profiles are used for comparisons. Figure 2 shows a sketch of electron density variation along with height and altitude coverage of used reference datasets. The GNSS, COSMIC, DMSP, Swarm-A, -B, -C and TOPEX/Poseidon orbits and the altitude coverages of spaceborne in-situ and TEC, and ground-based TEC data are indicated. In the following, we briefly describe the used databases and related missions.
The European Space Agency (ESA) launched the Swarm satellite mission on 22nd November 2013, aiming to measure the strength and direction of the Earth's magnetic field to new levels of precision. The mission consists of three identical satellites, namely Swarm-A, Swarm-B, and Swarm-C. Each Swarm satellite carries the Electrical Field Instrument that consists of the Langmuir Probe (LP) and the Thermal Ion Imager (TII). The LPs can measure electron temperature, electron density, and the electric potential of the ambient plasma, whereas the TII can characterize the ion drift and velocity in high resolution to determine the electric field around the Earth. Knudsen et al. (2017) highlighted a few novel scientific applications of Electric Field Instrument (EFI) data. For model validation, we used Swarm LP data for selected Low Solar Activity (LSA) and High Solar Activity (HSA) time periods which were downloaded from the ESA portal (see Table 1). Lomidze et al. (2018) assessed the accuracy and reliability of the Swarm LP data by using nearly coincident measurements from low-and middle-latitude incoherent scatter radars (ISRs), low-latitude ionosondes, COSMIC satellites. They found that the Swarm LPs systematically underestimate plasma frequency by about 10% (0.5-0.6 MHz). Considering this, the adjustment of LP data for each Swarm satellite is performed according to the corrections given in Lomidze et al. (2018). Note that the correction factor is given for the  Table 1).
The US Air Force has launched many DMSP spacecrafts since 1960's in sun-synchronous near-polar orbits of altitude about 800 km and inclination of 99°. The DMSP satellites are equipped with SSIES (Special Sensors for Ions, Electrons, and Scintillation) instrument package (Hairston & Heelis, 1996, Hairston et al., 1998 in order to monitor the thermal plasma in the topside ionosphere. The SSIES package includes a Faraday scintillation cup, a Retarding Potential Analyzer (RPA), a Langmuir probe, and an Ion Drift Meter (IDM) for measuring either the total ion density or the electron density in addition to a variety of other ionospheric parameters. The used DMSP data are originally processed by the Center for Space Sciences at the University of Texas (Dallas) and are obtained through the Madrigal portal (see Table 1). However, Madrigal files include the electron density data only from the scintillation cup instrument. Garner et al. (2010) conducted a comprehensive study on the quality of the DMSP topside electron density measurements, including data from the scintillation cup. They recommend undertaking a number of data quality control measures prior to the use of DMSP in situ electron density data. As Garner et al. (2010) recommended, in the present study, the scintillation cup data were rejected whenever the measurements were within the range of <200 el/cm 3 or >10 6 el/cm 3 . We found high variability in the data when satellites were flying over the Poles. To remove any possible erroneous measurements, we smooth the data along a continuous meridional arc by using a statistical filter that computes a median value at each measurement point, considering a region of 50 km surrounding that measurement point. That means at each measurement point the past and future data are considered within the same track for computing the median value. Our investigation shows that this approach successfully removes outliers detected, especially in the Polar region. We found that data from different combinations of DMSP satellite series were available during the investigation period. Data from the F13 to F17 satellite series were available during January and July of 2008. Data from the satellite series F15-F18 were available during January 2014, whereas data from F15-F19 were available during July 2014. It is worthy of mentioning that these satellites do not all have the same orbit local time, so changes in the availability of certain satellites can cause substantial changes in the mean electron density behavior for year-to-year or season-to-season.
However, this will not affect the model comparison since both NEDM and NeQuick2 are run for the same input conditions.
We used in-situ data from Van Allen Probes for model validation. NASA's Van Allen Probes were two robotic spacecraft used to study the Van Allen radiation belts surrounding the Earth. The spacecraft was launched on 30th August 2012 and deactivated in 2019 when they ran out of fuel. Both spacecraft flew at a geocentric, highly elliptical orbit of very low inclination of about 10°with a Perigee altitude of 618 km and an Apogee altitude of 30,414 km. Thus, Van Allen Probes cover a large altitude range, approximately 600 -30,500 km, for electron density comparisons. However, we used data up to 20,000 km altitude since our model is developed for use up to that height. The used Van Allen Probes electron density data are processed by Zhelavskaya et al. (2016) using the algorithm NURD (Neural-network-based Upper hybrid Resonance Determination). We used in-situ data from instruments RBSP-A and RBSP-B (Radiation Belt Storm Probes) for selected time periods during the LSA and HSA conditions (see Table 1) for the validation study.
Additionally, we used electron density in-situ data from the ICON (Ionospheric Connection Explorer) mission. NASA launched the ICON satellite on 11th October 2019, aiming to study the interaction between Earth's weather systems and space weather driven by the Sun (Immel et al., 2018). The spacecraft flew at a low inclination of about 27°with a Perigee altitude of 575 km. We used ionospheric plasma density obtained from the ion velocity meter (IVM) instrument at the location of the ICON satellite (mostly at low and equatorial latitudes). Heelis et al. (2017) described the configuration and operation of the IVM instrument as well as the procedures of determining the ion drift velocity. The data are processed by the University of California (UC) Berkeley's Space Sciences Laboratory, and data for selected time periods are downloaded from the UC Berkeley site (see Table 1). Although the used data were part of a preliminary release, UC Berkeley reported that the IVM was operating with low noise and producing highquality outputs throughout the period. Furthermore, good quality data were ensured by selecting times when ICON_L27_R-PA_FLAG and ICON_L27_DM_FLAG flags were zero giving the highest quality measurements.
COSMIC/FORMOSAT-3 is a joint Taiwan-U.S. mission launched during April 2006 with a constellation of six micro-satellites on the low Earth orbit (LEO) of 800 km and inclinations of 72°. Each COSMIC satellite is equipped with dual-frequency GPS receivers designed to collect topside GPS navigation data for precise orbit determination (POD) and bottom side RO measurements. The ionosphere is a frequency dispersive medium, and therefore, the dual-frequency navigation data can be used for computing TEC along GPS-LEO paths. Again, the dual-frequency RO data can be used for computing electron density profile up to LEO height using Abel inversion techniques. The COSMIC data are being processed by the University Corporation for Atmospheric Research (UCAR). The used POD TEC measurements and RO electron density profiles were obtained from the COSMIC Data Analysis and Archive Center (CDAAC, see Table 1). Yue et al. (2011) described the details of the differential code bias (DCB) estimation and LEO-based slant TEC derivation method, including the main error sources, the multipath calibration, the leveling of phase derived TEC to the pseudo-range derived TEC for ambiguity removal. They validated the slant TEC observations, comparing with empirical models and analyzing the TEC difference of COSMIC collocated clustered observations. Yue et al. (2011) found the accuracy of the LEO slant TEC to be within 1-3 TECU depending on the mission. To remove unrealistic data, the presence of negative TEC is checked in the complete phase-connected arc measurements. If there is any negative TEC value, the whole arc data are rejected. Slant TEC (STEC) is converted to vertical TEC (VTEC) using the mapping function given by Foelsche & Kirchengast (2002) with a single-shell height of 1400 km. Recently, Yuan et al. (2020) found that the performance of such a mapping function significantly improved for a single-shell height of 1400 km instead of a lower thin-shell height such as 1000 km. However, to minimize mapping function-related errors, those TEC observations with elevation angles lower than 45°were excluded. The location of Ionospheric Pierce Point (IPP) at the thin-shell height is considered the measurement location for the VTEC. We selected data from LSA and HSA periods during the summer and winter months, considering the natural variability of the ionosphere (for details, see Table 1). The COSMIC/FORMOSAT-3 mission is officially decommissioned on 1st May 2020 after 14 years of operation. Its successor, namely the COSMIC-2/FORMOSAT-7 constellation of six satellites, was launched on 25th June 2019. As a source of ground-based vertical TEC data, we used GIM from the Center for Orbit Determination in Europe (CODE) (Schaer et al., 1998a). GIM provides estimates of vertical TEC from ground up to GNSS satellite orbit height (~20,000 km). Daily IONEX (ionosphere exchange) files of GIM are downloaded from the Crustal Dynamics Data Information System (CDDIS) archive. Reference GIM VTECs are computed at the same IPP location as given by the COSMIC-GPS ray path using IGS-recommended spatial and temporal interpolation methods (Schaer et al., 1998b). The selected data period and data source are given in Table 1.
Dual-frequency altimeter missions such as TOPEX/ Poseidon and Jason become an excellent source of vertical total electron content data independent from GNSS-based measurements. The TEC is obtained by combining dual-frequency measurements over the oceans where the surface reflection of the signals is supported. The altimeter measurements are naturally vertical, so no mapping functions are required to convert them to the vertical. We used VTEC data from TOPEX/Poseidon mission for selected time periods during the LSA and HSA conditions (see Table 1) for NEDM2020 validation. Table 1 gives a summary of used data types, periods, satellite missions, and corresponding data sources.

Electron density model validations
In the following, we validated NEDM2020 results against the reference electron density and TEC data from spaceborne and ground sensors. Additionally, NEDM performance is compared with another 3D electron density model, namely NeQuick2. NeQuick2 is an ionospheric electron density model developed at the Aeronomy and Radiopropagation Laboratory of the Abdus Salam International Centre for Theoretical Physics (ICTP), Trieste, Italy . NeQuick2 is a three-dimensional model developed to calculate the electron concentration at any given location in the ionosphere for trans-ionospheric applications. It stores ITU-R (ITU Radio-Communication Sector) maps for the peak ionization and peak height information. The maps describe the spatial and temporal variability of the peak ionization parameters as a function of geographic latitude, longitude, universal time as well as the function of modified magneticdip latitude. The source code of the 3D electron density NeQuick2 model can be obtained from the ICTP site-https:// t-ict4d.ictp.it/nequick2/source-code. Figure 3 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference peak density values retrieved from COSMIC/FORMOSAT-3 RO data (indicated by FM). The left plot shows model comparisons during a medium solar activity year 2011, whereas the right plot shows comparisons during a high solar activity year 2014. We used full years of data that are available from all six satellites. The peak electron density Nm values are computed by NEDM2020 and NeQuick2 at the same time interval and geolocation as COSMIC RO data. For each year, all data are sorted into 2.5-degree latitude bins, and corresponding mean and Standard Deviation (STD) are computed and plotted in Figure 3, where error bars represent corresponding STD values. Comparing with COSMIC RO data, we see that both NeQuick2 and NEDM2020 underestimate the Nm values during the high solar activity year 2014.

Validation against RO peak density data
We computed the differences between the model outputs and reference COSMIC data (i.e., Nm model À Nm obs ) and estimated the mean, STD, and Root Mean Squared (RMS) differences. The values are given for each model in the plots. The RMS and STD estimates show that both models perform similarly in the global average. We find that the mean residual is about 0.07 À 0.17 Â 10 11 el/m 3 less for the NEDM2020, whereas the STD and RMS residuals are about 0.15 À 0.22 Â 10 11 el/m 3 less for the NeQuick2 model.
Abel inversion method is commonly used for retrieving electron density profiles from RO data. However, the assumption of spherical symmetry ignoring horizontal gradients used in the Abel inversion technique may degrade the reconstruction quality, especially at equatorial regions on both sides of the geomagnetic equator. Shaikh et al. (2018) found that the accuracy of the RO profiles changed with the azimuth of the occultation plane. They obtained a better accuracy for the azimuth range of 80°-100°, for which the effect of the horizontal gradient is minimal. It is well known that the horizontal electron density gradient is stronger along the north-south than the east-west direction. However, since the model values are compared with the monthly mean of RO data considering the full azimuth range of 0-360 of occultation plane, the data quality should not change the mean behavior.

Validation against Swarm in-situ observations
The left, middle and right panels of Figure 4 compare the performance of the NEDM2020 with the NeQuick2 model with respect to the reference in-situ data from Swarm-A, -B, and -C satellites, respectively. The top two panels show model comparisons during an HSA year 2014, whereas the bottom two panels show comparisons during an LSA year 2019. Again, for each year, we selected two months, January and July, representing winter and summer seasons in the Northern hemisphere. The seasons are opposite in the Southern hemisphere. The electron densities at the Swarm satellite height (approximately 430-520 km) are computed by NEDM2020 and NeQuick2 at the same time epoch and geolocation. For each month, all data are sorted into 2.5-degree latitude bins, and corresponding mean and STD are computed and drawn as error bar plots in Figure 4. Comparing with Swarm satellite data, we see that during the Northern winter and HSA condition (top panel), both NeQuick2 and NEDM2020 overestimate the monthly mean electron density, although the values lie within the error bars of Swarm data. However, the performance improves during the Northern summer for both models (2nd panel). The bottom two panels show that NEDM2020 mostly underestimates, whereas NeQuick2 follows the observations well during the LSA condition independent of seasons. Additionally, we computed the differences between the models and reference Swarm data and estimated the mean, STD, and RMS differences for each month. The values are given in each subplot for comparison. We find that during selected HSA and LSA periods, NeQuick2 performs consistently better than the NEDM2020. During the HSA period, the mean, STD, and RMS residuals are up to 0.12 Â 10 11 , 1 Â 10 11 , and 1 Â 10 11 el/m 3 less for the NeQuck2 compared to the NEDM2020 values. The corresponding residuals are up to 0.17 Â 10 11 , 0.21 Â 10 11 , and 0.21 Â 10 11 el/m 3 less for the NeQuck2 model during the LSA period. During the selected LSA period the residual error statistics are comparatively less than those during the HSA period.
The Swarm satellites undertook an orbital maneuver in early 2014 to transfer the satellites to their operational orbits, and the final constellation of the mission was achieved on 17th April 2014. The change in the orbit height due to orbital maneuver was investigated for all 3 satellites for the investigation period January and July in 2014 and 2019. Our investigation shows that the orbit height of Swarm A and C varies within~510-430 km, whereas the orbit height of Swarm B varies withiñ 520-490 km during the investigation period. More specifically, Swarm A and C measurements correspond to altitude range~480-510 and~460-480 km during January and July 2014, respectively, and~430-450 km during both months in 2019. Swarm C measurements correspond to altitude rangẽ 490-510 and~510-530 km during January and July 2014, respectively, and~500-520 km during 2019. It is worthy of mentioning that both NEDM2020 and NeQuick2 values are computed considering the actual measurement height. Thus, each subplot in Figure 4 represents the model validation in a slightly different altitude range as mentioned above. Figure 5 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference in-situ data from DMSP satellites. The left panel shows model comparisons during an LSA year 2008, whereas the right panel shows comparisons during an HSA year 2014. As before, each year, we selected two months January and July, representing the Northern winter and summer, respectively.

Validation against DMSP in situ observations
The electron densities are computed by NEDM2020 and NeQuick2 at the same time epoch and geolocation as DMSP satellites data (approximately at 850-870 km altitude). As before, for each month, all data are sorted into 2.5-degree latitude bins, and corresponding mean and STD are computed and plotted as error bars in Figure 5. Compared with DMSP satellite data, we see that NEDM2020 underestimates the electron density values during both seasons and years. The top panel shows that during Northern winter, NeQuick2 underestimates the electron density mainly in the Southern hemisphere. During Northern summer, NeQuick2 values are closer to observation values. We computed the differences between the models and reference DMSP data and estimated the mean, STD, and RMS of differences. The values are given for both models in  Figure 6 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference in-situ data from Van Allen Probes. Due to flying at a very low inclination orbit of about 10°, the Van Allen data only covers low latitude regions on both sides of the equator. However, the data covers a large altitude range, approximately 600-30,500 km, for electron density comparisons. The left panel shows model comparisons with the rbsp-a data, whereas the right panel shows comparisons with the rbsp-b data. Depending on data availability, we considered data from a medium solar activity year 2012 and a high solar activity year 2014. As before, we selected winter (Oct 2012, Jan 2014) months and summer months (Jul 2014).

Validation against Van Allen Probes data
The electron densities are computed by NEDM2020 and NeQuick2 at the same time epoch and geolocation as Van Allen Probes data (restricted to 600-20,000 km altitude). For each month, all data are sorted into 100 km altitude bins, and  Figure 6. In order to obtain representative results, only those bins which contain equal to or more than 50 data points are considered. Bins at a lower height (~<4200 km) contained a small number of reference data and are omitted from the plots and statistics. The left and right panels of Figure 6 show that NEDM2020 underestimates the electron density values at upper height (~>10,000 km), whereas it overestimates at lower height during all seasons and years. NeQuick2 underestimates the electron density at all altitudes during all seasons and years. We computed the differences between the models and reference Van Allen data and estimated the mean, STD, and RMS of differences. The values are given for both models in each subplots. Comparing residual statistics, we find that NEDM2020 gives smaller mean and RMS residuals compared to NeQuick2 in all cases. However, NeQuick2 gives smaller values for STD residuals in most cases. Figure 7 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference electron density data from the ICON mission. The reference data represent electron density estimates at a satellite location at about 513 km height. The left and right panels show model comparisons during December 2019 and May 2020, which correspond to the Northern winter and summer months, respectively. The ICON satellite is recently launched, and we do not have data from a high solar activity period. Again, due to low orbit inclination (27°), the data are limited to low latitude regions. The electron density values are computed with NEDM2020 and NeQuick2 at the same time epoch and geolocation as the ICON data. For each month, the data are sorted into 2.5-degree latitude bins, and corresponding mean and STD are computed and plotted as error bars in Figure 7. Statistical estimates of model differences (Ne model À Ne ICON ) are computed and given in each plot. Comparing with the ICON electron density data, we see that both NeQuick2 and NEDM2020 overestimates the Ne values during both winter and summer seasons. The global residual statistics show that NEDM2020 performs better than NeQuick2 in terms of mean and RMS residuals during both months. However, NeQuick2 performs better when considering STD residuals during both months. Figure 8 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference topside TEC data from COSMIC satellites. The left and right panels show model comparisons during a low and a high solar activity year 2008 and 2014, respectively. The top and bottom panels give comparisons during January and July, which corresponds to Northern winter and summer, respectively. The vertical TEC above COSMIC satellites (approximately at 760-870 km altitude) up to GPS satellite height of about 20,000 km are computed by NEDM2020 and NeQuick2 at the same time epoch and geolocation as COSMIC satellites. Monthly data are sorted into 2.5-degree latitude bins, and corresponding mean and STD are computed and plotted as error bars in Figure 8 as before. We see that in all four cases, both the NEDM2020 and NeQuick2 underestimate the topside TEC when compared with the reference COSMIC data. Comparing monthly mean TEC variation with latitude, we see that at low and middle latitudes on both Recently, Kashcheyev & Nava (2019) analyzed the difference between NeQuick2 and COSMIC TEC values and found NeQuick2 to underestimate the topside ionosphere and plasmasphere TEC. They found a global negative bias of À3.73 TECU with a spread of 2.58 TECU considering all data in the period 2006-2018. Although we considered a subset of periods, NeQuick2 shows similar performance as reported by Kashcheyev & Nava (2019).   Figure 9 compares the performance of the NEDM2020 with the NeQuick2 model with respect to the reference ground TEC data from CODE. The left and right panels show model comparisons during a low and a high solar activity year 2008 and 2014, respectively. The top and bottom panels give comparisons during January and July, which correspond to Northern winter and summer, respectively. The vertical TEC from ground up to GNSS satellite height of about 20,000 km are computed by NEDM2020 and NeQuick2 at the same time epoch and geolocation (2.5°Â 5°latitude and longitude grids) as the CODE GIM. Monthly data are sorted into 2.5-degree latitude bins and corresponding mean and STD are computed and plotted as error bars in Figure 9. Comparing monthly mean TEC variation with latitude, we see that at low latitude on both sides of the equator and Northern middle latitude, the NEDM2020 overestimates TEC during the Northern winter (see top panel). We see that NEDM2020 underestimates TEC at middle and high latitudes in the Northern hemisphere during the summer (see bottom panel). Global residual statistics show that NeQuick2 and NEDM2020 perform very similar during the LSA year 2008, whereas during the HSA year 2014, NEDM2020 performance is slightly better than NeQuick2 in terms of global mean and RMS residuals. For each year, all data are sorted into 2.5-degree latitude bins, and corresponding mean and STD are computed and plotted as error bars in the left panel of Figure 10. Statistical estimates of model differences (TEC model À TEC TOPEX ) are computed and given in each plot. Comparing with the altimeter TEC data, we see that both NeQuick2 and NEDM2020 underestimate the TEC values during high and low solar activity periods.

Validation against altimeter TEC data
The global residual statistics show that NEDM2020 performs better than NeQuick2 in terms of mean, whereas NeQuick2 shows less STD and RMS residuals during 2003. However, during 2008 NeQuick2 gives a smaller mean residual, whereas NEDM2020 gives a smaller STD residual. The underestimation of TEC by ionosphere models, when compared with the altimeter data, is also reported by Migoya Orué et al. (2008). They compared TEC values obtained from TOPEX/Poseidon satellite with the International Reference Ionosphere (IRI) 2001 and the NeQuick option for the IRI topside. They analyzed TOPEX/Poseidon TEC measurements from high and middle to low solar activity years (2000 and 2004) by binning the region every five degrees of modified dip latitude. They found that the IRI2001 and IRI (NeQuick topside) estimates are consistently lower than the TOPEX/Poseidon measurements during HSA. They found both model versions underestimate the magnitude of the equatorial anomaly crest by up to 20 TECU in some cases. However, it should be noted that the TOPEX/Poseidon TEC measurements contain systematic biases. The presence of systematic biases of about 2-3 TECU in the TEC data was first reported by Ho et al. (1995). Later Orús-Pérez et al. (2002) found the TOPEX VTEC measurements to have a systematic bias of about 2-5 TECU above the real ionospheric VTEC values when compared with other independent VTEC measurements. Similarly, Jee et al. (2004) assumed the TOPEX TEC bias up to 5 TECU. However, for the current analysis, we did not consider this bias since we confined our interest only to the general trends that the data and the models displayed.
5 Model scopes and scope of improvement NEDM2020 is particularly developed to calculate the electron concentration at any given location and time in the ionosphere for trans-ionospheric applications. The 2D key parameters models such as the NTCM, NPDM, and NPHM, as well as the combined 3D model NEDM2020, are optimized for the solar radio flux F10.7 <= 200 flux units. However, the model output is not restricted to F10.7 <= 200 flux units, rather the model still provides outputs with limited accuracy. Since the constituent models (e.g., NPSM, NTCM) are developed based on GNSS data covering altitude range up to about 20,000 km, the recommended use of the NEDM is the same altitude range that is about 20,000 km.
Our comparison with the DMSP in-situ data shows that NEDM2020 underestimates electron density in the high latitude and polar regions, especially during the summer months. When comparing with COSMIC topside TEC data, we found NEDM2020 to underestimate TEC in the middle and high latitudes, including the Polar region during all seasons. One reason may be that the used NTCM, NPDM, NPHM, and plasmasphere model perform poorly at the high latitude and polar regions. NTCM is based on GPS ground TEC data, NPDM and NPHM are based on GPS RO data and ground ionosonde data, and NPSM is based on GPS topside TEC data. On one hand the high latitude and polar regions are sparsely covered by GPS receivers, and on the other hand, the low inclination (55°in relation to the equator) of GPS satellites limits the number of satellites that can be tracked from those regions. In the end, we have a limited database at high latitude and Polar region   Clim. 2022, 12, 1 compared to that at middle and low latitude regions. However, during the last decade, many Galileo satellites have been launched and therefore, the chance of tracking more satellites at high latitudes is increased, although Galileo satellites have a similar inclination (56°) like GPS. Again, nowadays, many receivers are tracking GLONASS satellites with a slightly higher inclination of 64.8°.
Although, in general, comparable with the NeQuick2 model, the validation results of NEDM2020 against measurements from different instruments have still shown a number of deficiencies belonging to different altitude and geographical regions, partially referring to different model components. This information will help to improve subsequent modeling efforts in a synthesized approach of a compact 3D electron density model. Thus, this study has a twofold meaning: (a) it demonstrates the candidacy of NEDM2020 to be used as a background model for 3D electron density reconstructions, (b) it indicates specific deficiencies which should be mitigated in the separate models used and the synthesized approach we plan to develop. This gives justified hope that further studies allow further improvement of separate models like NTCM, NPDM, NPHM, and NPSM. In particular, it became clear that improving NPSM requires extended database including a full solar cycle to avoid additional corrections for low solar conditions as applied in this study.
Therefore, there are scopes of improvement for different key parameter models as well as for NEDM2020. Again, the used ionospheric key parameter models were developed based on GNSS space-and ground-based data and ionosonde measurements mainly from the solar cycle 23 during the 1998-2009 period. During the last decade, the availability of ionospheric data from different sensors has increased in many folds. On one hand, the ionosphere has been sounded by a large number of sensors providing a vast database. On the other hand, the accuracy of the measuring techniques has been improved significantly. Again, the used ionospheric key parameter models use a simple dipole model for geographic to geomagnetic latitude conversion. Our investigation shows that the dipole model has limited accuracy in reproducing ionospheric variability at and around the geomagnetic equator. The use of a more sophisticated geomagnetic field model such as the International Geomagnetic Reference Field (IGRF) model can significantly improve the model performance.

Conclusions
In the present study, we discussed a new electron density model, NEDM2020, for possible use in space weather services and mitigation of propagation errors for trans-ionospheric signals. The empirical model describes the ionosphere and plasmasphere up to GNSS orbit height of about 20,000 km. The electron density distribution is modeled by combining a Chapman layer representing the ionospheric F-layer, an E-layer layer model, and a plasmasphere model. Chapman layer parameters such as peak density, peak height, and vertical total electron content are derived from already existing NTCM based empirical models NPDM, NPHM, and NPSM. No separate merging algorithm is used for superposing ionosphere and plasmasphere models. Such an approach is justified because the peak density and scale height of E-layer, F-layer, and plasmasphere differ largely, and peak density heights are widely separated. Our studies show that the combined model gives continuous and smooth vertical electron density profiles over the entire globe during different geophysical conditions.
The new model has been comprehensively validated by using multi-instrument and multi-satellite space-and groundbased data. In order to scale/judge the NEDM2020 performance, we compared the model performance with the electron density model NeQuick2.
Comparing with reference RO data from COSMIC/FOR-MOSAT-3 missions, we found that both NEDM2020 and NeQuick2 perform similar on global average. Comparing with Swarm in-situ data, we found that during the selected HSA and LSA periods, NeQuick2 performed consistently better than the NEDM2020. We found that at DMSP satellite height of about 850 km, both models gave similar residuals in terms of RMS and STD values and thus performed very similar. Comparing with Van Allen's in-situ data, we found that NEDM2020 gave much smaller mean residuals compared to NeQuick2. Comparing with the ICON electron density data, we found that both NeQuick2 and NEDM2020 overestimated the reference Ne values. The global residual statistics show that NEDM2020 performs better than NeQuick2 in terms of mean and RMS residuals. However, NeQuick2 shows comparatively less STD values.
Comparing with COSMIC topside TEC, we found that at low and middle latitudes on both sides of the equator, the NEDM2020 values are closer to the reference data, whereas at high latitudes and Polar Regions NeQuick2 gives better results showing closer values to the reference data. The smaller values of monthly residual statistics indicate that NEDM2020 outperforms the NeQuick2. Global residual statistics with respect to CODE TEC data show that NeQuick2 and NEDM2020 perform very similar during the LSA year 2008, whereas during the HSA year 2014, NEDM2020 performance is slightly better than NeQuick2 in terms of global mean and RMS residuals.