Ionospheric Nowcasting Over Italy Through Data Assimilation: A Synergy Between IRI UP and IONORING

An accurate modeling of the ionosphere electron density is pivotal to guarantee the effective operation of communication and navigation systems, particularly during Space Weather events. Despite the crucial contribution of empirical models like the International Reference Ionosphere (IRI), their limitations in predicting ionospheric variability, especially under geomagnetically disturbed conditions, are acknowledged. The solution proposed in this work involves integrating real‐time, spatially distributed ionospheric measurements into climatological models through data assimilation. To enhance our predictive capabilities, we present an upgrade of the IRI UP data‐assimilation method, incorporating real‐time vertical total electron content (vTEC) maps from the IONORING algorithm for nowcasting ionospheric conditions over Italy. This approach involves updating the IRI F2‐layer peak electron density description through ionospheric indices, to finally produce real‐time maps over Italy of the ordinary critical frequency of the F2‐layer, foF2, which is crucial for radio‐propagation applications. The IRI UP–IONORING method performance has been evaluated against different climatological and nowcasting models, and under different Space Weather conditions, by showing promising outcomes which encourages its inclusion in the portfolio of ionospheric real‐time products available over Italy. The validation analysis highlighted also what are the current limitations of the IRI UP–IONORING method, particularly during nighttime for severely disturbed conditions, suggesting avenues for future enhancements.


