COSMIC GPS Ionospheric Sensing and Space Weather

As our civilization becomes more dependent on space based technologies, we become more vulnerable to conditions in space weather. Accurate space weather specification and forecasting require proper modeling which account for the coupling between the sun, the magnetosphere, the thermo sphere, the ionosphere and the mesosphere. In spite of the tremendous advances that have been made in understanding the physics behind different space weather phenomena, the ability to specify or predict space weather is limited due to the lack of continuous and extensive observations in these regions. Placing a constellation of GPS receivers in low-Earth orbit, such as the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC), provides an extremely powerful system for continuously and extensively measuring one of these regions, the ionosphere. COSMIC, by use of GPS occultations, will make it possible to obtain continuous and global3-dimensional images of electron density, irregularities and TIDs in the ionosphere and plasmasphere. COSMIC would provide nearly 5600 globally distributed occultations per day suitable for ionospheric sensing. Occultations can be processed individually to obtain vertical profiles of electron density, with vertical resolution of ~lkm, or collectively by means of tomography or data assimilation to obtain 3-D images of electron density or irregularity structure. In this paper we describe the GPS observables for ionospheric sensing and the occultation geometry. Our presentation evolves from discussing simple to more complicated inversion techniques starting with the Abel inversion, gradient-constrained Abel inversion, tomography, and finally data assimilation. In each of these techniques, the accuracy is assessed either via examination of real data from GPSIMET or via simulation. Brief discussions of measuring ionospheric irregularities and TIDs are given.

via simulation. Brief discussions of measuring ionospheric irregularities and TIDs are given.

INTRODUCTION
As we enter the twenty-first century, we are witnessing revolutions in several areas of technology. Remote sensing of the Earth environment is among those areas, where continu ous and global sensing of the Earth's atmosphere is becoming a technological necessity for our daily life. For instance, remote sensing of the neutral atmosphere has tremendously influ enced the way Numerical Weather Prediction (NWP) is done today and our reliance on it. And as our dependence on technology in space grows, a similar revolution will be witnessed in the manner we sense and model space weather. The awareness of the potentially hazardous effects of space weather on technological systems motivated the creation of the National Space Weather Program (NSWP) in the United States and similar programs elsewhere. NSWP places particular emphasis on the need to provide timely, accurate, and reliable space environment observations, specifications, and forecasts. These requirements are similar to what is already accomplished with a great measure of success in the NWP community.
Accurate specification and forecasting of space weather phenomena is very difficult be cause it requires ( 1) accurate modeling of the coupling between the sun, the magnetosphere, the thermosphere, the ionosphere, and the mesosphere, (2) continuous; reliable and accurate observations of all of these regions, (3) the ability to assimilate the data into the models in an optimal and self-consistent manner. These are great challenges for the space weather commu nity and require tremendous efforts and advances in several areas of research and technology.
In the next few years we will witness great advances in the way the space environment is sensed as several satellites or con$tellations of satellites are launched with space sensing in struments of various types. One such mission is the Constellation Observing System for Me teorology, Ionosphere, and Climate (COSMIC), which will consist of 8 microsatellites in 8 evenly spaced orbital planes with circular orbits at -800km altitude and 72 deg. inclination. Each COSMIC satellite is equipped with a receiver that tracks all Global Positioning System (GPS) satellites in view including those rising and setting behind the Earth's atmosphere and ionosphere. Limb sounding of the Earth's ionosphere with GPS-LEO links was first demon strated with the GPS/MET experiment and more recently with Oersted and SUNSAT. The COSMIC mission with its 8 microsatellites will provide nearly 5600 ionospheric occultations per day when counting rising, setting and "grazing" occultations. This relatively new ap proach of sensing th. e Earth's ionosphere with a GPS receiver in LEO will potentially change the manner by which we observe the ionosphere. GPS observables can be used to measure integrated electron density (known as total electron content-TEC) along the transmitter-re ceiver line-of-sight and scintillation. With a dense set of measurements, as would be obtained from COSMIC, it will possible to obtain instantaneous global 3-dimensional (3-D) maps of electron densities, ionospheric irregularities and Traveling Ionospheric Disturbances (TIDs) at all times. This paper has two purposes. First, to provide an overview of ionospheric sensing from GPS space measurements. Second, to discuss means of making use of the extremely rich set of measurements that could become available from COSMIC. We will describe the different methods of processing ionospheric data collected during an occultation in a progressive man ner starting from the simplest (Abel inversion) to the most sophisticated (data assimilation).
We first review the GPS observables and the various ways in which TEC, with different quali ties of noise, can be derived. This is done in section 2. In section 3 we discuss the coverage of COSMIC occultations in the ionosphere from which we will start to get a glimpse that a con stellation like COSMIC would provide nothing less than a revolution in ionospheric sensing.
In section 4 we review the Abel transform as a means of obtaining electron density profiles from individual occultations and we compare results obtained from GPS/MET with incoher ent scatter radar and ionosonde measurements. In section 5, we discuss a possible improve ment to the Abel inversion scheme by use of the "constrained-gradient" inversion and exam ine its performance through simulation experiments. In section 6, a tomographic scheme is discussed as a means of combining nearby occultations to provide a purely data driven speci fication of electron density. In section 7, we describe a full-fledged assimilation scheme using a recursive Kalman filter. An Observation System Simulation Experiment (OSSE) is de scribed in that section to assess the value of our assimilation scheme and the impact that a system such as COSMIC might have for space weather. In section 8 a brief discussion is given of the impact of GPS ionospheric measurements, including scintillation and detection of TID, on space weather. We conclude in section 9.

Occultation Geometry
A LEO-GPS occultation geometry is depicted in Fig. 1. An occultation takes place when the GPS satellite sets or rises behind the Earth's ionosphere or lower neutral atmosphere as seen by a receiver in LEO. Each occultation therefore consists of a set of links with tangent points ranging from the LEO satellite height to the surface. One occultation takes nearly 4-10 minutes depending on the relative geometry of the LEO and GPS satellites. When the signal passes through the Earth media, it experiences delay, bending and amplitude change (e.g., defocusing or scintillation). These changes can be used to extract information such as profiles of electron density in the ionosphere, or temperature, pressure and water vapor in the lower neutral atmosphere (0-50 km), in the vicinity of the tangent point.

