Over the years, sporadic E has been measured and monitored using a variety of different instrumentation and methods. This list includes sounding rockets (Hall et al., 1971;Mori & Oyama, 1998;Yamamoto et al., 1998), incoherent scatter radars (ISR) (Christakis et al., 2009;Mathews, 1998), ionosondes (Haldoupis et al., 2006;Merriman et al., 2021;Oikonomou et al., 2014), and GNSS radio occultations (RO) (Arras & Wickert, 2018;Arras et al., 2008;Chu et al., 2014;Gooch et al., 2019;Hocke et al., 2001;Hu et al., 2022;Wu et al., 2005;Zeng & Sokolovskiy, 2010). While sounding rockets are able to provide accurate in situ measurements, they are generally carried out in campaigns which limited spatial and temporal coverage. ISRs and ionosondes are both ground-based instruments that are able to provide E s measurement with good temporal resolution. However, ISR locations are very limited and ionosondes, while having much larger global coverage, are still lacking presence in some regions and over bodies of water. GNSS-RO, on the other hand, uses satellites in Low Earth Orbit (LEO) to receive signals from GNSS satellites in Medium Earth Orbit (MEO) to produce measurements that cover the globe with high vertical resolution. However, whereas ionosondes directly measure E s attributes like intensity (foEs) and height (hEs), these characteristics must be inferred from RO profiles.

In order to derive E s characteristics from RO measurements, data is needed which can be considered "true" and used to validate the inferred values. Since ionosondes can directly measure E s layer intensity and height, and are available at many locations across the globe, they are an ideal choice for the task. It should be noted that the heights given by the ionosondes are actually a virtual height (h′Es). This is a time-based measurement that assumes the signal from the ionosonde travels at the speed of light and reflects off a perfect electrical conductor at the virtual height. In reality, however, the ionization in the atmosphere causes the signal to bend and delay the return of the signal (i.e., h′Es ≤ hEs). At night and during the daytime below the E region maximum, the influence of ionization is fairly weak, and the virtual height is roughly equal to the actual height. However, during the daytime above the E region maximum, differences of a few km are possible (Haldoupis, 2019). Here, we assume that h′Es ≈ hEs.

With inputs (GNSS-RO measurements) and outputs on which to compare against (ionosonde measurements), all that is needed is to develop a model that is able to accurately estimate E s parameters from a given RO profile. In this work, we turn to a method that is currently absent from much of the literature in this area: machine learning (ML). Given a large enough training set under well-suited conditions, ML algorithms are capable of estimating a wide variety of functions very well. Although multiple ML algorithms including random forests and multi-layer perceptron networks were explored to develop models, convolutional neural networks (CNN) achieved the best performance and are the focus of this paper.

CNNs are a type of artificial neural network (ANN) that is traditionally composed of convolutional layers and pooling layers, followed by one or more fully connected layers. The convolutional layers act as feature extractors, where one or more kernels/filters are "convolved" with the input to learn data specific characteristics. The output from the convolution process, often called feature maps, are then sent through an activation function to increase non-linearity. These activated feature maps are then generally sent to a pooling layer to downsample the output in order to reduce the network computation costs while also make the model more robust to variations in the position of the features in the input. The cycle of convolution, activation, and pooling may be repeated several times. The learned feature maps generally represent lower (higher) level features in the earlier (later) stages. At the end of the convolution and pooling layers, the final feature maps are flattened and sent through one or more fully connected layers to learn the complex relationships between the features.

During training, the input is sent through the network and an output produced. This output is then compared with the known target variable and a loss, or cost, is calculated. For regression tasks, mean squared error or mean absolute error are common loss functions, while cross entropy loss is commonly used for classification tasks. The backpropagation algorithm is then used to calculate gradients for all the weights and biases in the network with respect to the loss function (Rumelhart et al., 1986). Optimizers such as Adam (Kingma & Lei Ba, 2015) can then be used to optimally update the weights in an effort to reduce the loss. During prediction, the input is sent through the network of optimized weights and biases to make a prediction.

In order to train a supervised learning model, input and target variable data sets are needed. For inputs, signal-to-noise ratio (SNR) and excess phase profiles from six Global Navigation Satellite System (GNSS) radio occultation (RO) missions during the years of 2008-2022 are used. These missions along with the years of data

