A Global Empirical Model of Electron Density Profile in the F Region Ionosphere Basing on COSMIC Measurements

The topside ionosphere accounts for a dominant part of the ionospheric total electron content, whereas accurate global modeling of topside ionospheric electron density (Ne) profile has not been fully achieved. In this study, a high precision Ne profile model, named α‐Chapman Based Electron Density Profile Model (α‐Chapman‐Based‐EDP), was built by using ∼4.5 million Ne profiles from the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC‐1) radio occultations. We first describe each of the profiles using five parameters of the α‐Chapman function, that is, peak density (NmF2) and height (hmF2) of F2 layer, scale height (Hm) as well as its altitude change rates, and then built a model for each of the parameters as a function of latitude, longitude, month, local time, and solar activity, through Empirical orthogonal function (EOF) analysis and Fourier expansion. Combining all the five models, we construct the α‐Chapman‐Based‐EDP. Compared with observations from COSMIC‐1 and ‐2, the model captures the ionospheric climatology well, such as solar activity dependence, seasonal variation, and spatial pattern, including the equatorial ionization anomaly and midlatitude trough as well as their variabilities. Our model can describe nearly 80% variability of Ne in F region. In contrast, the IRI2016 cannot well reproduce these characteristics, with errors higher than our model. The potential applications of our model were also discussed. A dense matrix data calculated by the model will be released in https://www.researchgate.net/profile/Qiaoling_Li5 with the permissions of COSMIC organizations.