GPS Observables
The GPS constellation consists of 24 satellites at -26500 km radius, -12hr period, orbit ing in 6 different planes inclined at -55°. Each P/ the pseudorange between transmitter i to receiver j; k =1 or 2 for Ll or L2 frequencies, respectively; (1) n/ an integer which is constant over a connected arc; A.. k the operating wavelength; ek phase measurement noise. Eqs.
(1 and 2) ignore (1) transmitting and receiving antenna phase center variations and "wind up" effects (Wu et al., 1993) which can be calibrated or modeled if needed, and (2) higher order ionospheric terms (order l/f 3 or higher), which result from the expansion of the Appleton Hartree formula (see, e.g., Bassiri and Hajj, 1993). These terms can be ignored for iono spheric retrievals. A subscript on any term in Eqs.
(1 and 2) implies that it is frequency dependent. The neutral atmosphere is non-dispersive at radio frequencies; however, because the L1 and L2 "split" in the ionosphere, at a given time they sense slightly different parts of the neutral atmosphere (as depicted in Fig. 1). This split is the reason for having subscripts on of and TEC/ in Eqs. (1and 2).

Ionospheric Observables
There are several methods for extracting the ionospheric delay along the LEO-GPS line of-sight. The selection of a specific method depends on many factors which include the data noise, desired retrieval accuracy and resolution, and the inversion technique used in retrieving geophysical quantities such as electron density profiles. For completeness, we summarize the different methods below.
2.3.1 Relative TEC from L1 and L2 (TECr,p) By ignoring the split of the L1 and L2 signals, we can calculate the relative TEC (subr,p script r refers to relative, and p refers to precise) along the LEO-GPS line-of-sight by taking the linear combination 1 J;2 J;2 TECr,p =- where B is an undetermined bias which is constant for a connected arc, and v is the TEC noise. The constant multiplying L1 -L2 is equal to 9.52X10 1 6 in IS units. Therefore, each 1 meter of differential delay between L 1 and L2 corresponds roughly to 10 TECU (1 TECU = 10 1 6 e/m3). The TEC noise is related to the L 1 and L2 noise by The precision of TEC when Anti-Spoofing (AS) is off or with enhanced codeless GPS rer, p ceivers such as the ones to be carried on COSMIC is better than 0.01 TECU.
where C is the sum of the transmitter and receiver differential code biases (DCB's) , which experience little variation and are usually assumed constant. Much work has been done on solving for these constants (Mannucci et al., 1998), and they are usually known to within 1-2 TECU. The TEC 1-CY noise ( r;-) is determined by the P1 and P 2 noise in a relation similar to Eq. (4). Since range measurements tend to be about 100 times noisier than those of the phase, a standard procedure is to solve for Bin Eq.
(3) by a least-square minimization of I, (TECa,n TEC )2, where the sum is over all measurements in a connected arc. This way, we obtain an r,p absolute and precise TEC • . P with a bias which is of order of 1-3 TECU.

2.3.3
Relative TEC fr om Ll and P1 (TEC. , n) Thus far, it has been assumed that the GPS receiver tracks both Ll and L2 signals. In many applications, the GPS receiver can be a single frequency receiver tracking Ll only. In this case the TEC can be derived by taking where B' is an undetermined bias which is constant over a connected arc. We note that TEC is a relative TEC measurement with noise ( � ) which is dominated by the Pl noise. If the range and clock terms ( ou, C;, Cj) in Eq.
(2) are known, then the sum of the atmospheric and ionospheric delays contributing to a phase measurement ( oijd X TECii I f 2 ) can be determined up to a constant bias. The range term can be determined from Precise Orbit Determination (POD) of the LEO and GPS satellites. GPS POD is routinely available by use of a ground network of receivers. LEO satellite POD is possible by use of GPS satellites in view. In addition, the clock terms can be solved relative to an accurate clock on the ground by a process similar to "double differencing" where four links (the occulting link and the dashed links in Fig. 1) are differenced such that all clocks cancel out. The details of this technique as it applies to neutral atmospheric sensing are described in (Kursinski et al., 1997). When the tangent point of the occulted signal is above 80 km, the neutral atmospheric delay is less than 2 mm and can be ignored, therefore a relative and precise TEC can be obtained.
This method has several advantages over the one described in 2.3.1: (1) it requires a single frequency receiver 7 , (2) it is more precise since the TEC L noise ( crv = ( + 1 2 I d)crv = r,p, I }J I 6.16x1016 crv1 in IS units) is significantly smaller than that ofEq. (4), (3) it does not ignore the effect of bending. (Based on GPS/MET data, Hajj and Romans (1998) showed that bending in 7 Sufficiently accurate POD and clock solutions can still be obtained with a single frequency receiver [see e.g., Bertiger et al., 1995].
the ionosphere during solar minimum is< 0.01 deg. If this value reaches up to 0. 1 deg. during solar maximum, ignoring bending will result in systematically over-estimating the HmF2since bending above the F-peak is upward-by up to 5 km). These advantages have to be weighed against the extra complications and time required in ensuring that precise orbits as well as double differencing geometries are available.
Ll and L2 amplitude measurements In addition to phase and range measurements, the receiver measures the voltage signal-to noise ratio (VSNR) for both Ll and L2 from which amplitude scintillation can be detected. Scintillation measurements, from which ionospheric irregularity maps can be obtained, are of vital importance to space weather because of their potentially detrimental effects on satellite and ground communications. Amplitude measurements can also be combined with the phase measurements to improve occultation retrievals in various ways including detecting and cor recting for atmospheric multipath effects, accoun ' ting for diffraction for the purpose of im proving the horizontal resolution (see e.g., Karayel and Hinson, 1997), and possibly localizing electron density irregularities along the GPS-LEO ray path (Vorob'ev et al., 1999;Sokolovskiy, 2000).