ELLIS ET AL.

3 of 19 used, number of measurements considered, and average maximum altitude of the occultation measurement are listed in Table 1.

In processing these profiles, only altitudes of 80-150 km are considered. The profile must also have a maximum altitude of at least 110 km. The profiles filtered based on these criteria are then interpolated to a constant spacing of 100 m using cubic interpolation. In order to normalize the SNR measurements, which may have different biases due to satellite geometries or sensors, each profile first has its respective mean subtracted and then the signal is divided by this mean. For excess phase measurements, the profiles are detrended by subtracting a linear least-squares fit of the data in 2 km windows. Any profiles that do not span 80-150 km are zero-padded to produce 701-point samples that are equally spaced by 100 m. This ensures a constant sized sample as input to the CNN. An example of the raw and processed RO measurements are shown in Figure 1. When E s is present (top panels), a large perturbation in the signal can be seen at the height of the E s layer. For target variables, E s intensity (foEs) and the height (hEs) values are obtained from the global network of ground-based Digisonde ionosondes running the latest Automatic Real-Time Ionogram Scaler with True Height (ARTIST) 5 scaling software. If manually scaled parameters are available, they are also considered. For auto-scaled parameters, only events with confidence scores of 40 or greater are used. If E s is not detected in a sample, the foEs and hEs values are set to zero. This allows for a binary classification (i.e., zero vs. non-zero foEs) that contributes to the overall model performance. In total, these criteria yielded 18.8 million ionospheric sounding events between the years of 2007-2022 across 61 locations. The distributions of these ionosonde samples are shown in Figure 2.

It should be noted that auto-scaled parameters may contain errors and biases due to the way that ARTIST5 extracts information. For example, in Stankov et al. (2023), an error analysis was carried out which compared ARTIST5 and manually scaled ionosonde parameters obtained from the DB049 station in Dourbes, Belgium during the period of 2011-2017. Here it was shown that for hourly foEs parameters at a standard sampling cadence, 67.74% of the time there was a perfect match between the scaling methods and 95% of the residual errors were found between -0.80 and +0.35 MHz, where a negative value indicates that ARTIST5 underestimated the true value and a positive value represents an overestimate by ARTIST5. For hourly h′Es values, the error bounds were -18 to +16 km with a perfect match 57.78% of the time. While we do look at model performance on a small set of manually scaled data, a more comprehensive analysis across all locations is recommended for future work. Additionally, since ARTIST5 confidence scores are primarily concerned with the quality of the overall profile inversion, thresholding based on these values may be of limited value for E s purposes. For this work, we note that the thresholding did not significantly change the E s parameter distributions (e.g., E s occurrence rate changes by less than 1%).

After filtering the RO and ionosonde data sets individually, the two sets are joined together using a 150-km and 30-min spatio-temporal window. If a RO profile is matched with multiple ionosonde soundings, only the closest in time is used. This process left a total of 36,521 samples available for training, validating, and testing the models. The breakdown of these samples based on ionosonde location and years is available in Supporting Information S1 (see Table S1).

Convolutional neural networks (CNNs) are used to develop both binary classification and regression models. The binary classification model is used to determine if E s is present (i.e., foEs > 0), while the regression models are used to estimate values for the intensity and height of the layers. CNNs were chosen as they are able to extract features contained in data. As the RO measurements inherently contain data points that are related to those around it, it is expected that CNNs will perform better than other model architectures.

As CNNs are commonly used with multiple input channels, both the 701-point L 1 SNR and excess phase profiles are included as inputs to the model. The L 2 profiles were excluded due to noise and more inconsistencies with the data. A flowchart describing the start to finish prediction of E s intensity and height from RO profiles is illustrated in Figure 3. The classification step here is used in order to remove some of the cases where there is no sporadic E to avoid a large class imbalance while performing the regression tasks.

For [3][4][5][6][7], the convolutional layer kernel size , the number of kernels for the first convolutional layer (doubles each consecutive layer) [8,16,32], the pooling layer kernel size [2,3], and the type of pooling layer [MaxPool1d,AvgPool1d]. The optimized parameters for each model can be found in Supporting Information S1.

