Space Weather Effects Observed in the Northern Hemisphere during November 2021 Geomagnetic Storm: The Impacts on Plasmasphere, Ionosphere and Thermosphere Systems

: On 3 November 2021, an interplanetary coronal mass ejection impacted the Earth’s magne-tosphere leading to a relevant geomagnetic storm (Kp = 8), the most intense event that occurred so far during the rising phase of solar cycle 25. This work presents the state of the solar wind before and during the geomagnetic storm, as well as the response of the plasmasphere–ionosphere–thermosphere sys-tem in the European sector. To investigate the longitudinal differences, the ionosphere–thermosphere response of the American sector was also analyzed. The plasmasphere dynamics was investigated through ﬁeld line resonances detected at the European quasi-Meridional Magnetometer Array, while the ionosphere was investigated through the combined use of ionospheric parameters (mainly the critical frequency of the F2 layer, foF2) from ionosondes and Total Electron Content (TEC) obtained from Global Navigation Satellite System receivers at four locations in the European sector, and at three locations in the American one. An original method was used to retrieve aeronomic parameters from observed electron concentration in the ionospheric F region. During the analyzed interval, the plasmasphere, originally in a state of saturation, was eroded up to two Earth’s radii, and only partially recovered after the main phase of the storm. The possible formation of a drainage plume is also observed. We observed variations in the ionospheric parameters with negative and positive phase and reported longitudinal and latitudinal dependence of storm features in the European sector. The relative behavior between foF2 and TEC data is also discussed in order to speculate about the possible role of the topside ionosphere and plasmasphere response at the investigated European site. The American sector analysis revealed negative storm signatures in electron concentration at the F2 region. Neutral composition and temperature changes are shown to be the main reason for the observed decrease of electron concentration in the American sector.