COSMIC COVERAGE AND RESOLUTION
During one LEO satellite revolution, each GPS satellite will in general set or rise once behind the Earth's atmosphere. Each COSMIC satellite, with nearly 14.6 revolutions per day, will see 14.6 x 24 x 2;;:700 daily occultations. The total number of daily occultations seen by the 8 COSMIC satellites is then over 5600. An example of the distribution of such a coverage based on actual GPS and representative COSMIC orbits is shown in Fig. 2 which has a total of 5616 occultations. These occultations are divided nearly equally between rising, setting, and grazing geometries, which are defined based on the angle between the occultation plane and the LEO orbital plane being between 0-45 deg., 45-135 deg., and 135-180 deg, respectively. Representing an occultation by one point, as in Fig. 2, does not reveal the complex struc ture of occultation links in the ionosphere. Figure 3 shows all the occultation links (at the rate of one per minute) obtained during 98 minutes (the LEO satellite period) as a function of local time. Links obtained from one COSMIC satellite are highlighted in order to see their distribu tion relatiye to the LEO plane. Lines that are nearly parallel correspond to links from the same occultation. Setting and rising occultations appear as a dense set of lines (in some cases nearly on the top of each other) parallel to the orbital plane. Grazing occultations appear as parallel lines that move considerably (in some cases half way across the globe) during an occultation span. This smearing of the occultation link makes the spherical symmetry approximation (used in the Abel inversion process) invalid for grazing occultations-which is the reason why these occultations are considered less useful for sensing the neutral atmosphere (0-50kmf However, because the scale height in the ionosphere is nearly 10 times larger than that of the lower neutral atmosphere, the horizontal scale over which the occultation links are sensitive in the ionosphere is nearly 10 times that of the neutral atmosphere (i.e., -2000 km in the iono sphere vs. 200 km in the lower neutral atmosphere). This fa ct allows us to combine links from nearby occultations in a tomographic manner when imaging the ionosphere. By examination of Fig. 3, it is possible to see that grazing occultations would be of extreme value to iono spheric sensing since they provide nearly orthogonal cuts to those of the setting and rising occultations. Even though the numbers of rising, setting and grazing occultations in Fig. 2 are nearly the same (1867, 1878 and 1868, respectively), their distribution as a function of latitude is considerably different as indicated in Fig. 4 which shows the average number of daily occulta tions per (I 000km)2 at different latitudes. Based on Fig. 4, we estimate that the distance between occultations for COSMIC (in one day) is on the order of 250-500 km depending on latitude, and on whether grazing occultations are used.

Vertical Resolution
The theoretical vertical resolution of ionospheric profiles obtained from GPS-LEO occul tations is set by the diffraction limit, which corresponds to the first Fresnel zone diameter given by Where Dcps and DLE0 are the distance from the tangent point ( Fig. 1) to the GPS and LEO, respectively. In the case of COSMIC, Dc P s "" 30,000 km and DLEo "" 3,000 resulting in D F "" 1.5 km. However, a more practical limit on the vertical resolution for the ionosphere is intro duced by horizontal structures which, if ignored, can alias into vertical structures. Therefore, the horizontal resolution depends on the retrieval scheme and the ability to capture the large horizontal scale. This issue will be addressed further later.

Topside Ionosphere and Plasmasphere coverage
In addition to the occultation links which probe the ionosphere below COSMIC height, COSMIC will collect positive elevation GPS data at 0.1 Hz rate which measure line-of-.s'ight TEC in the topside ionosphere and plasmasphere. Each COSMIC satellite will observe nearly 10 GPS satellites at positive elevations simultaneously in different directions of the sky pro viding a very dense coverage and mapping of electron density above the COSMIC altitude.

Complementary Ground Coverage
The continuously operating networks of GPS ground receivers provide a very rich set of GPS-TEC ground data complementary to the COSMIC space data. Data from over 200 GPS globally distributed ground receivers, each tracking 8-12 GPS satellites simultaneously ( Fig.  5), currently constitute a major network for sensing the ionosphere. These TEC data, for instance, were used to provide the first global-scale ionospheric maps of vertical TEC (Mannucci et al., 1998;Lu et al., 1998).