Class imbalance and anomaly detection is explored using cost-function weighting. In binary classification tasks where there are just two classes, a simple inverse balanced weighting is used. This makes the weight of each class proportional to the inverse of the number of its samples (i.e., the minority class will have larger weights). When these weights are multiplied by the cost function, the loss associated with minority class samples will increase.

As the network will attempt to learn parameters that reduce the loss, the performance of the minority cases should improve. A similar approach is taken with rare cases in regression. However, in regression the number of classes can be considered infinite since we are dealing with continuous values. To deal with this, the continuous target variable space is first divided into a discrete number of bins. This distribution is then smoothed using a symmetric kernel to help capture the relationship between nearby bins. Finally, that distribution is used on a sample basis to weight the loss function. Both an inverse of the distribution and the square root of the inverse are explored. An example of these weightings versus the foEs distribution is shown in Figure 4.

Regularization methods aimed to help generalize models are carried out in the form out dropout regularization and early stopping. In dropout regularization, nodes and their connections in fully connected layers are dropped from the network with a certain probability, p, during training. This forces the network to learn a more generalized and robust representation of the data, since it can not rely on any specific nodes to be available. In early stopping, the validation loss is monitored and training is stopped once the loss does not improve over a certain number of epochs (i.e., the number of iterations through a full cycle of the training data). This prevents the network from overfitting to the training data. A learning rate scheduler is also implemented to fine tune performance as the rate of model learning during training slows down.

The binary classification model is primarily evaluated by the F1-score, which is the harmonic mean of recall and precision and can be expressed as:

The F1-score is generally a good performance indicator for data sets with a class imbalance, as it considers both false positives and false negatives. However, it is not a perfect metric and the best metric will be dependent on where one would like the model to succeed (or fail) in a given situation. Accuracy, precision, and recall metrics are also included as additional data points.

For regression, the root mean squared error (RMSE) is used to determine the best model during hyperparameter tuning. The mean absolute error (MAE) is used when comparing models, as it is the most interpretable. Other metrics such as the relative mean absolute error (RMAE), bias, correlation coefficient (r), and coefficient of determination (R 2 ) are also provided as additional performance indicators.

In addition to analyzing the ML models, it is useful to compare the models with others found in literature. For this comparison, we use the methods from Yu et al. (2020), Gooch et al. (2019), andHu et al. (2022) to make predictions using the same testing set as the CNN models.

In Yu et al. (2020), S 4, max values are calculated from the COSMIC-I S 4 profiles over the years of 2006-2014. These S 4 values are calculated using the convention described in Briggs and Parkin (1963):

where I is the square of the SNR, 〈 〉 denotes a one-second time average, and 𝐴𝐴 ⟨𝐼𝐼⟩ is average intensity at each second calculated by applying a low-pass temporal filter in 〈I〉 (Yu et al., 2022). S 4, max is then the maximum S 4 value in the profile. The S 4, max values are then compared with foEs values at 25 ground-based ionosonde stations to develop the following relationship:

(3) 2022), the S 4 profiles are calculated from the L 1 SNR profiles retrieved from atmPhs and conPhs files instead of using the S 4 scintillation profiles provided directly from CDAAC (Syndergaard, 2006). Additionally, for this work, rolling variance and mean profiles of SNR 2 are calculated over a one-second window. S 4 is then taken as the square-root of the rolling variance profile divided by the rolling mean profile. These deviations may cause slight differences in the calculated S 4 values versus the profiles obtained from CDAAC.

In Gooch et al. (2019), the total electron content (TEC) profile is calculated using the L 1 and L 2 excess phase, which is then converted into an electron density assuming a constant E s layer thickness. The raw TEC, TEC r , is calculated using (Hocke et al., 2001):

where f 1 = 1.57542 GHz and f 2 = 1.22760 GHz are the L 1 and L 2 GPS carrier frequencies and ΔL 1 and ΔL 2 are the excess phase measurements in meters. The TEC r is then filtered using a 3rd order Savitzky-Golay filter (Savitzky & Golay, 1964) over a 25 km window. A detrended TEC, TEC d , is then calculated by subtracting the filtered TEC from TEC r and applying a 3rd order Savitzky-Golay filter with a 1 km window to reduce noise. The TEC perturbation of the layer, ΔTEC is then calculated from TEC r using the maximum and base altitudes of the TEC d profile. An electron density perturbation can then be calculated as:

where 𝐴𝐴 2 √ 2𝑅𝑅Δ𝑅𝑅 ≈ 176 km is the effective path length through the E s layer, assuming its geometry is that of a cylinder centered around the tangent point (Ahmad, 1999). foEs can then be calculated as

and hEs is the altitude of the max perturbation of TEC d .

Finally (Hu et al., 2022) uses the TEC profiles from the COSMIC-I and Fengyun-3C over the years of 2006-2019 to calculate foEs. This is done using the singular spectrum analysis (SSA) method to extract the TEC disturbed from the E s layer, TEC d . The S index is then calculated as the vertical gradient of the TEC d :

S max is then taken as the max of the S index, and the E s layer height is the altitude at which S max occurs.

In the remaining sections, the method from Yu et al. (2020) will be referred to as S 4, max , the method from Gooch et al. (2019) as TEC, and the method from Hu et al. (2022) as S max .

While foEs has generally been used as a measure of E s intensity, Haldoupis (2019) points out that the value provided by the Digisondes are in fact the sum of the normal E layer electron density along with the E s layer composed of metallic ions. A new value, foμEs, is then recommended as a measure of E s that better reflects the intensity of the metallic ion layer by removing the background contribution from the normal E region ionization. This is accomplished by first calculating the plasma density (NmE s ) from the plasma critical frequency (foEs) using

Next, the electron density of the background E region (N E ) is calculated from the International Reference Ionosphere model (IRI-2016) (Bilitza et al., 2017) by inputting the height (hEs) and time information associated with the foEs measurement. Finally, foμEs is calculated using the difference between the metallic ion density and the background E region density:

It should be noted that the use of values from the IRI may introduce additional uncertainties that aren't physical, as we are mixing measured values with those from an empirical model. Furthermore, the use of ionosonde derived virtual height as hEs can overestimate the actual height of the layer as discussed in Section 1, which introduces an additional altitude uncertainty into the conversion.

Using this method, foμEs was calculated for the entire ionosonde data set. Figure 5 shows the distribution of the foμEs values along with its associated height (which will be referred to as hμEs for naming convention consistency). Comparing with the foEs distribution from Figure 2, we see that there is no longer a gap between a zero frequency and the minimal detectable frequency by the ionosonde. Maybe not as apparent from the figures is that foμEs is non-zero only 22.7% of the time, roughly half of the time that foEs is non-zero. This creates a situation where the classes are more imbalanced, which could lead to increased difficulty for the CNN classification model to learn a good estimating function. Also, since the μ-variant parameters do have altered distributions, models that learn from these data sets may have some performance metrics that are consequently biased. This makes direct comparisons between models a little more difficult, and some metrics (e.g., correlations and relative errors) will likely be more useful than others.

In Figure 6, E s intensity is viewed as a function of hour (in local time) and solar elevation angle. In plots (a) and (c), we see that the global foμEs occurrence rates, and to a lesser extent foEs occurrence rates, using the ionosonde data set have peaks at around 0900 and 1600 local time (LT), which is encouraging since both diurnal and semi-diurnal behaviors are expected (Hodos et al., 2022;Yu et al., 2022). In plot (e), the geometry of the altitude-adjusted solar elevation angle (α) is depicted. If no E s is present in the ionosonde sample, α is calculated using an altitude of 110 km. In total, 59% of the ionosonde samples corresponded to α ≥ 0 (i.e., "day" condition at E s heights). As seen in plot (b), there is no foμEs seen at negative altitude adjusted solar elevation angles (i.e., "night"), whereas plot (d) shows roughly 18% of foEs occurring when α is negative. While much of this may be contributed to atmospheric tides, it also raises the question of the importance of photoionization in the enhancement of the layers. Alternatively, elevated background E-region electron density estimates from IRI may artificially remove the weaker nighttime E s layers, such that layers with foμEs ≥ 0 MHz are only present when illuminated by solar radiation.

