Space Weather Services for Civil Aviation—Challenges and Solutions

: This paper presents a review on the PECASUS service, which provides advisories on enhanced space weather activity for civil aviation. The advisories are tailored according to the Standards and Recommended Practices of the International Civil Aviation Organization (ICAO). Advisories are disseminated in three impact areas: radiation levels at flight altitudes, GNSS-based navigation and positioning, and HF communication. The review, which is based on the experiences of the authors from two years of running pilot ICAO services, describes empirical models behind PECASUS products and lists ground-and space-based sensors, providing inputs for the models and 24/7 manual monitoring activities. As a concrete example of PECASUS performance, its products for a post-storm ionospheric F2-layer depression event are analyzed in more detail. As PECASUS models are particularly tailored to describe F2-layer thinning, they reproduce observations more accurately than the International Reference Ionosphere model (IRI(STORM)), but, on the other hand, it is recognized that the service performance is much affected by the coverage of its input data. Therefore, more efforts will be directed toward systematic measuring of the availability, timeliness and quality of the data provision in the next steps of the service development.

the GNSS positioning performance [14]. The statutory procedures of aircrew radiation risk monitoring should be supported with estimates of effective dose rates at a multitude of different flight levels. Space weather-related problems in HF-communication can be triggered by four different types of events: absorption at auroral latitudes, polar cap absorption (PCA), solar X-ray flares and post-storm depression. The three first events are associated with electron density enhancements in the ionospheric D-layer, causing HF signal fadings. In the fourth case, the capability of the ionospheric F-layer to reflect radio signals decreases, due to significant drops in its electron density. The parameters to be followed are the level of global geomagnetic activity (for auroral absorption) as characterized by the Kp-index, attenuation in cosmic noise in the HF regime (for PCA), solar flux in the range of 0.1-0.8 nm (for X-ray flares), and maximum usable frequency (MUF, for post-storm depression).  30 50 Information on enhanced space weather activity causing potential problems in the above described impact areas is disseminated by advisories with a standardized structure. An example of such an advisory in its basic form is given in Figure 1. Besides some standard information (time stamp, service provider, advisory numbering, etc.), the advisory describes the space weather impact by its type, intensity and spatial extent. At the time of this writing, the extent is characterized by using a geographic grid with a resolution of 30°in latitude and 15°in longitude (hereafter called the ICAO standard grid). The advisory template includes also rows for forecasts with lead times of 6, 12, 18 and 24 h. Three global space weather centers were found to be compliant with the ICAO specifications and were asked to start official operations in November 2019. PECASUS is one of these three global centers, and is operated as a joint effort of research institutes from ten ICAO member countries. In this article, we provide a comprehensive review on the instrumentation and data analysis methodologies that are used in PECASUS in order to ensure timely issuance of validated ICAO advisories. The review is based on the authors' experiences gathered during years 2019-2021, when the PECASUS system was run in its pilot phase for ICAO. The core instruments and the data analysis methods behind PECASUS services are described in Section 2. Sections 3 and 4 provide short overviews of the PECASUS operational and self-evaluation activities. Section 5 showcases PECASUS performance during a post-storm depression situation, where the ICAO advisory thresholds were temporarily exceeded. Finally, Section 6 summarizes the lessons learned by the PECASUS team during its first years of operation and presents some prospects for future improvements in ICAO space weather services.

Observations and Methods
In the following subsections, the PECASUS approaches to monitor space weather activity in the three impact areas of ICAO interest-GNSS, radiation, and HF communicationare described. For each impact area, we first give some general information about space plasma processes that are behind the harmful impacts, and about some standard methods used in the space weather community to observe those processes. Subsequently, we describe the algorithms that PECASUS uses in its near-real-time services.

GNSS
In the impact area of GNSS, ICAO advisories are disseminated in two subcategories: in the cases of enhanced scintillation (Section 2.1.1) and of enhanced TEC (Section 2.1.2).

Scintillation
Small-scale irregularities in the ionospheric electron density cause scintillation in GNSS signals. During space weather storms, enhanced signal phase or amplitude scintillation can complicate the signal reception, e.g., by loss-of-lock events [15]. At high latitudes, phase scintillation is more typical than amplitude scintillation, and the disturbances are often related to auroras and so-called polar cap patches with localized gradients in electron density [16]. At equatorial latitudes, plasma bubbles building at local sunset cause amplitude scintillation [17]. Scintillation is monitored with specific GNSS receivers that have higher sampling rates than the standard receivers. As the coverage of scintillation, receivers are limited particularly in global scales; actual scintillation measurements are often supported with a proxy called rate of TEC (ROT) derived from networks of standard GNSS receivers with better coverage. The ROT index (ROTI) is retrieved per each receiversatellite link, according to the formula of [18]: ROTI = √ < ROT 2 > − < ROT > 2 in which ROT = dTEC/dt and indicates the temporal average calculated, e.g., over 5 min. As time derivatives are used, the challenging task of robust TEC calibration for absolute values can be avoided.
In PECASUS services, scintillation is monitored with high-cadence GNSS receivers (up to 50 Hz) whose carrier signal at L1 frequency (1575.42 MHz) is used to derive standard deviations in one-minute time windows. The amplitude scintillation index S 4 is normalized with the signal amplitude, and the phase scintillation index σ φ is derived from the detrended signal and given in radians. The geographic locations of the PECASUS scintillation receiver stations are shown in Figure 2.