Abel Inversion
When the ionospheric advance (for phase) or delay (for range) ( lk = +(d I h,2 )TEC, see Eqs. l and 2) is determined, it is differentiated (after proper smoothing) to obtain the extra Doppler shift induced by the medium. This extra Doppler shift can be used to derive the bending of the signal, a, as a function of the asymptotic miss distance, a (Fig. 1 ), by assuming a spherically symmetric atmosphere in the locality of the occµltation; the relationship between the direction of the signal's propagation and the extra Doppler shift, iJJ, is then given by where v and v are the transmitter and receiver velocity, respectively; k and k are the unit L r I r vectors in the direction of the transmitted and received signal, respectively; and k is the unit  where rt and r,. are the coordinates of the transmitter and the receiver, respectively; and n is the index of refraction. Eqs. (8) and (9) can be solved simultaneously in order to estimate a (a).
Solving Eqs. (8 and 9) requires knowledge of n at the satellites locations; however, the error acquired by setting n=l is negligible (Hajj and Romans, 1998;Schreiner et al., 1999).
The spherical symmetry assumption can also be used to relate the signal's bending to the medium's index of refraction, n, via the relation (Born and Wolf, 1980, p where a (r) = r x n(r) and r is the radius of the tangent point. This integral equation can then be inverted by using an Abel integral transform given by (see, e.g. , Tricomi, J 985, p. 39) Electron density can then be derived by use of the relation (-for phase advance, +for group delay) Eqs. (8)(9)(10)(11)(12) can be approximated and reduced to a single equation by ignoring bending. From Eq. (8) it is easy to establish the following relation between ionospheric bending and advance/ delay (I) where dldt is the time derivative, DcPs and Dwo are defined in Eq. (6), Vc P s.L and VLE O .L are the GPS and LEO velocity components perpendicular to the GPS-LEO link, r:,. is the radius of the tangent point of the straight line connecting the two satellites. By use of the approximation a = m ""' r1, Eq. (13), and expanding ln(n) in Eq. (11) we obtain Note that in both of these formulations, onl y the derivative of the ionospheric delay or TEC is used and therefore a constant bias is irrelevant. The upper limit of the integrals in Eq. (11 or 14) requires knowledge of a (a) or the derivative of TEC(r) up to the top of the ionosphere. While the GPS satellites are above the ionosphere, this is not true of the COSMIC instrument at -7170 km radius. Different methods have been used to overcome this for the GPS/MET experiment at 7110 km radius. Hajj and Romans (1998) have used an exponential extrapolation of a(a) based on information from a < 7110. Schreiner et al. (1999) have eliminated the contribution of the ionosphere above the LEO satellite b y assuming the occultation to be planar, the GPS satellite to be at infinity, and the ionosphere to be static over the duration of the occultation. In this way the contribution from the topside can be removed by subtracting a TEC measurement taken at point A ( Fig. 6) from that at point B.

Results from the GPS/MET Experiments
Many features and examples of GPS/MET retrievals can be found in (Hajj and Romans., 1998;Schreiner et al., 1999;Straus et al., 1999a and1999b). In this section we summarize some key features and present some new comparisons. Millstone Hill (42.6N & 2S8.5E) taken within 20 minutes of the time of the occultation. The general agreement is very good. Discrepancies between the ISR and the occultation profiles can be due to several factors including mis-collocation between measurements, spherical sym metry assumption of the Abel inversion and ISR measurement noise.

Examples of ionospheric retrievals
Other examples of GPS/MET electron density are shown in Figs. Sa, Sb and Sc which show one, several, and two sharp layers, respectively, at the bottom of the ionosphere 9 . The sensitivity to ionospheric structures of very short vertical scales, such as sporadic E, is a conse quence of the high vertical resolution charncteristic of the GPS occultation limb viewing ge ometry. (Data rate for retrievals in Figs. 7 and S is 0.1 and 50 Hz for points above and below 120 km, respectively, which explains the much finer resolution below 120 km.) It should be noted, however, that the gross structure in Fig. Sa below 200 km (i.e., the increase in density towards lower altitudes) is not realistic for the near-midnight local time of that occultation. This feature is most likely due to horizontal gradients which have a stronger effect on the retrieval at lower altitudes both in absolute (because the ray travels a longer distance in the ionosphere) and fractional (because the E-region densities are lower) terms.

NmF2 Comparisons with lonosonde
Ionospheric profile retrievals were derived at the University Corporation for Atmospheric Research (UCAR) using single frequency TEC data (as described in section 2.3.3) and the Abel transform for 40,000 GPS/MET occultations collected between Apr. 1995-Feb. 1997. Fi gure 9a shows a scatter plot of NmF2 values obtained from GPS/MET and nearby ionosonde. Figure 9b shows a histogram of the difference (NmF2(GPS/MET)-NmF2 (ionosonde)) with a standard deviation, CJ, of 0.15 x 1012e/m3• NmF2 values with differences larger than 4-a , (a total of 190 values or 4.3% of the entire data set) were excluded from the statistics.
We further examine these NmF2 differences as a function of local time (LT). Figure  Even though the differences are larger during the day in an absolute sense, the fractional difference is largest (-40%) during the pre-dawn sector and smallest (-20%) around local noon. The larger fractional error of NmF2 during the nightside reflects the fact that iono spheric horizontal gradients are generally larger on the nightside than on the dayside, which causes the spherical symmetry assumption used in the Abel inversion to be more in error.
These results are consistent with those found by Straus ( l 999a and 1999b ) . Differences of 20% or 40% in NmF2 implies differences of 10% or 20 % in foF2, respectively. It should be noted that these are average statistics for all seasons near solar minimum conditions (1995)(1996)(1997).

"HORIZONTALLY-CONSTRAINED" INVERSION
The dominant error in the Abel inversion is due to the spherical symmetry assumption imposed on the ionosphere. However, when information on horizontal structure is available, more general inversion schemes can be constructed. Hajj et al. (1994) presents one possible scheme where the ionosphere is divided into pixels on the basis of radius, latitude and longi tude. Pixels at the same layer (i.e., same radius) are then constrained to have a given functional  form with an undetermined scale factor. This functional form can be the same for all layers, or different for different layers. The problem is then reduced to solving for the single unknown scale factor for each layer, reducing the 3-D under-determined problem to a one-dimensional problem with a unique solution. Following the analysis by Straus (l 999b), in order to study the effects of horizontal gradi ents we replace the one dimensional electron density function, nJa), with a two-dimensional density n,.(a, </J ), where </J is the angle of a point along the ray path from the tangent point as measured from the Earth's center. We further assume that nJa, </J) can be separated into the product of two terms: nJa,O), the electron density profile at the tangent point, andf( ¢),a function which provides the ratio of electron densities along the line of sight to the density at the tangent point (note thatf( O) = I). It can be readily seen that only the portion off( </J) which has even symmetry about the tangent point (we denote by fJ ¢))contributes to gradient-in duced retrieval errors, since errors introduced by the odd terms will cancel upon integrating along the ray path.
To assess the utility of various types of ancillary data for constraining the horizontal structure and improving occultation retrievals, simulations of both occultation limb soundings and three different types of ancillary data have been performed using the RIEG ionospheric model (Reilly, 1993) to provide a "known" ionosphere based on climatology. The simulated ancillary data are in-situ plasma density, nadir-viewing extreme ultraviolet airglow, and total electron content maps derived from ground-based GPS receivers. The second type of mea surements will be available from an ionospheric photometer sensor on COSMIC. The third type of data is routinely available from Global Ionospheric Maps (GIM) derived at JPL. To focus on correction of errors strictly due to horizontal gradient effects along the line-of-sight for limb-viewing measurements, an idealized occultation geometry is imposed which is only susceptible to this particular error source. That is, we simulate occultations that include no motion of the tangent point throughout the measurement. Although the case of zero tangent point motion never actually occurs in practice due to the motion of the GPS satellite during an occultation, it is most closely approximated by an "in-track" occultation for which the direc tion to the GPS satellite is in the LEO orbit plane. Use of this idealized occultation geometry also allows us to perform the analysis on a regular spatial grid, avoiding the irregularly spaced spatial distribution of occultations which actually occurs under the realistic conditions gener ated by the orbital motions of various spacecraft involved.
The simulation proceeds as fol lows. A latitude/longitude grid, with 5° and 10' spacing, respectively, is constructed and a look direction for all occultations is selected. The driving parameters for the RIEG model (Fl0.7, Kp, day of year, universal time) are chosen and fixed. For each point on the grid, an in-track occultation (i.e., in the selected look direction) with zero tangent point motion is simulated. Simulated values of the various ancillary measurements are calculated along the plane containing the occultation lines-of-sight. For the cases in which these measurements are made by secondary sensors on the LEO satellite, this corresponds to the "in-track" assumption in which the satellite passes directly over the tangent point for the occultation. Electron density profiles are then determined using the Abel inversion and the gradient-corrected inversions for the three diffe rent types of ancillary measurements. The accuracies of the retrievals are judged in terms of NmF2 and HmF2, as compared to the "true" model values at the each of the grid points. The process is repeated fo r different look direc tions and ionospheric model drivers.
For this study, an additional simplifying assumption has been made in the case that the gradient information is obtained using the EUV technique. The specific emission simulated is the 91 J A emission associated with radiative recombination of atomic oxygen ions (Meier, 199 1;Dymond, et. al., 1997). The volume emission rate of this feature is proportional to the product of the oxygen ion and electron densities. For the purposes of this analysis, we have assumed that these densities are equal, so that the observed airglow is proportional to the integral of the square of the electron density below the satellite. An estimate of the horizontal gradients is then obtained from the variations in the square root of the simulated observations. To the extent that light or molecular ions are present, particularly at altitudes near the location of the maximum emission rate , this assumption will break down. However, since oxygen is known to be the dominant ion near the F2 peak under nearly all conditions, and since only the relative horizontal changes in the airglow magnitude are used to determine the gradients, our assumptions should introduce an insignificant error most of the time. Although the 91 1 A feature can be simulated globally, on the dayside it is difficult to isolate in practice from other dayglow emissions, limiting actual use of this correction technique to the nightside.
Error maps for NmF2 for one particular set of simulation results are shown in Figs. 11. This run corresponds to solar maximum, equinox, geomagnetically quiet conditions at 0 UT, with a northward occultation look direction (i.e., the direction in which the strongest average gradients are expected). The "true" ionospheric model values of NmF2 and HmF2, shown in Fig 11 a and 11 b, provide an indication of the horizontal structures in this simulation. Frac tional errors in NmF2 and actual errors in HmF2 obtained when using the simulated EUV measurements to constrain the horizontal structure are shown in Figs. I le and 1 ld, respec tively. A statistical summary of other simulations is given in Table 1. As might be antici pated, the largest NmF2 errors occur in the auroral zone and near the equatorial arcs, where gradients are strongest. Errors in HmF2 follow a similar phenomenology. On average, tech niques employing a gradient correction provide significant improvement over the spherically  Straus, 1999b) symmetric Abel transform, particularly for NmF2. In general, the TEC and EUV techniques for gradient correction outperform the in-situ based method. This is reasonable, since these data types provide information on gradients more localized to the vicinity of the peak of the layer, whereas the in-situ method relies on extrapolation of gradients measured at the satellite altitude to lower regions. However, improvement using any of the three gradient correction techniques is not guar anteed in individual cases. For example, consider the point near the peak of the equatorial anomaly near sunset at -10° latitude/270° longitude. Figure 12a shows ionospheric density contours in the plane of the occultation together with four representative lines of sight corre sponding to tangent altitudes of 300, 400, 500, and 600 km. The corresponding true profile at the tangent point is shown in Fig. l 2b together with Abel, in-situ, and TEC-based inversions. Despite the strong gradients in this region, the Abel inversion fortuitously performs best in reproducing the true NmF2. None of the techniques reproduce HmF2 or the bottomside par ticularly well. The situation can be understood by considering Fig. 12c, which shows the actual and estimated horizontal gradient function in the ionosphere, f ( ¢ ). Here the solid e curves represent the even asymmetries encountered by rays with the same four tangent altitudes as in Fig. 12a (the largest asymmetry corresponds to the lowest tangent altitude). At the highest tangent altitude, the asymmetry is modeled well by the in-situ based estimate (dashed curve), but the shape of the asymmetry function for lower altitudes is different. For tangent altitudes near HmF2 (between 400 and 500 km), neither of the two asymmetry estimates pro vides a good match to the actual asymmetries. This occurs primarily because the height of the layer is significantly varying within the region sensed by the occultation. If the height of the layer were not varying, then the gradients in the vicinity of the peak would be accurately  reflected in at least the TEC-based measure, and possibly the in-situ-based measure as well.
In addition to the equinox, north-looking case presented here, summer and winter solstice cases have been simulated. The overall statistics for these cases are shown in Table 1 and are similar to the equinox case. However, east-looking statistics, also shown in Table 1, are sig nificantly better for all methods. This is because ionospheric gradients are typically less strong for this direction of viewing. The most significant retrieval errors for the east-looking geom etry are near the sunrise terminator, in the auroral zone, and in the equatorial zone near the point where the magnetic equator bends away from the geographic equator, resulting in more significant gradients in the geographic east-west direction.
The simulation results presented in this section do not include ancillary data noise effects and are therefore optimistic. When realistic errors of the ancillary data are considered, the improvement seen in the ab ove simulations may not be as dramatic. For instance, Schreiner et al. (1999) have examined the use of vertical TEC ancillary data to constrain GPS/MET retriev als, but their derived NmF2 did not 11how any significant improvement over the Abel inver sion. The reasons for the discrepancies between the simulations presented here and statistical results presented by Schreiner et al. are not well understood, and further simulations and data analysis are needed to reach more definitive conclusions.