We first discuss the binary classification task that determines if E s is present (i.e., the model should return "1" if foEs or foμEs is greater than zero, and "0" otherwise). As the non-ML models only consider intensities greater than zero, here we look at just the CNN models. We analyze models with and without weighting of the cost function, as previously described. Table 2 lists the performance of the models on the final evaluation set that was unseen during training.

Here we see that the non-weighted foμEs has the highest accuracy, or fraction of predictions that are correct. However, this is to be expected, as a foμEs model would achieve almost 80% accuracy by always predicting the case of no E s . When comparing the weighted versus non-weighted models, it is apparent that the cost function weighting increases the F1-score and recall measures at the cost of decreasing the accuracy and precision. This is contributed to the model placing more emphasis on detecting the positive cases, which may introduce additional false positives. Finally, the foEs models performed better than foμEs on the F1-score metric, likely due to a more favorable class balance. Going forward, only the weighted versions of the models will be used, as they performed better on the F1-score.

It should be noted that even with a somewhat low recall for the models, roughly 90% of positive cases with intensities over 3 MHz were detected for every model except the non-weighted foμEs model, which only detected about 50% of those cases. This suggests that most of the errors came from low-intensity layers, which may not be as important in a practical setting.

For the analysis of the height of the E s layers, we start by comparing the performance of each model on the final evaluation set when it is known that E s is present from the ionosonde data (summarized in Table 3). In each of the comparisons, the prediction of the model is compared against the height value produced by the ionosonde. It can be seen that the CNN models perform similarly and do a better job of fitting to the ionosonde data, achieving lower RMSE, MAE, RMAE, and bias values along with better R 2 and r scores, compared to other models found in literature. Of particular note is that the non-ML models have a tendency to underpredict the ionosonde height. This may be at least partially explained by the difference in virtual height versus actual height of the E s layers, as these models use more of a physical indicator whereas the ML models learn what is present in the ionosonde data. Additionally, ionospheric irregularities encountered along the RO path before or after the tangent point are assigned to the tangent point altitude, which introduces an inherent geometry based altitude uncertainty from the RO measurements (Wickert et al., 2004) that generally underestimates the E s altitudes (Gooch et al., 2019). Furthermore, using ARTIST5 scaled parameters will also add some uncertainty in the model performance as there are errors in the way it extracts height values.

Figure 7 displays normalized heatmaps of the predicted versus true E s height values (as given by the ionosonde) for each of the models. Here we see that the CNN models have a tighter prediction window, but struggle to predict the heights of the E s layers at the extreme cases from 80 to 90 km and 120-140 km. This is likely due to the rarity of these high-hE s or high-hμE s cases. Class imbalance techniques could be implemented to improve prediction of these rare cases, but are not done here as most of the extreme height cases corresponded to lower intensity E s layers, which may not be as important. The non-ML models on the other hand have a wider prediction window but also tend to underestimate the ionosonde virtual height as previously mentioned.

For the intensity analysis, we again look at the performance of each model on the testing set when there is E s present in the ionosonde data (summarized in Table 4).

When comparing the CNN (foEs), S 4, max , TEC, and S max models, the foEs values retrieved from the ionosonde were used as the ground-truth target variables. For the CNN (foμEs) model, the foμEs values calculated using the method in Haldoupis (2019) were considered the true values. The CNN models generally performed better, with the foEs model performing better on the error metrics and the foμEs model on the bias, R 2 , and r metrics. This may be due to CNN models' ability to fit to a wide range of non-linearities in order to more accurately express the relationship between the RO and ionosonde data. For the non-ML models, it is interesting to note that while the correlation coefficients are reasonable, the R 2 and bias values are much worse. This likely points to the physically derived parameters such as S 4, max and S max being good indicators of E s intensity, but the model coefficients (e.g., slope and intercept) do not fit the data set used in this work well. This may lead to models with low variance, but high bias. Potential sources for these errors when making comparisons to these other models include calculation of the S 4 /TEC profiles versus retrieving them directly from CDAAC, ionosondes with different parameter scaling or locations, and different spatio-temporal windows when joining the RO and ionosonde data. Additionally, while there are errors in the ARTIST5 ionosonde parameter scaling, the error bounds as shown in Stankov et al. (2023) tend to be much tighter than the error bounds produced by the models. For example, the 95% residual error bounds of the CNN foEs model were -1.64 to +2.44 MHz, while the bounds of the residual error for the manual versus ARTIST5 scaling at the DB049 station was shown to be -0.80 to +0.35 MHz. However, it also should be noted that these ionosonde scaling errors may differ by location and is an area that still needs further exploration.