Total Electron Content
Besides scintillation, the ionosphere affects GNSS signal propagation by its bulk electron content. The signal gets delayed, which causes some error in its range used as input in the positioning procedures. In the first approximation, the range error (d) is directly proportional to the integrated electron density, along the signal path in the ionosphere, i.e., d = K f −2 N e (s)ds = K f −2 STEC, where K is a constant 40.3 m 3 s −2 , f is the signal frequency, N e (s) is the electron density, and s defines the signal path. STEC, i.e., slant total electron content, can be converted to the vertical total electron content (VTEC) with the following formula, where the ionosphere is described with a single layer approximation [14]: where R e is the Earth's radius, is the signal elevation angle, and h sp is the altitude of the single-layer ionosphere (typically 400 km). A common standard is to express VTEC values in TEC units (TECUs), with 1 TECU corresponding to 10 16 m −2 (i.e., 16.2 cm in the range error for GPS L1 signal). While ionospheric delay can be reduced from services using dual-frequency GNSS receivers, applications relying on single frequency information need additional support from ionospheric models for adequate positioning performance. Such applications are still used in aviation and, therefore, estimates on global N e are critical to enable safe traffic management and landing procedures. Statistical ionospheric models of volumetric N e , such as NeQuick [19] or VTEC, such as Klobuchar [20], provide convenient ways to provide such information as being embedded in the GNSS signals. It is obvious, however, that during space weather storms, range errors cannot be estimated reliably with statistical models. For improved error estimates, statistical models can be supported by TEC information from ground-based networks of dual-frequency GNSS receivers. For redundancy reasons, PECASUS hosts two parallel systems for global TEC monitoring: one developed under the German Ionosphere Monitoring and Prediction Center (IMPC; https://impc.dlr.de/, accessed on 27 August 2021) and the other developed under the Italian electronic space weather upper atmosphere system (eSWua; https: //doi.org/10.13127/eswua/tec, accessed on 16 September 2021). Both systems use as inputs freely accessible near-real-time (NRT) GNSS receiver data from global networks, such as EUREF (https://www.epncb.oma.be/, accessed on 27 August 2021) and the International GNSS Service (IGS, https://www.igs.org/, accessed on 27 August 2021). The resulting TEC maps have spatial resolutions of 2.5°in latitude and 5.0°in longitude and are updated once every 5 min by IMPC and once every 15 min by eSWua. Besides global maps, both systems also generate higher resolution TEC maps for the European region. These maps have a time resolution of 5 min, and the spatial resolution is 2°in latitude and longitude in the IMPC version and 0.5°in the eSWua version. For fine-scale investigations, eSWua is able to produce TEC maps over the Mediterranean territory with a spatial resolution of 0.1°and with a time resolution of 10 min. The exceptionally high resolution is achieved by the combination of IGS and RING (Rete Integrata Nazionale GPS; http://ring.gm.ingv.it/, accessed on 27 August 2021) GNSS receiver networks. Details about STEC calibration (removing hardware based biases and ambiguity in the carrier phase observable) for the IMPC TEC maps can be found in [14] and for the eSWua maps in [21][22][23].
In the IMPC approach, calibrated STEC values are converted to VTECs with Equation (1) and are combined with the empirical Neustrelitz TEC model (NTCM) [24] describing average TEC conditions. NTCM represents the spatio-temporal variations of TEC based on a polynomial consisting of linear terms corresponding to the diurnal and semi-diurnal variation, the annual and semi-annual variation, dependence on geomagnetic latitude, solar zenith angle and dependence on solar activity coefficients, multiplied by factors defined by linear least squares with high-precision data from the IGS long-term archives. In NRT operations, a streamlined version of the model is used that ignores seasonal and solar cycle dependencies, but allows rapid model adjustment according to the current day's GNSS data and thus, generation of TEC maps in near real time. For this purpose, the NTCM model coefficients and the corresponding TEC calibration of all input data streams are performed on the basis of a weighted least squares solution every 24 h, using the complete data set of the previous day, recursively adjusted to the TEC observations of the last 5 min.
Since ground-based GNSS data are often unevenly distributed, the inclusion of the NTCM in the TEC reconstruction helps to overcome such data gaps, which occur mainly over the oceans. Therefore, the assimilation uses observed VTEC values at their ionospheric piercing points to further fit the TEC values of the NTCM model with two-dimensional Gaussian weighting factors, depending on the distances between the piercing points and the model grid points. The final results represent measured TEC values near the measurement points, while slightly modified model values are provided at larger distances from the measurements.
In the eSWua approach, VTEC values at their piercing points are combined with the NeQuick model [19]. As the first step, data ingestion includes the calculation of differences between the observed VTECs and those by NeQuick at the corresponding piercing points. Subsequently, the differences are interpolated by kriging in order to create a global map of differences that is matching with the observations. Finally, the TEC nowcast is achieved by summing the map of differences with the TEC map by NeQuick. The second step of this process is under development, as standard kriging methods assume an isotropic distribution of data points, which obviously is not true for GNSS receiver networks. This deficiency can be partly overcome by handling differences instead of absolute TEC values, but as a further improvement the algorithm will, in its next release, manage the anisotropy by creating a variogram surface (for more details, see [25]).