TOMOGRAPHIC INVERSION
With the straight line propagation assumption, we can combine data from several occulta tions and obtain a linear system of equations given by where y is the m'th differential TEC measurement (TEC differences between consecutive "' measurements are used in this process to get rid of all TEC bias terms); x; is the electron density in a cell i (here we divide the ionosphere into pi xels and assume the density to be constant inside each pixel), A . is an element in the observation matrix relating the observation m to the density in pixel i; d"· ' is the standard deviation of the diffe�ential TEC measurement km) and horizontal (200-500 km) resolution will become available for the first time.

General Approach
The succes . : s of ionospheric modeling depends mainly on accurate knowledge of the forces (drivers) which enter into the hydrodynamic plasma equations for the ionosphere. These in clude solar EUV and UV radiation, magnetospheric electric fields, particle precipitation, dy namo electric fields, thermospheric neutral densities, temperatures and winds. Therefore, ide ally, one needs a coupled magnetospheric-thermospheric-ionospheric model that assimilates global and continuous data in all these regions. This is a goal of the space weather program which is likely to be reached and become operational Gust like operational weather models) in a few decades. However, a first step toward reaching that goal is to learn to do data assimila tion optimally in each region independently, before moving on to the more complicated coupled models. Furthermore, it is important in data assimilation to carefully study the impact of each new data type and its noise on the stability and accuracy of the model.
In order to assess the use of GPS-TEC data obtained from ground and space on iono spheric specification and forecasting, we start with a simplified low-and mid-latitude iono spheric model and use the Kalman filter in the context of an Observation System Simulation Experiment (OSSE). An OSSE is a controlled simulated experiment where data are synthe sized based on the "true" state of the medium being studied, and then assimilated after adding the proper kind of noise. OSSEs are commonly used in NWP to assess the impact of new data types and/or new assimilation schemes (Atlas, 1997).
In our OSSE scenario (Fig. 15) a "true" ionospheric state is generated based on model l (to be described later). This "truth run" is meant to represent the true state of the ionospheric electron density. Data are then generated by assuming a distribution of ground and space receivers tracking GPS. Measurement and representation error covariances are assigned to the synthesized data . A perturbation run is generated by a different model (model 2). The assimi lation run uses the apriori state produced by the perturbation run and model 2 for forward propagation in time, in conjunction with the data. The impact of the assimilated data is then evaluated by comparing the difference between the assimilation run from the truth (error 2) to the difference between the perturbation and the true runs (error I). More detailed descriptions of each of these steps are given below along with some results. But first we introduce some basic terminology and summarize the Kalman filter.