The normalized heatmaps in Figure 8 reinforce the results in 4 and show that the CNN models fit tighter to the line of optimum fit, whereas the non-ML models have more variability in their predictions. The large bias in the S 4, max model is apparent in the predicted foEs values, although the model does show a high density cluster, indicative of good precision. This bias is likely the result of using different ionosondes with different parameter scaling.

It can also be seen in the heatmaps that the CNN models tend to underpredict the high intensity E s layers. This is again likely due to the underrepresentation of these rare cases in the training set. Additionally, there is a physical restraint on S 4 amplitudes from E s induced scintillation (Emmons et al., 2022;Stambovsky et al., 2021), which may be a constraint in the CNN models for the stronger layers. In Figure 9, the inverse (INV) and the square root of the inverse (SQRTINV) cost function weighting schemes are applied to the foEs CNN model. The inverse weighting is much more aggressive and, while preforming the best on the highest frequencies, does so at the cost of the lower intensities. The SQRTINV weighting, on the other hand, strikes a balance between doing well on the rare cases while still maintaining respectable performance on the more common cases. The non-weighted model will be used in the next section, as the overall performance was better. However, the "best" model will be situationally dependent and based on if it is more important to get the rare cases correct, or perform better on average.

Section 4.1 discussed the errors of the binary classification of whether or not E s is present, whereas Section 4.3 quantified the errors for estimating the intensity. In a practical application, the errors of the two will compound since it is not generally known in advance whether E s is present, leading to the overall model architecture depicted in Figure 3. We now evaluate the overall error.

As can be seen in Figure 10, the general trend is that the error for every model increases with increased predicted intensities. Both CNN models perform well in the 0-2 MHz as expected due to the classification stage. Overall, the In Figure 11b, it isn't readily apparent if there is any correlation between the average predicted foE and the solar cycle (using the yearly averaged sunspot number as a measurement proxy). However, we do find a positive correlation of r = 0.39. Although not shown in the figure, when the data is split into day and night at the location of the E s layer, we find a positive correlation (r = 0.56) during the day, and almost no correlation at night (r = -0.01). This result is reasonable as during the daytime, foEs has a higher contribution from the background E layer which is highly solar driven. The number of samples in the daytime is also greater than those at night, so we expect to see a bias toward the daytime correlation value. These finding are in line with those reported in Zhang et al. (2015) and Zuo and Wan (2002), which noted a positive correlation for foEs during the daytime and a negative correlation at night for low and mid-latitude locations. However, when splitting the predicted foEs data into high (sunspot number ≥ 50) and low solar activity (sunspot number < 50), we find correlations coefficients of r = 0.46 and r = -0.32, respectively, for day and night. This could be due to the nature of RO profiles, which measure the signal perturbation from the background environment (Haldoupis et al., 2020). In times of high solar activity, background E layers may be enhanced and thus cause difficulties in accurately measuring and processing the relatively less intense perturbed signal (Hodos et al., 2022).

In Figures 11a and11c solar driven diurnal trend is apparent, with increased intensity in the daytime hours. In plot (d), we see that E s intensity peaks in the boreal summer (June-August) and reaches its minimum in the winter (December-February). Similarly, in plot (d), we see E s intensity peaking in the austral summer (December-February) and reaching its minimum in winter (June-August). Although in the Southern Hemisphere, the intensities do not get quite as high as the Northern Hemisphere, which has also been noticed in Hu et al. (2022) and Niu et al. (2019). Both the solar driven diurnal behavior and maximum intensity during summer months of the respective hemisphere agree with results reported in literature (Hodos et al., 2022;Yu et al., 2022). It should be noted that in the plots when the average intensity is low, the relative error is high. This can be attributed to factors such as limitations on the minimum detectable intensity of the ionosondes (∼1.5 MHz), the high number of cases of no E s , and errors introduced by incorrect classifications.