Introduction
Geospheric storms are of great interest not only from a scientific point of view but also for their impact on technological systems.Aviation is in many ways dependent on technology that is prone to space weather conditions causing strong spatiotemporal variations in the geospace environment, especially during storms, endangering the health of crews and passengers, also exposed to radiation issues.As a consequence, there is an increasing necessity to improve the nowcasting and forecasting capabilities related to space weather services, e.g., the one that has been recently set up to support the International Civil Aviation Organization (ICAO) [1,2].Although we can infer the general features, the intensity and peculiarities of each process impacting, modulating, and controlling the Solar Wind-Magnetosphere-Ionosphere-Thermosphere (SW-M-I-T) coupled system are not known a priori and depend on the particular geophysical situation.As a consequence, the need of studying the single event and the level of geo-effectiveness in the SW-M-I-T system and on the geomagnetic field is of paramount importance for the community.
During impulsive solar storm events, a great amount of coronal plasma is released outward in the interplanetary medium forming the Interplanetary Coronal Mass Ejection (ICME).ICMEs are the counterpart of CMEs that move out of the solar atmosphere, and can induce strong geomagnetic storms when they interact with the Earth's magnetosphere [3].CME events are typically observed during the solar-maximum phase, with most intense ones linked to active solar regions, but are not completely absent during the solar-minimum phase, with an occurrence (daily average) rate of ∼1 and ∼4, respectively [4,5].Experimental observations by using the LASCO experiment on the SoHO satellite proved that more than 70% of CMEs occurring during 1996-2002 was originated by prominent eruptions [6], while a smaller amount was associated with solar flare events strongly connected with active regions.Nowadays, it is well established that CMEs and solar flares can occur independently; nonetheless, the most intense ICMEs are associated with jointly CMEs/solar flare events.These transient energetic events strongly modify the interplanetary medium conditions, as testified by the earlier investigations conducted by Sheeley et al. [7] and by Bothmer and Schwenn [8] by using Helios 1 and P78-1 satellites.One of the main effects is an enhanced solar wind (SW) speed from hundreds to thousands of km/s with a consequent formation of an interplanetary (IP) shock with most IP shocks (∼72%) associated with CMEs [7].
Superimposed to the existing ordinary SW, the ICME can intercept the Earth's magnetosphere, transferring a small part of its energy through several mechanisms, with great efficiency occurring under negative (i.e., southward) vertical component B Z (<0) of the interplanetary magnetic field (IMF): consequently, the magnetosphere responds through activations or enhancements of its current systems.It is a matter of fact that the major effects of SW-magnetosphere interaction manifest at high latitudes as an intensification of geomagnetic activity through different mechanisms: one of them is represented by the superposition of SW-related directly driven and the tail-related (substorm) unloading components (for a review, Akasofu [9]).Another important effect of SW-magnetosphere coupling, mainly occurring during geomagnetic storms, is represented by the modification and intensification of Field-Aligned Currents (FACs): this effect leads to a sudden increase in the Joule heating effect in the upper atmosphere, through collisional processes, which in turn could drive traveling atmospheric disturbances (TADs), possibly leading to the formation of traveling ionospheric disturbances (TIDs) (see, e.g., Zakharenkova et al. [10]; Priyadarshi et al. [11] and references therein).In this regard, the magnetosphere system acts as a mediator converting and transferring a portion of the ICME energy to the Earth's ionosphere and thermosphere.As summarized by Di Mauro et al. [12], those mentioned mechanisms represent powerful manifestations of the SW-magnetosphere interaction, which could potentially affect the performance of space and ground infrastructures by disrupting operations and communications in multiple sectors of society.Economic consequences involved by those extreme events are reported by National Research Council [13] with an estimation of cost for damages or disruptions of tens or even hundreds of billions of Euros (see also MacAlester and Murtagh [14]).In addition, Vermicelli et al. [15] reported that the costs associated with the risks posed to critical space-borne and ground-based technologies by upper atmospheric events are high and comparable with those of other natural hazards, such as tsunamis, earthquakes, or floods.One of the results of the geomagnetic activity enhancements is the changing of the magnetospheric cold plasma distribution, concentration, and composition.Besides the direct impact on the magnetohydrodynamics (MHD), wave properties, and their propagation into the magnetosphere, the cold plasma characteristics can contribute to the growth of electromagnetic ion cyclotron (EMIC) waves and hence to the precipitation of energetic electrons into the atmosphere [16].
This plasma population (E∼ 1 eV) mainly originates in the ionosphere and usually extends from 1000 km to 4-6 Earth radii (R E ), surrounding the Earth with a torus-shaped co-rotating plasma known as the plasmasphere.The outer boundary, the plasmapause or plasmasphere boundary layer (PBL, [17]), changes its shape in response to the geomagnetic activity and is often observed as a sharp decrease of the electron number density.Due to the presence of heavy ions, and the different radial distribution of each ion species, the mass density decrease inside the PBL can be less dramatic or even absent [18].During prolonged quiet times, the PBL can actually vanish also in the electron number density, being replaced by a smooth decrease.Beyond the PBL, there is the plasmatrough, a region of tenuous plasma that streams towards the dayside magnetopause, controlled by the magnetospheric convection electric field.The plasmasphere and the ionosphere continuously exchange particles by ambipolar diffusion along the field lines [19].The outflow of plasma from the F2 region into the plasmaspheric magnetic tube of force is due to increased plasma pressure that depends on plasma density and temperature, during daytime at F2-region heights, while plasma pressure is decreased in the magnetic tube due to the outflow of plasma from the tube to the ionosphere during the previous night.During geomagnetic storms, the external part of the plasmasphere is eroded due to the enhancement of the convection electric field and the PBL can become highly structured (e.g., Sandel et al. [20]).During the recovery phase, the refilling process of the plasmasphere from the ionosphere is enhanced and the depleted flux tubes slowly recover to the pre-storm conditions.Traditionally, F2layer storm effects are considered in terms of "negative" and "positive" phases according to whether electron density in the F2-layer maximum is increased or decreased with respect to monthly median or quiet time pre-storm values.In particular, plasmaspheric depletions visible as density reduction by a factor of 2-3, have been reported concurrently with negative storms (e.g., Wang et al. [21]).Physical mechanisms forming both negative and positive F2-layer disturbances are known.They are related to global thermospheric circulation, neutral composition and temperature, electric fields, and plasmaspheric flux changes.The list of all related processes may be found in Rishbeth [22] and Prölss [23].The intensity of each process controlling the F2 region is not well known for each particular geomagnetic storm, so it is difficult to evaluate the ionospheric response.
Moreover, depending on prehistory and current state of the magnetosphere and thermosphere, the response will be different to the same forcing from above, but at present, a routine thermosphere and magnetosphere monitoring is absent.Morphological studies have shown that mid-latitude negative F2-layer perturbations are usually observed to follow magnetic storms started during the preceding night, while positive storm effects are generally associated with increased geomagnetic activity in the local daytime sector (see, e.g., Appleton and Piggot [24]; Jones [25]; Prölss [26]).Unlike negative storm commencements, which are the most frequent in the post-midnight-early-morning Local Time (LT) sector, and which are rare in the noon and afternoon hours [27], positive disturbances may start in any LT sector (e.g., Mikhailov et al. [28]).
During the year 2021, the Sun was in the ascending phase of its 25th activity cycle, therefore, CME events are less probable.Nevertheless, on 2 November 2021 at ∼03.00 UT, a halo CME originated from an M1.7-class solar flare in the sunspot region labeled AR2891 (https://spaceweather.com/archive.php?view=1&day=02&month=11&year=2021, https: //blogs.nasa.gov/solarcycle25/2021/11/(accessed on 7 November 2022).It is the first intense solar storm event observed in the current solar cycle 25 that led to a geomagnetic storm, with large ionospheric disturbances observed in the northern hemisphere.Based on the NOAA Space Weather Scales (https://www.swpc.noaa.gov/noaa-scales-explanation(accessed on 7 November 2022), this storm was classified as G3 (or Strong), with the geomagnetic Kp index reaching a maximum value of 8− on 4 November.The aims of our analysis may be summarized as follows: 1.
Investigating the main SW/IMF-ICME-related characteristics and the corresponding geomagnetic activity at high latitude to study different phases of the polar cap electrodynamic and related Joule heating effects; 2.
Investigating the level of geo-effectiveness of the storm, its impact on the plasmasphere and ionosphere dynamics, and their possible coupling through the joint use of ionosonde, GNSS receivers, and magnetometers; 3.
Retrieving neutral composition ), exospheric temperature T ex , and vertical plasma drift mainly related to thermospheric winds for the periods of interest using ground-based ionosonde observations, and to compare the retrieved parameters with modern empirical thermospheric models of neutral composition predicted for the period under study; 4.
Discussing the role of thermospheric parameters in the formation mechanism of the storm for the analyzed disturbed periods.

Data and Methods
The data and methods used to achieve the aims of this work are described in the following sections arranged by specific topics or studied environment.