.2 Basic Terminology and the Kalman Filter
We introduce the following definitions (commonly used in NWP) The true state: a discrete representation of the true underlying conti nuum of the ionospheric state (electron density) at time k.
The anal y sis: an estimate of x� given measurements at time k, and a forecast x{ .
The forecast: an estimate of xJ, given measurements up to time k-L The observations are related linearly to the true state through an ob servation o p erator Hk . £0 is an observational error. The observational error is composed of the measurement error, t:"' , and a re p resentativeness error, £,..
The true state at time k+ 1 is related to the state at time k via the for ward model 'Pk and a model error. Let Mk, Rk and Qk be the measurement, representativeness and model error covariances, re spectively. Then the Kalman filter can be summarized by the following set of Equations:

) an
Fk + l -TkFk T k + k' e K is known as the Kalman gain, P" and pl are the analysis and the forecast covariances.
Some of the terms defined here will be explained in detail below (see also Howe et al., l 998;Ruffini et al., 1999). During the data assimilation process, time is discretized (indexed by kin Eqs. 18) where the state during a time step (taken to be 15 minutes in our analysis) is assumed constant. At time k = 0, starting with a forecast (apriori) state, x , ; · , a forecast state covariance, P/ , and a set of observations, m2 , we can update the state by use of Eqs. (I 8a-c) to obtain the analysis and its covariance at time k = 0. Using the forward model 'Pk we can propagate the analysis and its covariance from time 0 to time 1 via Eqs. (18d and 18e) to obtain the forecast state and covariance to be used as apriori for the next time step, and so on.