Radiation
Earth's atmosphere and magnetic field protect us from the detrimental impacts of galactic cosmic rays (GCR) and solar proton events (SPEs). The shielding, however, is partial. According to the Lorentz force law, energetic charged particles propagating toward the Earth get deviated from their original trajectory by the geomagnetic field component parallel to the Earth's surface. As this component gets smaller at higher latitudes, also the shielding by the geomagnetic field becomes weaker there. Therefore, assessment of effective dose rates for cross polar flights is of particular interest to airliners. The shielding effect by the atmosphere occurs due to collisions that slow down penetrating particles, but on the other hand, they also generate showers of secondary particles that can enhance radiation levels at flight altitudes. Secondary showers are ubiquitous in Earth's atmosphere, due to the always-present GCR. Modeling the atmospheric radiation level due to GCR is easier than estimating the impact by SPEs because variations in GCR are slow, and the spectrum of primary particles is well known. On the other hand, additional radiation by SPEs is characterized by its abrupt appearance, short duration and energy spectra, varying significantly from case to case. In the context of very strong SPEs, secondary showers are detected as ground level enhancements (GLEs), i.e., enhanced count rates in ground-based neutron monitors.
In PECASUS, effective doses (µSievert/hour) at flight altitudes are estimated with the AVIDOS tool [26], which can handle both GCR and SPE impacts on radiation levels. The contribution of GCR on radiation levels at flight altitudes is described elsewhere [27]. For the impact of SPEs, the approach uses, as its primary data source, NRT count rates from the neutron monitor in Oulu (Finland, 65.05°N 25.37°E; [28]) with support by the ANeMoS GLE alerting service [29]. In the case of a GLE, AVIDOS first constructs two estimates for the primary proton energy spectrum-one soft and one hard version. In the method, the energy spectra of the primary protons are assumed to have a power law form in particle rigidity (particle momentum divided by its charge): where R is the primary proton rigidity expressed in GV, t is time, A 0 (t) is amplitude of the proton flux (in particles/(cm 2 · s · sr · GV)), and γ(t) is the spectral index (softness of the spectrum). A 0 (t) is derived from the neutron monitor data. The relationship between the primary spectra and neutron monitor count rates, N(R, t, z), can be described by the following formula [30]: where R min is the minimum considered rigidity (R min = 0.8 GV for Oulu neutron monitor station), z is the atmospheric depth of the station, and S(R, z) is the so-called yield function for the considered neutron monitor. For the yield function, AVIDOS uses the parametrization developed by [31]: where C mn are parametrization coefficients based on Monte Carlo simulations of neutron monitor responses (see [31] for more details and exact values of C mn ) and log stands for decimal logarithm. When the two estimates (hard and soft) of the proton spectrum are ready, the spectra are folded with results of Monte Carlo simulations of radiation transport through the atmosphere performed with the Geant4 application PLANETOCOSMICS ( [32,33]), which results in a large set of secondary particle fluences. Finally, to obtain the radiation protection quantity effective dose, the particle fluences are folded with appropriate fluence-to-dose conversion coefficients and are summed up. In AVIDOS, two types of such conversion coefficients can be used: that based on the older recommendations of the International Commission of Radiological Protection (ICRP) from 1990 [34] or that based on the more recent ones from 2007 [35].
AVIDOS performance was validated in an international comparison campaign performed by EURADOS' (the European Radiation Dosimetry Group, www.eurados.org, accessed on 27 August 2021) Working Group 11 on high energy radiation fields. In this exercise, nine codes estimating radiation exposure due to GLE were employed to calculate radiation doses for GLE 42 (September 1989) and GLE69 (January 2005) for three representative flight routes. The survey showed that estimated flight route doses may differ by up to a factor of ten [36]. Further research and verification by on-board aeroplane measurements are needed for more general level understanding on the model discrepancies. The NRT version of AVIDOS tackles this issue partly by the set-up of two estimates for the primary proton population spectra (hard and soft ones). Consequently, as output, the code yields maximum and minimum boundaries of the expected effective dose. By deducing the minimum-maximum range instead of a single value, the code gives implicit information about its performance (with the assumption that large ranges imply increased uncertainty in the estimates).
Complementary to ground observations, PECASUS also makes use of NRT space data to generate SPE forecasts, which are internally used to alert the operator for increased risks of radiation storms and possible impacts on aviation. Two predictor tools are adopted: HESPERIA UMASEP-500 ( [37,38]) and COMESEP SEP Forecast ( [39]). The HESPERIA UMASEP-500 tool produces real-time predictions of the intensity of the first hour of >500 MeV SPEs near-Earth on the basis of a time-delayed correlation analysis between the GOES soft X-ray flux and the GOES differential and integral proton fluxes in the energy channels above 500 MeV. When the correlation estimation surpasses a threshold, and the associated flare is greater than a specific X-ray peak flux, a >500 MeV SPE warning is issued before (up to tens of minutes) detection by neutron monitors. Particles of that high energy can easily reach the top of the atmosphere within minutes and create enhanced radiation levels at flight altitudes. The COMESEP alert system provides (in addition to geomagnetic storm alerts) SPE alerts in the lower energy range, i.e., >10 and >60 MeV. The SPE alerts are triggered by the automated detection of solar activity (>M1 solar X-ray flares). Based on a statistical analysis of historical data correlating flare and CME observations with SPE characteristics, a risk level estimate is derived from NRT data based on the expected occurrence probability and severity. The alerts are continuously updated along the timely evolution of event(s) in the case that new information on the solar eruption is received. As explained in the next paragraph on HF communication, these solar energetic particles mainly spiral-in into the polar regions, causing polar cap absorption events. In the PECASUS service, the COMESEP system is used as an advanced warning for increased risk of radiation impacts confined to the transpolar regions.

HF Communication
In the impact area of HF communication, ICAO advisories are disseminated in four subcategories: in the cases of shortwave fadeout, polar cap absorption, auroral absorption (all described below in Section 2.3.1) and of post-storm depression (Section 2.3.2).

