Interactive comment on “ Empirical model of the ionosphere based on COSMIC / FORMOSAT-3 for neutral atmosphere radio occultation processing

> This study discusses the use of an empirical electron density model with traditional > radio occultation processing packages. Also presented is a “scintillation” proxy index > for the identification of highly variable electron density profiles. Unfortunately, I see > this work as the combination of two separate studies that are each incomplete. Based > on the following reasoning, I believe that the present study is insufficient to warrant > publication and would require very significant modification to take it to a point where it > would be suitable for publication. I recommend that the authors undertake a significant > revision of the study and re-submit the study as two separate papers. >

EUMETSAT Polar System -Second Generation (EPS-SG) mission will contain a Radio Occultation payload whose primary objective will be to monitor the Earth's neutral atmosphere.The processing chain of the EPS-SG GNSS RO data will be based on a processor developed under EUMETSAT's lead for level 1B data (primarily bending angle over impact parameter), and on the ROMSAF's Radio Occultation Processing Package (ROPP) software for level 2 data (e.g.refractivity, temperature, see [Culverwell et al. (2015)]), which will allow to obtain the neutral-atmosphere-related products.In the current version of the ROPP software, the ionospheric effects (needed to estimate the bending angle of the neutral atmosphere) are modelled with an analytical formula.This analytical formula (called Zorro formula, see [Culverwell and Healy (2015)]) is the analytic integration of the ionospheric delay along the line of sight considering that the vertical profiles of electron density follow a Chapman layer model.
In this context, the EUMETSAT's ROPE study (Radio Occultation Profiling Evaluation, ITT number 15/210697) has, as one of the main objectives, to modify the ROPP Forward Model (FM) code so that it could accept an arbitrary electron density profile to model the ionospheric delay.While ROPE was primarily initiated to generate realistic ionospheric test data for processor testing, these modifications would also allow the retrieval of more accurate neutral atmospheric information.To this end, in order to test this new feature, an empirical model of the ionospheric electron density has been developed, also under the ROPE study.Such empirical model would provide electron density values for a representative set of ionospheric scenarios, where this new feature of the ROPP software would be then tested.
As an additional result of this effort, a new scintillation index based on the analysis of the vertical profile of electron density has been developed.Typically, when raw GNSS data is available, scintillation can be monitored using the C/N0 values (S 4 index, see [Kintner et al. (2007)]) or the carrier-phase measurements (σ φ or its lower-sampling rate version Rate of TEC Index, see [Prikryl et al. (2013)]).In this work, the Occultation Scintillation Proxy Index (OSPI) is being computed using the inverted profile, by analysing the variation in the topside of the vertical profile of electron density.
The paper is organized as follows: the following section includes a description of the ionospheric model and the data used for the generation of the model and different scenarios.The second section includes a description of the OSPI index, its generation and some comparison with current scintillation indices.The paper ends with a conclusions section.

Model generation
This section includes the description of the processing model to invert the observation profiles as well as post-processing tests and checks to automate the screening and discard unphysical or unrealistic profiles.

Profile inversion
The processing model described here uses the raw GNSS observables (i.e.dual-frequency pseudoranges and carrier phases), not the already processed profiles inverted through Abel transform by analysis centers such as e.g.UCAR ([Kuo et al. (2004)]).In order to obtain the vertical profiles of electron density from the GPS raw data, a modified peel onion mechanization of the Abel inversion has been used.This mechanization consists of: -Applying the Separability Hypothesis [Garcia-Fernandez et al. (2003)], which is a technique that uses Vertical Total Electron Content (VTEC) to account for the horizontal gradients of the ionosphere.Essentially, this hypothesis assumes that the electron density (3D field) can be separated into a horizontal component (the VTEC) and a vertical descriptor of the ionosphere (F(h) function or Shape Function) The shape functions F (h) can be understood as normalized N e profiles, with higher spatial and temporal correlation, as shown in Figure 1.
-Using Linear Mean Square (LMS) (instead of the iterative method) so that one can jointly process all occultation measurements with additional features such as phase bias estimation and bottomside constraints.See Figure 2 for an example of profiles inverted with LMS and Separability Hypothesis.The additional advantage of the LMS approach is the possibility to re-run various iteration, thus refining the estimation of ancillary parameters (such as phase bias or topside electron content).

Profile screening and selection
After a first COSMIC RO GNSS raw data check and editing, a (typically) large set of inverted profiles are available.However, not all these profiles are suited to build a profile database: some are incomplete, have artifacts or large noise or are outright  In order to perform this quality check in an automated way: -The profiles shall have a minimum height range, at least covering the main layers of the ionosphere (E, F1, F2, which translates in a typical height range between 150km to 500km) -Ideally, the integrated shape function (F (h)) along the vertical should in theory be 1 (i.e.inf 0 F (h)dh = 1) so that the integrated electron density (N e ) along the vertical profile yields the VTEC.In practice this does rarely hold because the profile does not usually account for the topside ionosphere.It should be, however, close to 1 (i.e.larger than 0.75).
-The maximum variation of the first order profile derivative should be less than 2000%.This threshold is wide enough to allow for ionospheric features with high variability but tight enough to discard unrealisitic profiles with large jitter due to measurement errors.
-The hmF 2 (maximum of the F2 layer peak) should be comprised between a physical and reasonable value (e.g.lower than a LEO satellite height but above the bottom E and D ionospheric layers).
-Due to the fact that negative N e values are unphysical, a strict positivity is required for all points of the profile.
One of the main advantages of this criteria is that it can be applied in an automated way.