Introduction
A comprehensive knowledge of the ionospheric channel is necessary to guarantee the effective operation, planning, and management of several radio-communication, telecommunication, and navigation systems.In particular, High-Frequency (HF) point-to-point communications, utilizing the Earth's ionosphere as the propagation medium, play an important role for many modern technological infrastructures and offer a valuable alternative to satellite communications (McNamara, 1991).However, the dynamic nature of the ionosphere poses challenges for HF communications.Indeed, Space Weather events have a notable impact on radio systems, even during relatively mild Space Weather events (see, e.g., Aa et al., 2021), with effects occurring immediately after the event and persisting in the following days (Cannon, 2013;Hapgood, 2019).Extreme solar events, such as flares and coronal mass ejections, result in significant alterations in solar wind and electromagnetic emissions, affecting the Earth's magnetosphere, ionosphere, and atmosphere (Messerotti et al., 2009).To determine the most suitable frequencies for HF communications in these conditions, timely and accurate information about the ionospheric channel is then of utmost importance (e.g., Pietrella & Pezzopane, 2020).
It is well-known that the critical frequencies of ionospheric layers are influenced by solar radiation and the composition and structure of the neutral atmosphere (Chapman, 1931;Ratcliffe, 1972;Rishbet & Garriott, 1969).At the same time, also the geomagnetic field configuration and the geomagnetic activity strongly affect the ionospheric layers (Buonsanto, 1999;Mendillo, 2006).Considering these factors, it becomes evident that a model capable of accurately capturing the ionospheric plasma variability due to solar and geomagnetic activity, and to changes in the background neutral atmosphere, is imperative.While empirical climatological models of the ionosphere, such as the International Reference Ionosphere (IRI) model (Bilitza et al., 2022), play a crucial role, they may fall short in predicting all ionospheric variability, particularly under disturbed geomagnetic conditions (Miro' Amarante et al., 2007;Pignalberi et al., 2016).The IRI model does incorporate a storm option (Araujo-Pradere et al., 2002;Fuller-Rowell et al., 2000) that attempts to simulate ionospheric plasma behavior during magnetically disturbed periods.However, despite ongoing efforts to enhance the model's performance under such conditions, accurately representing the ionosphere's response to severe geomagnetic storms remains a challenge.
To improve our capability of predicting ionospheric conditions, aligning with the prevailing trend in the ionospheric research community, it is essential to augment climatological models with real-time and spatially distributed ionospheric measurements through a data-assimilation process, as recently acknowledged by the IRI community (Bilitza et al., 2011(Bilitza et al., , 2022)).Assuming the availability of various spatially distributed ionospheric measurements, data assimilation entails the amalgamation of these measurements with a background model, to enhance the estimation of ionospheric conditions across the model's coverage area, even in regions lacking direct measurements.This approach allows for the spatial expansion of the model's effectiveness using limited measurements, while simultaneously increasing the accuracy of the model's predictions.Due to the current fast growth of ionospheric observations availability from large and spread ground-based networks (e.g., ionosondes and Global Navigation Satellite System (GNSS) receivers) and satellites' constellations (e.g., low-Earth-orbit in situ measurements and GNSS radio occultation), data-assimilation applications are rapidly growing in number and complexity.Many data-assimilation procedures, relying on IRI as background model, were developed in the past based on the assimilation of diverse ionospheric parameters from different ground-based and satellite-based facilities (Aa et al., 2016;An et al., 2019;Angling et al., 2009;Bilitza et al., 1997;Fridman et al., 2006;Galkin et al., 2012Galkin et al., , 2020;;Habarulema & Ssessanga, 2016;Hernandez-Pajares et al., 2002;Komjathy et al., 1998;Mengist et al., 2019;Migoya-Orué et al., 2015;Pezzopane et al., 2011;Pietrella et al., 2018;Pignalberi, Pezzopane, et al., 2018;Pignalberi et al., 2019;Schmidt et al., 2008;Ssessanga et al., 2015;Tang et al., 2022;Yue et al., 2012).Among these, we focus on the IRI UPdate (IRI UP) method (Pignalberi, Pezzopane, et al., 2018;Pignalberi et al., 2019), which is applied in this work in conjunction with the IONORING vertical total electron content (vTEC) mapping tool (Cesaroni et al., 2021) for a real-time specification of ionospheric conditions over the Italian territory.
The IRI UP method (Pignalberi, Pezzopane, et al., 2018;Pignalberi, Pietrella, & Pezzopane, 2018;Pignalberi et al., 2019;Pignalberi, Pietrella, Pezzopane, & Habarulema, 2021) is a data-assimilation procedure aiming to update the climatological output of the ionospheric F2-layer peak produced by the IRI model.Since the F2-layer peak is an anchor point for the entire vertical electron density profile in IRI, updating the F2-layer peak will update the entire profile.The core of the IRI UP method is the calculation of the IG12 eff ionospheric effective index through the assimilation of ionospheric measurements.IG12 eff is the effective value of the ionospheric index IG12, that is, the 12-month running mean of the index IG (Liu et al., 1983).IG12 is the index used by IRI to describe how the solar activity level affects the ordinary critical frequency of the F2-layer, that is, foF2 (Bilitza et al., 2022;Pignalberi, Pezzopane, et al., 2018), which is related to the F2-layer maximum electron density, that is, NmF2 (Davies, 1990).In particular, IG12 eff is the IG12 value that given as input to IRI exactly reproduces the foF2, and then NmF2, measured at assimilated reference stations.Pignalberi, Pezzopane, et al. (2018) introduced the IRI UP method and applied it over the European region by assimilating foF2 values routinely recorded by a network of ionosondes.IG12 eff values are first calculated for each assimilated ionosonde station, for the specific assimilation time (Pignalberi, Pezzopane, et al., 2018), and then spatially interpolated through the Kriging method (Kitanidis, 1997).Then, the calculated two-dimensional map of IG12 eff is ingested by the IRI model, thus updating the climatological description of foF2.The ingestion of spatially sparse ionospheric observations by ionosondes allows IRI UP to describe weather features at a smaller spatial and temporal scale compared to the underlying IRI climatological pattern.This was testified by the validation study performed by Pignalberi, Pietrella, and Pezzopane (2018) over 30 periods characterized by geomagnetically disturbed conditions.Thereafter, the IRI UP method was further developed by allowing the ingestion of vTEC values from GNSS ground-based networks (Pignalberi et al., 2019).To ingest vTEC values in IRI UP, Pignalberi et al. (2019) developed a statistical procedure to infer NmF2 from vTEC, which is essential to derive IG12 eff values at GNSS receivers' location.The procedure was further refined in Pignalberi, Pietrella, Pezzopane, and Habarulema (2021) by implementing and validating different vTEC calibration procedures in IRI UP.
In this work, the IRI UP data-assimilation method based on the ingestion of vTEC values will take advantage of vTEC maps provided in real-time by the IONORING algorithm (Cesaroni et al., 2021) over the Italian region.The main aim is the production of real-time foF2 maps over the Italian region for nowcasting purposes.To accomplish this task, the IRI UP method has been properly modified and adapted to the IONORING mapping to meet the requirements of a real-time application.Section 2 describes the data and models used in this paper.Specifically, data coming from ionosondes and ground-based GNSS receivers, along with the description of the IONORING algorithm.The same section introduces the IRI background model.Section 3 deals with the IRI UP-IONORING method by carefully describing how IONORING real-time vTEC maps are ingested by IRI UP, how vTEC values are then used to derive modeled NmF2 values at assimilated locations, and how modeled NmF2 values are used to calculate the IG12 eff map to be ingested by IRI to update the F2-layer peak electron density description.In Section 4, the performance of the IRI UP-IONORING method has been assessed for two validation periods encompassing geomagnetically disturbed days, and compared against the background IRI model and other different data-assimilation models.Section 5 discusses the results shown in Section 4 and traces the path for future developments and applications.

F2-Layer Peak Density From Ionosondes
Ionospheric vertical radio soundings are routinely performed by ground-based ionosondes through reflection of vertically propagated radio waves in the HF range (Budden, 1988).The sounding of the ionosphere in the HF range allows obtaining the frequencies reflected by the ionosphere as a function of the virtual height, which is the so-called ionogram (Hunsucker, 1991).Among the reflected frequencies, the most important are the ordinary critical frequencies, which are associated with the maximum electron density of the different ionospheric layers.For radio propagation purposes, foF2, which is associated to the absolute maximum of electron density in the whole ionospheric vertical profile, is the most significant parameter.foF2 is related to NmF2 through the wellknown relation NmF2 = (1.24•10 10 ) foF2 2 , with foF2 in MHz and NmF2 in electrons/m 3 (Davies, 1990).
Over the Italian territory, three ionosondes are installed and routinely perform ionospheric vertical incidence radio soundings.Table 1 collects their main information, and yellow circles in Figure 1 show their spatial distribution.Rome and San Vito ionospheric stations are equipped with DPS Digisondes (Bibl & Reinisch, 1978); foF2 data come from ionograms autoscaled by the Automatic Real-Time Ionogram Scaler with True height analysis (ARTIST) software (Galkin & Reinisch, 2008), and were downloaded from the Digital Ionogram DataBASE (Reinisch & Galkin, 2011) through the SAO explorer software (http://ulcar.uml.edu/SAO-X/SAO-X.html).The ionospheric station of Gibilmanna is equipped with an AIS INGV ionosonde (Zuccheretti et al., 2003); foF2 data come from ionograms autoscaled by the Autoscala system (Pezzopane & Scotto, 2005, 2007;Scotto, 2009) and were downloaded from the electronic Upper Space Weather atmosphere (eSWua, http://eswua.ingv.it/)portal (Upper atmosphere physics and radiopropagation Working Group, 2020a).Rome and Gibilmanna ionosondes are managed by the Istituto Nazionale di Geofisica e Vulcanologia (INGV).The considered foF2 time series have a 15-min time sampling in accordance with the sounding repetition rate of the ionosondes at minutes 0, 15, 30, and 45 of each Universal Time (UT).2022).According to Section 3.3 of Pignalberi et al. (2019), only ionograms with C-Score ≥75 have been considered to establish the statistical relation between NmF2 and vTEC.Differently, for validation purposes, ionograms recorded by the three Italian ionosondes have been manually scaled by an expert operator.

vTEC From Ground-Based GNSS Receivers and IONORING
Due to the ionospheric plasma's dispersive properties, multi-frequency GNSS receivers allow TEC estimates.TEC represents the integral of electron density along the path between the satellite and receiver at ground level, typically expressed in TEC units (TECU), which equate to 10 16 electrons/m 2 .Today, the presence of extensive networks of global, regional, and local GNSS receivers makes it easier to leverage TEC data for mapping and modeling the ionosphere across various spatial and temporal scales, facilitating the development of services tailored to different categories of end-users (Johnston et al., 2017;Kauristie et al., 2021;Vadakke Veettil et al., 2019).Among these services, IONORING (IONOspheric-RING) utilizes data from the National Integrated GNSS Network (RING, http://ring.gm.ingv.it/),which is maintained by INGV (INGV RING Working Group, 2016), to provide real-time TEC mapping over Italy.RING is a geodetic research infrastructure consisting of more than 200 GNSS sensors distributed across Italy (red circles in Figure 1).Each station is equipped with a professional GNSS receiver and a choke ring antenna.These characteristics make RING receivers well-suited for ionospheric studies, as evidenced in recent literature (Cesaroni et al., 2017;Musicò et al., 2018;Tornatore et al., 2021).Some of the RING receivers have been specifically chosen to provide input data for the IONORING algorithm for real-time TEC nowcasting (Cesaroni et al., 2021).This selection considers the balance between achieving comprehensive coverage of the ionospheric region and minimizing the time it takes for IONORING products to become available.To constrain the boundaries of the designated area, receivers in Ajaccio (AJAC), Graz (GRAZ), Marseille (MARS), Sarajevo (SRJV) and Cagliari (CAGL) (blue circles in Figure 1), which are available through the EUREF Permanent GNSS Network (https://epncb.oma.be/), are also included in the data collection process.
The input data for IONORING consists of GPS carrier phase observations on L1 and L2 frequencies (1,575.42 MHz and 1,227.60 MHz,respectively), Coarse Acquisition (CA) code observations on L1, and Y-code observations on L2.These observations are broadcasted at a 1-s interval through the Networked Transport of RTCM via Internet Protocol (NTRIP) by the selected receivers.Additionally, navigational information, such as orbital parameters for GPS satellites, is obtained through NTRIP streaming from the mentioned receivers or, alternatively, from IGS navigational information streaming.We obtain ionospheric slant TEC (sTEC) by analyzing code and phase measurements at L1 and L2 frequencies, which represent the integrated electron density along the slant path between the satellite and the receiver.From sTEC evaluations for each satellite-receiver couple, we derive maps of vertical TEC (vTEC).The procedure to retrieve sTEC and the corresponding mapping to vTEC is detailed in Cesaroni et al. (2021), specifically in their Section 2.
The vTEC maps are presented as the smoothed vTEC curve, which is calculated across a standardized grid with a spatial resolution of 0.1°, both in latitude and longitude.Specifically, IONORING interpolates the vTEC values estimated at every ionospheric pierce point (IPP), from all the considered satellite-receiver pairs, through a non- parametric local regression technique named LOWESS (Cleveland, 1979).The LOWESS interpolation technique mitigates the impact of unknown parameters but at the same time smooths out small-scale spatial variations in the final vTEC map.The IONORING grid encompasses the geographic area between latitudes 35°N and 48°N, and between longitudes 5°E and 20°E.The refresh rate of the map is 10 min and the entire processing chain is characterized by an average latency of under 1 min.The reliability of IONORING has been established through a comprehensive evaluation involving a 3-year timeframe from May 2017 to April 2020 (Cesaroni et al., 2021).During this period, vTEC maps generated by IONORING were compared to independent vTEC measurements using various techniques, revealing the high level of accuracy of IONORING products.The real-time vTEC mapping service is available in the eSWua portal (http://eswua.ingv.it/)managed by INGV as part of the provided data collections (Upper Atmosphere Physics And Radiopropagation Working Group, 2020b).An example map of nowcasted vTEC over Italy provided by IONORING for 21 September 2023 at 15:00 UT is shown in Figure 2.
Real-time vTEC maps over Italy nowcasted by IONORING feed the IRI UP data-assimilation method, as detailed in Section 3.1.Differently, vTEC values needed to establish the statistical relation between NmF2 and vTEC have been obtained directly by specific ground-based GNSS receivers co-located with Rome and San Vito ionosondes.Specifically, from the RING and EUREF networks, we selected the GNSS receivers nearest to the two ionosondes, and whose data set covers the ionosondes' data set extension.These constraints are satisfied by the INGR receiver of the RING network, which is exactly co-located with the Rome ionosonde, and by the MATE receiver of the EUREF network, which is 93 km apart from San Vito ionosonde.Detailed information on these GNSS receivers is provided in Table 2. vTEC time series from INGR and MATE GNSS receivers have been obtained by applying the Ciraolo et al. (2007) calibration procedure to daily RINEX files, containing L1 and L2 code and carrier phase data at a 30s time rate.In the calibration procedure, the equivalent vTEC at the IPP of 350 km is calculated assuming a thin spherical shell model of the ionosphere and a cut-off elevation angle of 50°.As a consequence, GNSS satellites whose IPP is within a distance of about 300 km from the vertical of the ground-based GNSS receiver are considered in the vTEC calculation.Considering this and the horizontal correlation distance of the ionosphere (Forsythe et al., 2020), MATE GNSS receiver and San Vito ionosonde can be safely considered as co-located.The reliability and robustness of the Ciraolo et al. (2007) vTEC calibration procedure in the IRI UP method has been carefully assessed by Pignalberi, Pietrella, Pezzopane, and Habarulema (2021).
Since the ionosondes' sampling rate is 15 min, vTEC time series at a 30-s time resolution were averaged and weighted by using a sliding window centered at minutes 0, 15, 30, and 45 of each UT, as detailed in Section 3.4 of Pignalberi et al. (2019).Details on the data outliers filtering is provided in Section 3.3 of Pignalberi et al. (2019).Finally, vTEC time series at 15-min time rate, simultaneous and co-located to ionosondes' observations, have been obtained for the years covered by ionosondes (see Tables 1 and 2).These vTEC time series have been used in conjunction with corresponding NmF2 time series to establish the statistical relation between NmF2 and vTEC, as detailed in Section 3.2.

The International Reference Ionosphere (IRI) as Background Model
IRI is an empirical climatological ionospheric model relying on ground-based and space-based observations of the ionosphere collected over the last decades on a global basis (Bilitza et al. (2022) and references therein).IRI project was initiated in 1968 sponsored by the Committee on Space Research (COSPAR) and the International Union of Radio Science (URSI).As of April 2014, IRI has become the official International Standardization Organization (ISO) standard for the ionosphere (ISO, 2009).The last version of the model is IRI-2020, which is the version used as background model in IRI UP.
Among the several ionospheric parameters modeled by IRI, here we focus on the electron density and, specifically, on the F2-layer peak, which is an anchor point for the whole vertical electron density profile.The F2-layer peak is characterized by the highest electron density recorded in the ionosphere (NmF2), which is associated with the maximum ordinary critical frequency reflected by the ionosphere at vertical incidence (foF2), at the altitude hmF2.Both F2-layer peak characteristics are empirically modeled by IRI through different options which are carefully described in Bilitza et al. (2022).In this work, for the modeling of foF2 we use the CCIR option (CCIR, 1967), which is recommended over land regions; while for hmF2 we use the default Shubin et al. (2015) option.Starting from the IRI-2001 version, an additional feature has been incorporated to estimate foF2 during periods of geomagnetic disturbance, known as the storm option (Araujo-Pradere et al., 2002;Fuller-Rowell et al., 2000).
The solar activity variation in the foF2 modeling is described by IG12, which is the 12-month running mean of the IG ionospheric index (Liu et al., 1983).IG is calculated from linear regression equations based on monthly median foF2 values recorded at noon at reference ionosonde stations (originally 11 stations, currently four) and corresponding CCIR-modeled values.The resulting index is the median of IG index values calculated in those reference locations, and is therefore global by definition.Recently, Brown et al. (2018) showed that hemispheric and monthly IG indices would improve the description of foF2 by IRI.The benefit of using the IG index is that it encompasses dynamic effects derived from ionospheric measurements, which are not accounted for by solar indices.IRI UP takes advantage of the foF2 dependence on IG12 by IRI to calculate effective indices through data assimilation (Pignalberi, Pezzopane, et al., 2018;Pignalberi et al., 2019).Section 3.3 details the calculation of ionospheric effective indices by IRI UP.
Over the years, IRI model performance was assessed in many validation studies performed by different groups.For a detailed validation of IRI F2-layer peak characteristics refer to Pignalberi, Pietrella, & Pezzopane (2021).Moreover, IRI was compared with other empirical and physics-based models, resulting in most of the cases the best, or at least on par with the best theoretical and assimilative models, in predicting NmF2 and hmF2 (Shim et al., 2011(Shim et al., , 2012(Shim et al., , 2017(Shim et al., , 2018)).IRI-2020 model Fortran code is available at http://irimodel.org/.

On the Application of the IRI UP Methodology to IONORING Real-Time vTEC Maps
The core of the IRI UP method was originally presented by Pignalberi, Pezzopane, et al. (2018) and Pignalberi, Pietrella, and Pezzopane (2018), and it was based on the assimilation of F2-layer peak characteristics from European ionosondes.Later on, the procedure to assimilate vTEC observations by a network of ground-based GNSS receivers in IRI UP was introduced by Pignalberi et al. ( 2019), focusing on the South African region.Pignalberi, Pietrella, Pezzopane, and Habarulema (2021) also investigated the impact of the use of different vTEC calibration procedures in the IRI UP method.The latter is the IRI UP methodology applied in this work, with the crucial difference that vTEC observations are assimilated from real-time IONORING maps, thus allowing a nowcasting of foF2 over the Italian territory.In this section, the main steps about the application of the IRI UP methodology to IONORING real-time vTEC maps are described.The focus is on the differences and new achievements compared to the previous IRI UP implementations; for a detailed description of some procedures, the reader will be redirected to the previous papers.

Assimilation of IONORING Real-Time vTEC Maps in IRI UP
The first step of the IRI UP method concerns the assimilation of vTEC values from real-time maps produced by IONORING over the Italian region.IRI UP runs over the same region covered by IONORING maps, that is, latitudes from 35°N to 48°N, and longitudes from 5°E to 20°E, but with a different spatial resolution.Indeed, after a preliminary testing phase, it resulted that assimilating vTEC values with the original spatial resolution of IONORING maps, that is, 0.1°in both latitude and longitude, was too costly in terms of the time needed to perform the entire data-assimilation process.By running at the native IONORING maps spatial resolution, it would have been difficult to meet the requirements of a real-time procedure providing an output in less than a minute.On the other hand, from our preliminary testing phase, it came out that no advantage, in terms of performance in reproducing the foF2 at ionosonde validation sites, would have been obtained from the use of such a high spatial resolution.This is probably due to the LOWESS interpolation technique applied by IONORING to spatially interpolate vTEC values.Despite being produced at 0.1°spatial resolution, both in latitude and longitude, IONORING vTEC maps do not show appreciable variations at this small spatial scale due to the smoothing applied by LOWESS.As a consequence, we decided to sample the IONORING map with a coarser spatial resolution.By considering the balance between the performance in predicting foF2 at ionosonde validation sites and the procedure running time, we found the best compromise at 0.5°spatial resolution in both latitude and longitude.
IONORING maps assimilated by IRI UP are slightly different from those delivered through the eSWua portal (http://eswua.ingv.it/).In fact, for data assimilation in IRI UP, vTEC values at the thin shell height surface are interpolated by IONORING up to a spherical distance of 1°from the nearest IPP.This is to avoid an excessive extrapolation of vTEC values in regions not covered by observations by the IONORING spatial interpolation procedure (Cesaroni et al., 2021).Moreover, since IRI UP already performs the mapping of ionospheric effective indices, no advantage would have been obtained by assimilating vTEC values too distant from real observations.An example of such IONORING vTEC map is shown in Figure 3a for 28 March 2022 at 00:30 UT, which clearly shows the effect of considering only observations within 1°from the nearest IPP when compared to the map in Figure 2. From Figure 3a, vTEC values sampled at 0.5°spatial resolution are considered for data assimilation in IRI UP, as shown in Figure 3b.

From vTEC to NmF2 Through Reference Ionospheric Stations
To calculate the ionospheric effective indices within the IRI UP approach, Pignalberi et al. ( 2019) developed a statistical procedure to derive a mathematical relation between vTEC and NmF2 based on long time series from selected co-located ionosondes and ground-based GNSS receivers.This approach is based on the well-established correlation between vTEC and NmF2 (Gerzen et al., 2013;Kouris et al., 2004;Krankowski et al., 2007;Spalla & Ciraolo, 1994;Ssessanga et al., 2014), which is also at the base of the ionospheric equivalent slab thickness definition (Pignalberi, Nava, et   To properly represent the diurnal and seasonal variations embedded in both NmF2 and vTEC, the data set from colocated stations has been sorted as a function of the Local Time (LT) and of the month of the year.Specifically, for the diurnal variation, 15-min wide bins have been used, then 96 bins in total, and, for the seasonal variation, the 12 months of the year have been considered.As a consequence, the original data set was split in 1,152 bins (96 LTs × 12 months), as a function of LT and month.Equation 1 has been applied to each of these subdatasets, then providing linear coefficients as a function of LT and month of the year, that is, a NmF2 (LT, Month) and b NmF2 (LT, Month), that are reported in Figure 4 in a matrix representation.
Figure 4 shows how both a NmF2 and b NmF2 exhibit a diurnal and seasonal variation, with large variations between day and night associated to notable gradients at solar terminator hours.By describing the relation between NmF2 and vTEC through Equation 1, the solar activity level variation is embedded in the slope coefficient a NmF2 , once the diurnal and seasonal variations have already been made explicit.Differently, the geomagnetic activity variation has not been sorted due to the lack of enough data.As a consequence, coefficients in Figure 4 have to be considered as climatological and representative of the median behavior between NmF2 and vTEC over the Italian region.
Pignalberi et al. ( 2019) and Pignalberi, Pietrella, Pezzopane, and Habarulema (2021) calculated a NmF2 and b NmF2 coefficients as a function of UT with one-hour time resolution.As a step forward, here we described the diurnal variation as a function of LT with a 15-min time resolution, that is, the same as the NmF2 and vTEC data sets.The description of the LTs around the solar terminator (both at sunrise and sunset), where large electron density gradients (both in space and time) lie, will benefit from these improvements.To evaluate the consistency of the proposed relation between NmF2 and vTEC, and its ability in modeling NmF2 as a function of vTEC with the coefficients provided in Figure 4, we performed a self-validation based on data from co-located ionosonde and GNSS stations (see Sections 2.1 and 2.2).Specifically, given the vTEC values measured by the GNSS receivers, NmF2 values were modeled through Equation 1 with a NmF2 and b NmF2 coefficients provided by Figure 4 as a function of LT and month of the year.NmF2 values were then converted into foF2, and were finally compared with foF2 values measured by co-located ionosondes.For each of the LT and month of the year bins, the following statistical metrics were calculated between the modeled (foF2 model ) and measured (foF2 ionosonde ) foF2 values: Mean of residuals (Equation 2), Root Mean Square Error (RMSE, Equation 3), Normalized RMSE (NRMSE, Equation 4), and Pearson correlation coefficient (R, Equation 5).

Mean of residuals
where N is the number of data of the time series, f oF2 ionosonde is the arithmetic mean of foF2 ionosonde , cov is the covariance, and σ the standard deviation.NRMSE ranges between about 8% and 16%, with maxima at sunrise and sunset for winter months, and at sunset for summer months.Hence, percentually, the lowest error in reproducing the measured foF2 is made during daytime when the correlation between foF2 and vTEC is known to be more consistent.To demonstrate that, panel (d) reports Pearson's correlation coefficient values.Correlation coefficients range between 0.5 and 1.0 with most of the values above 0.8.Specifically, the correlation between foF2 and vTEC is remarkably high during daytime for all seasons, and during nighttime for summer and equinoctial months.This testifies the suitability of Equation 1 linear relation between NmF2 and vTEC.The lowest correlation values characterize the nighttime hours in winter months, which are the conditions for which the contribution of the plasmasphere to the vTEC is not negligible compared to that provided by the ionospheric F2-layer (Klimenko et al., 2015;Yizengaw et al., 2008).
By averaging on all the LTs and months, from Figure 5 matrices we obtain an average mean of residuals = 0.027 MHz, an average RMSE = 0.57 MHz, an average NRMSE = 10.64%, and an average R = 0.87.These average values give a hint on the overall performance of the statistical procedure developed to infer NmF2 from vTEC observations.Equation 1, with coefficients provided in Figure 4, is used operatively to infer NmF2 whenever a vTEC value is available.Taking as example the case in Figure 3 for 28 March 2022 at 00:30 UT, assimilated vTEC values shown in Figure 3b are used in Equation 1, with a NmF2 and b NmF2 coefficients taken from matrices in Figure 4 for the

Calculation of Effective Index Map and Ingestion in IRI
The core of the IRI UP method is the calculation of ionospheric effective indices based on F2-layer peak characteristics assimilated at selected ionospheric stations (Pignalberi, Pezzopane, et al., 2018).In this work, we focus on the IG12 ionospheric index (Liu et al., 1983) which describes how the solar activity level affects the foF2 (and then NmF2) modeling in IRI (Bilitza et al., 2022;Pignalberi, Pezzopane, et al., 2018).Through the assimilation of foF2 observations by a network of ionosondes (Pignalberi, Pezzopane, et al., 2018) or vTEC by a network of ground-based GNSS receivers (Pignalberi et al., 2019), it is possible to calculate effective values of IG12, namely IG12 eff , at assimilated stations.IG12 eff is the value that, given as input to IRI, exactly reproduces the observed foF2 at the assimilated station.Spatially scattered IG12 eff values are then used to produce a regional map of IG12 eff to be ingested in IRI, and finally produce updated foF2 (and then NmF2) values over the considered region (Pignalberi, Pezzopane, et al., 2018).
IG12 eff values are calculated through observed foF2 values at assimilated stations through the procedure developed by Pignalberi, Pezzopane, et al. (2018).Briefly, for each assimilated station, we calculate the squared difference between the observed foF2 value (foF2 obs ) and those modeled by IRI as a function of IG12 (foF2 IRI (IG12)): Specifically, foF2 IRI (IG12) in Equation 6are foF2 values modeled by IRI at the location of the assimilated station, at the assimilation time, and as a function of IG12.IG12 eff is the value which minimizes Δ IG12 , that is, is the IG12 value that, given as input to IRI, exactly reproduces the observed foF2.
For a specific assimilation time, Equation 6 is applied to each assimilated point to calculate corresponding IG12 eff values.An example is shown in Figure 6b based both on assimilated vTEC values shown in Figure 3b and on NmF2 values modeled as described in Section 3.2.In Figure 6b, calculated IG12 eff values are shown as colored circles over assimilated points, which are overplotted to the IG12 eff map, and that have to be thought of as virtual stations.IG12 eff values at assimilated points are then used to produce a spatial map of IG12 eff over the considered region through the Kriging geostatistical mapping procedure (Kitanidis, 1997).Pignalberi, Pezzopane, et al. (2018) and Pignalberi et al. (2019) applied the Universal Kriging method to map IG12 eff over the European and South African regions, respectively.Since the limited spatial extension of the region covered by the Italian territory, we decided to apply the Ordinary Kriging method.In the Ordinary Kriging, the deterministic part, which describes the large-scale spatial variability, is taken as a constant; then, spatial correlations are described only by the stochastic part through the variogram model (for more details see Section 3.3 of Pignalberi, Pezzopane, et al. (2018)).
The core of the Kriging interpolation method is the building of the experimental variogram and the following empirical modeling.Calculated IG12 eff values at assimilated points are used to calculate semivariance values as a function of the distance between pairs of assimilated points, as detailed in Section 3.3 of Pignalberi, Pezzopane, et al. (2018).By definition, the semivariance γ is half the squared difference between a pair of measurements (Kitanidis, 1997;Pignalberi, Pezzopane, et al., 2018), namely: where x 1 and x 2 are the vector positions of a pair of IG12 eff measurements.Semivariance values are calculated for each pair of assimilated points, which for N stations are N(N 1) 2 in total.Since the large number of assimilated points provided by IONORING maps, semivariance values calculated through Equation 7 are binned as a function of the distance between pairs of assimilated stations, as described in Pignalberi, Pietrella, and Pezzopane (2018).The binned semivariance values plotted versus the distance between pairs of assimilated points is the experimental variogram, an example of which is represented in Figure 6a, where the red points are IG12 eff semivariance binned values.
The experimental variogram gives information about the scale of fluctuations of the represented quantity.In particular, the behavior of the experimental variogram at small distances between points describes the variability at small scales; instead, the behavior of the experimental variogram at large distances describes the variability at large scales.In a nutshell, the experimental variogram describes the spatial correlations between measurements, which are in turn related to the covariance and variance between measurements (Kitanidis, 1997).Experimental semivariance values are modeled through the fitting of simple mathematical functions to obtain the variogram model to be used in the Kriging interpolation method (Pignalberi, Pezzopane, et al., 2018).In this work, two semivariogram models are applied: the non-stationary Power variogram model (Equation 8), and the stationary Spherical model (Equation 9): where h is the distance between pair of measurements, c 0 is the nugget value, θ > 0 is the scale, 0 < s < 2 is the exponent, σ is the sill, and ⍺ the range.
The Power variogram model shows a non-stationary behavior with semivariance increasing with distance in the entire range of considered values (as in Figure 6a).Differently, the Spherical variogram model, after increasing at small scales, reaches a stationary value around the sill at the range distance.For a thorough description of the variogram models refer to Section 3.3 of Pignalberi, Pezzopane, et al. (2018).
The choice of the best variogram model fitting the experimental semivariance values has been automated by Pignalberi, Pietrella, and Pezzopane (2018) based on the Q 1 , Q 2 , and cR statistical metrics build on residuals between experimental and modeled semivariance values (for more details refer to Pignalberi, Pietrella, and Pezzopane ( 2018)).This automatic procedure has been implemented here for the choice of the best variogram model between the Power and the Spherical ones.Figure 6a  Since the nugget describes the small-scale spatial variations in assimilated data (Pignalberi, Pezzopane, et al., 2018), a null nugget testifies that assimilated data show no spatial variations at a scale lower than the minimum distance between pairs of assimilated points, that is, 0.5°both in latitude and longitude.This is probably due to the LOWESS interpolation technique used by ION-ORING which suppresses most of the spatial variations at a scale below 0.5°.In IRI UP-IONORING we obtained a null nugget in the variogram model in a very high percentage of cases (>95%).
IG12 eff maps are produced with a spatial resolution of 0.25°, both in latitude and longitude; this is also the spatial resolution of foF2 maps given as output by IRI UP.
The final step of the IRI UP method is the ingestion of the calculated IG12 eff map in the IRI model to update the IRI description of the F2-layer peak density.
For the case under study, that is, 28 March 2022 at 00:30 UT, the ingestion of the IG12 eff map (Figure 6b) in IRI produces the foF2 nowcasted map shown in Figure 7b.For comparison, Figure 7a shows the background foF2 map as modeled by the climatological IRI (IRI-2020 version) through the IG12 index.
Figure 7c shows the residuals between IRI UP (Figure 7b) and IRI (Figure 7a) maps.The ingestion of IG12 eff values changes both the magnitude and spatial distribution of modeled foF2 values.Differences between IRI UP and IRI reach about 1.0 MHz in specific sectors, which is significant considering that, for these conditions, foF2 lies between 3.8 and 5.0 MHz.

Assessment of the IRI UP-IONORING Performance in a Nowcasting Scenario
We assessed the performance of the IRI UP-IONORING nowcasting procedure in two periods characterized by scenarios spanning from quiet to moderately and severely disturbed geomagnetic conditions, and also including the occurrence of two strong solar flares.The two validation periods are: Here, we show and discuss the validation analysis results for the severely disturbed period 21-29 March 2023, with corresponding plots of Figures 8-14 and summary statistics of Table 3. Results for the moderately disturbed  4, while corresponding plots are provided in Supporting Information S1 (Figures S1-S6).
In the validation analysis, we consider as reference the foF2 observations from the three Italian ionosondes (Table 1).In particular, foF2 values have been manually scaled by an expert operator to ensure the highest reliability.
IRI UP-IONORING nowcasting procedure has been compared against different climatological and nowcasting procedures.Since the first goal of a data-assimilation procedure is to augment the performance of the underlying background model, the climatological IRI-2020 model, with the storm option on, has been considered in the validation analysis.To evaluate the accuracy and reliability of the use of nowcasting IONORING vTEC maps in IRI UP, in the validation analysis we considered also the IRI UP method but fed by vTEC values calibrated through the et al. ( 2007) procedure (see Section 2.2), which are obtained from each ground-based GNSS receiver of the RING network.Hence, in this case, assimilated vTEC values have the spatial distribution of the RING network Figure 1), which differ from the spatial distribution of vTEC values sampled from ION-ORING maps (see for example Figure 3b).It is worth mentioning that, compared to IONORING nowcasting maps, the Ciraolo et al. (2007) vTEC calibration procedure cannot be applied in a nowcasting scenario since it requires a complete arc of observations for each satellite in view from all the receivers at ground.Finally, the validation analysis includes also the Simplified Ionospheric Regional Model UPdate (SIRMUP, Pietrella et al. (2020) and references therein) nowcasting model.SIRMUP has been applied to generate foF2 nowcasting predictions over Rome, San Vito, and Gibilmanna ionosonde locations.SIRMUP is based on the background given by the SIRM model, which is fed with an effective value of the R12 (the 12-month running mean of the monthly sunspot number) index, that is, R12 eff .R12 eff is calculated on the basis of foF2 values autoscaled at reference ionosonde stations, and then used to generate a nowcasting foF2 map extending in longitude from 5°W to 40°E and in latitude from 34°N to 60°N.In this study, the reference stations are the permutation of available Italian ionosondes excluding the validation station under consideration.For example, to obtain foF2 SIRMUP nowcasting predictions over Rome, only the foF2 values autoscaled at San Vito and Gibilmanna have been used to calculate R12 eff .It is worth mentioning that SIRMUP, compared to the two IRI UP implementations, assimilates directly foF2 from ionosondes; as a consequence, it is considerably advantaged in the description of foF2.
Moreover, the underlying SIRM model was specifically developed for the European region, while the IRI model is instead global.Nevertheless, SIRMUP was added to the validation to serve as the benchmark.
Since IONORING vTEC maps are produced every 10 min at minutes 00, 10, …, and 50 UT, while ionosondes foF2 values are recorded every 15 min at minutes 00, 15, 30, and 45 UT, for a fair comparison, in the validation analysis only the values at minutes 00 and 30 UT have been considered.To summarize, the validation analysis is based on foF2 values observed by Italian ionosondes and those modeled by IRI-2020, IRI UP-IONORING, IRI UP-Ciraolo, and SIRMUP over the ionosondes location, for the two validation periods before presented, at 30min time resolution.
The validation analysis performed on the three Italian ionosondes' locations allows evaluating if the performance of IRI UP-IONORING is dependent on the location.In fact, considering the spatial distribution of RING GNSS receivers (Figure 1), data used by IONORING are mostly concentrated on the Italian land along the Apennines.As a consequence, Rome station is surrounded by many GNSS receivers compared to Gibilmanna and especially to San Vito.Moreover, Rome station is in the middle of the IONORING spatial domain, while Gibilmanna and San Vito are toward the border of such a domain.In the light of these considerations, we would expect better performance for Rome than for Gibilmanna and San Vito.5) and on the counts, that is, the number of points in the scatter plot.
Statistical metrics values reported in Figures 9, 11, and 13, for the 21-29 March 2023, are collected in Table 3 for an easy comparison.Table 4 contains corresponding values concerning the period 28 March 2022-4 April 2022.
Time series of foF2 at Rome for the period 21-29 March 2023, Figure 8a, clearly show the effect of the severe geomagnetic storm started on 23 March (Day Of the Year, DOY, 82), on the ionospheric F2-layer.Specifically, compared to previous relatively quiet days, measured values show an increase on DOY 82 at dayside, during the main phase of the storm.The initial increase in foF2 is then followed by a decrease on DOY 83 at dayside, and between DOYs 83 and 84 at nightside, during the recovery phase.Afterward, the foF2 slowly transitions to the climatological pattern.The IRI-2020 model shows a consistent underestimation at dayside during quiet days, while during disturbed days it correctly describes the negative phase of the storm albeit overestimating for both calibrated values show a consistently similar behavior for quiet conditions before the storm commencement, while they differ under disturbed conditions.The largest differences between the two are visible between DOYs 83 and 84 at nightside, where IG12 eff IONORING shows much lower values than Ciraolo ones.A similar behavior is visible also for the subsequent nights under disturbed conditions.This remarkable difference explains the fact that, under these conditions, the IRI UP-IONORING procedure obtains lower foF2 values than those obtained by the IRI UP-Ciraolo procedure, the latter showing a remarkable agreement with ionosonde observations.SIRMUP shows very consistent performance for both quiet and disturbed conditions with the exception of the overestimation at dayside for DOYs 83 and 84.
A general overview of the performances is provided by the histograms of residuals and scatter plots in Figure 9. IRI-2020 shows the largest dispersion and the poorest agreement with ionosonde observations.IRI-2020 mean of residuals is relatively good but this is due to the partial compensation between dayside underestimation under quiet conditions and a general overestimation under disturbed conditions.Both IRI UP implementations show a noticeable improvement for both the precision and accuracy in the description of foF2.However, while they show similar dispersion, IRI UP-IONORING shows a lower accuracy due to the consistent underestimation, which was already highlighted in Figure 8. SIRMUP shows very consistent performance but with some outliers when compared to IRI UP implementations.Overall, IRI UP-Ciraolo performances are comparable to SIRMUP, while IRI UP-IONORING suffers from the underestimation issue at nightside under disturbed conditions.However, by looking at slope coefficients from the linear fits, it comes out that IRI UP-IONORING has the best slope coefficient (the nearest to 1), and then a remarkable linear dependence from ionosonde observations.As a consequence, once the underestimation issue has been addressed, IRI UP-IONORING would show the best agreement, both in accuracy and dispersion, with observations.3).This is particularly important for IRI UP-IONORING because it is an indirect confirmation that the IONORING vTEC mapping procedure does not introduce any artifact in the spatial interpolation algorithm, except for the smoothing of small-scale variations by the LOWESS interpolation technique, and that the IRI UP Kriging spatial interpolation for IG12 eff mapping is robust.
The statistical analysis based on ionosondes' observations as reference gives a clear picture of how IRI UP can describe the foF2 variations as a function of time and geomagnetic activity, but limited to the ionosondes' site location.It is also interesting to see how IRI UP describes the spatial variations of foF2 over the Italian region by looking at foF2 maps.In Figure 7, we showed examples of foF2 maps generated by the IRI model and the IRI UP-IONORING method, and corresponding residuals, for a geomagnetically quiet period.In Figure 14 we show similar plots but for the 24 March 2023 (DOY 83), that is, the most disturbed day of our validation data set, for three different hours, specifically at 02:00 UT, 14:00 UT, and 22:00 UT.At 02:00 UT, the Dst reached its minimum value and Kp its maximum as highlighted in Figure 8c, while the hours 14:00 UT and 22:00 UT are in the recovery phase of the geomagnetic storm.However, the ionosphere shows a time lag of many hours, and effects persisting for several days, compared to the disturbance in the geomagnetic field, as evidenced by foF2 time series in Figures 8a, 10a, and 12a, showing the largest departures from the climatological background in the recovery phase of the storm.Specifically, at 02:00 UT foF2 shows little variations compared to previous quiet days, while at 14:00 UT, and even more at 22:00 UT, foF2 shows a remarkable negative phase.foF2 maps in Figure 14 quantify this negative phase effect and show corresponding spatial variations.Indeed, IRI UP maps obtained through the assimilation of vTEC values from IONORING maps show spatial variations of foF2 different from the underlying IRI model as testified by residuals' plots.Interestingly, at 02:00 UT and 14:00 UT the lowest residuals are in the northernmost sector, while at 22:00 UT the lowest residuals are in the southernmost sector although in a general context of lower foF2 values compared to IRI.This example highlights how IRI UP can improve not only the quantification of the geomagnetic disturbance effect on foF2, but also the corresponding spatial variations which can be different from the climatological pattern described by IRI under severe geomagnetic activity conditions.
Figures S1-S6 in Supporting Information S1 show the validation analysis for the period 28 March-4 April 2022, with statistical metrics reported in Table 4.This period is characterized by quiet to moderately disturbed geomagnetic conditions as shown by Dst and Kp indices in panel (c) in Figure S1 in Supporting Information S1.Moreover, this period is characterized by two powerful solar flares.The combination of solar flares and a moderate geomagnetic disturbance enhances foF2 above the quiet-time background pattern at dayside on DOY 90.Both IRI UP implementations and also SIRMUP catch this enhancement which instead is not described by IRI-2020 due to the moderate magnitude of the geomagnetic storm, and also because IRI does not include any short-term solar driver to account for solar flares effects.Overall, IRI UP-IONORING and IRI UP-Ciraolo well describe the foF2 behavior at the three ionosonde validation sites, for both the diurnal variation and magnitude values, under the different conditions experienced by the ionosphere in this period.It is worth highlighting that, in this case, IRI UP-IONORING and IRI UP-Ciraolo performances are very similar both in accuracy and precision.Specifically, the IRI UP-IONORING underestimation under severely disturbed conditions at nightside is not present under moderately disturbed conditions, which confirms that this issue is limited to very disturbed conditions.As a consequence, we verified that, excluding the severely disturbed conditions at nightside, IRI UP-IONORING and IRI UP-Ciraolo show similar performance across the different validation sites.This is very promising considering that only the IRI UP-IONORING procedure can work in real-time, with a performance comparable to that of the methodology using vTEC data calibrated through the Ciraolo's algorithm.

Discussion and Conclusions
In this work, we have proposed a new nowcasting procedure of the ionospheric parameter foF2 over the Italian territory based on both IONORING and IRI UP methods.The study showed that the new procedure is promising, being by far better than the background IRI model, and with corresponding results that are similar to those obtained by procedures that in principle should be more powerful, like the IRI UP-Ciraolo and SIRMUP ones.which set up small-scale time and spatial variations in the ionosphere, the fact that noticeable differences are present also for quiet conditions provides information on possible improvements to be brought to the IRI foF2 modeling.Specifically, panels (b) in Figures 8, 10, and 12 highlight a persistent diurnal variation of IG12 eff for quiet conditions, with higher values during daytime and lower during nighttime.Moreover, IG12 eff values also show variations at the scale of some hours associated to foF2 variations in this time range, which are not described by the climatological IRI.This indicates that there is room for improvement in the description of the foF2 diurnal variation in IRI.
The validation analysis highlighted that the results obtained applying the SIRMUP method are clearly the best ones, even when the model does not work at its best, namely, assimilating only one foF2 autoscaled value, instead of the two possibly available.This would suggest preferring it to the procedure proposed here.Nevertheless, it is worth highlighting some aspects for which in the end it is the IRI UP-IONORING method to be preferred for Space Weather purposes.The best performance of the SIRMUP model has to be ascribed to the fact that it predicts foF2 after ingesting the same kind of data, namely, autoscaled foF2 values.This is undoubtedly a great advantage for SIRMUP compared to the IRI UP-IONORING method, which instead nowcasts foF2 through Equation 1, based on vTEC values coming from the IONORING maps.There are two significant reasons for which the IRI UP-IONORING procedure should be chosen as the operational one for the Italian territory.First, differently from SIRMUP, besides providing nowcasting foF2 maps, it can also provide nowcasting three-dimensional electron density matrices.This is made possible thanks to the underlying IRI model for which the F2-layer peak is an anchor point for the entire ionosphere vertical profile.In particular, by feeding IRI with the IG12 eff map provided by IRI UP-IONORING, the F2-layer peak description by IRI will be updated in terms of foF2, with effect on the entire vertical profile.This kind of output is valuable as input for ray tracing tools (e.g., Pietrella et al., 2023) to evaluate the reliability of point-to-point communications, which is so crucial especially in a Space Weather context.Second, the IRI UP-IONORING procedure relies on a number of assimilation points much larger than those on which the SIRMUP model relies.If the two ionospheric stations of Rome and Gibilmanna simultaneously do not provide an autoscaled value of foF2, the SIRMUP model cannot generate any output, which is undoubtedly a limitation, especially during the occurrence of Space Weather events.In fact, it is well known that in conjunction with significant Space Weather events, autoscaled values of foF2 might be unavailable, due for instance to a large absorption suffered by the signal sent by the ionosonde (e.g., Fagundes et al., 2020).This feature cannot be ignored when thinking about a future Italian Space Weather service or, more in general, about a future Italian natural hazard alert service.
The results shown in the previous section highlight also that assimilating Ciraolo calibrated values guarantees slightly better performance than that of the IRI UP-IONORING procedure under peculiar conditions, that is, nighttime stormy ones.Unfortunately, currently the Ciraolo vTEC calibration method cannot be applied in real time because, to perform a reliable calibration, it needs complete arcs of GNSS observations.At the same time, the results also showed that the IRI UP-IONORING procedure and IRI UP-Ciraolo one, although based on two different calibration methods, are more or less comparable, except for the underestimation characterizing IRI UP-IONORING for disturbed nighttime periods.We delved into this by considering a number of different disturbed periods, to understand how to mitigate its occurrence.From a preliminary analysis (not shown here), it seems that a promising way is to fine tune the estimation of the vTEC maps, which led to a slight underestimation of the nighttime ionospheric behavior during severe Space Weather events.This feature is currently under investigation and will be addressed in the future development of the operational IRI UP-IONORING algorithm.The validation at different sites highlighted that both IRI UP implementations show similar performance in spite of the location.This testifies the reliability and robustness of both the IONORING vTEC mapping procedure and the IRI UP Kriging spatial interpolation for IG12 eff mapping.
A careful analysis of semivariance plots and corresponding variogram models, for the two considered validation periods, highlighted that IRI UP-IONORING shows null nugget values for more than 95% of cases.Since IONORING vTEC maps are spatially sampled at 0.5°, a null nugget value evidences how assimilated vTEC values show no marked spatial variations below about 0.5°.This is probably to be ascribed to the LOWESS interpolation technique used by IONORING because corresponding variograms obtained from vTEC values calibrated by any single GNSS receiver of the RING network, as in IRI UP-Ciraolo, show variations at a scale below 0.5°.Although ionospheric variations at such small scales are generally less relevant than large-scale variations for radio propagation purposes, they can intensify under severe Space Weather events.This is why further developments are needed to improve the description of ionospheric variations at small scales through the IONORING vTEC mapping procedure.
A feature of the IRI UP-IONORING procedure is that it relies on statistical relationships between vTEC and NmF2, which work very well during daytime but show some limitations during nighttime and at the solar terminator hours.This is not surprising, as it was recently highlighted by Pignalberi, Nava, et al. (2021) and

Space Weather
10.1029/2023SW003838 PIGNALBERI ET AL.Pignalberi et al. (2022) in their works focused on the ionospheric equivalent slab thickness, which is defined as the ratio between vTEC and NmF2.In fact, it is well known that remarkable differences arise between NmF2 and vTEC during nighttime due to the plasma interchange between the ionosphere and the plasmasphere along the geomagnetic field lines (Klimenko et al., 2015;Pezzopane et al., 2019;Yizengaw et al., 2008).Moreover, at solar terminator hours, due to a different solar illumination at different heights, NmF2 and vTEC show a time lag which is well reflected in the slab thickness.At dawn, higher values of the slab thickness mean that NmF2 is low, while vTEC is increasing due to the solar radiation affecting mainly the topside ionosphere; the reverse happens at dusk hours.Since the relation (Equation 1) between NmF2 and vTEC has a similar meaning to the slab thickness, many of the considerations made by Pignalberi, Nava, et al. (2021) and Pignalberi et al. (2022) on the slab thickness apply also to Equation 1.A possible future upgrade for the IRI UP-IONORING algorithm concerns the improvement of the estimation of coefficients in Equation 1, especially during nighttime and solar terminator hours, aimed at improving the relationship between vTEC and NmF2.
This study evidenced the reliability of the IRI UP-IONORING procedure, and a corresponding alert service concerning the actual ionospheric conditions over the Italian territory, based on corresponding foF2 values, could be set up, increasing the portfolio of ionospheric real-time products available in the eSWua web portal.However, the IRI UP method implementation here presented can be profitably applied to regions different from Italy.In fact, the statistical relationships between vTEC and NmF2 can be derived from any co-located ionosonde and groundbased GNSS receiver with enough simultaneously measured observations (at least one solar cycle long time series).Pignalberi et al. (2022) showed that, so far, at least 32 co-located ionospheric stations are available for different latitudes and longitudes.This allows applying Equation 1 at different locations.By considering the current large diffusion of many vTEC mapping procedures and vTEC map providers, both on a regional and global scale, it is understood how the procedures developed for IRI UP-IONORING can be applied everywhere.Moreover, the use of effective indices (IG12 eff ) with corresponding mapping through Ordinary Kriging with stationary variograms, along with the ingestion of effective indices maps in the IRI background model, makes the application of IRI UP potentially scalable on a global basis.
Based on the investigations by Pignalberi et al. (2019), the following relationship log 10 NmF2 = a NmF2 log 10 vTEC + b NmF2 , (1) turned out to be a robust indicator of the relation between vTEC and NmF2 at a selected location.In such formulation, the base-10 logarithm of both NmF2 and vTEC are linearly related, and a NmF2 and b NmF2 are the slope and intercept coefficients, respectively.In Equation 1, NmF2 is in electrons/m 3 and vTEC is in electrons/m 2 .In this study, NmF2 and vTEC values are those recorded by co-located ionosondes and ground-based GNSS receivers over the Italian territory, that is, the pairs Rome-INGR and San Vito-MATE (see Sections 2.1 and 2.2).Given the limited spatial extension of the Italian territory, Equation 1 has been applied to the data set obtained by joining the two pairs of co-located ionospheric stations.As a consequence, a NmF2 and b NmF2 coefficients obtained in this study are representative of the ionosphere behavior over the Italian region.

Figure 3 .
Figure 3. (Panel (a)) Real-time vTEC map from IONORING for the period 28 March 2022 at 00:30 UT. (Panel (b)) vTEC values assimilated by IRI UP from the IONORING map in panel (a).

Figure 5
Figure5collects the values of the statistical metrics described by Equations 2-5, calculated for each LT and month of the year, in the same matrix layout as Figure4.Mean of residuals (panel (a)) range between about 0.015 and 0.050 MHz, testifying a remarkable accuracy in reproducing foF2 through vTEC and Equation 1 supplied by coefficients of Figure4.Since the mean of residuals is almost null, RMSE (panel (b)) provides information on the dispersion of residuals, that is, the precision in reproducing the measured foF2.RMSE values range between about 0.3 and 1.0 MHz with the lowest values during nighttime in winter and highest around noon in winter.Overall, RMSE is characterized by a diurnal pattern with low values during nighttime and high values during daytime, with the notable exception of the summer months where the diurnal variation is less marked.Since RMSE is critically affected by the magnitude of foF2, we calculated the NRMSE which provides a normalized parameter to evaluate the precision of the method.Panel (c) shows the matrix of NRMSE values as percentages.NRMSE ranges between about 8% and 16%, with maxima at sunrise and sunset for winter months, and at sunset for summer months.Hence, percentually, the lowest error in reproducing the measured foF2 is made during

Figure 4 .
Figure 4. (Panel (a)) Matrix of the slope coefficient (a NmF2 ) of the linear relation (1) between log 10 (NmF2) and log 10 (vTEC), as a function of LT (x-axis) and month of the year (y-axis).(Panel (b)) Matrix of the intercept coefficient (b NmF2 ) of the linear relation (1) between log 10 (NmF2) and log 10 (vTEC), as a function of LT (x-axis) and month of the year (y-axis).

Figure 5 .
Figure 5. (Panel (a)) Mean of residuals between modeled (foF2 model ) and measured (foF2 ionosonde ) values of foF2, as a function of LT (x-axis) and month of the year (yaxis).(Panel (b)) Same as panel (a) but for RMSE.(Panel (c)) Same as panel (a) but for NRMSE.(Panel (d)) Same as panel (a) but for Pearson correlation coefficient values.Measured foF2 values are those from Rome and San Vito ionosondes (see Section 2.1).Modeled values are those obtained by using the vTEC values measured at co-located GNSS stations (see Section 2.2) in Equation 1 with coefficients given by Figure 4.

Figure 6 .
Figure 6.(Panel (a)) Experimental variogram plot for IG12 eff for the assimilation time 28 March 2022 at 00:30 UT.Red points represent the calculated experimental semivariance values (y-axis) as a function of the distance between pairs of assimilated stations (x-axis).Blue line is the best-fitted variogram model, which is the Power one in this case.(Panel (b)) Map of IG12 eff values obtained from the variogram in panel (a) after applying the Ordinary Kriging interpolation.The color inside each circle represents the value of IG12 eff at the corresponding point.

1. 28
March 2022-4 April 2022: quiet to moderately disturbed period, the Dst index reached the minimum value of 46 nT, the Kp index reached the maximum value of 5 , with the F10.7 ranging between 128.0 and 239.0 sfu.The latter extremely high F10.7 value is associated to two powerful M9.67 and X1.3 solar flares occurred on 30 and 31 March 2022, respectively; 2. 21 March 2023-29 March 2023: quiet to severely disturbed period, the Dst index reached the minimum value of 163 nT, the Kp index reached the maximum value of 8, with the F10.7 ranging between 147.3 and 159.4 sfu.Dst, Kp, and F10.7 indices were downloaded from the OMNIWeb Data Explorer NASA portal at https:// omniweb.gsfc.nasa.gov/form/dx1.html.

Figure 7 .
Figure 7. (Panel (a)) Map of foF2 background values by IRI.(Panel (b)) Map of updated foF2 values by IRI UP after assimilating the IG12 eff map of Figure 6b.(Panel (c)) Map of residuals between IRI UP (panel (b)) and IRI (panel (a)) foF2 values.The assimilation time is 28 March 2022 at 00:30 UT.

Figure 8 .
Figure 8. (Panel (a)) Time series of foF2 values measured at the Rome ionospheric station (in black) for the period 21-29 March 2023, from DOY 80 to 89 of the year 2023.foF2 values as modeled by IRI-2020, by IRI UP-Ciraolo, by IRI UP-IONORING, and by SIRMUP, are represented in green, red, blue, and cyan, respectively.(Panel (b)) Time series of the IG12 value (green), used by climatological IRI, and of IG12 eff indices as provided by IRI UP-Ciraolo (in red) and IRI UP-IONORING (in blue) procedures, respectively.(Panel (c)) Time series of Kp (in black) and Dst (in magenta) geomagnetic activity indices.

Figure 9 .
Figure 9. (Top panels, a-d) Histograms of residuals between modeled and measured foF2 values at the Rome ionospheric station for the period 21-29 March 2023.Panel (a) for the climatological IRI-2020 model, in green; panel (b) for the IRI UP-Ciraolo procedure, in red; panel (c) for the IRI UP-IONORING nowcasting procedure, in blue; panel (d) for the SIRMUP nowcasting model, in cyan.Statistical parameters are provided in the upper left textbox.(Bottom panels, e-h) Corresponding scatterplots of measured and modeled foF2 values reported in top panels.The solid line represents the linear fit with corresponding coefficients readable in the upper left textbox.The black dashed line represents the first-third quadrant bisector.

Figure 10 .
Figure 10.(Panel (a)) Time series of foF2 values measured at the San Vito ionospheric station (in black) for the period 21-29 March 2023, from DOY 80 to 89 of the year 2023.foF2 values as modeled by IRI-2020, by IRI UP-Ciraolo, by IRI UP-IONORING, and by SIRMUP, are represented in green, red, blue, and cyan, respectively.(Panel (b)) Time series of the IG12 value (green), used by climatological IRI, and of IG12 eff indices as provided by IRI UP-Ciraolo (in red) and IRI UP-IONORING (in blue) procedures, respectively.(Panel (c)) Time series of Kp (in black) and Dst (in magenta) geomagnetic activity indices.

Figure 11 .
Figure 11.(Top panels, a-d) Histograms of residuals between modeled and measured foF2 values at the San Vito ionospheric station for the period 21-29 March 2023.Panel (a) for the climatological IRI-2020 model, in green; panel (b) for the IRI UP-Ciraolo procedure, in red; panel (c) for the IRI UP-IONORING nowcasting procedure, in blue; panel (d) for the SIRMUP nowcasting model, in cyan.Statistical parameters are provided in the upper left textbox.(Bottom panels, e-h) Corresponding scatterplots of measured and modeled foF2 values reported in top panels.The solid line represents the linear fit with corresponding coefficients readable in the upper left textbox.The black dashed line represents the first-third quadrant bisector.

Figure 12 .
Figure 12. (Panel (a)) Time series of foF2 values measured at the Gibilmanna ionospheric station (in black) for the period 21-29 March 2023, from DOY 80 to 89 of the year 2023.foF2 values as modeled by IRI-2020, by IRI UP-Ciraolo, by IRI UP-IONORING, and by SIRMUP, are represented in green, red, blue, and cyan, respectively.(Panel (b)) Time series of the IG12 value (green), used by climatological IRI, and of IG12 eff indices as provided by IRI UP-Ciraolo (in red) and IRI UP-IONORING (in blue) procedures, respectively.(Panel (c)) Time series of Kp (in black) and Dst (in magenta) geomagnetic activity indices.

Figure 13 .
Figure 13.(Top panels, a-d) Histograms of residuals between modeled and measured foF2 values at the Gibilmanna ionospheric station for the period 21-29 March 2023.Panel (a) for the climatological IRI-2020 model, in green; panel (b) for the IRI UP-Ciraolo procedure, in red; panel (c) for the IRI UP-IONORING nowcasting procedure, in blue; panel (d) for the SIRMUP nowcasting model, in cyan.Statistical parameters are provided in the upper left textbox.(Bottom panels, e-h) Corresponding scatterplots of measured and modeled foF2 values reported in top panels.The solid line represents the linear fit with corresponding coefficients readable in the upper left textbox.The black dashed line represents the first-third quadrant bisector.

Figure 14 .
Figure 14.(First row panels, (a), (d), and (g)) Maps of foF2 background values by IRI.(Second row panels, (b), (e), and (h)) Maps of updated foF2 values by IRI UP-IONORING.(Third row panels, (c), (f), and (i)) Maps of residuals between IRI UP-IONORING and IRI foF2 values.The assimilation time is 24 March 2023 at 02:00 UT for the first column panels, 24 March 2023 at 14:00 UT for the second column panels, and 24 March 2023 at 22:00 UT for the third column panels.Color bar scale changes across different columns.

Table 1
Ionosondes Over the Italian Territory The data set extension is specified only for the ionospheric stations used to derive the relations between vTEC and NmF2 (see Section 3.2).
PIGNALBERI ET AL.In this work, foF2 values recorded by Italian ionosondes were used twice.First, foF2 time series spanning many years have been used to study the statistical relation between NmF2 and vTEC, as detailed in Section 3.2.Second, real-time observations by ionosondes, which are not assimilated by IRI UP, have been used for validation purposes, as detailed in Section 4.

Table 2
GNSS Receivers Co-Located With Ionosondes and Used to Derive the Relations Between vTEC and NmF2 PIGNALBERI ET AL.
shows the variogram model in blue, which in this specific case the Power one; variogram parameters and statistics are also reported in the plot.Semivariance values modeled by the variogram are then used to infer the covariance values needed for the Kriging interpolation through the Lagrangian multiplier method, as detailed in Section 3.3 ofPignalberi, Pezzopane, et al. (2018).The final output of the Kriging interpolation method is the IG12 eff map over the considered region, as shown in Figure6b.It is worth highlighting that, in Figure6b, interpolated IG12 eff values exactly reproduce the IG12 eff values calculated at assimilated points (shown as circles) due to the very reliable fit of the Power variogram at small distances and the null nugget.
PIGNALBERI ET AL. is

Table 3
Statistical Results of the Validation Against foF2 Values Recorded by Ionosondes for the Testing Period 21-29 March 2023