Solar Wind Parameters and Geomagnetic Field Measurements
The solar wind parameters and the interplanetary magnetic field are studied through upstream satellite datasets provided by CDAWeb.In this work, we used the 1-min OMNI data, time-shifted to the bow shock nose.On the ground, the geomagnetic activity was monitored by using Kp, Ap, and SymH indices, all derived by ground-based measurements.In the northern hemisphere, we estimated the polar cap electrodynamic activity through the Weimer [29,30] empirical model.In this work, we used 45-min moving-average SW and IMF data for estimating the high latitude electrodynamic conditions.This has been proven to be of support in identifying the possible driving processes leading to changes in the neutral wind circulation and possibly triggered by ionospheric heating at high-latitude [31].At low latitudes in the European sector, the geomagnetic activity is monitored by using the 1 min magnetic field variations collected at the Italian geomagnetic observatories of Duronia (IAGA code DUR, coordinates: 41.65 • N, 14.47 • E) and Lampedusa (LMP, coordinates: 35.52 • N, 12.55 • E), both run by Istituto Nazionale di Geofisica e Vulcanologia (INGV), represented as green dots in Figure 1.Hereafter, where not specified, the coordinates of any sites are reported in the geographic reference frame.

Plasma Mass Density Sounding
For investigating the plasmasphere dynamics, we used equatorial plasma mass density estimates (ρ eq ), remotely inferred by using magnetic observations recorded by the European quasi-Meridional Magnetometer Network (EMMA, [32]).The procedure used to derive ρ eq is extensively described by Del Corpo et al. [33,34].The core of the procedure is the determination of the geomagnetic field line resonance (FLR) frequencies by applying the gradient method [35,36] to 21 EMMA station pairs (the stations used are represented as black squares in Figure 1).The frequency is determined for the field line passing through the midpoint of each station pair by performing a 2-h cross-spectral analysis of the northsouth magnetic component detected simultaneously by the two stations.The mass density is then evaluated by solving numerically the governing MHD equation [37], using the TS05 magnetic field model [38] and assuming a power law dependence of the mass density along the field line given by: where r is the geocentric distance, and the suffix "eq" refers to equatorial values.The analysis was performed with a time step of half an hour.The equatorial distance determined for each station pair can vary with time, especially at high latitudes.To determine the plasma densities at fixed radial distances, the radial profiles with at least 3 data points were fitted with a smoothing spline [33].Fitted portions related to intervals larger than 2 R E in which experimental points were missing, were excluded in the analysis.In solving the MHD equation, the footprints of the field lines are considered fixed, that is, the ionospheric conductivity is assumed infinite.This approximation, valid during daytime hours, is not necessarily adequate during nighttime, when the conductivity can decrease by several orders of magnitude and free-end mode configuration of the field lines could be a better solution of the MHD equation [39].This makes the estimates obtained during nighttime somewhat questionable.Nonetheless, Thakahasi et al. [40] reported FLR detections across the midnight sector compatible with a fixed-end mode configuration; so, we decided to keep also the density estimates during nighttime, being aware that they are less reliable.