Profile regularization
For the purpose of processing of EPS-SG radio occultation data, and in particular for the applicability of the proposed empirical model into the ROPP-FM package of the ROMSAF software, it is necessary to regularize the profile and avoid sudden changes in it.In particular the following points were enforced in every profile that passed the screening test: -The topside of the profile has been extrapolated up to the ESP-SG orbit height (about 820km or beyond if need be) using an exponential model: where the coefficients α, β and h 0 have been computed by fitting the available topside of the profile, beyond the hmF2 (maximum of the electron density peak) plus a certain margin (typically 50km to 100km upwards).
-In order to avoid artifact errors when using the profiles in ROPP-FM, a smooth transition to N e = 0 in the bottomside was enforced.This was done by means of a masking function based on a sigmoid (see Figure 3).This sigmoid function defined a transition zone (from 1 to 0) of 20km that covered the last available samples of the profile bottomside.This approach allows guaranteeing a smooth transition to 0 while reducing at a minimum the modifications to the original profile.
-Removal of pseudo D layer.As it is known, using the Abel inversion implies that the error in the N e estimation increases with decreasing height.This is due to the fact that the retrieval of the lower layers need the estimation of the upper layers, thus the error accumulates in the bottomside.This can cause some artifact errors that in some cases might seem a fictitious sporadic layer in the D layer, which is not realistic (see example in Figure 4).To mitigate this effect, if a profile showed a peak under 90km (the upper boundary of the D-layer, see [Mitra (1951)]), this was considered a false peak in the D layer and thus the samples from 90km downward were removed (preserving, however, the rest of the profile).

Scenarios
As mentioned earlier in the paper, the effort of building a database of profiles was intended to give an empirical model for the ionosphere to correct its delay in the retrieval of neutral atmosphere using the Abel inversion.So far, the ROPP software 5

Scintillation and wave affected profiles
In order to provide a model for the scintillation, as present in the vertical profiles of electron density (or shape functions), the first step is to evaluate its morphology.In order to do this, the COSMIC occultations for a single day of Scenario 1 (2011, day of year 264) have been visually inspected (mainly looking at the shape of the F2 layer) in order to identify those profiles apparently affected by either scintillation (i.e.jitter-like noise in the profile) or wave-like structures.

5
The results of this manual selection is shown in Figure 6.The upper panel of this plot shows the distribution of all inverted profiles from the COSMIC satellites as well as the selected events with scintillation (red) and selected events with waves (green).Over all 757 occultations for this day 1 , 119 of them have been labeled as profiles affected with scintillations and 7 have been labeled as profiles affected by wave-like structures.In order to illustrate the morphology of scintillating as well as wave-affected profiles, the figure include examples of both in the bottom panels.The distribution of the scintillation events

