Available online at www.sciencedirect.com Advances in Space Research 43 (2009) 1825–1834 www.elsevier.com/locate/asr Proposal of new models of the bottom-side B0 and B1 parameters for IRI D. Altadill a,b, J.M. Torta a, E. Blanch a,* b a Observatorio del Ebro, Universitat Ramon Llull – CSIC, Crta. Observatori No. 8, E43520 Roquetes, Spain Center for Atmospheric Research, University of Massachusetts Lowell, 600 Suffolk Str. 3rd Floor, Lowell, MA 01854, USA Received 30 October 2007; received in revised form 18 August 2008; accepted 18 August 2008 Abstract The time series of hourly electron density profiles N(h) obtained from 27 ionosonde stations distributed world-wide have been used to obtain N(h) average profiles on a monthly basis and to extract the expected bottom-side parameters that define the IRI profile under quiet conditions. The time series embrace the time interval from 1998 to 2006, which practically contains the entire solar cycle 23. The Spherical Harmonic Analysis (SHA) has been used as an analytical technique for modeling globally the B0 and B1 parameters as general functions on a spherical surface. Due to the irregular longitudinal distribution of the stations over the globe, it has been assumed that the ionosphere remains approximately constant in form for a given day under quiet conditions for a particular coordinate system. Since the Earth rotates under a Sun-fixed system, the time differences have been considered to be equivalent to longitude differences. The time dependence has been represented by a two-degree Fourier expansion to model the annual and semiannual variations and the year-by-year analyses of the B0 and B1 have furnished nine sets of spherical harmonic coefficients for each parameter. The spatial– temporal yearly coefficients have been further expressed as linear functions of Rz12 to model the solar cycle dependence. The resultant analytical model provides a tool to predict B0 and B1 at any location distributed among the used range of latitudes (70°N–50°S) and at any time that improves the fit to the observed data with respect to IRI prediction. Ó 2008 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Bottom-side Ionosphere; IRI model; Ionospheric thickness and shape parameters; Spherical Harmonic Analysis; Solar activity dependence 1. Introduction The Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI) promoted the International Reference Ionosphere model (IRI) to produce an empirical standard model of the ionosphere. The IRI model provides averages of the ionospheric parameters for given location and time (Bilitza, 2001). IRI describes the electron density profile N(h) of the bottom-side of the F2 layer as follows (Bilitza, 1990): N ðhÞ ¼ NmF 2 * expðX B1 Þ ; coshðX Þ X ¼ ðhmF 2  hÞ : B0 Corresponding author. E-mail address: eblanch@obsebre.es (E. Blanch). ð1Þ NmF2 in Eq. (1) means the electron density at the F2 peak, hmF2 means the F2 peak height and B0 and B1 are the parameters that brought about the thickness and shape of the bottom-side of the F2 layer, respectively (Bilitza et al., 2000). B0 equals to the difference between hmF2 and the height where the electron density equals to 0.24 times NmF2 in absence of the F1 layer or the F1 peak height (hmF1) if the latter occurs. B1 describes the shape of the profile between the two heights from which the B0 is estimated (e.g., Reinisch and Huang, 1998; Bilitza, 1998). The current IRI version has two options to model the parameters B0 and B1, the standard and Gulyaeva options. The standard option for B0 is based on an empirical table of B0 values deduced from the measured bottom-side profiles as result of the analysis of a large amount of ionosonde data (Bilitza et al., 2000). This table contains the solar 0273-1177/$36.00 Ó 2008 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2008.08.014 1826 D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 cycle, seasonal, diurnal and geographical variations of B0. The standard option for B1 represents the diurnal variation by the following Epstein transition function (e.g., Gulyaeva, 2007)   1 1  ; ð2Þ B1 ¼ 2:6  0:7 1 þ expðtr  tÞ 1 þ expðts  tÞ where t means local time, tr and ts mean sunrise and sunset times, respectively, and all times are expressed in hours. The Gulyaeva’s option for B0 is an analytical function in terms of hmF2 and the bottom-side semi-thickness (Gulyaeva, 1987). The Gulyaeva’s option for B1 uses a constant value B1 = 3. According to Eq. (1), B0 and B1 are two key parameters to specify the bottom-side profile. The modeled values of foF2 and hmF2 by IRI agree reasonably well with the average values and measured values for quiet periods (e.g., Fuller-Rowell et al., 2000; Blanch et al., 2007). Nonetheless, IRI predictions for B0 and B1 show significant discrepancies with the observed and average values (e.g., Sethi and Mahajan, 2002; Lei et al., 2004; Blanch et al., 2007; Zhang et al., 2008). Recent investigations showed that a Local Model (LM) of B0 and B1 at midlatitudes improves significantly the IRI prediction (Blanch et al., 2007; Altadill et al., 2008). These LMs have an analytical formulation in terms of Fourier harmonics whose coefficients depend on a single parameter, the solar activity. The LMs reproduce the diurnal variations under quiet conditions of the B0 and B1 in terms of season (month) and solar activity (Sunspot number). The natural step forward is to enlarge the usefulness of the above formulation for a wider region, better if global, and to asses whether such formulation can improve the global IRI prediction for B0 and B1. The IRI model is being improved and updated continuously with the results of the annual workshops and specific sessions at general meetings (e.g., Bilitza, 2003) and one of the long-term goals of IRI is ‘‘to replace the current tabular form of the bottom-side parameters of the IRI model with appropriate mathematical functions” (see the IRI 2006 workshop Report by Bilitza at http://modelweb.gsfc. nasa.gov/ionos/iri/iri_06_report.html). In order to achieve the above goal, the aim of this paper is to present a new tool to predict B0 and B1 at any location distributed among the latitude range from 70°N to 50°S and at any time. This tool is an analytical model for each parameter whose spherical surface coefficients are expressed in terms of season (month) and solar activity (Sunspot number). The following sections describe the data base used to build the models, the data analysis technique, the model obtained, and the improvements of the model with respect to the IRI. The paper ends with a discussion section along with the envisaged new developments. 2. Data and modeling approach Vertical incidence (VI) ionograms recorded by DGS or DPS systems were used for this study. They are available in the Digital Ionospheric Data Base (DIDB) of the Center for Atmospheric Research (CAR) of the University of Massachusetts, Lowell (Reinisch et al., 2004). Data covered the time interval from 1998 to 2006, practically the entire solar cycle 23. The stations of the DIDB having the largest time of records and seeking for the best geographical distribution possible were selected. The digisonde stations used in this study are listed in the Table 1. Fig. 1 depicts the time span of the data and the geographical distribution of the stations listed in Table 1. The ionogram’s traces were inverted into ‘‘true” height electron density profiles N(h) with the True Height Profile Inversion Tool (NHPC) included on the Digisonde Ionogram Data Visualization/Editing Tool (SAO-X). This software is available at the website of the CAR (http://ulcar.uml.edu/), and a brief description is available at Reinisch et al. (2005) and references therein. In order to obtain the averages of the ionospheric parameters defining a quiet ionosphere to be modeled, the Monthly Averaged Representative Profile (MARP) technique (Huang and Reinisch, 1996) was applied. The MARPs were computed for a given station, month and hour with a percentage of exclusion of 25%, i.e., 25% of the individual N(h) profiles having the largest deviations compared to the average N(h) were excluded for computing the MARP. Therefore, the ‘‘extreme” N(h) profiles, which are most likely related to disturbed ionospheric conditions, were not used to obtain the average profiles and the MARPs were accepted to be the expected N(h) profiles for quiet ionospheric conditions (see Huang and Reinisch (1996) for details). The experimental values of the B0 and B1 were obtained from the best fit of Eq. (1) to the computed average N(h) profiles (Reinisch and Huang, 1998). The IRI-2001 model (ftp://nssdcftp.gsfc.nasa.gov/models/ionospheric/iri/iri2001/) was also used to obtain the B0 and B1 parameters for the same locations and time intervals from which experimental data was gotten to compare whether the new models predict the experimental values better than IRI does. Both options of IRI were considered for computing B0, standard and Gulyaeva, and only the standard option for computing B1. The foF2 storm model was turned off in IRI-2001 because this study deals for quiet ionospheric conditions. The Spherical Harmonic Analysis (SHA) was the analytical technique chosen for the modeling. The intrinsic ability to solve Laplace’s equation in spherical coordinates converted SHA into a powerful tool to analyze spherically distributed data from potential fields since the epoch of Gauss. SHA can also distinguish the fields from their external and internal sources making it especially useful for modeling the different constituents of the geomagnetic field variations (e.g., Sabaka et al., 2002). Nonetheless, SHA has been widely used in geophysics in general and in aeronomy in particular (e.g., Lazo et al., 2004) because it is capable to analytically model a general time-dependent function on a spherical surface as follows: D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 Table 1 Catalogue of ionospheric stations used in the analysis listed from high to low geographic latitude. The station codes, geographical and geomagnetic coordinates and the dip latitudes are also given for each station. Geomagnetic coordinates and dip latitudes are computed from IGRF (Macmillan and Maus, 2005) evaluated at 1998.5, and the latter at an altitude of 300 km. Station URSI code Qaanaaq Sondrestrom Goose Bay Chilton Ebro, Roquetes Wallops Is. Anyang El Arenosillo Osan AB Point Arguello Wuhan Eglin AFB Okinawa Hainan Ramey Kwajalein Sao Luis Fortaleza Ascension Is. Jicamarca Learmonth Madimbo Cachoeira Pau. Louisvale Grahamstown Bundoora Port Stanley THJ77 SMJ67 GSJ53 RL052 EB040 WP937 AN438 EA036 SN437 PA836 f ðh; u; T Þ ¼ Geog. lon. Geom. lat Geom. lon. 77.5 66.98 53.3 51.5 40.8 290.6 309.06 299.7 359.4 0.5 87.98 76.35 63.55 53.77 43.24 13.14 35.63 15.14 84.46 81.29 82.31 70.89 60.48 48.55 36.32 37.9 37.39 37.1 37.1 34.8 284.5 126.95 353.3 127.0 239.5 48.39 27.35 40.92 27.06 41.27 355.27 196.47 72.80 196.53 304.49 48.83 33.78 31.83 33.43 40.22 WU430 EG931 OK426 HA419 PRJ18 KJ609 SAA0K FZA0M AS00Q JI91J LM42B MU12K CAJ2M 30.5 30.4 26.7 19.4 18.5 9.0 2.6 3.8 7.95 12.0 21.8 22.39 23.2 114.4 273.2 128.2 109.0 292.9 167.2 315.8 322.0 345.6 283.2 114.1 30.88 314.2 20.03 40.5 16.75 8.88 28.99 3.45 6.75 4.98 2.17 1.52 32.27 24.25 13.66 185.43 342.61 198.35 180.50 4.8 237.74 27.50 33.58 56.34 354.84 186.18 97.92 24.24 26.74 41.87 21.09 13.36 27.98 3.78 0.48 5.05 20.42 0.6 36.08 39.20 16.99 LV12P GR13L BV53Q PSJ5J 28.5 33.3 37.7 51.6 21.2 26.5 145.05 302.1 28.47 34.13 45.83 41.32 86.99 91.0 222.58 11.23 45.10 44.97 52.13 30.8 N X n X Geog. lat. Dip lat. 1827 not available for the data base we got and a rather good distribution of them in latitude can only be guaranteed, we made use of a ‘‘zero order” approximation. We assumed that the ionosphere under quiet conditions remains approximately constant in form over a given day, and since the Earth rotates 360° under a Sun-fixed system, we considered the local time differences to be equivalent to longitude differences. An important issue for this approach is the choice of the coordinates defining the parallels that pass over the location of each original station. These parallels should define the geographical latitudinal position of the 24 fictitious stations, 15° in longitude apart from each other, arising from our ‘‘zero order” approximation. Geographic, geomagnetic, magnetic dip, modified dip and magnetic apex coordinates were tested for this study. Among the physical sources that cause diurnal variation of the ionospheric parameters, it is well known that the photo-ionization and heating are mainly determined by geographic latitude and the location of the Sun, but electric fields and neutral wind effects, playing significant role on P mn ðcos hÞfgmn ðT Þ cosðmuÞ n¼0 m¼0 þ hmn ðT Þ sinðmuÞg: ð3Þ h and u in Eq. (3) mean the geocentric spherical coordinates, colatitude and longitude, respectively, and T is time in any convenient unit. Pnm are the associated Legendre functions of the first kind, or simply Lengendre polynomials, of n and m integer degree and order, respectively. gnm and hnm are the spherical harmonic or Gauss coefficients that define the model, which can be functions of T for a time-varying solution. 3. Data analyses and results The inverse problem to infer the model given a set of observations is reasonably straightforward when the observations are available over the whole sphere, or when the observations are distributed on a solid network to produce an over-determined system of equations that can be solved by least squares. Since a global distribution of stations is Fig. 1. Top: Time span covered by the data base for each station used in the analysis. The stations are distributed from North to South geographic latitude. Black lines indicate the availability of the experimental data and grey lines indicate the time and stations with synthesized data as described in Section 3. Bottom: Geographical distribution of the stations under analysis. 1828 D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 the hmF2 variation, are sensitive to the direction of the main geomagnetic field (e.g., Rishbeth and Mendillo, 2001). So that it seems more natural to work with the symmetry provided by magnetic coordinates, either magnetic dip (which take into account the ‘‘true” magnetic inclination, or dip, of any location) or geomagnetic coordinates (given by the simplification provided by assuming a dipolar field only). Both, in any case, are in practice obtained from a model of the geomagnetic main field at a given epoch, such as IGRF (Macmillan and Maus, 2005). The modified dip latitudes (modip) introduced by Rawer (1963) are those used by IRI (see Bilitza, 2001). The modip latitudes come near the geomagnetic inclination at low latitudes and get closer to the geodetic latitude as latitude increases. Magnetic apex coordinates (Richmond, 1995) are aligned with the geomagnetic field lines. That is why they are useful for describing the horizontal ionospheric currents, equivalent currents and ground-level magnetic perturbations, which are strongly horizontally organized by the ambient magnetic field, and have been also applied to map electric fields and plasma-drift velocities at ionospheric heights. Although not shown here, different coordinate systems have been tested to define the parallels that pass over the location of the original stations to distribute the 24 hourly values in longitude, concluding that magnetic dip coordi- nates are among the best choices for B0 and that geographic coordinates are among the best choices for B1. The physical reasons that can explain why the local time variations of B0 and B1 accommodate best to the magnetic dip and geographic symmetry, respectively, are beyond the scope of this paper. The rest of coordinate systems have provided worst results from the point of view of the expected seasonal inter-hemispheric symmetries. Accordingly, magnetic dip and geographic coordinates were chosen for grouping the stations into different latitudinal belts for the inspection of the spatial–temporal behavior of the B0 and B1, respectively, that will be shown next. 3.1. Experimental time–space pattern The local time variation of B0 (Fig. 2) and of B1 (Fig. 3) has a distinct pattern for different stations located within different latitudinal ranges, indicating a systematic behavior. The plots in Fig. 2 depict B0 as function of local time at indicated latitudinal belts, seasons and solar activity. The same apply for the plots in Fig. 3 but for B1. Summer season means May–July at the Northern Hemisphere and November–January at the Southern Hemisphere. Winter season means November–January at the Northern Hemisphere and May–July at the Southern Hemisphere. The Fig. 2. Local time variation of B0 for different stations located within the indicated latitudinal belts and for indicated seasons and solar activity levels. From left to right: B0 during summer for high solar activity (HSA), B0 during winter for HSA, B0 during summer for low solar activity (LSA) and B0 during winter for LSA. Note that stations are grouped within dip latitudes. D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 1829 Fig. 3. As in Fig. 2 but for B1. Note that stations are grouped within geographical latitudes. grey dots of Figs. 2 and 3 correspond to the experimental values of B0 and B1, respectively. Thick lines represent the polynomial functions that best fit the given parameter as function of local time, showing the dominant daily pattern. The experimental behavior of B0 can be summarized as follows (Fig. 2). The diurnal variation dominates the local time pattern during summer for all magnetic dip latitude ranges and solar activity levels, noon values of B0 are larger than midnight values. The diurnal amplitude during summer decreases as latitude increases and it increases as solar activity increases. Moreover, summer B0 values for high solar activity are larger than for low solar activity. The local time pattern of B0 in winter is significantly different from that in summer; it is clearly latitude dependent. The winter local time pattern at low latitudes shows a clear diurnal variation for all solar activity levels, with larger noon values than midnight ones. In addition, the winter B0 values and the diurnal amplitude at low latitudes for high solar activity are larger than for low solar activity. The winter local time pattern at mid-low latitudes shows a clear semidiurnal variation for high solar activity, but it shows a weak diurnal variation for low solar activity. Moreover, the winter B0 values at mid-low latitudes for high solar activity are larger than for low solar activity. Finally, the winter local time pattern at mid-high latitudes shows a dominant semidiurnal variation, there being better developed for high solar activity than for low solar activity. The winter B0 values at mid-high latitudes for high solar activity are also larger than for low solar activity. Fig. 3 shows the analysis results for B1 and a first inspection indicates that B1 is clearly noisier that B0, especially for low solar activity levels. Moreover, B1 during winter is noisier than during summer. The diurnal variation dominates the local time pattern of B1 during summer for all latitude ranges and solar activity levels, noon values of B1 are smaller than midnight values. The diurnal amplitude during summer decreases as latitude increases and summer values of B1 for high solar activity are smaller than for low solar activity. There is no clear daily pattern of B1 during winter for low solar activity, B1 is very noisy and only a weak diurnal variation manifests at low latitudes. However, diurnal variation of B1 dominates the daily pattern at low geographic latitudes during winter at high solar activity, whereas semidiurnal variation becomes dominant for mid-low and mid-high latitudes. In addition, B1 values during winter tend to be larger than those during summer. The above experimental results are consistent with those reported in previous studies (e.g., Adeniyi and Radicella, 1998; Sethi and Mahajan, 2002; Lei et al., 2004; Blanch et al., 2007), indicating that the zero order approximation is reasonable. 1830 D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 Fig. 5. Cross section plots of B0 as function of geographical longitude and latitude for different seasons and solar activity levels. (a) and (b) July and January for high solar activity (Rz = 119.6), respectively; (c) and (d) July and January for low solar activity (Rz = 15.2), respectively. The plots show the results at 0 UT. Fig. 4. Results of the yearly model of 1998 obtained without synthesized data (grey lines) and with synthesized data (black lines). The results are compared with the experimental data for given stations (grey dots). 3.2. Models parameterization and results Slightly different parameterizations were chosen for the models of B0 and B1. The maximum degree of the expansion (N in Eq. (3)) was set to 6 for both B0 and B1. The latter is a reasonable compromise between the allowance for different behaviors at high, middle and equatorial latitudes (see Figs. 2 and 3), and the number of stations distributed in colatitude. The maximum order of the expansion for the longitudinal (local time in this case) dependence (m in Eq. (3)) was limited to 4 for B0, but only to 2 for B1. Dealing with hourly values does not justify the allowance for more than four diurnal harmonics. Moreover, the large noisy daily variations observed in the B1, especially for low solar activity (Fig. 3), suggested us limiting the representation of this parameter to just the mean, the diurnal and the semi- diurnal terms. The time dependence over the year was represented by a two-degree Fourier expansion to model the annual and semiannual variations, according to the experimental results (Figs. 2 and 3) and the results from previous investigations (e.g., Altadill et al., 2008 and references therein). Therefore, the Gauss coefficients of Eq. (3) were expressed as ( )  m  X 2 gamn;q gn ðT Þ ¼ cosð2pqT =12Þ ð4Þ hamn;q hmn ðT Þ q¼0 ( m ) gbn;q sinð2pqT =12Þ: þ hbmn;q The year-by-year spatial–temporal analyses as described above furnished nine sets of spherical harmonic coefficients for each parameter B0 and B1. However, some of the yearly models suffered from the lack of data in some criti- D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 cally located stations, especially from those at the northernmost or southernmost latitudes that constrain the model behavior or from those located at regions with a distinct pattern and without any neighbor station. Fig. 4 shows how the lack of data leads to an unrealistic yearly model (e.g., SMJ67 and GR13L). In order to overcome this lack of data disturbing some of the yearly models, synthetic data were added at critical latitudes as follows. We isolated those years with enough data in those strategically located stations which allow obtaining realistic yearly models for both parameters. These years are 2002, 2004, 2005 and 2006, practically covering the whole range of solar activity levels. We noticed that the spatial–temporal yearly coefficients are solar activity dependent in a similar way as obtained for the local models of B0 and B1 at midlatitudes (Blanch et al., 2007; Altadill et al., 2008). Thus, the Gauss coefficients from these four years were expressed as linear functions of the solar activity index Rz12. Finally, these initial (or preliminary) models for B0 and B1 were used to create synthetic data for the following critical stations when experimental data were not available; SMJ67, PRJ18, JI91J, AS00Q and GR13L. Fig. 1 shows the data gaps corresponding to the aforementioned stations that were filled with synthetic data. In addition, we created a couple of fictitious stations in the Southern Hemisphere located at 60° and 70° dip magnetic South (SYN60, SYN70 in Fig. 1) for B0, and a single one at 70° geographic South for B1, to constrain the behavior of the model at high latitudes since there is a clear lack of stations there. We synthesized the data corresponding to these fictitious stations from those initial models at 60° and 70° dip magnetic North for B0, and at 70° geographic North for B1 (where they are well constrained from existing real data) but appropriately ‘‘mirrored” by shifting them by 6 months to accommodate seasonal differences. The last step, previous to build the global models, was to reanalyze the new data set containing real and synthetic data. These reanalyses brought the definitive nine sets of spherical harmonic coefficients for each parameter B0 and B1. Fig. 4 also shows that current yearly sets of spherical harmonic coefficients provide realistic yearly models. Once the definitive sets of Gauss coefficients was gotten by this simple iterative approach, the analytical functions determining the linear dependence of these coefficients with the solar activity were obtained. The yearly average of the Sunspot activity number Rz12 was again selected as proxy of the solar activity for the modeling purposes. According to the different parameterization described above, the spherical harmonic (SH) model for B0 is fixed by 430 coefficients (solar activity dependent) and the SH model for B1 is fixed by 230. Figs. 5 and 6 depict examples of the model results for B0 and B1, respectively, at different seasons and solar activity levels. The contour plots of Fig. 5 (Fig. 6) show B0 (B1) as function of latitude and longitude. Despite the models allow obtaining the parameters into whole latitude range, the plots of Figs. 5 and 6 were restricted to the latitude range from 70°N to 50°S. This geographical 1831 Fig. 6. As Fig. 5 but for B1. range is the one where data was gotten to feed the model and where the model is believed to be valid. A detailed analysis comparing the model results (Figs. 5 and 6) with the experimental results (Figs. 2 and 3) shows that the empirical model reproduces the dominant patterns noticed on Section 3. 4. Models goodness and validity Direct comparisons between experimental data and modeled values indicate that the spherical harmonic (SH) models improve the IRI modeling. Fig. 7 shows an example of the above comparisons for B0 parameter using both IRI options, standard and Gulyaeva, at a particular solar activity level. The SH model reproduces the diurnal and semidiurnal variations as well as the seasonal variations noticed in Section 3 for all latitudes better than IRI does. The standard option of IRI underestimates the diurnal and seasonal variations of B0 for high latitude stations (e.g., SMJ67), it underestimates the diurnal variation for midlatitude stations (e.g., EB040, GR13L) and it provides good results for equatorial latitudes (e.g., JI91J). The Gulyaeva option 1832 D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 Fig. 7. Comparison of the experimental values of B0 for a particular year (grey dots) with model results from SHA (dashed lines), with the IRI standard option (grey lines) and with the IRI Gulyaeva option (black lines). of IRI gives good results for midlatitude stations, but it fails for high and equatorial latitude stations (e.g., SMJ67 and JI91J). The latter agrees with previous results; the Gulyaeva option of IRI behaves better than the standard one for B0 at midlatitudes while the standard option gets better results at equatorial latitudes (Bilitza, 2001). Although not shown here direct comparisons as described above confirm that the results of the SH model fits the representation of both B0 and B1 better than IRI does for all latitudes and solar activity levels. The root mean square errors (RMSE) of the models were compared to asses the goodness and validity of the SH models (Fig. 8). These comparisons were done on a yearly basis (i.e., for different solar activity levels) and on a global basis taking into account both IRI options for B0 and the standard option for B1. The results show that SH models give less error than current IRI for both B0 and B1 and for all solar activity levels. Whereas the RMSE of the SH model over the nine years for B0 is 17.77, the RMSE of the standard option is 24.43 and that of the Gulyaeva option is 25.84. Therefore, SH prediction of B0 gives 37.5% less RMSE than the standard option and 45.4% less RMSE than Gulyaeva option. These results indicate a significant improvement of the SH model compared with current IRI for B0. As refers to B1, the RMSE of the SH model over the nine years is 0.67 and that of the standard option is 0.80. Thus, SH prediction of B1 provides 19.4% less RMSE than the current IRI. The improvement of the prediction of the SH model for B1 with respect to the IRI is less than the improvement for B0 but still significant. This fact can be result of the noisier pattern of B1 compared to the clear pattern of B0 (Figs. 2 and 3). 5. Summary and discussion Many papers have shown significant disagreements between the modeling of B0 and B1 by the IRI and their observed behavior. Consequently, the IRI community is decided to replace the current form of B0 and B1 with appropriate analytical functions. This work has broached the above goal aiming at improving the current IRI prediction of B0 and B1 under quiet ionospheric conditions developing new models for both parameters. D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 1833 Although the SH models can represent B0 and B1 at global scale, the validity of the model should be limited between 70°N and 50°S, where empirical data were available. Another weakness of the SH models arises from the ‘‘zero order” approximation. This could hide some longitudinal effects on the variation of the parameters. The above limitations are the result of restricted availability of the empirical data to build the models. Therefore, further works for developing better models would seek for better data coverage in both time and space. Efforts obtaining data proxies of B0 and B1 from classical ionosondes (e.g., Zhang et al., 2008) and of satellite data (e.g., Gulyaeva, 2007) would help on this task. Acknowledgements Fig. 8. Comparisons of the root mean square errors (RMSE) obtained with the models involved in this study. The plot at the top shows the solar activity levels corresponding to each year, the middle plot shows the RSME for B0 and the bottom plot shows the RSME for B1. The electron density profiles measured over 27 stations distributed world-wide for the time interval 1998–2006 have been analyzed to obtain the average profiles representing quiet ionospheric conditions. The experimental values of B0 and B1 have been deduced from the average profiles, serving for obtaining the time–space pattern of B0 and B1 and for modeling parameterization. The local time behavior of the parameters has been considered as longitude variation to avoid the poor density of measurements over the globe. Moreover, preliminary models for both parameters have been used to synthesize data at critical latitudes when data gaps appear and to constrain the behavior of the models at high latitudes of the Southern Hemisphere. The Spherical Harmonic Analysis has been used for modeling purposes, allowing the analytical modeling of B0 and B1 by time-dependent functions on a spherical surface. The SH models have been parameterized according to the time–space pattern of both parameters. The analyses have provided nine sets of spherical harmonic coefficients for each parameter B0 and B1. Finally, the coefficients corresponding to each parameter have been further expressed as linear functions of the solar activity providing two models, one for B0 and another for B1. The yearly average of Rz12 has been selected as proxy of the solar activity. The SH model for B0 has been fixed by 430 coefficients and the SH model for B1 by 230. The two empirical models proposed in this investigation have been tested against the IRI resulting that SH models improve the prediction of B0 by 40% and that of B1 by 20%. Therefore, it is worthwhile considering the SH models as potential options into further IRI versions. This research has been supported by Spanish projects CGL2006-12437-C02-02/ANT of MEC, and 2006BE00112 of AGAUR, and also by USAF Grant FA8718-L-0072 of the AF Research Laboratory. The authors wish to express their gratitude to the DIDB team and to the ionospheric stations contributing to it for making the data available, to the IRI team for making the IRI model available and to Arthur Richmond for providing the code to compute the magnetic apex coordinates. References Adeniyi, J.O., Radicella, S.M. Variation of bottomside profile parameters B0 and B1 at high solar activity for an equatorial station. J. Atmos. Sol. Terr. Phys. 60, 1123–1127, 1998. Altadill, D., Arrazola, D., Blanch, E., Buresova, D. Solar activity variations of ionosonde measurements and modeling results. Adv. Space Res. 42, 610–616, 2008. Bilitza, D. The International Reference Ionosphere 1990. National Space Science Data Center, NSSDC/WDC-A-R&S Reports 90-22, Greenbelt, Maryland, November 1990. Bilitza, D. Improving the standard IRIB0 model, in: Radicella, S.M. (Ed.), Proceedings of the IRI Task Force Activity 1997, International Center for Theoretical Physics. Report IC/IR/98/9, Trieste, Italy, pp. 6–14, 1998. Bilitza, D. International Reference Ionosphere 2000. Radio Sci. 36 (2), 261–275, 2001. Bilitza, D. International Reference Ionosphere 2000 – examples of improvements and new features. Adv. Space Res. 31 (3), 757–767, 2003. Bilitza, D., Radicella, S., Reinisch, B.W., Adeniyi, I.O., Gonzales, M.M., Zhang, S.R., Obrou, O. New B0 and B1 models for IRI. Adv. Space Res. 25 (1), 89–96, 2000. Blanch, E., Arrazola, D., Altadill, D., Buresova, D., Mosert, M. Improvement of IRI B0, B1 and D1 at mid-latitudes using MARP. Adv. Space Res. 39 (5), 701–710, doi:10.1016/j.asr.2006.08.007, 2007. Fuller-Rowell, T.J., Codrescu, M.C., Wilkinson, P. Quantitative modeling of the ionospheric response to geomagnetic activity. Ann. Geophys. 18, 766–781, 2000. Gulyaeva, T.L. Progress in ionospheric informatics based on electron density profile analysis of ionograms. Adv. Space Res. 7 (6), 39–48, 1987. Gulyaeva, T.L. Variable coupling between the bottomside and topside thickness of the ionosphere. J. Atmos. Sol. Terr. Phys. 69, 528–536, doi:10.1016/j.jastp.2006.10.015, 2007. Huang, X., Reinisch, B.W. Vertical electron density profiles from digisonde ionograms. The average representative profile. Ann. Geofis. 39 (4), 751–756, 1996. 1834 D. Altadill et al. / Advances in Space Research 43 (2009) 1825–1834 Lazo, B., Calzadilla, A., Alazo, K., Rodrı́guez, M., González, J.S. Regional mapping of F2 peak plasma frequency by spherical harmonic expansion. Adv. Space Res. 33, 880–883, 2004. Lei, J., Liu, L., Wan, W., Zhang, S.R., Holt, J.M. A statistical study of ionospheric profile parameters derived from Millstone Hill incoherent scatter radar measurements. Geophys. Res. Lett. 31, L14804, doi:10.1029/2004GL020578, 2004. Macmillan, S., Maus, S. International Geomagnetic Reference Field – the tenth generation. Earth Planets Space 57, 1135–1140, 2005. Rawer, K. Propagation of decameter waves (HF Band), in: Landmark, B. (Ed.), Meteorological and Astronomical Influences on Radio Wave Propagation. Academic Press, New York, pp. 221– 250, 1963. Reinisch, B.W., Huang, X. Finding better B0 and B1 parameters for the IRI F2-profile function. Adv. Space Res. 22 (6), 741–747, 1998. Reinisch, B.W., Galkin, I.A., Khmyrov, G., Kozlov, A., Kitrosser, D.F. Automated collection and dissemination of ionospheric data from the digisonde network. Adv. Radio Sci. 2, 241–247, 2004. Reinisch, B.W., Huang, X., Galkin, I.A., Paznukhov, V., Kozlov, A. Recent advances in real-time analysis of ionograms and ionospheric drift measurements with digisondes. J. Atmos. Sol. Terr. Phys. 67, 1054–1062, doi:10.1016/j.jastp.2005.01.009, 2005. Richmond, A.D. Ionospheric electrodynamics using magnetic apex coordinates. J. Geomag. Geoelectr. 47, 191–212, 1995. Rishbeth, H., Mendillo, M. Patterns of the F2-layer variability. J. Atmos. Sol. Terr. Phys. 63 (15), 1661–1680, 2001. Sabaka, T.J., Olsen, N., Langel, R.A. A comprehensive model of the quiet-time, near Earth magnetic field: phase 3. Geophys. J. Int. 151, 32–68, 2002. Sethi, N.K., Mahajan, K.K. The bottomside parameters B0, B1 obtained from incoherent scatter measurements during a solar maximum and their comparisons with the IRI-2001 model. Ann. Geophys. 20 (6), 817–822, 2002. Zhang, M.-L., Wan, W., Liu, L., Shi, J.K. Variability of the behavior of the bottomside (B0, B1) parameters obtained from the ground-based ionograms at China’s low latitude station. Adv. Space Res. 42, 695– 702, 2008.