In Figure 12, examples of the predicted intensity and height from the CNN model are illustrated using ionograms along with corresponding radio occultation inputs. In (a), a strong E s layer is present on the ionogram while a large perturbation is also present in the input RO signal. In (b), a less intense layer is present on the ionogram while the RO inputs are much weaker and feature more noise-like perturbations. Finally, in (c), there is no E s layer seen on the ionogram while again the RO inputs are weak and lack any well-defined large perturbations. In all of these situations, the CNN model was able to predict the presence, intensity, and height of the E s layer with good accuracy.

Figure 13 shows similar trends for the foμEs model, though as expected, the average predicted values are lower since the foμEs occurrence rates and average intensities are lower than foEs. However, when calculating the correlation coefficient between the predicted foμEs and yearly sunspot numbers, we see only a slight correlation of r = 0.09. This is reasonable, as the foμEs values attempt to remove the intensity contributions from the solar driven background E. In general, the error for the foμEs and foEs models tend to increase (decrease) as the predicted E s intensity increases (decreases).

Although six different RO missions were used, most prediction errors were similar across the constellations. Table 5 lists the average MAE for foEs and hEs predictions across all models. The relative performance for each constellation when compared with the average total error across all constellations is listed in parentheses. In both intensity and height predictions, it can be seen that the GeoOptics satellites fared the worst with relative errors of 151% and 118% for foEs and hEs predictions, respectively.

Additionally, in Section 2 it was noted that the ionosonde data set is primarily obtained ARTIST5 auto-scaling software, with a small portion of the data set scaled manually. Again, ideally all the data would be manually scaled as there may be errors and biases in the way ARTIST5 extracts the parameters (Stankov et al., 2023). However, with large amounts of data, that can become impracticable. In Table 6, the relative errors for both manual and ARTIST5 scaled data are listed for each model. Here we find that on average, the relative predicted foEs errors on the manually scaled data are slightly higher than those of the ARTIST5 scaled data, while the hEs values are very close (note that both foEs manual and ARTIST5 scaled height values may overpredict the actual heights as they represent a virtual height). These differences may in part be due to the fact that the manually scaled foEs values are slightly higher than those that are scaled by ARTIST5, while the hEs values are very close (<1 km on average). Since each model tends to perform worse at predicting foEs on more intense layers, an increase in error is reasonable. However, there are many variables that come into play and a more comprehensive study of manual versus auto-scaled parameters across all ionosonde locations is still recommended for future work.

In this work, CNN models were developed to detect E s along with estimating its height and intensity. These models were trained using L 1 SNR and excess phase profiles from GNSS-RO missions as inputs, along with foEs and hEs values provided by ionosondes as the target variables. In total, 36,521 total samples across 45 locations and ranging from the years 2008-2022 were available after pre-processing for training and testing. Models were also explored using E s intensity values that attempted to remove the background E layer ionization contribution from foEs values (i.e., foμEs). Each of the models performed well, although the foμEs CNN classification model struggled more due to a class imbalance and additional uncertainties introduced by the foμEs conversion process, which relies on IRI background E-region electron density estimates. However, it was shown that with different loss function weighting schemes, improvements can be made. The CNN approach also performed favorably when compared with other models found in literature, achieving lower errors when predicting both the E s intensity and height. The E s intensity error attained by the model generally grew with the strength of the layers, as these cases occurred less often during training. This may be compensated for using techniques to combat class imbalance, but at a potential cost of worsening performance of the commonly occurring cases. It is also noted that while the CNN models performed well on height prediction, there are uncertainties associated with the virtual heights calculated by the ionosondes and the actual E s layer height. This may cause the ionosonde to overestimate the height of the E s layer. Additionally, errors in parameters extracted by ARTIST5 may impact model performance, and a more comprehensive analysis of scaling methods across multiple ionosonde locations is recommended. In the end, however, the relatively low error in the models makes them a good choice for applying to large RO data sets for E s climatology purposes, or even near real-time detection if enough RO sensors become available.