Shortwave Fadeout, Polar Cap and Auroral Absorption
Solar X-ray flares and energetic particle precipitation increase the electron density in the lower part of the ionosphere (D-layer). While the shortwave fadeout impact by X-ray flares lasts from some minutes to a few hours, an absorption event due to energetic particle precipitation can last for days. In both cases, the interaction between radio waves and D-layer plasma enhances when compared to normal conditions. As this interaction is affected strongly by the ionization of neutral atmospheric constituents and by electron collisions with atmospheric neutral particles, radio waves are absorbed instead of reflected, like in the case of F-layer interaction. D-layer absorption can be a problem particularly in the low end of the HF-spectrum because electron density in the D-layer is too low to affect higher frequencies, even during disturbed conditions. A solar X-ray flare causing D-layer absorption is known as a shortwave fadeout, a legacy term from radio wave applications using frequencies lower than HF. Airliners need to be aware of shortwave fadeouts because they challenge HF-communications in the sunlit ionosphere. At the time of this writing, PECASUS advisories on fadeouts are based on X-ray observations by the XRS instrument (0.1-0.8 nm; [40]) of the geostationary GOES-16 satellite-maintained by the U.S. National Oceanic and Atmospheric Administration (NOAA).
The D-RAP model (https://www.swpc.noaa.gov/content/global-d-region-absorptionprediction-documentation, accessed on 27 August 2021) is widely used in space weather services-including PECASUS-to estimate the spatial extent and frequency impact of D-layer absorption, due to solar X-ray flares and SPEs. As inputs, the model uses X-ray and proton fluxes by the GOES-16 satellite and the Kp index (more details about Kp are given later in this section). For X-ray flares, the model provides an estimate of the highest affected frequency in the dayside ionosphere and spatial extent of the impacted area, the latter of which is used as supporting information in ICAO services. However, the intensity of impact (MOD/SEV) is determined from direct X-ray flux measurements (c.f. Table 1). For SPE-driven polar cap absorption, D-RAP gives estimates of 30 MHz signal attenuation in decibels (dB) based partly on empirical relationships as described by [41]. The relationship between proton flux and day-and nightside absorption levels (in dB) can be expressed as follows: and where A d and A n are the dayside and nightside absorption levels (with sun elevation angles >10°and <−10°, respectively), J(E>Ec(InvLat)) is the proton flux (in (cm 2 · s · sr) −1 ) for energies larger than E c . E c (in MeV) is the cut-off energy above which protons are able to overcome the geomagnetic shield and, thus, precipitate into the ionosphere at the given location, whose position is characterized with the McIlwain L-parameter [42]. In the twilight zone, absorption levels are estimated with linear bi-interpolation from A d and A n according to the sun elevation angle and InvLat. D-RAP performance in the polar cap absorption assessment was evaluated with riometer data from six high-latitude stations at different longitudes for 13 space weather storms during the years 1998-2005 [43]. The survey reveals that the model performs reasonably well for high-latitude stations and can provide useful qualitative information also during very disturbed periods. However, for more accurate quantitative estimates, more work is still needed, as D-RAP seems to under-estimate absorption events, particularly in the European sector.
The D-RAP approach cannot capture energetic electron precipitation of magnetospheric origin, which can cause HF-signal absorption at auroral latitudes. Energy spectra and spatial distribution of electrons is linked with magnetospheric dynamics in such a complicated way that they cannot be characterized using the single-point measurement approach (at geostationary satellites) used for proton fluxes. This complication can be partly overcome with the knowledge that the energetic electron precipitation being at high levels is particularly common during geomagnetic storms. Therefore, ICAO has chosen to use the Kp-index as the indicator of enhanced auroral absorption in its advisories. Kp-index is composed from measurements of 13 well-established geomagnetic observatories located at sub-auroral latitudes. It describes geomagnetic activity as the deviation from the regular daily variation with a three-hour resolution [44] and with a scale from 0 to 9 [45]. Statistical analysis of Kp-archives since 1933 reveals that the occurrence of high enough Kp-values for ICAO advisories on auroral absorption varies, according to the solar cycle being on average seven per year during years around the solar maxima and once per year around the solar minima. Typical durations of geomagnetic storms with peak Kp ≥ 8 are in the range of 1-3 days.

Post-Storm Depression
Geomagnetic storms are associated with post-storm depression of electron density in the ionospheric F2 layer. Depressions are usually more intense, longer-living and widespread during the main and recovery phases of intense geomagnetic storms. According to a common recommendation between the ICAO space weather centers, post-storm depression advisories may be issued up to 96 h after Kp ≥ 6 has been reached. With decreased electron density, the ionosphere's capability to reflect radio waves lowers, in particular, at the higher end of the HF-regime. Post-storm depression is caused by changes in the thermosphere-ionosphere interaction as a consequence of enhanced Joule heating and particle precipitation at auroral latitudes [46]. Heating at high latitudes causes up-welling, where neutral molecules are lifted to typical F2 layer altitudes (300-400 km) and cause the loss of electrons due enhanced recombination. Heating generates also horizontal pressure gradients and advection of molecule-rich plasma from auroral latitudes to lower latitudes. The spatial distribution and intensity of a post-storm depression depends strongly on season and local time. During summer time, the background neutral wind pattern has an equatorial component, which supports equatorward propagation of molecule-rich plasma fronts. Therefore, post-storm depression events penetrate to lower latitudes and are stronger there during summer than during winter months. Furthermore, observers in early local time morning hours at mid-latitudes experience more severe depression effects than those in the afternoon sector. The recovery of ionosphere-thermosphere conditions back to normal may take 1-2 days after a strong post-storm depression event.

Ionosonde Measurements
Post-storm depressions can be detected with ionosondes, which are vertically sounding pulsed radars operated in a range of HF frequencies (typically 1-20 MHz). According to this measurement concept, the ionosphere is probed with gradually increasing frequencies.
Signals transmitted from the ground are reflected back from the ionosphere at the altitude where the ionospheric plasma frequency ( f p , f p ≈ 9 √ N e ) is greater than or equal to the signal frequency. The highest frequency that is reflected back gives a good estimate for the maximum N e in the ionosphere. Approximation of the reflection altitude can be obtained from the signal travel times, with the assumption of the signal propagation speed being that of light. The maximal electron density is typically in the F2 layer. URSI (International Union of Radio Science) specifies the ionospheric parameter codes. f oF2 is used for the F2 layer plasma frequency peak, and h F2 is used for the virtual altitude of f oF2. The latter is the height obtained directly from the measured radio wave travel times; it can be converted to true height hpF2 or hmF2 only by using a theoretical fitting curve [47].
The maximum usable frequency listed in the ICAO advisory threshold parameters (c.f. Table 1) can be related with URSI parameters MUF3000F2 and M3000F2, the former of which is the F2 layer maximum usable frequency for a 3000 km path and the latter of which is MUF3000F2/ f oF2. Obtaining f oF2 from ionosonde measurements is straightforward, as it can be easily scaled directly from the ionogram. Estimating MUF3000F2 or M3000F2 requires either fitting the ionogram traces with a semi-empirical curve [47] or knowledge of the true peak height hmF2. The latter, in turn, needs to be calculated from the measured h F2, using some theoretical model curve to take into account the delay of the signal due to ionization below the reflection altitude.

The PECASUS Approach
In PECASUS, both f oF2 and MUF3000F2 are used to monitor post-storm depression. As pointed out above, f oF2 is available as a standard parameter for all ionosonde types, which makes it a robust source for global monitoring. The network of ionosondes that is used in the PECASUS approach at the time of this writing is presented in Figure 3. f oF2 depression, δ f oF2, is estimated for each available ionosonde as follows: where < f oF2 > is the 30-days' median value of f oF2. Positive values of δ f oF2 correspond to depression of f oF2. f oF2 values from each available ionosonde are checked and compared with the 30-day median. The outputs of the f oF2 algorithm-δ f oF2 values and excesses of the ICAO MOD/SEV thresholds are given in the resolution of the ICAO standard grid (c.f. Section 1) and with a time resolution of 15 min. If, for a specific ionosonde, there are fewer than 20 values for the past 30 days, the 30-day median is estimated, using the NeQuick model [19]. The δ f oF2 value from a particular ionosonde is naturally considered to be representative for its own cell in the ICAO standard grid, but also for the neighboring cells if the central angle between their boundaries and the ionosonde location is less than 10°. As stormy conditions can, in some cases, challenge reliable ionosonde measurements, the strategy here presented is to apply some conservative rules in the δ f oF2 algorithm in order to avoid false warnings. If several ionosondes contribute to δ f oF2 for a given cell and all those values are below the ICAO advisory thresholds, then the median of the individual values is used to represent the cell conditions. If only one of the contributing ionosondes is suggesting excess of a ICAO threshold, that measurement is ignored, and the cell obtains the median value from the remaining ionosonde measurements. If more than one of the contributing ionosondes are giving higher δ f oF2 than the thresholds, then the median is determined from those ionosondes, and the other ionosondes are ignored. A similar rule is applied also when creating the maps of ICAO threshold excesses. A cell is denoted with yellow or red shading only if two or more ionosondes contributing to its content are showing strong enough depressions for MOD or SEV advisory issuances. The color is determined according to the stronger activity level. Examples of δ f oF2 maps created with these principles are shown in Figure 4. The method used to produce MUF3000F2 global maps is composed of two models for f oF2: global dynamic model (GDMF2) [48] and EUROMAP [49]. MUF3000F2 maps are created from f oF2 maps by using propagation factor M3000F2 estimates from the empirical international reference model (IRI; [50]). The space weather impact on MUF is characterized by the ratio between MUF3000F2 and background conditions, which come from the empirical model by [51]. Maps are created with a one hour time resolution, based on NRT data and indices. As many of the indices are available also as forecasts, the methodologies described below are suitable for forecasts with some moderate modifications.
GDMF2 is a global empirical model particularly for mid-latitude post-storm depression conditions, which relates ∆ f oF2 with a neutral composition and temperature taken from any thermospheric model, for instance, NRLMSISE-00 [52]. The impact of solar and geomagnetic activity is taken into account with the solar radio flux at the wavelength of 10.7 cm (F10.7) and with the Kp-index converted to ap (for conversion; see, for example, ref. [48]), respectively. ap is a linear variant of Kp [53] with values in the range of 0-400. In order to describe the inertia in thermospheric response to enhanced Joule heating, GDMF2 uses instead of ap its time-weighted accumulation as defined by [54]: where t is time and τ is a parameter characterizing the thermospheric response time.
Typical values for τ are around 0.7-0.8, which correspond to response times of 8-13 h [48]. The EUROMAP approach is applied to well-established ionosonde stations that can provide NRT data and have long-term archives of polished products available for the analysis. EUROMAP includes regression models to describe strong negative disturbances (ap(τ, t) > 30) and training models over the previous 28 days. Both these components are specifically tailored for each ionosonde station used in the model. Constructing a storm regression model typically requires a data archive from the station of interest of more than 20 years in order to ensure good enough statistics on storm conditions. The F10.7 is used to characterize the impact of seasonal and 11-year cycle variations of solar activity on f oF2, while T-index is used to construct the f oF2 background for each ionospheric station (for more details on T-index derivation, see [55]). The training model is based on the knowledge that ∆ f oF2 correlates well with its previous value during quiet conditions. Therefore, ∆ f oF2(t) for the ongoing hour can be estimated with the following formula: where ∆ f oF2(t − 1) is the ∆ f oF2(t) of previous hour and coefficients C 0 , C 1 , and C 2 are defined with the least squares fit to the data of the previous 28 days (one solar rotation period).

PECASUS Operational Activity
The first and central task for PECASUS is to operate a framework capable of monitoring all space weather activity that is relevant for aviation, and to issue and disseminate advisories as the ICAO requirements specify. The PECASUS framework, to which all PECASUS consortium members contribute, consists of a number of central hubs and expert groups.
The main Advisory Dissemination Hub, responsible for the provision of advisories to the aviation community, is hosted by the Finnish Meteorological Institute, Finland (FMI), which is also the consortium lead. The Advisory Consolidation and Production Hub is operated by the Solar Terrestrial Center of Excellence, Belgium (STCE). The STCE operators are responsible for the monitoring of aviation-relevant space weather through a specifically developed web-application tool, which supports them to produce the space weather advisories when the pre-defined conditions are met (see Table 1). The content and format of the draft advisories are checked in the Advisory Dissemination Hub before they are pushed to the Aeronautical Fixed Telecommunication Network. The Met Office, United Kingdom (UKMO) has the role of the Resilience Hub to support the other two hubs in standby mode with the capability of taking over either role on short notice. This ensures the continuity of service in the case of technical problems with the advisory production or with the dissemination processes. The three hubs are manned with over a dozen operators, working in shifts, each equipped with a dedicated training in space weather forecasting and PECASUS operations. The provision of data, observations and models are supported by several expert groups, members of the PECASUS consortium: Seibersdorf Laboratories, Austria (SL), German Aerospace Center, Germany (DLR), Space Research Centre of Polish Academy of Sciences, Poland (SRCPAS), Istituto Nazionale di Geofisica e Vulcanologia, Italy (INGV), Frederick University, Cyprus (FU) and South African National Space Agency (SANSA), UKMO, and STCE. Self-evaluation is conducted by the Royal Netherlands Meteorological Institute, the Netherlands (KNMI).
The official ICAO space weather monitoring that was initiated with three global centers was organized with a bi-weekly rotational principle, where a center serves either as an On Duty Center (ODC), a Primary Backup Center (PBC) or a Secondary Backup Center (SBC). The official handovers between the centers is set on Tuesdays, every two weeks, at 08:00 UTC. While PECASUS serves as the official ODC only for two weeks out of every six weeks, the PECASUS operational services are run continuously, with the only difference being that, during the four weeks of the back-up role, the PECASUS advisories are circulated only among the consortium. At the time of this writing, PECASUS issued only a small number of advisories, as the operational services started close to the solar minimum of activity. Nevertheless, several short-lived events had already been observed and a sequence of advisories had been issued for GNSS scintillation, HF-COM PCA, HF-COM post-storm depression, and HF-COM short-wave fadeout. An example post-storm depression event is discussed in Section 5.
Data and models directly related to the ICAO requirements for issuing advisories are continuously monitored. Such critical inputs (described in Section 2) are provided by the PECASUS expert groups and are centralized at STCE and FMI. For resilience of the continuous operations, redundancy is implemented in the data streams when possible: PECASUS is using and monitoring multiple vertical TEC global maps built independently by DLR and INGV, two D-RAP models run independently by UKMO and the NOAA Space Weather Prediction Center, as well as various MUF and f oF2 maps provided by SRCPAS and INGV. PECASUS also monitors the following K-indices: Kp-index from NOAA and Potsdam, K-index measured in the local station in Dourbes (Belgium), as well as NRT Kpindex from NOAA. As a proxy for GOES X-ray measurements, X-ray fluxes at wavelengths of 6-20 nm and 17-80 nm from the LYRA instrument on-board PROBA2 ( [56]) can be used. This solution was tested successfully on 1 June 2021, when data interruption in the GOES X-ray flux was observed around 15:00 UTC and 20:00 UTC. However, it is good to recognize that the PROBA2 back-up arrangement is available only for exceptional cases. LYRA data are usually downloaded with several hours of delay and measurements may contain gaps due to missing down-links or calibration campaigns.
The incoming real-time space weather data are monitored partially automatically. Special routines are implemented to alert the operator on duty when one or more of the critical data or model outputs exceeds specific thresholds. Once the operator is alerted, she/he analyzes the data manually to decide whether an advisory needs to be issued or not. The driving strategy in the automatic routine is to minimize the number of missed events, thereby accepting the potential triggering of false alerts that the operator on duty needs to disregard. Table 2 provides the set of criteria used to alert the operator on duty, who, later on, performs further manual checks to analyze the event. Shortwave Fadeout GOES-16 X-ray flux reached the moderated threshold.
Post-storm depression One or several areas in δ f oF2 maps remain above the moderated threshold for 45 min and one of the K and Kp index was greater or equal to 6 over the past 5 days.
Indirectly related data and models are also intensively used during the operations activities. For example, increases in the GOES-16 proton flux are monitored as possible precursors of HF-COM polar cap absorption and used for surveillance of radiation events in addition to the ground based neutron monitors (see [27,57,58]). Additionally, the solar wind speed, density, and Interplanetary Magnetic Field (IMF) measurements of CMEs and HSS are closely monitored since they are the most relevant parameters in terms of their impact at the Earth, or geoeffectiveness ( [59,60]). STCE forecasters use several methods for obtaining the CME parameters. One of them is CACTus ( [61]), developed at STCE, which detects coronal mass ejections in images sequences from coronagraphs (such as SOHO/LASCO and STEREO/COR2). Furthermore, it allows a real time estimation of the initial speed and angular width of CMEs.

Self Evaluation, User Feedback Assessment, and Training
PECASUS self evaluation includes the verification and validation of products, using independent data sources, analyses of anomalous situations, and annual post-processing of PECASUS data archives and error logs. The work contains also coordinated reporting procedures for off-nominal situations. These are evaluated as part of the ISO quality management system, and the total performance is reported to the Meteorological Operations Group (MOG) that is part of the Meteorological Panel of ICAO. Based on the lessons learned, performance will be a subject of continuous attention.
During the forthcoming years, there will be the important but difficult task of postprocessing the issued advisories. This includes the validation of nowcasts and verification of forecasts. This will provide insight for the correctness and timeliness of advisories, and, in some cases, may also reveal cases where an advisory has incorrectly not been issued. The challenges for validation are caused by the fact that many space weather data (observations and models) are not yet available in a strictly standardized, quality-controlled way, like meteorological data. An important task for self-evaluation is to find potential systematic discrepancies between actual observations and their estimates based on the methods described in Section 2.
As the aviation world is multifaceted, organizing user feedback is not trivial. We hope, however, to organize user engagement events or present space weather at aviation meetings. This will be more effective in the forthcoming years when the Sun's new cycle progresses and the number of advisories issued are increased.
Finally, a space weather service needs trained staff, especially operators. PECASUS duty officers receive dedicated training provided by the Space Weather Education Center (SWEC, https://www.stce.be/SWEC, accessed on 27 August 2021) of the STCE. PECASUS duty officers learn to analyze and interpret the PECASUS data and to navigate through the PECASUS dashboard. As a standard, the STCE organizes 'Space Weather Introductory Courses-SWIC' several times a year. The SWIC includes live lectures and exercises, written course material with references and an evaluation with certificate issued by the STCE. New PECASUS duty officers follow this course to become acquainted with space weather and the STCE space weather bulletins. Tailored courses for aviation end-users are organized on request by the STCE. An online User's Guide to Space Weather Advisories for Pilots is also available (https://www.stce.be/PECASUS_guide4pilots, accessed on 27 August 2021 or http://pecasus.eu/?page_id=593, accessed on 27 August 2021).

PECASUS Advisories on Post-Storm Depression during 28-30 September 2020
The first space weather event with official ICAO advisories took place from 28 September to 1 October 2020. Activity was triggered by high speed solar wind (≈ 600 km/s) emanating from a solar coronal hole. High-wind speed combined with a southward component of IMF enhanced the energy input from the solar wind to the magnetosphere, which launched a geomagnetic storm with provisional Kp = 6− during the three last UT-hours of 27 September. High Kp-values were observed also earlier, around noon of 24 September. The first signs of post-storm depression were observed around 05:30 UT on 28 September by the Australia-Canada-France-Japan consortium (ACFJ), which served as the on-duty center at that time. During the event, ACFJ sent a total of 13 advisories about subsequent depressions in different regions of the globe over a period of 3.5 days. The analysis below is concentrated on the period of 28-30 September UT, when depression was observed in the European sector.
As Figure 5 shows, enhanced geomagnetic activity (with Kp-values 4-5) continued still throughout the day on 28 September (UT hours 0-24 in Figure 5), but activity ceased over the next two days (UT hours 24-72 in Figure 5). Over the next three days (28)(29)(30), moderate post storm depression took place during daytime at all stations. In Nicosia, which is the lowest latitude ionosonde station among those shown, the effect was smallest, with a clear depression only between 5:00 and 11:00 UT on 28 September (UT hours 5-11 in Figure 5). During night times, however, clear depression was observed at all the other stations. The day-night difference in depression intensity is due to the different thermospheric circulation during daytime and nighttime hours. The poleward component in the daytime neutral wind restricts depression to higher latitudes during daytime hours. A comparison between EUROMAP and IRI(STORM) models and observations at Pruhonice and Rome is shown in Figure 6. It is encouraging to notice that during this moderate negative MUF depression event, EUROMAP follows better the observations than the IRI(STORM) model [62].  Figure 7 shows the spatial extent of night time depression on 28 September as characterized by the ratio of MUF3000F2 to its background, suggesting ICAO threshold excesses in the Northern mid-latitude band at longitudes ca. from 30°W to 105°E. This result can be compared with the advisory recommendations by the f oF2 algorithm, which are shown for hours between 19:00 and 24:00 UT in Figure 8. The f oF2 service recommends advisories on MOD level for the European region, but besides mid-latitudes, high latitudes are impacted. On the other hand, the longitudinal extent of the impacted area in the f oF2 map is narrower than that in the MUF3000F2 map. These differences between the two products can be explained by their different approaches: f oF2 maps are mostly controlled by the availability of NRT ionosonde data, while MUF3000F2 maps reflect information also from long-term statistics. Both aspects are considered to be valuable assets in PECASUS operations.

Conclusions and Future Prospects
The present work describes the PECASUS system of data analysis methods and models, which provides advisories on enhanced space weather, compliant with the Standards and Recommended Practices of the International Civil Aviation Organization (ICAO). The present review is based on the authors' experiences gathered during the years 2019-2021, when the PECASUS system was run in its pilot phase for ICAO. In the ICAO framework, advisories are disseminated in three impact areas: radiation levels at flight altitudes, GNSS based navigation and positioning, and HF communication. Space weather activity levels and thresholds for advisory issuances are characterized by a set of geophysical parameters, which are widely used among the space weather community both in research and services. The PECASUS machinery with its core models, their data inputs and augmenting models are summarized in Table 3.

Radiation
Avidos ANeMos Neutron Monitors COMESEP [39] GOES protons HESPERIA GOES X-rays UMASEP-500 [37]; [38] HF EUROMAP [49] IRI [50] ionosondes GDMF2 [48] MUF background [51] F10.7 foF2 warnings NeQuick [19] Kp T-index D-RAP [41] GOES protons GOES X-rays Overall 24/7 solar wind Space Weather Monitoring density, CACTus [61] velocity, IMF Kp NRT Kp (NOAA) Local K GOES X-rays PROBA2 X-rays PECASUS space weather services rely on software solutions and NRT observations that were originally developed for research purposes, and often the development work originated outside the PECASUS consortium. Both these factors naturally pose a risk for the reliability, availability and maintainability of the PECASUS system, which, according to the ICAO requirements, must be above 95%. This challenge is partly tackled by duplicating some of the critical products-such as GNSS TEC, scintillation indices and post-storm depression estimates-with more than one independent methodology. The critical elements are typically assets of extensive use by their developers and by other well-established national and international space weather services, such as that in the NOAA Space Weather Prediction Center (SWPC) and in the Space Safety program of the European Space Agency. For example, GOES-16 X-ray flux and GOES-16 Proton flux are crucial external data taken from SWPC that are directly monitored by the PECASUS operators and are also used as input for models, such as the D-RAP model, and AVIDOS maps. Wide use in operational services motivates the continued maintenance and development of these critical products, although in some cases, sustained funding mechanisms supporting them are missing. ICAO's joining to the stake holders of space weather services will hopefully improve this situation in future.
An important part of PECASUS development is to combine the actual service provision with careful quality management, which includes a clear division of responsibilities, systematic self-evaluation with user feedback collection, and education activities. The responsibility to maintain automated data and product provision routines is shared in PECASUS among the original developers of the algorithms. The 24/7 semi-automatic and manual operations are conducted in the PECASUS Advisory Consolidation Hub, Dissemination Hub and Resilience Hub. Product verification and validation and service performance are part of the self-evaluation activity of PECASUS, as well as the development and setting up of an education program.
The next steps in the PECASUS development will include improvements in its spatial coverage. Particularly in scintillation and post-storm depression, the service performance is constrained by the coverage of input data used in the analysis. In scintillation measurements, information by the station network of 50 Hz receivers (c.f. Figure 2) can be supported with the ROTI proxy (c.f. Section 2.1.1). PECASUS already has routines to create global ROTI maps with time and space resolutions of 15 min and 2.5°in latitude and 5°in longitude. More work is still needed to gain a better understanding of the relationship between ROTI values and ICAO scintillation thresholds, which most likely will vary, according to latitude and level of activity. As improvements in the global TEC mapping algorithms, the PECASUS research teams are investigating possibilities to use direction-based variograms for kriging interpolation and novel methods for differential code bias estimation [63].
For the case of post-storm depressions, no proxies exist to replace actual ionosonde measurements. A way forward would be to add more stations with historical archives to the process that is currently used to create MUF3000F2 maps for the European region [64]. The algorithm combines NRT measurements and historical information from different locations by the ordinary kriging method for improved spatial coverage and better resilience against occasional gaps in observations.
The template of ICAO advisories include also entries for space weather forecasts. The forecast lead times (6, 12, 18 and 24 h, c.f. Figure 1), which come as heritage from meteorological services, are, in most cases, too challenging for space weather services with the current physical knowledge of the phenomena that are of interest for ICAO. However, the common understanding among ICAO, airliners and service providers is that, even with nowcasting information alone, the space weather advisories are a welcome step forward in enhancing safety in civil aviation against space weather related risks. It is clear, however, that demand for forecasts will become stronger once the service users become more familiar with the space weather discipline.
PECASUS partners' work for forecasts of post-storm depressions is a natural starting point to tackle the challenges on space weather forecasts. EUROMAP is already now a working f oF2 forecast method for the European region based on historical and real-time f oF2 observations. Severe and moderate MUF depressions related to high geomagnetic activity with forecasts of 1-24 h lead times can be achieved with higher accuracy than those by IRI(STORM) [49], while no advantage over IRI(STORM) is obtained for forecasts with lead times greater than 3 h during low geomagnetic activity. The next step will be to achieve forecasts of moderate and severe MUF depressions during low geomagnetic activity for all lead times with better accuracy, and to use this model globally for all ionosonde stations providing NRT data. NRT MUF3000F2 maps over Europe are also generated as a PECASUS service [64] and planned to be improved in several ways-by means of data assimilation for a wider network of ionosondes, by better definition of the background ionosphere, and by an interpolation procedure aimed to reduce undesired local and border effects.
PECASUS will continue also the work for improved TEC forecasts. The eSWua TEC services (https://doi.org/10.13127/eswua/tec, accessed on 27 August 2021) include a novel empirical model for global TEC forecasts with a 24 h lead time [65]. The approach applies a nonlinear autoregressive neural network and is able to reproduce ionospheric conditions especially during disturbed periods. This product by [66] is implemented as a real-time service in the Ionosphere Prediction Service (IPS) prototype (http://ips.gsceuropa.eu, accessed on 27 August 2021). Machine learning and data mining covering observations from several solar cycles may be also more generally an attractive way to make progress in space weather forecasts. Long-term archives enable also probabilistic forecasts, whose interpretation is familiar for airliners from meteorological services.
As a step towards physics based modeling, the European Heliospheric Forecasting Information Asset (EUHFORIA, ref. [67]) was implemented at the PECASUS Advisory Consolidation and Production Hub with the goal to support space weather forecasting. EUHFORIA is a 3D Magneto-Hydro-Dynamic solar and heliospheric model that simulates solar wind (including HSS) and the evolution of CMEs from 0.1 AU to 2 AU. The code was specifically designed to estimate arrival times of CMEs and HSSs, and their impact in space weather operational activity. For the estimation of the potential impact of a CME, EUHFORIA needs input data of CME initiation in order to launch a simulation. Such input is available from the CACTus tool [61] which is already in routine use at the Advisory Consolidation and Production Hub. After its validation in an operational environment, EUHFORIA can be used by PECASUS as a supporting asset when considering inputs for the forecast entries of ICAO advisories.
Finally, some effort was directed toward measuring the availability, timeliness and quality of the data provision used in PECASUS and the issued advisories. However, this is an ongoing process that needs to be developed in the future in order for PECASUS to be able to evaluate and quantify its performance and the validity of the issued advisories. For example, validation campaigns will be organized, and statistical studies involving the issued advisories (which will become more numerous as the solar activity increases) will be carried out. While the PECASUS system is still in its infancy, it has already demonstrated that a self-organized community of research institutes is capable of providing a 24/7 operational service for a world organization as ICAO. We expect that-as more space weather-affected industries formulate their requirements-similar concepts and ways of working can be expanded to an increasing number of services.

Funding:
The ICAO member countries of the PECASUS consortium have agreed to maintain the above described space weather services at their own cost during the period from November 2019 to November 2022. Kapanen and Theresa Hoppe for their valuable assistance in preparing the figures for this paper and Petri Koskimaa for the statics on Kp-index.

Conflicts of Interest:
The authors declare no conflict of interest.