Dynamical Modeling
Our physical modeling is based on a simplified version of the University of Sheffield Plasmasphere/Ionosphere Model (USPIM) (Bailey and Shellek, 1990;Bailey et al., 1993) for middle and low-latitude regions. In our model we solve the conservation of mass and momen tum equations for a plasma accounting for production, loss, transport of the maj or ion species in the F-region (O + ), neutral wind, and electric forcing. These equations can be written as an at +V•(nV) =(P-L), n is the ion number density; V is the velocity; P and L are production and loss rates, respec tively ; k11 is Boltzmann's constant; Tis temperature, Mis the ion mass, g is gravitational accel eration, E and B are electric and magnetic fields, respectively; v is the collision frequency for moment transfer between the atomic oxygen ion and the neutral particles; U is neutral wind. An equation similar to Eq. (20) can be obtained for electrons, and after ignoring terms that are multiplied by the electron's mass, we obtain In addition we also have -\7(nek8TJ -ene (E + "'., X B) = 0.
(2 1) We solve for the ion and electron densities by solving the above equations and making use of the empirical or parameterized models of the thermosphere (MSIS), thermospheric winds from the Horizontal Wind Model (HWM), solar EUV (SERF2, (Tobiska, 1991)), and electric fields (e.g., Heppner and Maynard, 1987;Fejer and Scherliess, 1997).
Traditionally, the equations of many ionospheric models are rewritten in a moving Lagrangian coordinate frame. The motion of this coordinate frame is dictated by the plasma drift perpendiculru· to the geomagnetic field lines. That approach introduces significant com putational efficiency by transfor ming a time-dependent partial differential equations in 3-D into a family of time-dependent differential equations in 1-D. However, this 1-D approach introduces a significant difficulty to data assimilation, since the measurements are taken in 3-D space from ground-to-satellite and satellite-to-satellite links, making the mapping between data and the model parameter space very difficult to construct. In addition, the moving frame creates technical difficulties in applying some well-established data assimilation algorithms such as the adj oint approach.
Because of these considerations, we solve Eqs. (19-22) in a fixed Eulerian frame. Figure   16 shows the region being modeled and its division into pixels. These pixels are bounded by constant geomagnetic field lines, constant geomagnetic potential lines and local time, where the B field is that defined by an eccentric tilted dipole.

"Truth" Run
A "true" ionospheric state was created as a function of latitude, height and local time. The 2-D region modeled and the p-q coordinates are depicted in Fig. 16. The region modeled is bounded by (-30N, + 30N) geomagnetic latitude, 100-1000 km altitude, with 2 deg. horizontal and 40 km vertical resolution (a total of -1100 pixels). The run was done for 16 hours between 700 and 2300 LT (after an initialization run of 72 hours) based on EUV, EXE zonal drift, electron temperature, and neutral densities, temperatures and winds, derived from climato logical models for solar maximum summer solstice conditions. The EX B zonal drift used is shown in Fig. 17 (red curve).

Perturbation Run
A perturbation run (to be used as an apriori) was also generated for the same region and time period by use of the same physical model described in Sec. 7.3 except that some driving forces were modified. In particular, we use a different vertical EX B drift (shown in blue in Fig. 17), and reduce the neutral density of atomic oxygen by a factor of 2. The true run, perturbed run and the difference are shown at 700 LT in Fig. 18. The electron densities for the perturbed run are nearly half of those of the "truth" run and they remain so at all local times, mainly due to reducing the neutral atomic oxygen density by a factor of 2.

Generation of Data
For the purpose of our study we assumed 4 GPS transmitters and 6 GPS receivers: 5 ground receivers evenly separated by 20 deg. tracking GPS at the rate of one data point every 3 minutes, and one LEO receiver at 700 km altitude with fore and aft antennas tracking the GPS satellites at the rate of 0. 1 Hz. All stations and satellites are assumed to stay in the same 2-D magnetic longitudinal plane at different local times of the runs, which implies that the satellites are precessing at the same rate as the Earth. This is clearly not the case in real life, but with a constellation of LEO satellites such as COSMIC, it is not far fe tched since one can combine data from different LEO and GPS satellites at various local times. In practice, the problem is inherently 3-D, but our aim here is to simplify the problem sufficiently in order to demonstrate the use of ionospheric data assimilation. A further simplification is made by using absolute TEC (as suppose to differential TEC) as our measurements and ignoring any biases (which are normally -1-3 TECU as described in Sec. 2.3.2). Figure 19 shows the data links for the first two hours of our OSSE. The color indicates the time when the data were collected. The dense samples correspond to LEO data acquisition starting from 30 deg. eleva tion relative to the LEO local horizon plane down to the surface. Some of the links do not pass through the region being modeled and were therefore discarded.
Given the receiver-transmitter geometry, synthetic TEC values are generated based on the "underlying continuum of the true" electron density (which will be explained more later). Line-of-sight TEC values are shown in Fig. 20a   shapes in Fig. 20a corre spond to two occultations, one rising and one setting as seen in Fig.  20b. Since our OSSE did not include a plasmasphere, TEC is nearly 0 at the beginning or the end of an occultation. On average, there are nearly 200 ground data links and 2 occultations every 1 hour.

.Observational Errors
A major part of data assimilation is to assign appropriate covariances to the measure ments and representativeness errors. Measurement errors are due to instrumental noise. Rep resentativeness error is due to the fi nite representation of the model to the underlying con tinuum both in time and space. To understand the latter type of error better we consider the difference between TEC data generated based on a very fine space and time gridding (which is the actual data being assimilated) and TEC data for the same links and same driving forces but for the coarser spatial (shown in Fig. 16) and temporal grid (15 minute). The difference in TEC for each link is a measure of the representativeness error and is shown in Fig. 2la and  21 Based on the histogram of the representativeness error, and by excluding differences larger than 20 TECU, we estimate a 1a representativeness error of 2 TECU. Therefore, in our assimilation we take R in Eq. (18b) to be ' xa� where l is the unit matrix and <Jp= 2 TECU. We also take Min the same equation to be j! XO'� where a 0 = 0. 1 TECU.

Assimilation Runs
In order to carry out the assimilation process according to Eqs. ( 18), we still need to define the covariance for the apriori state, P 0, and the model error covariance, Qk. We choose P0 and Qk to be diagonal matrices with elements given by Pa = Diag(n;1 (O),n;2 (0), n;3 (0),...  where ne (k) is the electron density in pixel i at time k. Fig. 22a shows the "true" electron density at 1900 LT. Fig. 22b shows the electron density "analysis" obtained by using the recursive formulas of Eqs. (18). Fig. 22c shows the difference between the truth and analysis. The "true" electron density state, including the equatorial anomaly, is very well captured by the analysis, with fractional error of -5% near the maximum densities. This is encouraging because it implies that when there is a sufficiently dense set of TEC links, as would be the case with COSMIC, it is still possible to get accurate electron densities even if the driving forces are poorly known (e.g., substantially reduced neutral density and no prereversal enhancement in vertical drift). However, as demonstrated below, the physics still has a significant impact on the inversion. For operational space weather purposes, some practical considerations such as computer memory and time limitations ought to be considered. Even though we will not discuss these issues in detail here, we do point out that no existing operational weather model (for the neu tral atmosphere) uses the full Kalman filter (i.e., keeps and updates the full covariance of the state, which is usually of order 10 6 parameters). For the recursive Kalman filter 16-hour run presented above, it takes nearly 8 hours of CPU time on a dual processor DEC Alpha to as similate -4000 links. As a first attempt to reduce this processing time, we tested an approxi mate solution to the Kalman filter where Eq. 18c is replaced by while the other equations in (18) are kept the same. This saves us the expensive operation of computing Kk HkP /; the state covariance is still updated through Eq. (18e). This reduces the computer run time from 8 hours to 1 hour without any significant impact on the solution as shown in Fig. 22d and f. The reason for the similarity in the two solutions is that the rank of Hk  (the observation matrix for data taken within J 5 minutes) is small. Therefore, the impact Kk Hk P/ has on updating P/ through Eq. (18c) is much smaller than the impact of the for ward propagation of P; ' and the addition of the model error covariance through Eq. (l Se).
This is the first evidence that the physical model plays an important role in the assimilation.
To further assess the impact of the physical model we test a further approximation where the recursive inversion is done using Eqs. (1 Sa, l Sb, 24) and (25a) and P/� 1 =Pi:", (25b) where no forward model is used. We refer to this approximation as tomography, although this is not strictly accurate, since in tomography Eq. (18c) can still be used. The analysis of this approximation and the associated error are shown in Figs. 22f and g. The degradation of the reconstruction relative to the previous solution is apparent. The benefit of data assimilation can be seen in Fig. 23, which is a measure of the performance of each technique: no data assimilation (perturbed run), tomography, approximate, and full Kalman.

BENEFITS TO SPACE WEATHER
There are two important aspects to space weather, one scientific and one operational. On the scientific side the focus is on trying to improve our understanding and modeling of space weather, while on the operational side the focus is on the users and systems that might be affected by the space weather environment. COSMIC will make a tremendous contribution to both these aspects in three main areas: ( 1) Electron density specification and forecasting which have been the main focus of our discussion thus far, (2) ionospheric scintillation and irregularities, and (3) TIDs. These three areas are not entirely decoupled, but for the purpose of clarity we divide our discussion below into the contribution of COSMIC to each of these areas.

Electron Density Specification and Forecasting
Accurate electron density specification allows us to better understand and model several  very important space weather related phenomena. Most importantly, it will allow us to under stand the response of the ionosphere to geomagnetic storms and substorms during their growth, main and recovery phases. Measurements of electron density will improve our understanding of the time scale over which the ionosphere responds to certain changes in the driving forces.
Continuous and global measurements of the ionosphere and plasmasphere are needed for bet ter characterization of diurnal, monthly, seasonal, and inter-annual cycles and for improving climatological models. In addition, electron density measurements will improve modeling of several space weather phenomena (see e.g., Schunk and Sojka, 1995) such as the production and evolution of boundary and auroral blobs, the generation of spread-F, sporadic E and inter mediate layers. In addition to direct measurement of electron density, several inferences con cerning the thermosphere and the magnetosphere can be made from know ledge of the time evolution of electron density. For instance, NmF2 and HmF2 and their variation in time can be used to infer dynamo electric field and thermospheric wind. On the operational side, specification of electron density can be used to estimate iono spheric delay along any line-of-sight which can be used for single frequency GPS users, deep space observations in the microwave band (such as from the Deep Space Network), single frequency altimeters (such as Geosat Follow On-GFO) and military surveillance applica tions.

Ionospheric Scintillations and Irregularities
Most of our discussions thus far has focused on electron density specification and fore casting from space TEC data. Equally important to space weather are GPS phase and ampli tude scintillation measurements. Similarly to TEC measurements, scintillation measurements integrated over different lines-of-sight can be combined in a tomographic manner to produce 3-D maps of irregularities (Kunitsyn and Tereshchenko, 1992). These measurements can also be used to localize electron density irregularities along the GPS-LEO ray path if the data quantity is insufficient for a tomographic localization (Vorob' ev et al., 1999;Sokolovskiy, 2000).
Furthermore, by assimilating scintillation measurements into an ionospheric scintillation model, much can be learned about the mechanism of the generation and dissipation of these irregularities. For instance, ground-based GPS TEC measurements, in conjunction with NCAR' s TIE-GCM, revealed that an irregularity event occurred across the northern U.S. continent during a major geomagnetic storm in October 1995 (Pi et al., 2000). The observed irregulari ties are believed to be due to the penetration of plasma convection and strongly sheared flow down to the middle latitudes, which not only deepened the middle latitude ionospheric trough but also may have produced plasma undulations and Kelvin-Helmholtz instabilities. An ex ample of ionospheric irregularities observed during a geomagnetic storm through the global GPS network using the ROTI index (Pi et al., 1997) is shown in Fig. 24.

Traveling Ionospheric Disturbances
Observations of the ionospheric E-and F-regions have shown the presence of oscillations resembling atmospheric planetary waves, typically observed in the mesosphere and lower ther-  [Pi et al., 1997] measures dual-frequency GPS phase fluctuations. One ROTI cor responds to 1 TEC Unit change per minute.
mosphere. These have indicated a 2-day oscillation in foF2, in addition to longer period waves of 5, 10 and 16 days, suggesting a causal connection between events in the neutral atmosphere and the ionosphere (Forbes and Zhang, 1996;Harris and Vincent, 1993). In their analysis, Forbes and Zhang use a database of foF2 and foE from ionosonde measurements extending over a long period of time, while they approximate HmF2 by an empirical relation. Their database is limited to 20 geographical locations corresponding to 20 ionosonde stations. Continuous, global and accurate recovery of electron density at all heights will facilitate a detailed study of the quasi two-day (QTD) wind oscillation as well as longer period oscilla tions at different altitude ranges. Because the GPS occultation technique provides accurate measurements of neutral atmospheric density below 60 km and electron density above that height, these measurements are particularly suitable to establishing a causal connection, if any, between waves generated at the stratosphere and lower mesosphere and those in the iono sphere.
As an example, Fig. 25a illustrates the sensitivity of radio occultation to atmospheric waves. Close agreement in both ampl itude and phase with a radiosonde sounding taken some 300 km away (the profile shown in blue is in almost exact agreement with the GPS/MET profile) implies that the wave has a horizontal wavelength that is much larger than 300 km. The horizontal and vertical scales of this wave are consistent with a 5-day westward propagat ing Rossby-gravity wave (Kursinski et al., 1996). On the other hand, electron density derived in the region between 60-120 (25b) shows an oscillation that might be indicative of a possible causal connection between the stratospheric wave and that in the D-and E-region. It should be noted that even though the electron density estimated in Fig. 25b could be biased (due to the spherical symmetry assumption used in the retrieval), the point-to-point structure of the pro file should be accurate.
COSMIC will provide the possibility to detect and analyze individual wave events in  68E). Also shown is a radiosonde profile obtained from a ship at 350 km from the occultation confirming the presence of what might be a Rossby-gravity wave. (Fig. (a) is after Kursinski et al., 1996) order to quantify and catalog the wave parameters, and to study and identify the relation be tween the source event and the resulting gravity wave and TIO.

CONCLUSION
COSMIC OPS occultations in the ionosphere will yield accurate· 3-D specification of electron density, ionospheric irregularities and global maps of TID's. The -5600 global daily occultations from COSMIC can be processed individually to produce high-resolution vertical profiles of electron density, or collectively by properly assimilating them into a physically based ionospheric model. This kind of continuous global 3-D specification of the ionosphere will have tremendous payoffs for improving our understanding and the modeling of various space weather phenomena, as well as serving operational needs such as communication and navigation. Just as the number of OPS ground receivers has grown substantially over the past decade, the number of space-based OPS receivers collecting occultations will al so witness substantial growth in the next decade, to the point where OPS ionospheric occultation mea surements will truly revolutionize the manner in which we sense the ionosphere and our un derstanding of it.