Researchers mainly built the topside Ne profile with or without separating ion constitutes (Zhu et al., 2015). From hmF2 to the upper transition height, the dominant ion species shifts from O + to H + and He + . One possible option is to build the topside profile by the sum of ion components (González-Casado et al., 2013;Kutiev & Marinov, 2007;Kutiev et al., 2010Kutiev et al., , 2009Stankov et al., 2003;Zhu et al., 2015). Another more common way is direct modeling of topside profile without separating ion components (Bilitza & Reinisch, 2008;Di Giovanni & Radicella, 1990;Nava et al., 2008;Radicella & Zhang, 1995;Rawer et al., 1978). Describing a profile of Ne or ion components at least needs two kinds of information: the key anchor points and the altitude gradient. Different mathematical functions have been developed to describe the profile of Ne or ion component, including Chapman-type, exponential, parabolic, and Epstein functions (Chapman, 1931;Rawer, 1988;Reinisch et al., 2004;Stankov et al., 2003;Venkatesh et al., 2011). The topside Ne profile is modeled using key anchor points (peak density [NmF2] and hmF2; upper transition height) and the altitude gradient of Ne or ion component (scale height, or other thickness parameters) in the profile functions (e.g., Bilitza et al., 2018;Kutiev et al., 2009;Nava et al., 2008;Stankov et al., 2003;Zhu et al., 2015). In addition, several researchers also built the "topside model" by modeling of Ne at different heights (e.g., Huang et al., 2015;She et al., 2017;M.-L. Zhang et al., 2020).
Unlike the bottom Ne profiles, which can easily be measured by the widespread ground-based ionosondes, the topside Ne profiles are mainly measured by incoherent scatter radars (ISR), and satellite-based technique. ISR runs in discontinuous time at few specific locations due to expensive cost, and the in-situ measurements provide Ne at the height of the satellite, so they are used to develop station-specific or altitude-limited models of the topside Ne profile (e.g., Huang et al., 2015;. The topside sounders onboard Alouette and ISIS recorded topside ionograms covering from altitude of hmF2 to ∼1,000 km during the 1960-1980s, and previous studies mainly built the global topside Ne profiles based on these data (e.g., International Reference Ionosphere [IRI] model). The emergence of the low Earth orbiting radio occultation (RO) technique, especially the launch of the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), has largely increased the quantity of the topside Ne profile (Yue et al., 2016). COSMIC's Ne profiles have the advantages of high accuracy, high vertical resolution, and full global coverage (Anthes, 2011), and the Ne profiles have been accumulated over one solar cycle. Thus, the accumulated data provides unprecedented Ne profiles with good spatiotemporal resolution and coverage in the ionospheric F region ).
Due to the limited spatial coverage of the data from topside sounders (Bilitza et al., 2009), previous models displayed non-trivial shortcomings in presenting the topside ionospheric Ne profile. IRI is the most widely used empirical climate model of Earth's ionosphere, and it was recognized as the official standard model for Earth's ionosphere by many international organizations such as International Standardization Organization (Bilitza et al., 2018). It nearly stands for the highest level of the empirical models of the Earth's ionosphere. Previous studies showed that IRI provides reliable NmF2 and hmF2, but large discrepancies are reported for the Ne above hmF2, even the best performed NeQuick choice (Bilitza et al., 2009(Bilitza et al., , 2012(Bilitza et al., , 2006Mengist, et al., 2020;Lühr & Xiong, 2010). The discrepancies reflect the defects of the topside scale height in IRI. An outstanding shortcoming of IRI is the representation of equatorial ionization anomaly (EIA) (Appleton, 1946;Liang, 1947). EIA is a distinct feature in the low latitude ionosphere, covering about half of the global area in 24 h and varies with diverse geophysical conditions . It presents two Ne peaks at about ±15° magnetic latitude and a Ne trough near the magnetic equator. With increasing height, EIA's double peaks merge into a single peak (Eccles & King, 1969;Tulasi et al., 2009). Unfortunately, IRI fails in representing this merging signature (Bilitza et al., 2006). The failure reflects the poor performance of the altitude change rate of the scale height in IRI. Moreover, Liu et al. (2008),  and Li et al. (2019) disclosed the significant longitudinal variation of topside scale height in mid and low latitudes ionosphere based on COSMIC data. In comparison, the topside sounders' data have insufficient and inhomogeneous distribution at longitude (Figures 2 and 3 in Kutiev et al., 2006), which indicates the poor longitudinal performance of scale height in IRI. IRI's above-mentioned shortcomings indicate its poor presentation of the topside scale height and the large influence due to the limited topside sounder data.
The empirical model's performance critically depends on the spatiotemporal resolution and the coverage of the measured data. COSMIC has better spatiotemporal coverage than topside sounders in the F region ionosphere, allowing to develop the substitute Ne profile model of IRI within 180-850 km. Many researchers have built models of NmF2, hmF2, and topside scale height of Ne based on COSMIC data (e.g., Gowtam et al., 2017;He, Liu, Wan, & Wei, 2011;Liu et al., 2019;Pignalberi et al., 2020;Prol et al., 2018;Shu, 2015;Tulasi et al., 2018;Themens et al., 2017;Wu et al., 2016;M.-L. Zhang et al., 2010M.-L. Zhang et al., , 2014, but fewer studies were taken to build the models of Ne profile (Gowtam et al., 2019;M.-L. Zhang et al., 2020). M.-L. Zhang et al. (2020) built a global Ne model in terms of geomagnetic latitude, local time (LT), height, month, and solar activity based on RO Ne profiles from COSMIC, CHAMP, and GRACE. However, Zhang's model did not include longitude. Gowtam (2019) built a global three-dimensional (3-D) ionospheric Ne model based on ionosonde, topside sounders, and RO measurements using an artificial neural network. Gowtam's model used the constant scale height retrieved between hmF2 and hmF2+200 km in the ionospheric F region, and therefore might be subject to large errors at high altitudes given that the scale height is altitude dependent (Liu et al., 2008). Therefore, it is necessary to develop a COSMIC-based Ne profile model considering longitudinal and altitudinal variations of the topside scale height.
Using the first 2.7 million COSMIC profiles and through empirical orthogonal function (EOF) analysis, He, Liu, Wan, and Wei (2011) constructed an empirical model called F2COSMIC, which represents NmF2 and hmF2 as functions of latitude, longitude, month and LT at the low solar activity level. F2COSMIC had great success in extracting the details of non-migrating tidal signatures in the equatorial ionosphere. Also using the COSMIC profiles through EOF analysis, Li et al. (2019) built a global model of the α-Chapman scale height of the lower topside ionospheric Ne profile (within 200-km altitudes above hmF2). The current study tailors and expands the F2COSMIC and the scale height model to delineate the ionospheric variabilities associated with altitude and solar activity. We chose the α-Chapman function with height-varying scale height (Lei et al., 2005;Liu et al., 2008) to describe the Ne profile in the ionospheric F region. The function only takes five independent parameters, and each Ne profile in the ionospheric F region could be described by the five parameters (Lei et al., 2005), which indicates that we can concisely build the global model of Ne profile in the F region by simply building the five parameters' models.
In this study, using ∼12-years about 4.5 million Ne profiles from COSMIC-1, we built a high precision global empirical model of COSMIC-1 Ne profile covering 180-850 km in the ionospheric F region. Our model can well capture the typical Ne climate phenomena in the ionospheric F region (including the merging of the EIA), and its performance is better than IRI in the F region. The model can be a substitute of IRI in the F region.