Scintillation index
In order to automate the identification of the profiles with scintillations, the Rate of Total Electron Content Index (ROTI), used to identify scintillating environments in the ionosphere, is being adapted in the case of electron density profiles.This proxy index (named Occultation Scintillation Proxy Index, OSPI) is defined as the standard deviation of the variation between consecutive vertical profiles samples (of the topside ionosphere sampled at ca. 1 to 3km) normalized by the value of the profile at the maximum (i.e.N mF 2).This is expressed by the following formula: where σ denotes the standard deviation operator -∆N e is the difference between consecutive samples of the electron density vertical profile.Only the samples of the profile comprised between h topside,min and h topside,max are being considered (typically 550km and 650km above Earth's surface respectively).
-N mF 2 is the maximum electron density of the profile.
As an example, Figure 7    -The OSPI at 600km (i.e. between 550km and 650km) is the one that provides the best agreement relative to the profiles labelled as scintillating using the naked eye approach, compared to the OSPI at other heights (such as e.g.F2 layer or E layer, at lower height intervals) -The OSPI index (i.e. standard deviation normalized with the electron density) outperforms other similar index such as 10 the standard deviation of the electron index normalized by the RMS of the electron density along all the profile or the standard deviation of the electron index without normalization.
-The right panel of Figure 9 shows the consistency between the proposed indices (taking as reference the OSPI at 600km), which is larger compared to the case when taking as reference the naked eye.The correlation check between the scintillation parameters given by COSMIC and the proposed OSPI index is provided in Figure 10.Even though the correlation seems weak between OSPI and S 4 or SN R L1 , the upper pictures show that there are some dependency that can be exploited.Note however that a perfect agreement of the S 4 and SN R L1 relative to the naked eye profiles or with the OSPI index are not expected because these indices refer to the scintillation in amplitude, rather than the scintillation in phase, which is the one affecting the vertical profiles.

Conclusions
As part of the activities of the EUMETSAT's ROPE study the ROMSAF's ROPP software has been updated to add a new feature consisting of the processing of inverted profiles of electron density in order to retrieve neutral atmospheric bending angle more accurately.This is in particular valuable for scenarios where the ionosphere deviates from its mean behaviour (high geomagnetic activity, sporadic events in the E-layer, travelling ionospheric disturbances,...).In this context, an empirical model of the ionospheric electron density has been developed for this study.This model has been built using raw GNSS data from COSMIC constellation, using an Abel inversion mechanization based on Linear Mean Square (rather than the classic peel onion approach) and the Separability Hypothesis (which overcomes the spherical symmetry assumption of the Abel inversion).
This dataset has been developed for 4 characteristic scenarios, covering various states of the ionosphere.
In addition to the ionospheric model, a proxy index for scintillation monitoring based on the retrieved profiles of electron density (OSPI) has been proposed.Results have shown that the OSPI has relevant correlation with the S 4 amplitude scintillation index.

Figure 1 .
Figure 1.Separability hypothesis generates shape functions (F (h)) instead of electron density profiles.The plot shows two inverted occultations that are ca.2000km apart.The similarity between both Shape Functions is higher than the corresponding electron density profiles (obtained assuming spherical symmetry).

Figure 2 .
Figure 2. Example of LMS inversion for a COSMIC occultation in combination of Separability Hypothesis.
unphysical.Therefore they have to undergo a quality check, as already suggested by several authors (e.g.[Uma et al. (2016)]).

Figure 3 .
Figure 3.A transition of the bottomside to reach Ne = 0 at the bottomside was enforced by means of a masking function.

Figure 4 .
Figure 4. Example of a pseudo D layer in an inverted profile.In this case, the samples below 90km are removed.

Figure 5 .
Figure 5. Scenarios of the EUMETSAT's ROPE study.Each period (grey rectangle) is shown against the Solar Flux and the Kp index, to illustrate the geomagnetic and Solar activity conditions of each scenario.

10 shown
in top panel of Figure6(mostly lying at southern hemisphere and top latitudes) is consistent with Figure1ofBasu and Basu (1989) (in page 134).1This number corresponds to the complate data set of profiles for this day without the quality check based on the solving strategy, described in the paper 7 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-217Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 18 July 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 6 .
Figure 6.Distribution of occultations (in local time and latitude) for the COSMIC satellites on 2011 doy 264.Manually detected occultations that are scintillation or wave affected are also included in the plot.The bottom panels include two examples of scintillation-and wave-affected profiles.
includes the OSPI values for a clean profile (without scintillation after visual inspection) and a profile affected with scintillation.Note that the OSPI values can be one order of magnitude larger depending on the scintillation effect.More generally, the histogram of the OSPI values for clean and scintillating profiles (obtained from COSMIC on 2011 / 264) are shown in Figure8.Among the 119 profiles that have been manually labeled 'with scintillation', 68 of them yielded an OSPI value greater than 0.003141.This represents a 57% of coincidence.There is a dubious area with OSPI values between 0.001 and 0.0031 with no clear distinction between clean or scintillating profiles.Despite this uncertainty region, this proposed index can constitute a useful tool for automatizing the detection and selection of scintillation-affected occultations.Further tuning can also improve the performance.

Figure 7 .
Figure 7. Example of OSPI values for a clean profile (left) and a profile affected by scintillation (right).

Figure 8 .
Figure 8. Normalized histograms of the OSPI values for the COSMIC profiles of 2011 / 264.The histograms have been split into OSPI values from clean profiles and OSPI values from profiles affected with scintillation.

Figure 9 .
Figure 9.Comparison of OSPI obtained at different height intervals with amplitude scintillation indices provided by COSMIC.Left panelshows the agreement relative to the profiles labelled with "naked eye" while right panel shows the agreement relative to OSPI between 550km and 650km (i.e.OPSI @ 600km).

Figure 10 .
Figure 10.One-to-one comparison between the different scintillation indices provided by COSMIC and OSPI.Upper left shows the OSPI vs. S4, upper right shows OSPI vs. SNR of L1 and lower panels show the correlation between the SNR of L1 and S4 in linear (left) and log (right) scale.

Table 1 .
List of scenarios defined for the ionospheric data set, prepared for EUMETSAT's ROPE study The characteristics of each scenario are summarized in Table1and shown in Figure5.Each scenario comprised several days to include a sufficiently high number of COSMIC occultation events.The range of each scenario has been defined so that it Culverwell and Healy (2015)lay using an analytical integration of the Chapman function (i.e. the so-called Zorro equation, seeCulverwell and Healy (2015)).In order to test this new feature, implemented in the context of the ROPE study, different representative scenarios in terms of ionosphere (geomagnetic conditions, solar activity, season,...) have been prepared.10 used to define the scenario is also indicated in Table 1.6 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-217Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 18 July 2017 c Author(s) 2017.CC BY 4.0 License.