Ionosonde Measurements
In this analysis, ionosonde data from five stations are used: Juliusruh (54.6 • N, 13.  [41], and provide data in real-time (Upper atmosphere physics and radiopropagation Working Group, 2020a) available in the electronic Space Weather upper atmosphere (eSWua) data portal (http: //eswua.ingv.it).The ionosonde in Moscow is an AIS PARUS ionosonde, while those in Pruhonice and Millstone Hill are DPS-4D ones (Reinisch et al., 2008).Moscow (geomagnetic latitude Φ = 52.5 • N), Juliusruh (Φ = 51.4• N), and Millstone Hill (Φ = 51.0 • N) have close geomagnetic latitudes and should manifest similar variations both in ionospheric and thermospheric parameters, while Rome (Φ = 36.3• N) is a lower latitude station, and storm variations in all parameters are expected to be different.As a reference to identify the quiet behavior of the F2-layer peak density (foF2) during the storm, we have also evaluated foF2med, i.e., the median foF2 in the 27 days preceding 3 November 2021.

Total Electron Content Measurements
Four ground-based GNSS receivers (red dots in Figure 1) have been used to evaluate the TEC.The receivers have been chosen to be as close as possible to the ionosondes.The receiver in Rome (id: INGR, coordinates: 41.83 • N, 12.51 • E) is located in the headquarter of the INGV; it is part on the Italian "Rete Integrata Nazionale GNSS" (INGV RING Working Group [42]), and it is co-located with the Rome AIS-ionosonde [43,44].The INGR receiver is part of the IONORING system [45], able to provide real-time monitoring and modeling of the ionospheric Total Electron Content (TEC) over Italy (Upper atmosphere physics and radiopropagation Working Group [44]), which is again available in the above-mentioned eSWua website.The receivers in Mendeleevo (id: MDVJ, coordinates: 56.02 • N, 37.21 • E), in the Geodetic Observatory Pecny (id: GOPE, coordinates: 49.54 • N, 14.47 • E), and in Sassnitz Island of Ruegen (id: SAS2, coordinates: 54.41 • N, 13.64 • E) are part of the EUREF Permanent GNSS Network.Those geodetic receivers provide data in RINEX format with a cadence 30 s maximum.By applying the TEC calibration by Ciraolo et al. [46], it is possible to derive the slant ionospheric TEC for each receiver-satellite couple pair.Additionally, the technique is able to estimate the vertical TEC (vTEC) over the given receiver.The details of its retrieval are given in Cesaroni et al. [45] and Tornatore et al. [47].An elevation mask of 20 • is considered for the GNSS observation.As a reference to identify the vTEC behavior during the storm, we have also evaluated vTECmed, i.e., the median vTEC in the 27 days preceding 3 November 2021.

Method to Retrieve Thermospheric Parameters
The state of the mid-latitude ionosphere reflects the state of the surrounding neutral atmosphere and the intensity of the solar ionizing radiation.Therefore, in principle, it is possible to solve the inverse problem of aeronomy to retrieve basic thermospheric parameters (neutral composition, temperature, wind) and solar EUV ionizing flux from ionospheric observations.This has opened a way to monitor the state of the thermosphere using the world-wide ionosonde network observations [48].Recently, a new method has been proposed which uses plasma frequency f p at 180 km and foF2 ionosonde observations [49]; it is simpler than the basic version [28] and can be used routinely for automatic operations.The method is based on solving an inverse problem of aeronomy using observed noontime foF2 and five time instants at (10,11,12,13,14) LT values of plasma frequency fp180 at 180 km height as input information along with standard indices of solar (F 10.7 ) and geomagnetic (Ap) activity.An advanced version of our method additionally uses hmF2 autoscaled from ionograms as fitted parameters.The list of the retrieved parameters includes: neutral composition (O, O 2 , N 2 concentrations), exospheric temperature T ex , the total solar EUV flux with λ ≤ 1050Å, and vertical plasma drift W mainly related to the meridional thermospheric wind, V nx .The empirical model Mass-Spectrometer-Incoherent-Scatter (MSISE00) [50] is used for comparison.

Experimental Results
This section presents the main results of our investigations.Section 3.1 provides an overview of impacting ICME characteristics, together with geomagnetic activity indices, polar cap electrodynamic effects and geomagnetic field variations at the Italian territory.The dynamics of the plasmasphere and the ionosphere are described in Sections 3.2 and 3.3, respectively.Section 3.4 presents the results of the thermospheric parameter analysis.

The SW, IMF Conditions, and Geomagnetic Activity from High to Low Latitudes during 2-5 November 2021
Figure 2 shows (panels a-d) the IMF total strength B and B X , B Y , B Z components in the geocentric solar ecliptic (GSE) reference frame, SW velocity V SW , plasma (mainly protons) density n, and temperature T during 2-5 November 2021.On 3 November 2021 at ∼20:23 UT, an IP shock was registered at Earth's orbit (it is marked with a vertical orange stripe in panels a-d): it can be classified as a Fast Forward (FF) IP shock since all four above-mentioned parameters increased around their respective maximum rate.Following Paschmann and Daly [51], the minimum variance analysis (MVA) of magnetic field data method was applied for calculating shock properties: we found, in the GSE reference frame, that the IP shock impacted on the magnetopause with an angle of ∼30 • with respect to the Sun-Earth line with a velocity of ∼620 km/s; the shock normal angle of ∼19 • , i.e., the angle between the average B upstream and IP shock front normal direction n = (−0.85,0.52, 0.04), is also computed (bold characters indicate vectors).These results suggest, according to Sheeley et al. [7], that this kind of IP shock is probably related to the CME emitted by the Sun: indeed, on 2 November 2021, a halo CME was jointly originated with an M1.7-class solar flare occurring at ∼3 UT in the sunspot region labeled AR2891 (https://spaceweather.com/archive.php?view=1&day=02&month=11&year=2021, https://blogs.nasa.gov/solarcycle25/2021/11/(accessed on 7 November 2022).
With the arrival of the IP shock at Earth's orbit, the geomagnetic activity quickly enhanced, as it can be seen from both Kp and SymH indices (Figure 2g-h).In the northern hemisphere, we estimated the average Joule heating J (Figure 2f) and electric potential E (Figure 2e) by using the Weimer [29,30] polar cap electrodynamic model, hereafter indicated as W05 model.J was calculated as in Spogli et al. [31] and Alfonsi et al. [52] but with some updates in order to provide a physical meaning to the resulting J: it was evaluated on a regular grid with a discrete cell resolution of ∼2 and 4.5 degrees along latitude (in the range 0-90 • ) and longitude (in the range 0-360 • ), respectively, each one characterized by a discrete surface δS; then, the resulting J was averaged over the whole spherical cap surface S. With the same just mentioned procedure E was computed and, in addition, it was divided into positive (E + ) and negative (E − ) electric potential and averaged over the corresponding surface.E can be regarded as indirectly (model-derived) related with AL (E − ) and AU (E + ) auroral geomagnetic indices, respectively, (at present, not freely available from the World Data Center for Geomagnetism Kyoto official website https://wdc.kugi.kyoto-u.ac.jp (accessed on 7 November 2022)), while J can be regarded as a possible driver of TADs.It can be seen that J attained high values on 3 November 2021 at ∼21:15 UT just after the IP arrival and in correspondence with B Z < 0 that triggered the enhanced E − .This stage also corresponds to the Sudden Storm Commencement (SSC) whose knee occurred at ∼19:40 UT (vertical dashed lines) well identified by SymH and |∆H| (less by ∆D) at DUR and LMP (Figure 2i-l, solid and dashed lines) geomagnetic observatories; H and D are the magnetic field components in the horizontal local plane, along the local magnetic meridian (H) and perpendicular to it (D), respectively; ∆H and ∆D are obtained by subtracting from H and D their respective average values computed for the whole analyzed time interval.After this almost isolated peak, J increases with increasing both E + and E − , attaining high values around 07:00 UT and 12:00 UT, which corresponds to the main phase of the geomagnetic storm well identified by maximum values of Kp and |SymH|, as well as by maximum values in |∆H| at DUR and LMP.During the next hours of recovery phase, the B Z was positive for several hours till ∼03:30 UT on 5 November, then reached negative values at ∼14:40 UT on 5 November whose oscillations are well followed by ∆H at both observatories.After this stage of the recovery phase, both Kp and SymH approached their respective pre-storm values.Therefore, the geomagnetic storm can be considered as concluded.

Dynamics of the Magnetospheric Plasma
For describing the variations of the magnetospheric plasma density, it is useful to refer to quiet time conditions (see, e.g., Pezzopane et al. [53], Vellante et al. [54]) when the plasmasphere reaches a state of saturation and is in equilibrium with the underlying ionosphere [55].The period under investigation was preceded by several days characterized by a moderate geomagnetic activity in which the Kp index reached a maximum of 5− with an average of 2.3, so we extended our analysis back to 29 October, which is the last day of a very long period with low geomagnetic activity.The average Kp in the period 24-29 October was ∼0.6, and it dropped to ∼0.3 during 27-29 October.This makes the configuration of the plasmasphere on 29 October ideal to serve as a reference for the saturation conditions in the European local time sector.
Figure 3 shows, from top to bottom, the SymH index, the Kp index, and the equatorial plasma mass density variations inferred at 5.5, 4.5, 3.5, and 2.5 R E , respectively, from 29 October to 9 November.Red dots are the plasma mass densities inferred as described in Section 2.2, while the blue dotted curves are smoothing splines obtained from the observation points, drawn to guide the eye of the reader through the long-term variations of the plasma density.Shadowed areas indicate times in which one (light-gray) or both (dark-gray) footprints of the field line passing to the observation point, were not lighted by the Sun.On 29 October, the plasma density is relatively high at all the radial distances and changes very slowly during the day, suggesting a condition of equilibrium with regard to the ionosphere.This idea is confirmed by the radial profiles' inspection.Figure 4a shows the equatorial density estimates derived from the EMMA station pairs (open circles) and the related fitted radial profiles (solid lines), at 07:30 UT for selected days.On 29 October, (in black) the density profile decreases smoothly with the radial distance assuming the typical pattern of a saturated plasmasphere (e.g., Carpenter and Anderson [55]).The average plasma mass densities observed on 29 October at the four selected radial distances are reported as blue dashed lines in Figure 3. Starting from 30 October, the enhanced convection contributes to the erosion of the plasmasphere, visible in Figure 3 as a progressive reduction of the daily average level of the density (see blue dotted lines).The erosion is more pronounced at higher geocentric distances, and is less evident at 2.5 R E , where the density oscillates around 2 × 10 3 amu cm −3 .
Figure 4a shows more clearly how the erosion takes place in the 09:30 LT sector (07:30 UT).A plasmapause appears around 4 R E and progressively moves earthward until it reaches ∼ 2 R E during the main phase of the storm (blue radial profile).Back to Figure 3, we can see the effects of the ionospheric refilling superimposed to the plasmasphere erosion and visible as a monotonic increase of the density during daytime hours.Typical examples are the diurnal patterns of 2 and 3 November.Even if the coverage is poor, we can see traces of the plasmaspheric density reduction during the night between 2 and 3 November, mostly at 4.5 and 5.5 R E .That decrease could be due to the downstream of plasma towards the ionosphere or to the crossing of the magnetospheric regions still dominated by the convection electric field and characterized by a lower density.
During the main phase of the storm, we observe a decreasing trend at all the radial distances considered in Figure 3, with the exception of the lowest one.At the same time, if compared with the preceding day, the average density level decreases by a factor of ∼2 at 2.5 and 3.5 R E while it remains similar or even increases at outer magnetic shells.This could result from the combined effects of the plasmasphere depletion more evident in the inner shells, and the sunward surge, an effect that moves a large amount of plasma from the nightside to the dayside, which is more visible at the external shells.At 08:30 UT, near the end of the main phase, a sudden increase of the plasma density is clearly visible at all geocentric distances.Such an increase cannot be explained with a plasma exchange with the ionosphere, which generally takes place on a time scale of a day or more.It is most probably due to the crossing of a drainage plume which enters in the line of sight of EMMA as the network rotates with the Earth.Figure 4b shows the radial profiles at 07:30 (blue) and 10:00 UT (red) highlighting a consistent plasma density increase in the trough region between 2.5 and 7 R E , with a maximum around 4-5 R E .The slope of the radial profile is very similar to the saturated conditions, although the density assumes smaller values.
The erosion apparently continues until 6 November when the typical diurnal refilling pattern is visible at all geocentric distances and is followed by a slow recovery in the following days.The slow recovery is also confirmed by the radial profiles shown in Figure 4c.On 8 November, the plasma density at 3 and 5 R E was still at about 50% and 20% of the saturated conditions, respectively.

Morphology of the Ionospheric Storm
This section reports the relative behavior of the bottomside ionosphere, as probed by ionosonde data, and the integrated ionospheric behavior, as probed by GNSS receivers.As mentioned in Section 2 and reported in Figure 1, ionosondes and GNSS receivers have been selected with a criterion of co-location (or maximum proximity), to ensure the reference to the same ionospheric sector.Time series of foF2 (red dots) and of vTEC (black curve) are reported in Figure 5, separately for each ionosonde-GNSS receiver couple.In the same figure, the red and black dashed lines indicate the corresponding values of foF2med and vTECmed, respectively, and as defined in Sections 2.3 and 2.4.The time series focuses on the period 3-5 November 2021, during which J reached the highest values, i.e., more than 2.2 mW/m 2 and 3-h Kp index up to 8 (Figure 2f,g).
To further highlight the deviations from the quiet reference resulting in positive and negative phases of the ionospheric storm, Figure 6 reports the time profilesof the ratio vTEC/vTECmed (black) and of the same quantity evaluated for foF2 (red dots).The black shaded area indicates the [0.8; 1.2] variability range, corresponding to a reference quiet time variability.
The first splash of auroral activity at 20-21 UT on 3 November with the J increasing up to 1.8 mW/m 2 and Kp equal 6-between 21-23 UT.
The most striking feature of the relative foF2-vTEC behavior is the fact that, at the beginning of the ionospheric storm, foF2 at Moscow/MDVJ and Juliusruh/SAS2 differs from vTEC.Specifically, the double-peaked positive storm signature in vTEC at MDVJ has not the same clear signature in Moscow.The vTEC minimum between the two peaks at about 06:00 UT on 4 November fits with the quiet time behavior, while the corresponding foF2 values indicate a negative phase of the storm.The foF2 negative phase found in the range 00:00-14:00 UT at the Juliusruh ionosonde corresponds with a normal behavior of the vTEC.Again, for the northernmost couples, a clear signature of a negative storm (for both vTEC and foF2) is found starting from 14:00 UT on 4 November and lasting until the early hours on 5 November.Being located at almost the same latitude, but in a different longitudinal sector, the behavior of the northernmost ionosonde-GNSS receiver couples indicates a local effect of the ionospheric storm development and a possible role of the topside ionosphere and plasmasphere in the found mismatches between the vTEC and foF2 behavior.
By focusing on the 2 couples Juliusruh/SAS2 and Pruhonice/GOPE, a similar behavior is found, with less pronounced and less lasting negative phases, whenever found, with respect to Moscow/MDVJ.This is expected from ionospheric sectors at the same longitude but at slightly different latitudes in the mid-latitude sector.The behavior of the Rome/INGR couple significantly differs from the one of the couples in almost the same longitudinal sector but at higher latitudes.In fact, positive storm signatures on vTEC time series are found after 06:00 UT on 4 November and lasting almost all the rest of the day, while little deviations from the quiet behavior is found on foF2 data.
On 3 November, around 19:40 UT, SSC took place.After this, on 4 November at 02 UT (Figure 7) begins the negative phase until 10 UT of 5 November; after this, the conditions come back to median behavior.Interesting is the comparison between the foF2 variations at Millstone Hill and Juliuruh, Moscow stations located all at the same geomagnetic latitude.

Thermospheric Retrieved Parameters
Let us check if our method to extract thermospheric parameters confirms the presentday mechanism of midlatitude F2-layer storms.Unfortunately, we cannot follow the development of this process in time as we can analyze only daily (at 12 LT) variations.The main thermospheric parameters retrieved from the observed ionospheric observations are listed in the Table 1.Following the same criteria adopted for plasmasphere and geomagnetic conditions, we selected 29 October as reference day also for the ionosphere-thermosphere system (see Section 3.2).
The obtained results (Table 1) confirm the upsurge of the equatorward wind V nx (positive W = V nx sin I cos I, where I is the magnetic inclination) related to the increase of auroral heating on 4 November.According to the F2-layer storm mechanism (e.g., Prölss [56]), the enhanced equatorward wind is followed by the disturbed thermospheric composition with low [O]/[N 2 ] ratio.At the three higher latitude stations, our results indicate a [O]/[N 2 ] 300 decrease on 4 November, particularly evident at Juliusruh.Similar behavior is obtained using the MSISE00 empirical model [50].For all stations, the exospheric temperature (T ex ) increased, while [O] 300 is close to the reference day.
On 5 November, at all stations, strong poleward thermospheric winds can explain the quick recovery of thermospheric parameters.
For European stations, the difference between the [O]/[N 2 ] retrieved and predicted by MSISE00 on the disturbed day (4 November) is 3 and 2 times less with respect to the reference day.For the other days, the difference is less.
Looking at the thermospheric parameters on 3 November in MH, an increase of T ex , [O] and [N 2 ] with respect to the reference day (29 October) is seen.On 4 November, as for the stations in the European sector, there was an upsurge of the equatorward wind V nx , an increase of T ex , [N 2 ], and a decrease of [O].The ratio [O]/[N 2 ] 300 decreases strongly and this is seen in the negative ionospheric phase observed during all the day (Figure 7).On 5 November, a strong poleward thermospheric wind transferred neutral composition from lower latitudes to restore the normal (undisturbed) state of the thermosphere at higher latitudes.
The values of T ex and [N 2 ], calculated with the empirical model MSISE00, are underestimated with respect to the retrieved parameters, while [O] and [O]/[N 2 ] 300 values are overestimated.It should be stressed that an observed decrease in electron density cannot be explained with the thermospheric parameters predicted by the MSISE00 model.

Discussion
We analyzed the thermosphere-ionosphere-plasmasphere system response to the geomagnetic storm that occurred on 3-4 November 2021, during the rising phase of the current solar cycle 25 [57].A documented aspect of the ionosphere-plasmasphere coupling during disturbed times is the formation of storm-enhanced density (SED) in the sub-auroral ionosphere region that maps into plasmaspheric drainage plumes in the equatorial plane [58][59][60].From our plasmasphere analysis, we found evidence of a plume crossing around 10 UT (∼12 LT) on 4 November (Figures 3 and 4b).Before the plume interception, the equatorial plasma mass density profile suggests that the plasmapause was at approximately 2R E .This means that the ionosonde GNSS receiver at couple Rome/INGR were sounding regions magnetically connected with the plasmasphere, while the magnetic filed lines, linked to the other ionosonde receiver pairs, passed from the plasma trough to the plume.It is challenging to find a connection to the plume in the ionospheric data presented.However, MDVJ TEC time series shows an enhancement compatible with a SED (top panels of Figures 5 and 6) between 08:00 and 11:30 UT.The corresponding foF2 behavior does not report any significant deviation from the quiet behavior, indicating that the enhancement is likely occurring in the plasmasphere.Of course, we have no information on the topside ionosphere, which can contribute to the found differences between vTEC and foF2 time series, because unfortunately, no satellite observations from Swarm and GRACE over Rome were available on this day and during these hours.The same effect is not clearly visible in the SAS2 data (second panels of Figures 5 and 6), which is approximately at the same latitude of MDVJ, but separated in LT by about one hour.A possible explanation is that the drainage plume is fading during such time or co-rotating towards dusk exiting from the line of sight of SAS2.This idea would be supported by the absence of clear TEC enhancements also at GOPE that was in the same LT sector of SAS2.The positive storm visible in INGR data and in foF2 variations is more probably related to equatorward thermospheric wind.
Of particular interest is the plasmaspheric density fall off that occurred during the night between 5 and 6 November and is clearly visible in Figure 3 for all the magnetic shells presented.That decrease is very fast and hardly attributable to a downward flux towards the ionosphere.It is more likely due to a boundary crossing of the flux tubes monitored by EMMA which progressively moves to a region dominated by the convection electric field [16].This interpretation is compatible with the foF2/vTEC observations in Moscow/MDVJ, Juliusruh/SAS2, and Pruhonice/GOPE.In fact, starting from 18:00 UT on 5 November, the negative storm signatures are present mainly on vTEC data rather than on foF2.This indicates that the vTEC depletion in that part of the recovery phase is ascribed to phenomena occurring above the peak of the F2 layer.
One of the aims of this paper is the discussion of the role of neutral composition.In this regard, we remind the reader that the mechanism of ionospheric storms is related to changes in the thermospheric circulation with corresponding changes in neutral composition [23,56] and to the appearance of electric fields of magnetospheric origin [61].Daytime NmF2, in accordance with Rishbeth and Barron [62], closely follows the q(O + )/β variation at the hmF2 height.The q(O + )/β ratio is controlled mostly by the [O]/[N 2 ] ratio, and the relationship of NmF2 with [O]/[N 2 ] during geomagnetic storms is well established [23].
The main aeronomic parameters, responsible for the formation of daytime mid-latitude F layer, have been extracted from ionospheric observations, using the earlier developed method confined to middle latitudes and noontime hours [49].The method allows us to highlight how the main (controlling) aeronomic parameters have been changed for a particular geophysical situation using the observed changes in the electron concentration at F1 and F2 layers.
In particular, the event of 4 November clearly manifests variations in the Ionosphere-Thermosphere system in both longitudinal sectors during this storm period (Table 1).The thermospheric circulation was inverted during daytime hours (positive W corresponds to the equatorward V nx ).Exospheric temperature, T ex , was strongly elevated, more than the temperature predicted by MSISE00.The [O]/[N 2 ] 300 ratio decrease is mainly due to an increase of [N 2 ] related to the T ex increase.[O] 300 undergoes small variations.In fact, the atomic oxygen abundance in the thermosphere was strongly decreased but, due to a larger T ex , [O] small variations at 300 km are found.The found variations of thermospheric parameters confirm the expected response of the ionosphere and thermosphere to geomagnetic storms.It should be stressed that MSISE00 does not reproduce the retrieved thermospheric parameter variations during this storm event.
The effect of SSC on the American and European ionosonde stations is clearly seen looking at Figure 8.In this figure, the green lines represent hmF2 measurements during the quietest day of the month, which corresponds to 1 November 2021.On 3 November after the occurrence of SSC (∼19:40 UT), an hmF2 uplift at different stations in the American sector can be seen (Figure 8, left panels), the same behavior can be observed in the European sector (Figure 8, right panels).A different behavior of the Ionosphere-Thermosphere system can be observed analyzing the ionosonde data located in the American sector.On 4 November, there was a splash of Joule heating at 11:30 UT (Figure 2), as a reaction to this increased heating a TAD was formed which moved to the equator (Figure 9) and we see corresponding variations of hmF2 and NmF2.After TAD, we observe a decrease of electron density and The impact of the storm in the Ionosphere-Thermosphere system is very different in the two longitudinal sectors.Looking at Table 1, on 3 November, we can notice that the vertical plasma drift in Europe is between (−18.7,−23.8) m/s with respect to MH −16.5 m/s; this means that the poleward V nx is stronger in Europe than in America so that a storm-induced equatorward V nx is blocked to higher latitudes in Europe than in America.The consequence is that the disturbed neutral composition in the American sector can reach lower latitudes than in Europe.In addition, though MH, Juliusruh, and Moscow have the same geomagnetic latitude, MH is a "near-to-pole" station manifesting larger storm variation with respect to "far from pole stations" as in Europe [64].

Conclusions
The main results can be summarized as follows: 1.
The possible interplay between plasmasphere and ionosphere during the storm has been highlighted.This has been done by integrating the plasmaspheric information provided through EMMA, the bottom-side ionosphere information provided through ionosondes and the integrated ionospheric information through GNSS receivers.Specifically, we identified the possible presence of a SED in the sub-auroral ionosphere spanned by the Moscow/MDVJ ionosonde-GNSS receiver couple between 08:00 and 11:30 UT on 4 November.2.
Again, for the cross-talk between ionosphere and plasmasphere, during the recovery phase of the storm, in the range 18:00-24:00 UT on 5 November, we identified a vTEC depletion ascribed to phenomena occurring above the F2-layer peak in correspondence with a plasmaspheric density decrease, probably due to boundary crossing of the flux tubes monitored by EMMA moving into a region dominated by the convection electric field.

3.
In the European longitudinal sector, we found on 4 November an increase of T ex and a decrease of [O]/[N 2 ] at 300 km in the higher latitude stations (Juliusruh and Moscow) with upsurge of the equatorward wind V nx related to the increase of auroral heating.4.
Negative foF2 disturbances at Moscow and Juliusruh also take place during daytime hours (Figures 5 and 6), and this may be attributed to a low [O]/[N 2 ] 300 ratio.In Rome, the [O]/[N 2 ] 300 decrease is not as strong as at the other two stations, and also T ex and [O] are close to quiet time values, although there is a thermospheric equatorward wind that probably determines the positive F2-layer storm phase in foF2.

5.
At MH, the occurrence of a TAD is observed, which could be related to enhanced Joule heating occurring during the SSC around 19:40 UT, which determines NmF2 The stronger effect of the geomagnetic storm on the American sector with respect to the European one is due to the stronger poleward meridional thermospheric wind in Europe that blocks the disturbed neutral composition at higher latitude and also because MH is a "near-to-pole" station.7.
On the storm onset, the plasmasphere was already partially eroded due to the moderate geomagnetic activity in previous days.During the main phase, the enhanced convection electric field contributed to the formation of a drainage plume detected by EMMA.In the aftermath, the plasmasphere remains far from the saturated conditions, although a partial recovery is visible between 2.5 and 6R E .The procedure used to infer the plasma mass density provided some estimates also during nighttime, showing in some cases a trend compatible with the daytime profiles.These results could confirm that the fixed-end mode of toroidal standing Alfvèn waves can occur also during nighttime.Further and more accurate analysis, also using a finite ionospheric conductivity as boundary conditions for solving the MHD governing equation, could clarify the reliability of those points.

Figure 1 .
Figure 1.Position of the ground-based instrumentation used to perform the analysis in the European sector.Black squares are EMMA stations, blue dots are ionosonde stations, red dots are GNSS receivers and green dots are geomagnetic observatories.See text for further details.

Figure 2 .
Figure 2. The SW, IMF parameters and geomagnetic activity during 2-5 November 2021.(a) Interplanetary magnetic field components and strength; (b) radial solar wind speed; (c) solar wind plasma density and (d) temperature; (e) positive and negative average horizontal electric potentials in the northern hemisphere and (f) average Joule heating J; geomagnetic activity indices (g) Kp and (h) SymH; geomagnetic field fluctuations at 1 min time resolution at Italian observatories of (i) DUR and (l) LMP.The vertical red strip (a-d) refers to the FF IP shock.The vertical black dotted line (h,i,l) refers to the SSC.

Figure 3 .
Figure 3. From the top: the SymH index, the Kp index, and the equatorial plasma mass density inferred at 5.5, 4.5, 3.5, and 2.5 R E , respectively.Red dots are the observation points, while the blue dotted curves are smoothing splines highlighting the long-term variations.Blue dashed lines are the average plasma mass density values observed on the quiet day 29 October, and used as a reference for saturated conditions.Shadowed areas indicate times in which one (light-gray) or both (dark-gray) footprints of the field line passing to the observation point are not lighted by the Sun.

Figure 4 .
Figure 4. Radial plasma mass density profiles on selected days.(a) Plasmasphere erosion from 29 October to 4 November.(b) Possible plume crossing during the main phase of the storm on 4 November.(c) Plasmasphere partial recovery in the interval 4-8 November.All the profiles but the red one on panel b refer to 07:30 UT.The red profile is taken at 10:00 UT.Solid lines are the fit obtained from EMMA density estimates (open circles).Dashed curves are the fit portions obtained between EMMA points separated by a radial distance larger than 2 R E .Those parts of the fit are not used to evaluate the time series in Figure 3.

Figure 5 .
Figure 5.Time profiles of vTEC (black curve) and foF2 (red dots) for the considered couples ionosondes-GNSS receiver, during 3-5 November 2021.The black and red dashed lines represent the corresponding median values evaluated in the 27 days preceding 3 November 2021, and here used as the quiet reference.

Figure 6 .
Figure 6.Time profiles of the ratio between vTEC, during 3-5 November 2021, and the median values evaluated in the 27 days preceding 3 November 2021 (black curve) and of the same quantity evaluated for foF2 (red dots).Black shaded area indicates the [0.8; 1.2] variability range during quiet times.

Figure 7 .
Figure 7. foF2 observations at Millstone Hill in the American sector, during 3-5 November 2021.Upper panel: Observed hourly (black line) and median values evaluated in the 27 days preceding 3 November 2021 (green line) foF2 variations; bottom panel: foF2/median variation (black line), for the 3-5 November 2021 geomagnetic storm.Grey shaded area indicates the [0.8; 1.2] variability range during quiet times.

Figure 8 .
Figure 8. hmF2 observations, on 3 November 2021, in the American sector (left panels), at Millstone Hill, Boulder and Eglin, and in the European sector (right panels), at Moscow, Juliusruh, Pruhonice, and Rome.The reference values are shown as green lines (see text for details).
[O] both 1.7 times less with respect to the reference day while [N 2 ] increases.[O]/[N 2 ] is around six times less at 300 km.Perturbed neutral composition and temperature variations may have been transferred from high latitudes by equatorward wind during post midnight [63].Therefore, the main reason for the observed decrease in electron concentration at F2-layer heights on 4 November appears to be a strong decrease in the ionization rate due to [O] decrease and an increase of the O + .

Figure 9 .
Figure 9. foF2 and hmF2 variations at Millstone Hill between 08 and 16 LT on 4 November 2021.The blue arrows indicate the hmF2 and foF2 increase due to the TAD.
and hmF2 variations followed by a decrease of electron concentration at the F2-layer heights on 4 November 2021 and the [O] decrease with respect to the pre-storm quiet time level.A comparison with the MSISE00 model predictions shows that the model does not describe properly the [O]/[N 2 ] relative variations during the period under consideration.MSISE00 predicts, at MH, a [O]/[N 2 ] ratio decrease on 4 November, 1.8 times lower with respect to the reference day, while the [O]/[N 2 ] calculated from retrieved parameters is six times lower with respect to the reference day.6.

Table 1 .
Geomagnetic Ap index, vertical plasma drift W, exospheric temperature T ex , and neutral composition ([O], [N 2 ] and [O]/[N 2 ] ratio) retrieved at Moscow, Juliusruh, Pruhonice, Rome, and Millstone Hill at 12 LT.Retrieved exospheric temperature and neutral composition are given in comparison with MSISE00 model values at 300 km.Data for the reference quiet geomagnetic day are given in bold.