Data and Processing
As the first constellation of satellites dedicated primarily to RO, COSMIC-1 comprises six identical micro-satellites and is launched into a circular low-Earth orbit on April 15, 2006 (Anthes et al., 2008). Each satellite carried a GPS RO receiver to measure the retardation and bent angle of the radio signal from GPS. The TEC along the receiver-GPS ray could be derived from the retardations of the dual-frequency GPS signals. The successive TEC observations during an occultation event can be inverted to the ionospheric Ne profile along the tangent points by Abel inversion (Cheng et al., 2006;Rocken et al., 2000). Abel inversion performs very well above 180 km altitude and could provide reliable Ne profiles as well as hmF2 and NmF2 (Hu et al., 2014;Kumar et al., 2016;Lei et al., 2007;Schreiner et al., 2007;Yue et al., 2016;Yue et al., 2011. COSMIC-2, an operational follow-on mission to the COSMIC-1, was successfully launched on June 25, 2019 (Anthes & Schreiner, 2019). It is an equatorial constellation of six satellites and each satellite carries an advanced Tri GNSS RO system instrument (Schreiner et al., 2020). COSMIC-2 has been providing high precision Ne profiles within ±45° latitude.
In the current study, Ne profiles from COSMIC-1 (COSMIC-2) were used to build (evaluate) the empirical model. Figure 1 presents the daily profiles of COSMIC and the solar and geomagnetic indices for the analyzed periods. We collected 4,487,884 pieces of Ne profiles from COSMIC-1 during the period from the day of year (DOY) 194 in 2006 to DOY 242 in 2018. The maximum, minimum and mean F10.7 (10.7-cm solar flux index, unit: sfu, 1 sfu equals 10 −22 W·m −2 ·Hz −1 ) are 255, 64.6, and 95.2 sfu, respectively. The data with daily geomagnetic activity index Ap larger than 30 nT are excluded in the modeling. We also did not use the Ne below 180 km altitude. Considering the Ne profile from COSMIC is slanted, the profiles whose longitude span (∆Longitude) within 45° and latitude span (∆Latitude) within 15° are adopted in this study (e.g., Figures 3b and 3c; examples of latitudinal and longitudinal distributions of a Ne profile). To make the utmost of the data, for profile does not satisfy the geographical spans, we only use the lower part of the Ne profile within geographical spans, with a cutoff altitude at 450, 500, ..., or 850 km. We took the COSMIC-2 Ne profiles from October 1, 2019 to December 31, 2020 (1,117,551 pieces profiles) as the independent data set to quantitative evaluate our model and IRI2016. The chosen period is during the low solar activity. The maximum and minimum F10.7 index are 113.1 and 63.4 sfu, respectively.

Solution
As specified in equations (1-4) in Figure 2, the α-Chapman function can well describe the ionospheric Ne profile in the F region in terms of only five independent parameters (Chapman et al., 1931;Lei et al., 2005;Liu et al., 2008). The five parameters include key anchor points of NmF2 and hmF2, altitudinal gradient parameter of scale height (Hm) and its altitude change rates (a 1 : above hmF2, and a 2 : below and equals hmF2). Figure 3 shows an example of α-Chapman function fitting on a Ne profile and the corresponding latitude and longitude distributions. Actually, overwhelmed profiles present well matches as Figure 3a. We took statistical analysis of the correlation coefficients (R) between COSMIC (4,487,884 pieces) and α-Chapman function fitted Ne profiles, results show that 98.03% R is larger than 0.97 ( Figure 1S). The Ne profile model in the F region ionosphere is developed after separately building empirical models of the five parameters of the α-Chapman function. The current study extended F2COSMIC built by He, Liu, Wan, and Wei (2011) and 2PCAFourier-Hm built by Li et al. (2019) to model the five parameters. The data were first fitted to predefined grids through spherical harmonic (SH) analysis and then modeled through EOF analysis and Fourier expansions.

Algorithms
The EOF analysis, also known as Principle Component Analysis, was first invented by Pearson (1901). It decomposes a matrix into a linear combination of a sequence of EOFs determined from the matrix (Jolliffe, 2002). As instanced by Wan et al. (2012) and reviewed by He and Vogt (2018), the EOF analysis grasps typically the inherent characteristics of data with limited EOFs, and it could often separate the interlaced physical contributions into different EOFs. EOF analysis is popular in analyzing and modeling the ionosphere (e.g., Chen et al., 2016a;He, Liu, Wan, & Wei, 2011;Le et al., 2017;Lei et al., 2012;She et al., 2017;Wan et al., 2012;Yu et al., 2015;M.-L Zhang et al., 2014;S.-R. Zhang et al., 2013;Zhao et al., 2005).
The process flow diagram in Figure 2 summarizes our modeling processes which are similar to the construction of the F2COSMIC (as detailed in He & Vogt, 2018). Briefly, our modeling comprises the following steps: (1) Each profile was fitted using the α-Chapman function according to equations (1-4), which is different from the cubic spline regression used in F2COSMIC. Then we built a model for each of the five parameters individually using the same algorithms. In the following steps, we take hmF2 as an example to introduce the detailed algorithms LI ET AL.
10.1029/2020SW002642 5 of 23 (2) We binned hmF2 by month and longitude (glon) using the 45-days and 45° windows with 1 month and 15° resolutions (12×25 bins in total). In each bin, according to equations (5-7), we carried SH regression to fit hmF2 in terms of latitude (glat), LT, and solar activity (F10.7). P n m are associated Legendre polynomials of degree n and order m. C nm and S nm are SH coefficients. C nm and S nm were further denoted as the quadratic regression of F10.7, given that ionospheric parameters exhibit complex dependence on solar activity (Chen & Liu, 2010;Lei et al., 2005;Liu et al, 2006). The expansion coefficients, including h nm , h nm ′,h nm ″, j nm , j nm ′, and j nm ″, were obtained from the least squares fit to equations (5-7). In each bin, the regular 3-D matrix of hmF2 are calculated by bringing the expansion coefficients into equations (5-7) on the 181 × 49 × 5 mesh grids of glat (90°, 89°, …, −90°), LT (0:00, 0:30, …, 24:00), and F10.7 (70, 90, …, 150 sfu). Thus, hmF2 is represented by a five-dimensional (5-D) matrix (12 × 25 × 181×49×5) (3) The 5-D matrix of hmF2 was reshaped into two-dimensional (2-D) matrixes by separating variables into two groups in terms of (glat, LT, and glon) and (month, F10.7). Then we performed EOF analysis to the 2-D matrix of hmF2 as equation (8). Where E k is the kth order EOFs representing the variations of glat, LT, and glon; A k is the corresponding coefficient of E k , and it denotes the variations of the season (month) and solar activity. The order (1-k) of EOF analysis are sorted by their relative contributions to the total variances of the 2-D matrix. A typical method to obtain E k and A k is to calculate the eigenvalues and eigenvectors of the covariance matrix of the 2-D matrix. Further details about EOF analysis can be found in the tutorial by Smith (2002) (4) A k was further decomposed through the three-order Fourier regression specified in equations (9-10) based on the laws of the prevalent seasonal and solar activity variations (e.g., Liu et al., 2009). b km , c km , e km , and f km are the Fourier expansion coefficients obtained from the least squares fit. A k in any month within F10.7 in the intervals of 70-150 sfu could be reconstructed by the Fourier harmonics and its coefficients. Further, a global empirical model of hmF2 was built in terms of E k combined with the reconstructed A k from the three-order Fourier regression The algorithms adopted to hmF2 were also applied to the other four parameters of the α-Chapman function. Finally, the Ne profile model was built by taking the five parameters into the α-Chapman function. We named it as α-Chapman Based Electron Density Profile Model (α-Chapman-Based-EDP). The model can LI ET AL.
10.1029/2020SW002642 6 of 23 give six-dimensional (month, glon, glat, LT, F10.7, and altitude) Ne of F region ionosphere. Among the six dimensions, only the first four are included in the F2COSMIC (He, Liu, Wan, & Wei, 2011). Figure 4 presents the relative contribution of the first ten EOFs to the total variance of each parameter. The chocolate, orange, purple, green, and blue colors represent NmF2, hmF2, Hm, a 1 , and a 2 , respectively. We adopted the first four orders to build models, and the total percentages of the five parameters (shown in the brackets of Figure 4) are larger than 98%.

Results of EOF Analysis
Figures 5-9 presents the E k of NmF2, hmF2, Hm, a 1 , and a 2 obtained from EOF analysis, respectively. They describe the different spatial variation modes at different local times. Figure 10 displays the (A k ) coefficients of different E k , which manifests the variations of season and solar activity. All E k and A k were normalized by its maxima to show the spatiotemporal variations better.
As shown in the top panels in Figures 5-9, the first EOF component (namely A 1 × E 1 ) generally denotes the ionospheric F region's primary spatiotemporal variations. E 1 in Figures 5-7 shows the general latitude and longitude patterns at different local times, which is largely consistent with previous climate studies (e.g., Liu et al., 2009;Li et al., 2019;Le et al., 2017;He et al., 2009). For example, hmF2 displays peaks at midlatitudes during nighttime because of the lifting by neutral winds; the daytime Hm reaches a huge peak in the equatorial region with two troughs at low latitudes due to the electrodynamics (Liu et al., 2008;. Many typical phenomena are clearly shown in E 1 . The typical daytime EIA (Appleton, 1946), Weddell Sea Anomaly , and nighttime midlatitude trough (Yang et al., 2018;He, Liu, Wan, & Zhao, 2011) could be clearly seen in the E 1 of NmF2 (top panel of Figure 5). The midlatitude longitudinal variations associated with neutral winds dynamics are also clearly shown in E 1 of NmF2, hmF2, Hm and a 1 (Chen et al., 2016b;He et al., 2009). The low latitude wavenumber 3 and 4 longitudinal structures also could be found in hmF2, Hm, and hmF2, especially during the daytime. The first column of Figure 10 shows the amplitude A 1 as a function of F10.7 and month. A 1 increases with F10.7 for NmF2, hmF2, and Hm, which manifests the positive solar activity dependence of Ne in the F region. A 1 decreases with F10.7 for a 1 and a 2 . It should be noted that E 1 of a 2 is negative, so the anticorrelation between A 1 of a 2 and F10.7 manifests the positive solar activity dependence of the altitude change rate of scale height in the bottom F region. At the same time, the A 1 of a 1 manifests that the altitude change rate of scale height in the topside ionosphere decreases with solar activity. The seasonal variation modes of A 1 (the first column of Figure 10  a 2 display an annual variation with a weak semi-annual component at 70 sfu, and the semi-annual variation becomes a distinct seasonal variation during high solar activity; Hm displays a weak seasonal variation. For NmF2, hmF2, and Hm, the second EOF component mainly manifests the hemispheric/seasonal asymmetry during solstice seasons. We take NmF2 as an example to explain it. The E 2 of NmF2 displays opposite signs in the northern and southern hemispheres, and the corresponding A 2 is positive (negative) during June (December) solstice. Multiplying E 2 with A 2 , the second EOF component increases (decreases) NmF2 in the northern (southern) hemisphere during April-August, while the situation is reversal in other months. In addition, the second EOF component of a 1 denotes the latitudinal variation modulated by a complex seasonal variation.
For NmF2, hmF2, and Hm, the third EOF component is essential in describing equatorial nighttime enhancement and contributes to the midlatitude longitudinal variation. The equatorial E 3 is positive during nighttime, and the corresponding A 3 is dominant by a semi-annual variation and gradually becomes LI ET AL.
10.1029/2020SW002642 8 of 23 positive from negative with increasing solar activity. Multiplying E 2 with A 2 , their negative effect becomes weak from low to high solar activity, so the equatorial third EOF component reflects the equatorial nighttime enhancement. In addition, the third EOF component of a 1 mainly modulates the latitude variation with an annual variation.
The fourth EOF component presents complex latitude and longitude variation with a regular semi-annual variation.

Model Error
The 4,487,884 pieces of Ne profiles from COSMIC-1 after denoising (but not griding) were compared with the five α-Chapman parameters' model. We discarded the Ne profile, if either of NmF2, hmF2, Hm, a 1 , LI ET AL.
10.1029/2020SW002642 9 of 23 and a 2 is out of three standard deviations from the mean, these profiles account about 3.6%. It is worth noting that the referenced observations of model's evaluation include the geomagnetic disturbed periods. Figure 11a displays the histogram of the relative errors between observations and model. The solid red line gives the results from the normal fitting, the μ and σ denote the mean and variance of the relative error, respectively. The relative errors mainly distribute within 40%, 20%, 75%, and 50% for Hm, hmF2, NmF2, and a 1 . The relative errors of a 2 have a wide distribution. Figure 11b presents the Bivariate histogram comparing the models and the observations. The details of the comparison are measured in root-mean-square-errors (RMSE), relative errors, and R in Formulas S1-S3. The black solid line shows that the values of the model are equal to the observation. More than 90% (marked at the top of each panel) data distributes in bins with samples higher than the limits of 100, and these scatters mainly tightly spread near the two sides of the black lines. The RMSE are 9.33 km, 22.69 km, 1.64 × 10 5 cm -3 , 0.04, and 0.2 for hmF2, Hm, NmF2, a 1 , and a 2 , respectively. The corresponding R is 0.75, 0.87, 0.92, 0.56, and 0.51, respectively. R 2 describes the ratio between the variances of model and observation, therefore measuring the goodness of the model. Our model LI ET AL.

10.1029/2020SW002642
10 of 23 can describe 57%, 75%, 84%, 31%, and 26% of variances of Hm, hmF2, NmF2, a 1 , and a 2 . It is the first time the models of scale height and its altitude gradient are measured with R 2 .
Besides the internal validation, we also check the predictive performance of our model and IRI2016 based on COSMIC-2 Ne profiles. We collected the parameters of Ne covering 180-850 km from COSMIC-2 (samples are 318,688,592 in total), including year, day of year or month, F10.7, latitude, longitude, and height. The predicted Ne was obtained by bringing these parameters into α-Chapman-Based-EDP and IRI2016.
Then we compared the predicted and observed Ne. Figure 12 shows overestimate the COSMIC-2 Ne at different extents. The distribution of relative errors has a long tail toward positive, considering α-Chapman-Based-EDP and IRI2016 both describe the climatological scale. Thus, the highly asymmetric distribution of relative errors indicates the complex and changeable day-to-day variability of Ne in ionospheric F region. There still a long way to well understand/capture/predict the day-to-day variability. The R in (a1) displays narrower distributions than (b1); the μ, σ, and RMSE from our model are less than IRI2016. View from R, our model (IRI2016) can predict 80% (71%) of Ne variations. These results, indicates that our model has less errors and perform better in prediction than IRI2016. Previous studies showed that good consistencies of Ne between COSMIC and ISR and ionosondes (Hu et al., 2014;Kumar et al., 2016;Lei et al., 2007), so we did not validate our model against the ISR and ionosondes.
View from Ne profiles of COSMIC in Figure 1, it is easy to find that the profiles are more during the ascending phase of the 23rd solar cycle. Thus, our model may perform better during the solar rising phase LI ET AL.

10.1029/2020SW002642
12 of 23 Figure 9. Normalized first four base functions (E 1 , E 2 , E 3 , and E 4 ) of a 2 (scale height' altitude change rates below and equals the peak height of F2 layer) as a function of latitude and longitude at different local times. The dotted lines denote the zero contours. than declining phase. Earth's ionosphere displays evident hysteresis phenomenon. That is, the same solar input produces different ionization in the rising and declining phases of solar cycle (Huang, 1963;Huang et al., 2020;Kane, 1992). How to effectively reproduce the ionospheric hysteresis is still a challenge problem in models . Due to the data limitation, we did not consider the hysteresis effect in current study. As a continuation of this study, we plan to improve our model under different solar phases in future.

Validity at Climatological Description
The validity in the physical description is also important to check the performance of a model, and it is actually more important than the mathematical errors. First, we validated the models of the five parameters of the α-Chapman profile. Panels of Figure 13 from top to bottom display the global variations of NmF2, hmF2, Hm, a 1 , and a 2 obtained from the α-Chapman-Based-EDP in December with F10.7 of 70 sfu. Our model well captures many typical ionospheric phenomena. NmF2 shows (1) latitude structure: EIA, nighttime midlatitude trough around the subauroral region in the northern hemisphere (Yang et al., 2018), and midlatitude band structure (Zhong et al., 2019); (2) longitudinal structure: wave-like longitudinal structure at midlatitudes with mag. lat. within 45°-60° (Q. L. Wang et al., 2016;S.-R. Zhang et al., 2012;Zhao et al., 2013), nighttime enhancements, including the ones during dusk-to-midnight (Chen et al., 2016b) and near midnight ( Q. H. . Hm and hmF2 display an evident hemispheric asymmetry, and longitudinal variations in equatorial and midlatitude regions (Li et al., 2019). It is worth noting that a 1 nearly displays an opposite spatial pattern with Hm. Prol et al. (2018) disclosed that the opposite responses between a 1 and Hm result from low latitude electrodynamics. Further studies should be taken to disclose these opposite responses at midlatitudes. Generally speaking, our model could well capture the typical temporal-spatial variations of NmF2, hmF2, Hm, and a 1 .
LI ET AL.

10.1029/2020SW002642
13 of 23  How about the validity of the α-Chapman-Based-EDP for the Ne in F region? EIA nearly is the most important phenomenon at low latitudes; it presents two Ne peaks at about ±15° mag. latitude and a Ne trough near mag. equator in F region, the two peaks merge into a single peak at high altitudes and display interhemispheric asymmetries during solstices (e.g., Lin et al., 2007;Tulasi et al., 2009). A stronger crest first presents in the winter hemisphere and would transit to the summer hemisphere at noon/afternoon local times during solstices. Huang et al. (2018)     transitions appear at 12 and 15 LTs in June and December, respectively. The transition time is earlier at 95°W (0°E) in June (December). The double-peak merges into a single peak lower than 400 km altitude. The representation of EIA shows good consistency with previous studies Lin et al., 2007;Tulasi et al., 2009), indicating that our model could well capture the ionospheric variation in F region at low latitudes. Figure 15 displays the corresponding results from IRI2016. IRI could not well capture the merging signature of EIA, the double-peak structure does still not merge into a single peak at 700 km. For 0°E, the transitions happen at 12 (17)  observations has been adopted in IRI2016, while the ionospheric Ne profile information, especially the topside ones, are still not included in IRI2016. The truthful representation of EIA is still a challenge for IRI.
The midlatitude trough is strongest in local winter (He, Liu, Wan, & Zhao, 2011;Lee et al., 2011). Figure 16 shows the midlatitude trough in local winter by NmF2 calculated by (a) α-Chapman-Based-EDP and (b) IRI2016 as well as (c) COSMIC-1 observations. Figure 16 a1 and c1 (a2 and c2) show the results in the northern (southern) hemisphere in December (June), while Figure 16 b1 (b2) displays results in DOY of 360 (180). As shown in Figure 16c, midlatitude trough distributes widely at longitudes during the local winter, and its minimum aligns well with geomagnetic latitudes at ±60°. Figure 16a and 16c indicate that α-Chapman-Based-EDP model well represented the midlatitude trough including width, longitudinal pattern, and temporal variation (Lee et al., 2011;Yang et al., 2018). However, the midlatitude trough calculated by IRI2016 appears at higher latitudes and has wider latitude distribution than observations. Its longitudinal LI ET AL.

10.1029/2020SW002642
19 of 23 pattern also displays a significant difference from that obtained from the observations, especially for the southern hemisphere. What's more, Figures 16a and 16c indicate that our model well represents the longitudinal variation at midlatitudes (with mag. lat. lower than 60°) (Q. H. Wang et al., 2016;S.-R. Zhang et al., 2013;Zhao et al., 2013), while IRI2016 shows different longitudinal patterns compared with observations. Figure S2 shows the corresponding relative error of Figure 16. Our model represents NmF2 with less error (<50%) than IRI2016 in the midlatitude trough. These comparisons disclosed that our model performs better than IRI2016 in capturing the trough and longitudinal variation at midlatitudes.
The good representations of EIA and midlatitude trough indicates that the α-Chapman-Based-EDP could well capture the typical ionospheric climatology of Ne in the F region, and its performance is better than IRI in the ionospheric F region.
(1) Our model organizes the 12-years Ne profiles of COSMIC more efficiently; the gridded Ne can be used in ionospheric climate research, such as seasonal variation, solar activity dependence, spatial structures. And the gridded data also can act as a good background for ionospheric assimilation in F region.
(2) The Model can be a substitute of IRI in the ionospheric F region and could help to improve the validity of IRI in the ionospheric F region. For example, the Hm obtained from our model could help to improve the performance of IRI in the topside ionosphere.
(3) The model can provide a useful reference of Ne in the F region for radio communication. The Earth's ionosphere is the transportation channel for radio wave, the characteristics of Ne in F region is important for engineering applications related to a radio wave. For example, the NmF2 limits the maximum usable frequency for terrestrial signal propagation on Earth and limits the lowest usable frequency for a signal transmitted to Earth's outer space (Hoque & Jakowski, 2011).

Summary
In this study, using ∼4.5 million Ne profiles from COSMIC, a high precision global Ne profile model in the F region ionosphere was built by EOF analysis and Fourier regression analysis of its coefficients to the five parameters of the α-Chapman function. The model has the following advantages: (1) The model can be a substitute of IRI in the ionospheric F region. It well captures the typical Ne climate phenomena in the ionospheric F region, its performance is better than IRI. The model could well describe EIA (the merging of the double-peak signature into a single peak, the asymmetry between the two crests and their transition), trough and longitudinal variation at midlatitudes (2) Our model displays less error than IRI2016, providing an alternative ionospheric Ne model in the F region. The model can give six-dimensional (latitude, longitude, altitude, local time, month, and solar activity) Ne around the world covering 180-850 km altitude and 70-150 sfu F10.7 ranges, as well as the key parameters of Ne profile including NmF2, hmF2, Hm, and its growth rates at height. Our model describes 57%, 75%, 84%, 31%, 26%, and 80% of variability of hmF2, Hm, NmF2, a 1 , a 2 , and Ne (3) The model has wide applications, and it could be used for climate studies, ionospheric assimilation, as well as some related engineering applications in radio communication. A dense matrix data calculated by the model will be released in https://www.researchgate.net/profile/Qiaoling_Li5 with the permissions of COSMIC organizations