The Ionospheric Equivalent Slab Thickness: A Review Supported by a Global Climatological Study Over Two Solar Cycles

The ionospheric equivalent slab thickness (τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document}) is a parameter characterizing both the distribution of the plasma in the ionosphere and the shape of the corresponding vertical electron density profile. It is calculated as the ratio of the vertical total electron content (vTEC) to the ionospheric F2-layer electron density maximum (NmF2). Since its definition dated back in the 60s, a lot of information on the behavior of τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document} for different helio-geophysical conditions has been cumulated and the connection with several plasma properties has been also demonstrated. The beginning of the Global Positioning System (GPS) era in the 90s had a strong effect on the studies about τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document} because GPS signals allow to obtain the vTEC up to about 20000 km of altitude. Recently, τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document} has also found application in many data-assimilation methodologies, especially for the improvement of empirical ionospheric models based on near real-time data. All of these topics are reviewed and discussed in this paper based on the literature published in the last sixty years. Moreover, to highlight and summarize the main global climatological features of τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document}, in this work we selected thirty-two ionospheric stations globally distributed and co-located with ground-based Global Navigation Satellite System (GNSS) receivers, for the last two solar cycles. This allowed to collect a dataset of NmF2 and vTEC that represents the largest and most complete ever analyzed for studies concerning τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\tau $\end{document}, which gave the chance to deeply investigate its spatial, diurnal, seasonal, and solar activity variations. The corresponding results are presented and discussed in the light of the existing literature.


Introduction
The ionospheric equivalent slab thickness (τ ) is an important parameter characterizing the shape of the vertical electron density profile and, as a consequence, the plasma distribution in the dynamical system composed by the ionosphere and plasmasphere. Its definition, in terms of ratio of the vertical total electron content (vTEC) to the ionospheric F2-layer electron density maximum (NmF2), and early corresponding applications, date back to the Extended author information available on the last page of the article 60s when the first measurements of vTEC through Faraday rotation or Doppler frequency shift from geostationary satellites became available (Bauer 1960;Garriott 1960;Ross 1960;Roger 1964;Taylor 1965;Titheridge 1972). Since then, a bulk of evidences regarding the τ behavior for different geophysical conditions have been accumulated. In addition to this, the ionospheric equivalent slab thickness found applications in various modeling approaches. Its relation to other important ionospheric parameters, like the plasma scale height and the neutral and plasma temperature, was also investigated. All these matters are reviewed and discussed in this work.
This review is complemented with the most remarkable dataset of τ compiled so far. Specifically, thirty-two ionosonde stations globally distributed and co-located with groundbased Global Navigation Satellite System (GNSS) receivers, for the last two solar cycles (from 1997 to 2019), allowed to collect a valuable dataset of NmF2 and vTEC to get the main τ climatological features at equatorial, low, middle, high, and polar cap latitudes for both hemispheres. To perform the climatological analysis, the same methodology and procedures adopted by Pignalberi et al. (2021a) are applied. The analysis is performed for each of the thirty-two ionospheric stations considered in terms of: a) two-dimensional grids of binned medians of measured τ as a function of the local time (LT) and the month of the year, for three different levels of solar activity; b) boxplots representing the median and dispersion of τ as a function of the LT, for four months representative of different seasons, and for three different levels of solar activity. In order to highlight the most important characteristics of τ , the paper focuses on selected ionospheric stations representative of specific latitudinal sectors in terms of the Quasi-Dipole (QD, Laundal and Richmond 2017) magnetic latitude. The reader can find the results for all the thirty-two stations in the Supplementary Material.
The paper is organized in three parts: • Part I (Sect. 2 to 4) pertains the review section. In Sect. 2, the ionospheric equivalent slab thickness is defined and its physical and geometrical meanings are described along with the implications that the different methods applied over the years to measure the vTEC have on it. In the same section, the slab thickness sensitivity to both NmF2 and vTEC is analyzed to highlight the corresponding strengths and weaknesses. Section 3 presents the historical background and a review of the current knowledge on τ . The most relevant past studies on τ are here described to highlight the corresponding main spatial, diurnal, seasonal, solar and magnetic activities variations. The topics of this section, except the magnetic activity dependence, will be further discussed and deepened in Sect. 6 considering the results of the original climatological analysis reported in Sect. 5. In Sect. 4 a list of theoretical and modeling applications of τ developed by different authors are reported and discussed; • Part II includes the climatological analysis based on the new and original τ dataset (Sect. 5), and the corresponding discussion of the main patterns (Sect. 6). Section 5 first describes both the dataset of thirty-two co-located ionosonde stations and ground-based GNSS receivers used to obtain the τ time series, and the corresponding data filtering, binning, and statistical analysis. Then, the results of the statistical analysis are presented for five different latitudinal ranges by highlighting the main spatial, diurnal, seasonal, and solar activity variations of τ . Section 6 discusses the results of Sect. 5 in the light of the current knowledge on τ as outlined in Sect. 3. The most relevant features are here outlined and critically discussed by complementing and expanding the current knowledge on τ ; • Part III aims to bring to light what are the main gaps in the current knowledge on τ both from a theoretical and modeling point of view, and which future developments (Sect. 7) might be of help to this regard. Conclusive remarks are provided in Sect. 8. Fig. 1 Schematic representation of an ionospheric vertical electron density profile with highlighted the F2-layer electron density maximum NmF2 (red circle), the vertical total electron content vTEC (green double arrow segment), and the equivalent slab thickness τ (blue shadowed rectangle)

Definition and Calculation
The ionospheric equivalent slab thickness is defined as the ratio of vTEC to NmF2: where NmF2 is expressed in electrons per m 3 (el/m 3 ), and vTEC is measured in electrons per m 2 (el/m 2 ). As a consequence, τ is expressed in meters (m) and represents the thickness of an ideal ionospheric slab of constant electron density equal to NmF2 with a vTEC value equivalent to that of the entire vertical electron density profile. Figure 1 is a schematic representation of an ionospheric vertical electron density profile with highlighted NmF2 = 10 12 el/m 3 (red circle) at F2-layer peak, vTEC = 30 TECU (1 TECU = 10 16 el/m 2 ) (green double arrow segment), and τ = 3 · 10 5 m = 300 km (blue shadowed rectangle) as an ideal ionospheric slab with thickness 300 km around the F2-layer peak. NmF2 is obtained from the ordinary critical frequency of the F2 ionospheric layer (foF2) through the following relation NmF2 = (1.24 · 10 10 ) foF2 2 , with foF2 in MHz and NmF2 in el/m 3 (e.g., Davies 1990). foF2 is scaled from ionograms recorded by ionosondes either manually by an expert operator or automatically by robust autoscaling software :00 UT and autoscaled by the Autoscala algorithm developed at the Istituto Nazionale di Geofisica e Vulcanologia. In this case, the autoscaling software outputs a correct value for foF2 equal to 8.8 MHz, identified in correspondence of the asymptotical trend (highlighted by the dashed black vertical line) of the F2-layer ordinary ray. The figure shows also the vertical electron density profile obtained by applying the inversion procedure proposed by Scotto (2009) on the identified trace represented in blue. The ionogram has been downloaded from the Electronic Space Weather upper atmosphere database system (eSWua) owned by the Istituto Nazionale di Geofisica e Vulcanologia and managed by the Upper Atmosphere Physics and Radiopropagation group (www. eswua.ingv.it -HF data, version 1.0, https://doi.org/10.13127/eswua/hf) (Galkin and Reinisch 2008;Pezzopane and Scotto 2007). Of course, values validated by an expert operator are more reliable than those output by the autoscaling software. Nevertheless, the manual validation of the ionograms is really demanding and time consuming and, most of all, cannot be performed in real time, which is limiting for Space Weather purposes. This is why since the 80s much work has been done to develop algorithms able to automatically scale the ionograms giving as output the standard ionospheric characteristics and an estimation of the vertical electron density profile (e.g. Reinisch and Huang 1983;Fox and Blundell 1989). Figure 2 shows an example of ionogram recorded at the ionospheric station of Rome (41.8 • N, 12.5 • E; QD Lat. 35.9 • N) on 07 April 2022 at 19:00 UT and autoscaled by the Autoscala algorithm developed at the Istituto Nazionale di Geofisica e Vulcanologia Scotto et al. 2012). Figure 2 shows a well-defined ionogram, the F2-layer trace is very clear up to the critical frequency and therefore foF2 (and consequently NmF2) can be easily autoscaled by Autoscala from the vertical asymptote of the trace. Unfortunately, the ionograms are not always like that shown in Fig. 2, and often the F2-layer trace near the critical frequency is not clearly recorded due to multiple reasons: interference, absorption, scattering, multiple reflections, presence of sporadic E layers. In these cases, the autoscaling analysis that follows might be wrong and the output for the foF2 value unreliable (Reinisch et al. 2005;Galkin and Reinisch 2008;Scotto and Pezzopane 2008). Bullett et al. (2010) tested thoroughly the accuracy of the Autoscala software on the ionograms recorded by the Vertical Incidence Pulsed Ionospheric Radar ionosonde and found that, when the ionogram trace is clear, foF2 is automatically scaled with an accuracy of 0.1 MHz for 75% of the dataset and with an accuracy of 0.5 MHz for 99% of the dataset. Instead, when the ionogram trace is noisy or characterized by spread-F phenomena they found that the percentage of the autoscaling accuracy of 0.1 MHz drops to 50% and that of the autoscaling accuracy of 0.5 MHz drops to 94%. Therefore, when carrying out studies based on autoscaled data, one must be sure that these are largely reliable, as on the other hand it will be highlighted in Sect. 5.1.
The methodologies and techniques applied to calibrate vTEC values have evolved over the years and different approaches exist with different performances (Pignalberi et al. 2021b). Before the advent of the Global Positioning System (GPS) era in the 90s, vTEC was retrieved through Faraday rotation or Doppler frequency shift of geostationary satellites radio signals Titheridge 1972). However, Faraday rotation measurements are affected by the assumed hypothesis of a constant magnetic factor along the ray path, and this limits the retrieval of vTEC up to about 2000 km of height (Titheridge 1972;Davies 1990). The current availability of GPS and GNSS satellites have had a deep impact on the vTEC calculation, both for the ever increasing number of GNSS satellites constellations and for the possibility of using different vTEC retrieval techniques based on single-frequency and multi-frequency signals. vTEC retrieved from GNSS signals has the advantage of representing the contribution of the ionosphere-plasmasphere system up to about 20000 km of height. For the purpose of calculating τ , it has to be considered that the seminal works dated back to the 60s-80s were based on vTEC values retrieved through Faraday rotation measurements, while the most recent works are based on vTEC from GNSS signals. As a consequence, when comparing results obtained with different vTEC retrieval techniques, differences are expected because of the plasmaspheric contribution in the range 2000-20000 km of height (Pignalberi et al. 2021a). vTEC values retrieved from dual-frequency GNSS signals, as explained below, are used in the climatological analysis described in Sect. 5.
In recent times, TEC is usually measured by combining GNSS observables on, at least, two different frequencies (e.g., GPS L1 and L2). In fact, the ionospheric effects on GNSS signal propagation include a phase advance and a range delay, being proportional to the slant TEC (sTEC, i.e., the electron density integrated along the ray path connecting a given satellite-receiver couple) according to: in which sTEC is expressed in el/m 2 , f in Hz, and I is the delay in length unit induced by the ionosphere in the pseudo-range. Given Eq.
(2), sTEC can be obtained by calculating the so-called geometry-free linear combination, that is defined as: in which subscripts 1 and 2 refer to two different GNSS frequencies (e.g., L 1 and L 2 GPS frequencies), respectively, while λ refers to the corresponding wavelength, N is the ambiguity on carrier phase measurements, "arc" refers to continuous carrier phase observations (i.e., for which the product λN can be considered constant), t is the frequency-dependent delay due to the receiver (R) and satellite (S), c is the speed of light in vacuum, and ε is the noise on the carrier phase observation combination in length units. From (2) and (3), the following can be derived: in which B R = c(t 1,R − t 2,R ) f 2 40.3 and B S = c(t 1,S − t 2,S ) f 2 40.3 are the "inter-frequency biases" (IFBs) for carrier phase observables induced by the receiver and the transmitter, respectively, C arc = f 2 40.3 (λ 1 N 1 + λ 2 N 2 ) is the ambiguity for L 4 , and ε L = f 2 40.3 ε is the carrier phase noise in TECU.
A similar linear combination for the code measurements can be obtained as: where b R and b S are the so called "differential code biases" (DCBs), due to the receiver and satellite, and ε P is the noise on code measurements. In Eq. (5) there are no ambiguity terms since, differently from phase observables, code measurements are "absolute", i.e., they do not suffer from cycle ambiguities. Considering a single arc of observations, by subtracting the two new observables, one for code measurements and the other one for phase measurements, the average difference on a single arc (brackets in Eq. (6)) between carrier phase and code observables can be written: in which ε L can be neglected (Braasch 1996). By subtracting (6) from (4), the following can be derived: being the carrier-phase ionospheric observable "levelled" to the code-delay ionospheric observable. The term ε P arc is the effect of the noise and the multipath along the arc on carrierphase observations, as described in Ciraolo et al. (2007).
There are several approaches that can be used to estimate the b R , b S , ε P arc "bias" terms based on different assumptions (e.g., ε P arc ≈ 0) estimating the biases for a single station or for multiple stations at once (Seemala and Delay 2010;Conte et al. 2011). One possible approach to estimate sTEC without neglecting ε P arc consists in estimating a single bias for each arc, instead of separate biases for the receiver and satellites. In this case, Eq. (7) can be written as:L in which β arc is the sum of all the errors with non-zero mean. This is the equation defining the TEC calibration, following the method described by Ciraolo et al. (2007). In the specific, sTEC values are mapped as a two-dimensional surface by means of the classical thin shell method as: where vTEC(φ 1 , φ 2 ) is the unknown describing a surface in the reference frame defined by a coordinates couple (φ 1 , φ 2 ) over the thin shell (bi-dimensional), and χ is the angle formed by the ray-path and the perpendicular to the shell at the ionospheric pierce point (IPP) (Pi et al. 1997;Mannucci et al. 1998). The geometry of the thin shell mapping approach to retrieve vTEC from sTEC is outlined in Fig. 3. Once defined the functional form of vTEC(φ 1 , φ 2 ), Eq. (9) can be solved via the standard least squares method. Thin shell geometry representation. R e is the radius of the Earth, χ is the angle formed by ray path and the vertical crossing the ionosphere in the ionospheric pierce point (IPP), H IPP is the height of the thin shell ionosphere. Red line represents the sTEC while the blue line represents the vTEC 2.2 Meaning and Implications for the Ionosphere τ is a significant ionospheric parameter because depending on both vTEC and NmF2 observations can provide information on both the topside and bottomside ionosphere, and on how the plasma is distributed between these regions and the overlying plasmasphere. To highlight that, top panels of Fig. 4 represent τ curves from 100 to 500 km in steps of 100 km for NmF2 and vTEC values usually observed in the ionosphere-plasmasphere system. Top right panel of Fig. 4 is a zoom of the left one for low NmF2 and vTEC values that are usually observed at night or at the solar terminator passage. From its mathematical definition and the observation of Fig. 4 top panels it is clear how τ is very sensitive to NmF2 or vTEC changes for very low NmF2 and vTEC values, while it is more stable for high NmF2 and vTEC values. At very low NmF2 and vTEC values, a little change in vTEC of the order of few TECUs (i.e., of few 10 16 el/m 2 ), or a corresponding one in NmF2 of the order of 10 11 el/m 3 , allows stepping from a τ curve to another. For example, at NmF2 = 10 11 el/m 3 (foF2 ≈ 2.8 MHz) stepping from 2 to 3 TECU, τ rises from 200 to 300 km (corresponding curves are identified respectively in orange and green in the top right panel of Fig. 4). On the contrary, at NmF2 = 20 · 10 11 el/m 3 (foF2 ≈ 12.7 MHz), to obtain the same τ change, vTEC has to step from 40 to 60 TECU. This fact has a lot of implications when comparing τ values retrieved by the oldest methodologies based on geostationary satellites data and the newest ones based on GNSS satellites data. In fact, the plasmaspheric contribution between 2000 and 20000 km is of the order of magnitude of units of TECU and, more importantly, during night the plasmaspheric contribution to vTEC is relevant while the ionospheric one is at its minimum. From these considerations, one might expect that most of the differences between the results shown by different authors in the past could be ascribed (especially when vTEC values are low) to the technique applied to retrieve vTEC data for the τ calculation, and to the errors associated to the vTEC calibration itself. On the contrary, the NmF2 retrieval from ionograms has not changed over the years and the error associated to this kind of remote sensing technique is negligible compared to the one affecting the vTEC retrieval and calibration. Moreover, while the NmF2 is a punctual measure associated to the F2-layer peak, the vTEC is an integral one from the base of the ionosphere (at about 50-60 km of height) to the GNSS satellites height (at about 20000 km of height). As a consequence, the vTEC includes the contribution of the bottomside ionosphere (up to the F2-layer peak height), of the topside ionosphere (from the F2-layer peak height to the plasmasphere) and of the plasmasphere itself. However, the separation in these three regions is a matter of nomenclature because there are continuous exchanges of plasma between them and such fluxes of charged particles are not stationary being dependent on the ongoing helio-geophysical conditions; thus, diurnal, seasonal, geographic, solar and magnetic activity variabilities are found. The only vTEC value cannot describe all these complications, providing only an overview of the ionosphere-plasmasphere system; as a consequence, also the information given by the slab thickness is partial and it is difficult to effectively separate the ionospheric and the plasmaspheric contributions. However, by making simple assumptions on the vertical electron density profile shape, it is possible to infer very useful information from the slab thickness value. According to this, in order to understand the geometrical meaning of the slab thickness, we used a prototype version of the NeQuick model (Nava et al. 2008) in which the τ value (and hmF2, where hmF2 is the height of the F2-layer peak) can be constrained, to show the effect that changing the τ value has on the shape of the vertical electron density profile, for fixed vTEC or NmF2 values. In the bottom panels of Fig. 4 the prototype NeQuick model was run for τ values from 100 to 500 km in steps of 100 km in two different ways: by fixing the NmF2 value at 10 12 el/m 3 (≈ 8.98 MHz) and leaving vTEC free; by fixing the vTEC value at 20 TECU and leaving NmF2 free. In both cases, the NeQuick model was run according to the following parameters: latitude = 45 • N, longitude = 15 • E, month = April, hour = 14 Universal Time (UT), hmF2 = 300 km. From the bottom panels of Fig. 4 it is evident the effect of the slab thickness on the shape of the profile: in a nutshell, small τ values are representative of "thinner" electron density profiles, while high τ values are associated to "thicker" profiles, at the F2-layer peak. Therefore, τ can describe the distribution of the plasma between the lower (bottomside) and the upper (topside) ionosphere by modifying the profile shape.
The τ sensitivity to variations of vTEC and NmF2 can also be analytically studied through a two-variables Taylor series expansion around specific vTEC and NmF2 values, as outlined below. From the mathematical theory of two-variables Taylor series expansion, at the first order, a two-variables function can be approximated by: with In our specific case, the two-variables function is τ , then: By applying Eqs. (10)-(11) to Eq. (12), we obtain: with Equation (13) describes the first-order variations of τ to variations of both vTEC and NmF2 around the point (NmF2 0 , vTEC 0 ). In case of only vTEC or NmF2 variations, the corresponding τ variations are described by Eq. (15): It is worth noting that Eqs. (13) and (15) are locally valid as a first-order approximation around the point (NmF2 0 , vTEC 0 ). Thus, they are useful to estimate τ variations due to possible measurement errors characterizing either vTEC or NmF2, or both of them.
(19) In this case, Eqs. (18)-(19) show the very strong impact that an error in the calibration on vTEC has on the τ estimation, much higher than that associated to NmF2 one. As before, the major impact is when vTEC and NmF2 are low (as during nighttime) with percentage relative error of about 8.3% and 33.3% for vTEC equal to 0.5 and 2 TECU, respectively.
These analytical reasonings, along with Fig. 4 results and corresponding discussion, point out the very high sensitivity of τ to vTEC variations, especially for low vTEC and NmF2 values, which are usual during nighttime. Then, even a relatively small error in the vTEC calibration can produce a quite high misrepresentation of τ . This has also many implications when comparing results obtained from vTEC values obtained from Faraday rotation measurements versus those obtained with modern GNSS ground-based receivers. In fact, the contribution to vTEC of the plasmaspheric region between about 2000 and 20000 km of altitude ranges from tenths to units of TECU, thus producing a relevant underestimation in the τ prediction when Faraday rotation measurements are used.

Historical Background and Current Knowledge on the Slab Thickness
Ionospheric and plasmaspheric processes affect concurrently both the bottomside and topside ionospheric variability making the τ behavior rather complex. The main diurnal, seasonal, latitudinal, solar and magnetic activities climatological features of τ have been investigated by many authors. From the wide existing literature, it emerges a very complicated picture that reflects the inherent complexity of this ionospheric parameter and, as a consequence, the need of further studies.
Since τ depends significantly on the latitude, the most relevant studies here cited are first grouped taking into account the location of the stations used by the authors, and the main diurnal and seasonal features are described. Thereafter, the main results about the τ variations connected to the solar and magnetic activity are shown.
This Section aims at briefly introducing the main climatological features of τ as outlined so far by many authors. The topics here presented will be further discussed, complemented, and deepened in Sect. 6 considering the results of the original climatological analysis reported in Sect. 5. The only exception is the magnetic activity dependence which is not faced by Sect. 5: as a consequence, a more detailed picture of this dependence is given in Sect. 3.6.

Slab Thickness at Low Latitudes
With regard to studies over the equatorial sector, the diurnal variation of τ over Kodaikanal (10.1 • N, 77.3 • E; QD Lat. 2.7 • N) during the winter of 1966-1967 shows the presence of a pre-sunrise peak with high values visible throughout the day (Das Gupta et al. 1975). Iyer and Rastogi (1978), analysing data again over Kodaikanal during half a solar cycle (1964)(1965)(1966)(1967)(1968)(1969), found a semi-annual variation of τ characterized by equinoctial maxima both during daytime and nighttime, with a pre-sunrise peak in winter and summer months for high solar activity. Values greater in summer and winter, unusual summer pre-sunrise peaks along with minima during post-sunrise and sunset hours were found by Kenpankho et al. (2011), who investigated the τ behavior for medium solar activity over Chumphon At low latitudes, the τ behavior at Hawaii (19.0 • N, 154.0 • W; QD Lat. 19.3 • N) relative to one year of high (1981) and low (1985) solar activity has revealed the occurrence of pronounced pre-sunrise peaks in winter, with daytime yearly mean significantly higher than nighttime ones for high solar activity conditions, while the reverse is true for solar minimum activity . Likewise, also studies conducted over Nicosia (35.0 • N, 33.0 • E; QD Lat. 29.2 • N) during minimum solar activity (from January 2009 to August 2010) by Haralambous (2011) confirmed the presence of pre-sunrise peaks in winter with a post-sunset increase and mean nighttime values higher than the mean daytime ones for all seasons except for summer. Venkatesh et al. (2011) -23 • N], with τ decreasing as the latitude increases, was observed for medium solar activity during daytime in winter by Das Gupta et al. (1975) and for all seasons by Venkatesh et al. (2011), who analyzed the hourly behavior of τ at some equatorial and low-latitude stations.

Slab Thickness at Middle Latitudes
The complexity of the τ behavior at Northern and Southern middle latitudes has been highlighted by several works (Bhonsle et al. 1965;Titheridge 1973;Essex 1978;Davies and Liu 1991;Minakoshi and Nishimuta 1994;Spalla and Ciraolo 1994;Goodwin et al. 1995aGoodwin et al. , 1995bBreed et al. 1997;Leitinger et al. 2004;Jin et al. 2007;Kouris et al. 2008Kouris et al. , 2009Stankov and Warnant 2009;Mosert et al. 2013;Huang and Yuan 2015;Sharifi and Farzaneh 2015) which highlighted: a) pre-sunrise peaks in all seasons for both low and high solar activity; b) post-sunset peaks in winter for low solar activity and in summer for both medium and high solar activity; c) the occurrence of pronounced peaks during post-midnight hours; d) daytime values smaller than nighttime ones for all seasons for low solar activity, except for summer where the reverse occurs; e) daytime values significantly larger in summer than in winter. a) pre-sunrise peaks in all seasons for low, medium, and high solar activity, which stand out clearly in winter; b) post-sunset peaks in winter for low and medium solar activity; c) high values for high solar activity for some equinoctial and summer months in the LT sector 19-22; d) daytime values markedly smaller than nighttime ones, regardless of solar activity (except in summer); e) a higher dispersion of measurements during nighttime and solar terminator hours.

Moreover
The studies carried out by , analyzing the τ behavior for fourteen world-wide stations located at different latitudes from the geomagnetic equator to the North pole for high solar activity, showed that τ can behave in a really complex manner with latitude; they found a clear decreasing trend of τ as the latitude increases only in winter during daytime.

Slab Thickness at High Latitudes
Compared to studies performed at low-and middle-latitude stations, the investigation of the τ behavior at high and polar cap latitudes are very rare.  found nighttime values higher than daytime ones at Goose Bay (53.3 • N, 60.3 • W; QD Lat. 60.2 • N) for both solar maximum and minimum conditions; while the occurrence of pronounced post-sunset peaks was observed in winter and equinoctial months for high solar activity. Analysing data over the Antarctica station of Casey (66.3 • S, 110.6 • E; QD Lat. 80.6 • S) during a period of high solar activity, Yadav and Bhawre (2020) found nighttime values significantly higher than daytime ones, mostly during equinoctial months, and presunrise peaks were clearly detected during summer and equinoctial months.

Global-Scale Variations of the Slab Thickness as Derived from vTEC Maps
The τ behavior has been analyzed also by means of vTEC maps derived from groundbased GPS receivers and Constellation Observing System for Meteorology, Ionosphere and Climate (COSMIC) Radio Occultation (RO) measurements (Anthes et al. 2008). To this regard, Guo et al. (2011)  a) during nighttime, values are higher than during daytime for non-auroral latitudes in all seasons; an opposite behavior is observed at auroral latitudes in summer; b) at high latitudes, higher values are observed in winter with a constant trend of around 350 km during daytime, while a decreasing trend was observed after the pre-sunrise peak for equinoxes; c) minimum values around 200 km were noticed in summer during daytime; d) for all seasons, the diurnal variation in the Southern hemisphere is smaller than that in the Northern one; e) the diurnal variations observed in the winter hemisphere are comparatively greater than those in the summer one; f) pre-sunrise peaks at middle latitudes are higher than those observed at equatorial and high latitudes; g) peaks near the magnetic equator around mid-noon are observed in all seasons; h) post-sunset peaks are observed especially in winter and equinoxes at middle latitudes.
In a next investigation, from the global diurnal maps of τ reconstructed from Global Ionospheric Map (GIM) TEC and COSMIC RO observations collected during the period from 2006 to 2014, Huang et al. (2016) found: a) a depletion in the Equatorial Ionization Anomaly (EIA) region, with a pre-sunrise enhancement observed at all latitudes and the occurrence of a post-sunset enhancement only in some areas at equinoxes for low and high solar activity; b) during June solstice, daytime values larger in the Northern hemisphere (summer) than in the Southern hemisphere (winter), while pre-sunrise and post-sunset peaks larger in the winter hemisphere than in the summer one; c) particularly large values at high and polar cap latitudes for all seasons and for both minimum and maximum solar activity; d) an enhancement at nighttime in the longitude sector 0 • -120 • E and latitude sector 30 • -75 • S during the Southern hemisphere summer.

Solar Activity Dependence
The τ variability in relation to the solar activity cycle has been also widely studied in order to find possible correlations. With this regard, there exists a large variety of results depending on the latitudinal sector and the location of the considered station. For example, in all seasons, no distinct dependence on the solar activity was observed over Kodaikanal near the magnetic equator (Iyer and Rastogi 1978). On the contrary, over Chumphon, again close to the magnetic equator, a diurnal variation clearly increasing with the solar activity was found (Kenpankho et al. 2011). With regard to low latitudes, a positive correlation with the solar index F10.7 data was found in winter during daytime by  at Hawaii; analogously, at Lunping (23.0 • N, 121.9 • E; QD Lat. 16.3 • N), positive correlations were observed in all seasons with solar indices EUV (170-190 Å), F10.7, and SSN (solar sunspot number), with the best one associated to the EUV flux (Dabas et al. 1993). Also Chuo (2007) over Chung-Li (24.9 • N, 121.0 • E; QD Lat. 18.3 • N) found a positive correlation with the solar index F10.7. Again, Iyer and Rastogi (1978) report that for low and middle latitudes τ increases with solar activity. At middle latitudes, a positive correlation over Sagamore Hill (42 • N, 288 • E; QD Lat. 51.3 • N) was observed for all seasons at noon (Davies and Liu 1991); while, more recently, an increase of τ with solar activity was observed over Dourbes (50.1 • N, 4.6 • E; QD Lat. 45.8 • N), with variations of up to 24% and 14% in the daytime and nighttime hours respectively (Stankov and Warnant 2009). At high latitude, an increase of τ with F10.7 was observed at Goose Bay but only for high solar activity . In a study conducted at Casey, located over the Southern polar cap, positive correlations were found between τ and solar radiation fluxes X ray (1-8 Å) and EUV (260-340 Å), while a poor correlation was detected with F10.7 (Yadav and Bhawre 2020). Finally, also a significant negative correlation during winter nighttime hours was observed in the Southern hemisphere at Salisbury (34.8 • S, 138.6 • E; QD Lat. 45.6 • S) during the years 1991-1994 (Breed et al. 1997).

Magnetic Activity Dependence
The effects of magnetic activity disturbances on τ is a topic that has received much attention since the first studies dating back to the 60s due to the very important implications that such impulsive and short-time events have on the F2 layer and on the plasma distribution in the ionosphere. Although not strictly connected to the climatological behavior of τ , the magnetic activity dependence is a relevant topic and the most remarkable findings of past studies are here reported. The difficulties to relate the magnetic activity to corresponding τ variations is testified by the results shown by different studies: a little daytime effect in summer (Ross 1960); a not clear or a very weak correlation (Taylor 1965;Dabas et al. 1979); no quiet-time and storm-time significant correlation (Fox et al. 1991); conflicting results, with τ increasing and decreasing according to the season (Bhonsle et al. 1965).
On the contrary, from several studies has emerged that considerable and clear increases of τ just occur in the course of a geomagnetic storm in all seasons and hours of the day and in both hemispheres. For example, in several cases analyzed by Garriott (1960), a remarkable increase of daytime τ values was observed during a geomagnetic storm occurred in winter; while during a geomagnetic storm occurred in summer in the northern hemisphere, τ was observed to increase by 30% (Bauer 1960), and it was seen even to double during the development of a geomagnetic storm occurred in the southern hemisphere in winter (Munro 1962).
By plotting the average τ value calculated for the period July 1961-October 1964 in three separate daytime intervals against the geomagnetic index Kp, Yeh and Flaherty (1966) from τ observations made at the middle latitude station of Urbana (40.1 • N, 88.1 • W; QD Lat. 50.4 • N) found a tendency of τ to increase as the magnetic activity increases. Hibberd and Ross (1967) analyzed storm-time variations of τ/τ (being τ the departure of τ measurements from climatological values) calculated during 18 geomagnetic storms occurred in the period July 1961-June 1962, from the beginning of the sudden storm commencement or when the geomagnetic index ap > 39 nT. They found a positive correlation as described by the positive τ/τ variations clearly observed up to three days later the storm commencement, with greater values occurring within the first ten hours. Titheridge and Andrews (1967) studied the behavior of τ at a southern hemisphere middlelatitude station near Auckland (36.8 • S, 174.7 • E; QD Lat. 42.6 • S) during the geomagnetic storm of 16 June 1965. They observed τ values above the monthly median level for almost one day (≈ 06:00 16th -06:00 17th ) and some hours after the storm commencement (23:00 15th LT), correspondingly to a clear increase of the Kp index. τ values below the monthly median value and a slow gradual decrease of Kp were observed in the next 42 hours (06:00 17th -24:00 18th LT). Mendillo et al. (1972) investigated the average storm time and local time variations of τ from the monthly median during 28 geomagnetic storms occurred in the period December 1967-November 1969 at the middle-latitude station of Hamilton (42.6 • N, 70.8 • W; QD Lat. 51.7 • N). During such storms, τ enhanced with peaks clearly observed during the first two days of the storm in the early morning (10:00-12:00 LT) and in the afternoon hours (12:00-15:00 LT). Kerseley and Hajeb-Hosseinieh (1976) found an increasing linear trend between τ and geomagnetic activity during moderate solar activity levels at the middle-latitude station of Aberystwyth (52.4 • N, 4.0 • W; QD Lat. 49.0 • N). Balan and Iyer (1978) selected storm periods characterized by Ap ≥ 30 nT and studied the departures of τ from the mean of the 7 days before the sudden storm commencement ( τ ) to find whether geomagnetic activity can influence the τ behavior observed at Hamilton station. They observed positive phase ( τ > 0) and negative phase ( τ < 0) occurring more frequently during daytime and nighttime, respectively. The effect of storm magnitude on τ was investigated plotting τ vs. Ap. The results showed that the increase of the geomagnetic activity leads to a clear increase of τ , while for negative values of τ no clear dependence on geomagnetic activity was observed. Chauhan and Gurm (1981) studied the effects of geomagnetic storms on τ during the solar minimum period of 1975-1976 for two Indian low-latitude stations. At Ahmedabad (23.0 • N, 72.6 • E; QD Lat. 17.0 • N) station, located near the crest of the EIA, τ exhibited values above its mean level during two moderate geomagnetic storms. Differently, at Delhi station, far from EIA, the results turned out to be in general inconclusive without any clear dependence on magnetic activity. Huang (1983) considered weak (Ap < 30 nT) and strong (Ap ≥ 30 nT) geomagnetic storms occurring during the period March 1977-February 1982 to investigate the geomagnetic storm effects on the variation of τ at the low-latitude station of Lunping. The results highlighted that, when strong storms occur, two hours after the storm commencement τ starts to increase reaching its maximum three hours later, and a positive phase around midnight or pre-sunrise hours and a negative phase in daytime are observed until the end of the second storm day. Instead, during weak storms, τ shows features that are similar to those observed for the strong ones, even though unusual negative phases during pre-storm hours are observed.  investigated the effect of magnetic disturbances on τ by comparing the mean diurnal variations of τ for magnetically quiet days (Ap < 10 nT) with the ones obtained for disturbed days. They found a significant enhancement of τ , both in daytime and nighttime hours, at the middle-latitude station of Boulder (40.0 • N, 105.3 • W; QD Lat. 48.3 • N) and at the high-latitude station of Goose Bay during the years of solar minimum (1985) and solar maximum (1981), respectively. Gulyaeva and Stanislawska (2005) studied the storm-time τ behavior through spatial maps of τ during five magnetic storms occurred in the years of minimum (1995)(1996) and maximum (2002) of solar activity using data from several stations spread throughout Europe. They performed a super-posed epoch analysis based on an index defined to describe how much the measured map differ from the climatological one. Their study highlighted an enhancement of τ for more than 2 days from the onset of the storm at solar maximum, while at solar minimum the increase of τ appears delayed by 12 to 24 hours with respect to the onset of the storm and with a shorter duration. The τ storm-time behavior was also investigated at the European middle-latitude station of Dourbes by Stankov and Warnant (2009) during two severe geomagnetic storms occurred on 15 July 2000 and on 8 November 2004, both characterized by a Dst geomagnetic index close to −400 nT. For the first storm, during both the main phase (≈ 19:30 15th -23:00 15th UT) and the early hours of the recovery phase (≈ 23:00 15th UT-06:00 16th UT), a peak of τ with a value greater than 600 km was observed. In the second storm, during the main phase (≈ 21:00 7th UT-06:00 8th UT) and the early stages of the recovery phase (≈ 06:00 8th UT-09:00 8th UT), a very significant increase of τ was observed. A peak of τ at 650 km, corresponding to 250% above the monthly mean, was found three hours after the beginning of the recovery phase; after that, τ started to decrease but always with values higher than the quiet reference level for all the storm period. Chuo et al. (2010) found an apparent increase of τ with respect to its monthly median behavior at Wuhan (30.6 • N, 114.4 • E; QD Lat. 24.5 • N), a middle-latitude station close to the northern crest of EIA, during the main and recovery phase of a moderate geomagnetic storm (minimum Dst = −238 nT) occurred on 22 October 1999. Studies analyzing τ during storm time at equatorial latitudes are relatively few. In a recent survey, Zhang et al. (2021) defined a disturbance index (DI) based on τ measurements recorded at Guam (13.6 • N, 144.8 • E; QD Lat. 6.1 • N) station, located near the magnetic equator (5.54 • N dip latitude), during the years 2012-2017. DI quantitatively describes the effects of geomagnetic storms on τ variations. The comparison between the storm-time and quiet-time DI diurnal variations leads to the conclusion that geomagnetic storms enhance τ during most of the storm period except for some cases when they have negative effects around sunrise.
The above mentioned studies highlighted a very cryptic picture with different results depending on location, solar activity level, local time hour, and the level of magnetic disturbance. In fact, from past studies it is difficult to draw a clear and definite pattern of the magnetic activity influence on τ variations. In the investigation of the τ variations triggered by magnetic activity disturbances, an issue that has to be taken into account is that it is not always easy to establish the quiet level of τ to which refer. Moreover, it may not be wholly correct to correlate the time series of τ with the simultaneous ones of a whatever geomagnetic index or with its daily average (Hibberd and Ross 1967). For example, time lags of one or two hours between magnetic data and τ time series should also be investigated to see whether better results may be achieved.

Main Applications and Models of the Slab Thickness
Over the years, the slab thickness has found application in several fields linked to both the explanation of physical processes and modeling purposes. We review here the most relevant of these applications, although some of them are now outdated.
Several studies have shown that τ is related to a variety of physical processes: a) it allows to estimate the scale height of the ionized constituents through the relationship linking τ to the neutral vertical scale height H n = k B T n /m i g (Hargreaves 1992;Rishbeth and Garriott 1969), where k B is the Boltzmann's constant, T n is the neutral temperature, m i is the ion mass and g is the acceleration due to the gravity field. Specifically, in the simplifying working hypothesis of a single α-Chapman layer with constant scale height H n , it can be demonstrated that τ = 4.133H n (Wright 1960;Hibberd and Ross 1966;Titheridge 1973). However, like discussed in Pignalberi et al. (2021a), this formula implies the assumption of a constant value of H n for the whole topside ionosphere, which is a strong assumption, as pointed out by recent studies (e.g., dos Santos Prol et al. 2019; Pignalberi et al. 2020a). Moreover, and more importantly, the α-Chapman function represents the F2 layer of the ionosphere only under very strict assumptions unrepresentative of the real ionosphere whose vertical distribution of plasma is driven by the plasma scale height H p = k B T p /m i g that, in turn, depends on the plasma temperature T p and other physical quantities (e.g., Pignalberi et al. 2018Pignalberi et al. , 2020a; b) it is related to the plasma temperature T p = (T e + T i ) even if it does not seem to be a good indicator of electron (T e ) and ion (T i ) temperatures (Amayenc et al. 1971;Furman and Prasad 1973;Stankov and Jakowski 2006). Also in this case, albeit a slight dependence of τ on T p is not excluded, since the electron density has a certain degree of correlation with the electron temperature, it is too simplistic to reduce the τ trends to the T p ones; c) it allows to estimate the thermospheric neutral temperatures (T n ) through the relationship τ = 0.216T n (Titheridge 1973) or T n = (τ + 250)/0.5 (Jakowski et al. 2017). As above, T n regulates many chemical reactions that are fundamental for the ionization of the particles constituting the ionosphere and affects the vertical extension of the atmosphere as a whole; then, a slight dependence of τ on T n has to be expected. However, no strict correlation between τ and T n has been ever demonstrated.
For what concerns modeling applications, the slab thickness found application in several data-assimilation methodologies as the parameter linking vTEC to NmF2 (the vice versa is also possible but not very useful). In fact, its knowledge is of practical interest to predict NmF2 at locations where only vTEC values are known from observations, a condition of easy occurrence given that the number of ground-based GNSS receivers is by far greater than that of ionosondes. For example, Krankowski et al. (2007) used τ as a factor to convert vTEC maps from GPS observations of the EUREF permanent GNSS network (https:// www.epncb.oma.be/) to foF2 values. Specifically, they found τ = 310 km as the optimal conversion factor on the basis of the analysis of foF2 observations recorded at five European ionospheric stations from 1996 to 2001 and corresponding vTEC values from co-located GPS receivers. Gerzen et al. (2013)  It is worth highlighting that the sensitivity of τ with respect to not climatological processes makes it a good candidate for Space Weather monitoring and alerts (Galkin et al. 2022).
τ is applied also to improve the modeling of the electron density profile shape when ingested by empirical models like NeQuick (Nava et al. 2011). Gulyaeva et al. (2021) explored the possibility to use the center of mass (Hc) of the ionosphere as a proxy for the effective ionospheric shell height (to be used in the conversion of sTEC to vTEC). For this purpose, Hc has been determined from the electron density profiles within the boundaries of τ in the height range (hmF2 − τ bot , hmF2 + τ top ) where τ bot and τ top are the bottomside and topside equivalent slab thickness. Fox et al. (1991) developed a model describing the τ climatological behavior for middle latitudes based on NmF2 values recorded by the ionosonde installed at Wallops Island and vTEC data measured at Hamilton using Faraday rotation of signals from geosynchronous satellites, during the period 1965 to 1986. Pignalberi et al. (2021a) compared the output of the Fox et al. (1991) model to actual τ values recorded at Wallops Islands for the last two solar cycles, highlighting several differences most of which are due to the different methodologies applied to derive vTEC data for the τ calculation, as discussed in Sect. 2.
Recently, Jakowski and Hoque (2021) developed an empirical global climatological model of τ which is the only available so far. They modeled τ based on a dataset of NmF2 values derived from CHAMP (CHAllenging Minisatellite Payload), GRACE (Gravity Recovery And Climate Experiment), and COSMIC/FORMOSAT-3 radio occultation measurements, and vTEC values from vTEC maps generated by CODE (Center for Orbit Determination in Europe) based on IGS data archive. Through a dataset spanning from 2001 to 2018, they described the main spatial, diurnal, seasonal, and solar activity variations in τ in a very compact way based on 12 coefficients via non-linear least square methods. Jakowski and Hoque (2021) model reproduces the main τ climatological variations and anomalies but the comparison with τ values derived by five co-located ionosondes-GNSS receivers at different latitude revealed the very-high variability and dispersion of τ (particularly for the nighttime and solar terminator hours) which the model can hardly represent. As a consequence, there is still a lot of room for improvement about τ global empirical modeling to better describe its behavior also under non-equilibrium conditions.
It is of course possible to calculate τ as a by-product of the NmF2 and vTEC values modeled by IRI and NeQuick (Nava et al. 2008) global models. However, such models are not designed for the τ calculation and the datasets on which they are based are heavily unbalanced towards the lowest part of the ionosphere. So, the description of the topside part of the ionosphere made by such models needs improvements (Pignalberi et al. 2016(Pignalberi et al. , 2018.

vTEC and NmF2 Datasets
Section 3 summarized the most relevant investigations of the τ behavior performed in the last 60 years, most of which are related to very limited datasets. This happens because it is not easy to get datasets of NmF2 and vTEC measurements that are simultaneously available for a long period of time over a given location (i.e., a specific ionospheric station). As a consequence, the selection criteria of an ionospheric station suitable for studying the behavior of τ heavily rely on the availability of an operational ionosonde and a co-located  Table 1 for the identification number (ID) specification GNSS receiver, both of them continuously working for a long period of time to allow for a climatologic study based on a robust statistics.
Specifically, in this study, the pairs ionosonde station-GNSS receiver routinely working for at least ten years, so as to cover at least a solar cycle, were selected. When available, we selected exactly co-located ionosonde station-GNSS receiver pairs; however, in some cases this was not possible due to the current spatial distributions of such observational facilities. Anyway, we selected pairs of ionosonde station-GNSS receiver whose distance is well below the ionospheric horizontal correlation distance (Forsythe et al. 2020); the maximum distance between pairs is the one of Jicamarca station (759 km, see Table 1). According to this criterion, we were able to select thirty-two ionospheric stations (i.e., pairs of co-located ionosondes-GNSS receivers) distributed all over the world (see Table 1 and Fig. 5). NmF2 and vTEC data, recorded respectively by co-located ionosondes and GNSS receivers, were used to calculate τ values during the last two solar cycles, from the beginning of 1997 to the end of 2019. Figure 6 represents the time spanned and yearly percentages of available τ values for each ionospheric station. As far as we know, the dataset considered in this study is the longest and most complete used so far for a climatological analysis of τ .
It is worth noting that most of the selected ionospheric stations are located at Northern middle latitudes (Fig. 5), while for both the equatorial and very high latitudes, there are only a few ionospheric stations satisfying our requirements. Overall, it is almost the same for the Southern hemisphere. This is due to the unevenly distribution of both ionosondes and GNSS receivers between the two hemispheres and to the difficulties in installing and maintaining such instrumentation at equatorial and high latitudes. With the growing number of GNSS receivers installed throughout the world and newly installed ionosondes, over the next few years it will be possible to collect enough data to extend and enlarge this investigation at equatorial and high latitudes.
As mentioned in the Introduction, this investigation represents an extension of a previous work recently published by Pignalberi et al. (2021a) concerning the τ climatology over three Northern middle-latitudes ionospheric stations. In this study, we applied the same methodology and validation analysis proposed by Pignalberi et al. (2021a). As a consequence, only Table 1 List of ionosonde stations and corresponding GNSS ground-based receivers, with the corresponding network of belonging, used in this study. The first column marks the identification number (ID) of each ionospheric station (see Fig. 5). The geographic coordinates of each ionosonde station and GNSS ground-based receiver are reported along with the distance between them in km (last column). The Quasi-Dipole magnetic latitude (Laundal and Richmond 2017) is also reported for each ionosonde station ID the essential information about NmF2 and vTEC datasets and corresponding data selection and binning are provided here, the full details and additional information can be found in Sect. 2 of Pignalberi et al. (2021a). NmF2 data were retrieved from ionograms recorded by DPS Digisondes (Bibl and Reinisch 1978) at the thirty-two ionospheric stations, and autoscaled by the Automatic Real-Time Ionogram Scaler with True height analysis (ARTIST) software (Galkin and Reinisch 2008). Only ionograms characterized by a suitable Confidence Score (C-Score) index were selected (C-Score ≥ 75) so as to guarantee the reliability of the ionogram itself and consequently of the autoscaled NmF2 value. According to the sounding repetition rate of the ionosondes, the NmF2 time series have a time sampling of fifteen minutes (at minutes 0, 15, 30, and 45 of each UT hour).
vTEC data were obtained from daily Receiver INdependent EXchange (RINEX) files, recorded by GNSS ground-based receivers selected as near as possible to the ionosonde stations, through the application of the Ciraolo's calibration method (Ciraolo et al. 2007). The ionospheric pierce point altitude was set to 350 km, and the minimum elevation angle to 20 • . Differently from the NmF2 time series, vTEC ones have a thirty-seconds time sampling. To make the two datasets consistent, a procedure making use of a Gaussian weight was adopted so as to provide the vTEC values actually used for the τ calculation as weighted mean of vTEC values sampled at a thirty-seconds rate (vTEC) in windows fifteen-minutes wide centered at minutes 0, 15, 30, and 45 of each UT hour. The details of this procedure are given in Pignalberi et al. ( , 2021a.

Slab Thickness Calculation and Binning
NmF2 and vTEC time series concurrently available were used to get the τ time series, for each ionospheric station listed in Table 1, through the application of Eq. (1). As a consequence, τ time series at a fifteen-minutes time interval were obtained for each ionospheric station of Table 1, according to the time span and availability shown in the bottom panel of Fig. 6.
Since our investigation concerns the climatological behavior of τ , we focused on magnetically quiet time periods. Similarly to Pignalberi et al. (2021a), quiet time periods were selected on the basis of the Sym-H (http://isgi.unistra.fr/Documents/References/Iyemori_ et_al_2010.pdf) and AE (Rostoker 1972) geomagnetic indices. Further details about the choice of Sym-H and AE geomagnetic indices are provided in the Sect. 2.2 of Pignalberi et al. (2021a). Geomagnetic activity indices, at one-minute interval, were downloaded from the OMNIWeb Data Explorer website (https://spdf.gsfc.nasa.gov/pub/data/omni/high_res_ omni/). τ values actually employed in this survey are those calculated for the following conditions: −25 nT ≤ Sym-H ≤ 5 nT AND AE ≤ 300 nT.
The solar activity variation of τ was investigated by using the F10.7 81 solar index, i.e., the 81-day running mean of the daily value of the solar radio flux at 10.7 cm (F10.7) wavelength (Tapping 2013). This index is particularly suitable for a climatological study since it smooths the short-time variability of F10.7. F10.7 daily data were downloaded from the OMNIWeb Data Explorer website (https:// omniweb.gsfc.nasa.gov/form/dx1.html). The solar activity dependence has been studied by binning the τ time series according to the following three ranges: The three solar activity ranges were selected to optimize and uniform the data distribution in each of the three considered bins taking into consideration the F10.7 81 values characterizing the last two solar cycles (from 1997 to 2019). The F10.7 81 time series is represented in the top panel of Fig. 6 along with the chosen thresholds.
For each of the thirty-two ionospheric stations considered, and for each level of solar activity, the diurnal and seasonal τ variability were investigated by binning data in bins fifteen-minutes wide in LT, for each month of the year. For each bin, the following statistical metrics were calculated: • 5th percentile, i.e., the lower whisker; • 25th percentile, i.e., the first quartile; • 50th percentile, i.e., the second quartile, representative of the median; • 75th percentile i.e., the third quartile; • 95th percentile, i.e., the upper whisker.
Note that the median is the more appropriate statistical parameter to be used when dealing with climatological studies because it cuts off the outliers of the distribution. The difference between 75th and 25th percentiles defines the interquartile range (IQR) which provides a measure of the data dispersion around the median value, while the whiskers highlight the tails of the distribution.

Diurnal, Seasonal, and Solar Activity Variations for Different Latitudinal Sectors
Due to the extension of the analyzed dataset, in this section we present the results relative to selected ionospheric stations (among those listed in Table 1

Equatorial Sector
The climatological behavior of τ at equatorial latitudes is investigated through the Jicamarca station that has a percentage of available data (40.4%) much greater than Fortaleza (22.2 %) and Sao Luis (31.8%), and it lies right above the magnetic equator (QD Lat. 0.2 • N). We remind to the reader that for Jicamarca station the GNSS considered as co-located is at more than 700 km south of the ionosonde. Especially in this geomagnetic sector this could have an impact on the results. Figure 7 shows the grids of fifteen-minutes binned medians of τ values measured at Jicamarca for LSA, MSA, and HSA. Figures 8-10 show the corresponding boxplots for March, June, September, and December.
Overall, the main results shown by Fig. 7 Figure SM1 in the Supplementary Material. With regard to the solar activity dependence, a positive correlation is found in September in daytime hours and in all seasons in nighttime hours ( Figure SM2 in the Supplementary Material).

Low Latitudes: Main Differences Between Northern and Southern Hemispheres
The results for Ramey (QD Lat. 27.5 • N) and Learmonth (QD Lat. 32.3 • S) stations are shown as representative of the low latitudes for the Northern and Southern hemispheres, respectively. The data availability is 28.8% and 35.7% for Ramey and Learmonth, respectively. Figure 11 shows the grids of fifteen-minutes binned medians of τ values measured at Ramey and Learmonth stations. Figures 12, 13 and 14 show the corresponding boxplots for March, June, September and December. For an easy seasonal comparison between the two stations, the sequence of the months of Learmonth grids and boxplots is from July to June with respect to Ramey's one that is from January to December. The same holds for all the other stations located in the Southern hemisphere.
For LSA the results of Fig. 11 show that the τ median behavior at Ramey is morphologically similar to that of Learmonth in terms of: pre-sunrise peaks observed in all seasons; values at nighttime and around solar terminator hours larger than daytime ones; an evident increase of values between 8-12 LT in summer. For MSA, Learmonth exhibits daytime values higher than Ramey in summer, and also at nighttime in equinoctial months. Conversely, Ramey values around the morning solar terminator hours in spring and winter are considerably higher than Learmonth ones.
Looking at the boxplots for LSA, pre-sunrise peaks are larger at Learmonth than at Ramey during equinoctial months and summer, while a more defined pre-sunrise peak is observed at Ramey in winter. On the contrary, for MSA and HSA, more evident pre-sunrise peaks are observed at Ramey in equinoctial months and winter, while no pre-sunrise peak is  The presence of post-sunset peaks is visible above all at Learmonth in winter, regardless of solar activity. With regard to the main differences between the two sites we can assert that during daytime Learmonth presents values that are systematically higher than those at Ramey and the most significant differences are found in the vernal equinox for LSA and MSA, in summer for MSA and HSA, and in winter for HSA. During nighttime, higher values are found at Ramey in winter, for LSA and MSA. Moreover, in the vernal equinox for HSA, Ramey shows values noticeably higher than those observed at Learmonth. In all the other cases, significantly higher values are observed at Learmonth. The average values of medians observed at Ramey and Learmonth for both daytime (08:00-16:00 LT) and nighttime (21:00-05:00 LT) hours, and for each level of solar activity, as well as the amount of the most significant variations (at least 20 km), are shown in Fig. 15. Finally, for both sites no significant correlation with solar activity is found.

Middle Latitudes: Main Differences Between Northern and Southern Hemispheres
Regarding  Fig. 16, while the corresponding boxplots are shown in Figs. 17, 18 and 19. Port Stanley is the only station whose data are binned at thirty-minutes because of its ionosonde time sampling.
The results of Fig. 16 show a τ median behavior at Boulder qualitatively similar to that of Port Stanley. At both sites, and for each solar activity level, in equinoctial and winter Looking at the boxplots, while post-sunset peaks are clearly visible in winter, the same cannot be said for pre-sunrise peaks. During daytime, values at Boulder are significantly higher than those at Port Stanley independently of solar activity level and season, except in summer for HSA. During nighttime, the situation appears more variegated: for LSA, except for the vernal equinox, higher values are observed at Boulder; for MSA, higher values are observed at Port Stanley for the vernal equinox and at Boulder in winter, while during summer and autumn equinox no significant difference is found; for HSA, except in winter, higher values are observed at Port Stanley. The average values of medians observed at Boulder and Port Stanley for both daytime and nighttime hours and for each solar activity level are reported in Fig. 20, along with their most significant differences (at least 20 km).
With regard to the solar activity dependence, a negative correlation is observed in equinoctial months at Boulder for nighttime hours. Port Stanley instead exhibits a positive correlation in summer for daytime hours, while a negative correlation is observed in the autumn equinox for nighttime hours ( Figure   The results of Fig. 21 show that the median behavior observed at Rome is qualitatively comparable to that observed at Hermanus. The climatological features previously described for Boulder and Port Stanley hold also for Rome and Hermanus. This means that, independently of solar activity, greater values occurring during nighttime and around the solar terminator hours in equinoctial and winter months, a decrease of values during daytime hours, and higher values during the central hours of the day in summer, characterize the climatological trend at middle latitudes. Looking at the boxplots, pre-sunrise and post-sunset peaks are clearly observed only in winter at both sites. With regard to the main differences between the two sites, it can be said that: in daytime hours, apart from the exceptional case observed for HSA in the autumn equinox, whose explanation would require further investigations, higher values are observed in Rome than Hermanus, but with differences not exceeding 33 km; during nighttime hours, except for the summer LSA case, Rome shows values systematically by far higher than Hermanus ones, independently of solar activity. The average values of medians observed at Rome and Hermanus for both daytime and nighttime hours and for each solar activity level are shown in Fig. 25 along with the extent of their most representative variations (at least 20 km).
Concerning the solar activity dependence, during daytime hours no clear dependence is observed for both sites, while for nighttime hours a negative correlation is found in equinoc-

High Latitudes
The Tromso (QD Lat. 66.5 • N) station, with a data availability of 34.9%, has been selected as representative of the τ climatological behavior at high latitudes. Figure 26 shows the grids of fifteen-minutes binned τ medians measured at Tromso for the three solar activity levels, while the corresponding boxplots are shown in Figs. 27, 28 and 29. Figure 26 shows that Tromso is characterized by exceptionally high τ values during nighttime and dawn hours (0-8 LT), occurring mainly at equinoxes and in winter. Except for the summer season, the examination of boxplots also highlights important pre-sunrise peaks more evident for LSA and MSA, while post-sunset peaks are clearly observed only in winter, except for HSA.
In principle, given that Tromso is situated at a very high geographic latitude we would expect a marked seasonal dependence because the intensity of the solar radiation has an important effect on both NmF2 and vTEC. Nevertheless, this is not the case because Tromso median grids offer an overall picture substantially similar to that found at middle latitudes.
During daytime hours, low values are observed in the vernal equinox for LSA and MSA and in the autumn equinox for HSA, while high values are observed in winter for LSA and HSA and in summer for MSA; during nighttime hours, low values are observed in summer, while high values are detected in the vernal equinox. The average values of medians observed for the four selected months, for both daytime and nighttime hours and for the three levels of solar activity, are shown in Figure SM5 in the Supplementary Material. Finally, no correlation with the solar activity cycle is found.

Polar Cap Sector
For the polar cap sector, the choice is limited to the only Thule (QD Lat. 84.5 • N) station, characterized by a relatively high percentage of available data (57.8%).  Figure 30 shows that the results obtained for Thule are by far different than those obtained for Tromso, and in general they are very different from those characterizing other latitudes. In fact, in this case, the trend is strongly driven by the intensity of the solar radiation, thus a very strong seasonal dependence is observed for such extremely high latitudes. For LSA and MSA, the highest values are observed mainly in summer, for both daytime and nighttime hours, when the solar radiation is always present, and the lowest values are observed mainly in winter when the solar radiation is practically absent.
Nevertheless, for HSA, the seasonal variability is less marked and a peculiar diurnal variability in February and March is observed. The explanation of these features needs further specific studies.
Due to the extremely high latitude, there are no effects ascribable to the solar terminator passage; as a consequence, distinct pre-sunrise and post-sunset peaks are never observed. This is also supported by boxplots results that show an almost flat diurnal trend for all seasons.
For both daytime and nighttime hours, minimum values are observed in winter for LSA and MSA and in the autumn equinox for HSA. During daytime hours, maximum values are visible in summer, while during nighttime hours they occur in summer for LSA and in the vernal equinox for MSA and HSA. The average values of τ medians observed for the four selected months, for both daytime and nighttime hours and for the three levels of solar activity, are shown in Figure SM6 in the Supplementary Material. The dependence of τ with solar activity shows a positive correlation, independently of the season and the hour of the day ( Figure SM7 in the Supplementary Material).

The Ionospheric Equivalent Slab Thickness Main Patterns
In the present paper, the τ climatological behavior has been investigated over thirty-two ionospheric stations spread throughout the Northern and Southern hemispheres. Results of Sect. 5 and those shown in the Supplementary Material highlighted very distinct diurnal, seasonal, and solar activity variations for different latitudinal sectors. Some of the outcomes here presented are consistent with the ones achieved in the past by several authors, while others are completely new.
To summarize the main results here presented and discuss them in the light of what is already known from the literature (Sect. 3), in the following the main points concerning the diurnal, seasonal, solar activity, and latitudinal τ variations will be discussed.

Diurnal and Seasonal Variations
Results of Sect. 5 highlighted how the diurnal variation of τ is strongly dependent on latitude and season, so that it is not possible to find a general diurnal trend common to all analyzed stations. However, there are some specific diurnal features shared by many stations over a wide range of latitudes, e.g., daytime/nighttime differences in the τ magnitude, and the presence of pre-sunrise and post-sunset peaks.
The most remarkable diurnal differences are due to the latitude with three main different behaviors for equatorial, low to high, and polar cap latitudes, respectively. Equatorial stations like Jicamarca (see Sect. 5.3.1) exhibit a diurnal trend characterized by a well-defined pre-sunrise peak and daytime relative maximum values, with slight seasonal variations. Stations in a wide range of latitude from low to high (see .4) share a common diurnal pattern characterized by highest values at nighttime and lowest at daytime (with the exception of the summer season), and the presence of pre-sunrise and post-sunset peaks, with local variations induced by seasonal and solar activity dependence. Instead, at polar cap latitudes (see Sect. 5.3.5) the diurnal variation is very faint compared to the seasonal and solar activity ones, and very different compared to other stations.
A remarkable diurnal feature exhibited by τ at many stations is the presence of presunrise and post-sunset peaks. In this regard, results of Sect. 5 highlighted that: • at equatorial latitudes, regardless of solar activity, a clear evidence of pre-sunrise peaks and an absence of post-sunset peaks is observed; • at low latitudes, for LSA, both hemispheres exhibit pre-sunrise peaks in each season, while post-sunset peaks are visible only in the autumn equinox and in winter; for MSA and HSA, pre-sunrise peaks are detected in the Northern hemisphere except in summer, while less distinct pre-sunrise peaks are observed in the Southern hemisphere; post-sunset peaks are detected only in winter in both hemispheres; • at middle latitudes, pre-sunrise peaks are barely visible at equinoxes, but clearly defined in winter, for different solar activity; regardless of solar activity, post-sunset peaks are found for both hemispheres in winter; • at high latitudes, for LSA and MSA, pre-sunrise peaks are mostly visible at equinoxes and in winter, while post-sunset peaks only in winter; • at polar cap latitudes, neither pre-sunrise nor post-sunset peaks are detected, independently of solar activity.
The results for the equatorial latitude station of Jicamarca (QD Lat. 0.2 • N) are fully in accordance with those achieved by Kenpankho et al. (2011) andVenkatesh et al. (2011) who detected pre-sunrise peaks in all seasons over Chumphon (QD Lat. 3.1 • N) and Trivandrum (QD Lat. 0.9 • N), respectively. Likewise, the occurrence of low-latitude pre-sunrise and post-sunset peaks observed in winter at Ramey (QD Lat. 27.5 • N) and Learmonth (QD Lat. 32.3 • S) are consistent with the winter pre-sunrise peaks observed over Hawaii (QD Lat. 19.3 • N) islands ) and equinox/winter pre-sunrise and post-sunset peaks over Nicosia (QD Lat. 29.2 • N) (Haralambous 2011). At the middlelatitude stations considered in this study, peaks are observed at equinoxes and more clearly in winter with a tendency to be less pronounced as the solar activity increases. These facts agree with the results obtained by several authors (Titheridge 1973;Goodwin et al. 1995a; A. Pignalberi et al.  The physical mechanisms that can concur in the explanation of peaks observed in this investigation are: • since the equatorial ionosphere is not connected to the plasmasphere by magnetic field lines, there is no plasma replenishment to balance the loss of electrons in the lower ionosphere due to the recombination (Kenpankho et al. 2011); this circumstance could play a role for the observed equatorial pre-sunrise peaks; • during the early evening hours, the equatorial pre-reversal electric fields produce an upward vertical E × B drift causing a plasma fountain rising to heights larger than 1200 km at the equator, with accompanying plasma density depletions in the F region (Balan and Bailey 1995;. After the pre-reversal strengthening, the E × B drift rapidly reverses causing a downward plasma flux which is forced by the neutral winds to move away from equator, thus increasing the topside electron content in the range of magnetic latitude between 10 • S and 20 • S, and between 10 • N and 20 • N. The plasma depletion in the F region created by the pre-reversal strengthening is instead still maintained in the range of magnetic latitude between 30 • S and 30 • N . This would explain why at Jicamarca, as well as at Fortaleza (QD Lat. 7.1 • S), an increasing trend of τ with a clear peak is practically never observed during post-sunset hours, while at Cachoeira Paulista (QD Lat. 18.8 • S), just located near the equatorial anomaly crest, an increase of τ with well-defined peaks is observed especially in the equinoctial and summer months around 20-21 LT (see Supplementary Material). This circumstance is also in agreement with the nighttime phenomenon of the equatorial electron temperature anomaly, occurring predominantly at equinoxes and clearly associated to the deposition of plasma in the range of latitude |10 • |-|20 • | resulting from the sunset electrodynamical processes ) strictly connected to the existence of strong electric fields over the equatorial latitudes (Modi and Iyer 1989); • due to the field-aligned plasma flow from the magnetosphere to the ionosphere, large downward fluxes of protons are injected in the topside at pre-sunrise hours, thus lowering the upper transition height, i.e., the height where [O + ] is equal to [H + ] (Mosert et al. 2013). Evans and Holt (1978) revealed that the duration of such fluxes in the winter hemisphere, in particular during the post-midnight hours, decreases as the solar activity increases. Belehaki et al. (2003) found an increased plasmaspheric content during nighttime. These findings indicate an increase of the topside ionosphere content, thus explaining the presence at the various latitudes of the sharp pre-sunrise and post-sunset peaks observed in the winter season especially for LSA; • at solar terminator hours, the ionosphere is sunlit differently at different altitudes; then, NmF2 and vTEC manifest a time lag which is well reflected in τ . For example, at sunrise, higher values of τ mean that NmF2 is low while vTEC is increasing due to the solar radiation impinging from above. The reverse is usually observed at sunset hours; • changes of the meridional neutral winds direction most pronounced in winter months, and enhanced zonal electric fields, are responsible for the downward movement of the ionosphere. Consequently, the abrupt increase in the loss rate and the corresponding rapid decrease of NmF2 before sunrise, while smaller changes in the topside ionosphere occur, would explain the winter peaks observed at middle latitudes (Titheridge 1973); • the so called morning overshoot phenomenon, namely, the rapid increase of the electron temperature observed around sunrise (Oyama et al. 1996), which increases the plasma scale height, thus justifying the observed pre-sunrise peaks. Specifically, as observed around the magnetic equator in the range of magnetic latitude between 10 • S and 10 • N for the period 05:30-06:30 LT, the enhancement of the morning overshoot following the combined action of the downward plasma vertical E × B drift and neutral winds (Oyama et al. 1996) could play a role in the development of pre-sunrise peaks observed in all seasons at Jicamarca.
A quite common diurnal feature shared by many stations from low to high latitudes is the daytime/nighttime differences in τ absolute values. Specifically, results of Sect. 5 highlighted how, from low to high latitudes, nighttime values are greater than daytime ones in equinoctial months and even more in winter months, while the reverse is true for summer months.
The highest daytime (10-16 LT) values, regardless of solar activity, are always observed at the magnetic equator (e.g., Jicamarca station), except in winter. Such feature constitutes a second maximum, in addition to the pre-sunrise one, that is not noticed in the other sites The fact that during daytime hours the largest values are observed near the magnetic equator is also confirmed by Kenpankho et al. (2011) who demonstrate that the slab thickness at Chumphon is much larger than that found at low, middle, and high latitudes.
The two distinct maxima observed during sunrise hours and around the midday at Jicamarca well agree with the results of Venkatesh et al. (2011) who found strong pre-sunrise peaks with a secondary maximum around noon-time hours during all seasons at different Indian equatorial stations. This fact confirms that the maximum around the central hours of the day is a typical feature of the equatorial latitudes. This unusual feature of the equatorial region can be explained in terms of the fountain effect: the equatorial vertical E × B drift during the morning hours is positive (Fejer et al. 1995;Oyama et al. 1996;, which means that the ionospheric plasma is lifted upward. The resulting effect is that while NmF2 does not show any significant increase (Huang 1974), an enhancement of the topside electron content occurs, thus leading to an enhancement of τ over the magnetic equator (Das Gupta et al. 1975).
The diurnal variations observed at low, middle, and high latitudes, characterized by presunrise and post-sunset peaks and generally higher nighttime values than daytime ones during the equinoxes and winter (while the opposite occurs in summer), are in agreement with the results of other authors (Chuo et al. 2010;Guo et al. 2011;Mosert et al. 2013). The explanation of the diurnal trends can be given regarding them as a superposition of two diurnal curves: a first one which maximizes at noon corresponding to the typical thermospheric variation, and a second one with a double peak at sunrise and at sunset which corresponds to the plasmaspheric variation (Belehaki et al. 2003) which is most prominent in winter and it is due to the lowering of the upper transition height (Mosert et al. 2013). Anyway, it is worth Fig. 26 Same as Fig. 7, but for Tromso noting that the unusual larger daytime values observed for HSA in the autumn equinox at Hermanus constitute a feature that has no clear explanation.
The extremely high values observed at Tromso for each level of solar activity, and the relatively high values noticed at Thule, for HSA, are in agreement with the outcomes of Huang et al. (2016) who, analyzing global maps of τ reconstructed from vTEC data obtained from GIM and NmF2 retrieved from COSMIC data, found particularly high values at high and polar cap latitudes for all seasons and for both minimum and maximum solar activity.
The really high values observed at Tromso, in the nighttime and solar terminator hours, could be due to the high-latitude ionospheric trough phenomenon, that is a limited region of depleted plasma density occurring in the ionosphere within the auroral oval or in the polar cap (Rodger et al. 1992;Voiculescu et al. 2016). Troughs were originally considered as phenomena occurring at night during winter (Moffet and Quegan 1983;Rodger et al. 1992). More recent investigations, based on satellite records, have shown that troughs form during all seasons and also during daytime (Aa et al. 2020;Karpachev 2019). Ma et al. (2000) using European Incoherent Scatter Scientific Association (EISCAT) observations showed that high-latitude ionospheric troughs occur in the dawn sector, extending from early morning to pre-noon, but sometimes also in the afternoon and pre-midnight. Winser et al. (1986) again using EISCAT observations found an ionospheric trough characterized by a drop of electron density of 54% on 6 March 1984 right over Tromso. These studies support the idea that the ionospheric trough can play a significant role in causing an enhancement of τ at high latitudes, as that observed at Tromso.

Solar Activity Variations
Results of Sect. 5 highlighted how the solar activity affects the τ magnitude. Except for the very high-latitude station of Thule, no general solar activity trend in τ can be found, with variations that are function of local time, season, and latitude. The most remarkable variations induced by the solar activity are: • at equatorial latitudes, for nighttime hours, a positive correlation is observed in all seasons; • during daytime hours, from equatorial to high-latitude sectors, no significant correlation is observed in both hemispheres, except for rare cases for which a positive correlation is found; • during nighttime hours, negative correlations are found essentially in equinoctial months for both Northern and Southern middle latitudes; • at polar cap latitudes, for both daytime and nighttime hours, a positive correlation is observed in all seasons.
This study showed that at equatorial, low, middle, and high latitudes a correlation of τ with solar activity is practically absent in daytime except for a few cases. The result obtained in the vernal equinox at Jicamarca (QD Lat. 0.2 • N), for which a positive correlation is observed, is consistent with that of Kenpankho et al. (2011) who, analyzing data over Chumphon (QD Lat. 3.2 • N), found in equinoctial months a diurnal variation of τ in the years 2004-2006 clearly increasing for years of higher solar activity. The absence of any correlation at low latitudes is in contrast with both the positive correlation observed at Hawaii (QD Lat. 19.3 • N) by  and with the results of Dabas et al. (1993) who, from the examination of data over Lunping (QD Lat. 13.3 • N), observed a positive correlation between τ and three different solar indices. However, our results are in accordance with Chuo (2007) who, analyzing data over Chung-Li (QD Lat. 18.3 • N), found a poor correlation with the solar activity index F10.7.
At middle latitudes, except for the winter case of Port Stanley (QD Lat. 38.7 • S) for which a positive correlation was noticed, our results are controversial when compared to those of other authors. In fact, the absence of correlation in daytime hours is in sharp contrast with the positive correlation found by Stankov and Warnant (2009) at Dourbes (QD Lat. 45.8 • N), and by Davies and Liu (1991) at Sagamore Hill (QD Lat. 51.3 • N). At high latitudes, the lack of any correlation between τ and solar activity observed at Tromso (QD Lat. 66.5 • N) is not in accordance with the positive correlation obtained by  at Goose Bay (QD Lat. 60.2 • N). On the contrary, the results found for Thule (QD Lat. 84.5 • N) indicate a positive correlation with solar activity in all seasons and match with the observations carried out at the Antarctic station of Casey (QD Lat. 80.6 • S) for which a positive correlation of τ was found with each of the three types of solar radiation considered: F10.7, X-ray flux, and EUV flux (Yadav and Bhawre 2020).
The situation is different during nighttime hours when negative correlations are observed at middle latitudes, especially in equinoctial months. Such results are in contrast with those of other authors. For example, both Stankov and Warnant (2009) at Dourbes and Davies and Liu (1991) at Sagamore Hill obtained a positive correlation.
Differently from the Northern hemisphere, in the Southern hemisphere the results found are somewhat consistent with what can be read in the literature. For instance, the results by Breed et al. (1997) agree with ours, even though in our case the negative correlation was found mostly in equinoctial months. The physical mechanism that could explain the observed middle-latitude negative correlations is the following: a lowering of the [O + ]-[H + ] transition level in the upper F2 region for low solar activity due to large downward fluxes of [H + ] from the magnetosphere, occurring generally in winter during post-midnight hours, whose duration decreases as the solar activity increases (Evans and Holt 1978), can increase the topside electron content, which results in higher τ values at a time of reduced solar activity (Breed et al. 1997).
The circumstances above mentioned tell us that the behavior of τ for different latitudinal sectors varies in a complicated way in relation to the solar activity cycle, depending on the location of the observing station. An easy classification is then not possible and this topic deserves attention for future studies.

Latitudinal and Inter-Hemispheric Variations
τ exhibits large variations with very different diurnal and seasonal patterns as a function of the latitude that show some characteristics features as: • a decreasing trend of τ from equatorial to low latitudes is generally observed during daytime hours; • in the Northern hemisphere, during daytime hours, a tendency of τ to increase from low to high latitudes is detected for LSA in the vernal equinox, for MSA at both equinoxes, and for HSA in summer; while in winter no significant variations are observed between low and middle latitudes; • in the Northern hemisphere, during nighttime hours, an increasing trend of τ from the equator to middle latitudes is observed for LSA in the autumn equinox and in winter, and for MSA in winter.
The decreasing trend of τ in daytime hours as the latitude increases is deduced comparing the values observed at Jicamarca with those observed at Ramey and Learmonth. A further confirmation comes from the fact that the daytime values of τ at Jicamarca are almost always higher than those detected at Cachoeira Paulista, a station located near the crest of the equatorial anomaly (see the Supplementary Material).
Our results agree with those by Das Gupta et al. (1975) and Venkatesh et al. (2011). The former, analyzing some low-latitude stations, found that during daytime in winter τ falls off rapidly as one moves away from the magnetic equator, attaining a minimum value at the stations situated below the crest of the equatorial anomaly; the latter, examining the hourly τ behavior of some equatorial and low-latitude stations, found τ values decreasing as the latitude increases, for all seasons.
The observed daytime decreasing latitudinal trend can be explained in terms of the fountain effect: the vertical electrodynamic drift lifts upward the equatorial ionospheric plasma, which diffuses along the magnetic field lines to be deposited at low latitudes. This causes an enhancement of plasma around the crests of the equatorial anomaly and then a reduction of τ , which maximizes around noon when the equatorial anomaly is fully developed (Das Gupta et al. 1975). Anyway, the analysis about the latitudinal variation revealed that τ can behave in a complex manner and a clear latitudinal signature does not emerge. This circumstance is also confirmed by the studies carried out by  who, analyzing the τ behavior at 14 world-wide stations located at different latitudes for high solar activity, observed a clear decreasing trend of τ as the latitude increases only in winter during daytime.
Another interesting spatial variation is the presence of hemispheric differences, the most remarkable ones are: • at low latitudes, nighttime and daytime values in the Northern hemisphere are smaller than those in the Southern one; • at middle latitudes, for both LSA and MSA, during daytime and nighttime hours, Northern values are clearly greater than Southern ones; while no clear pattern emerges for HSA.
The Northern-Southern differences observed at low latitudes by day and night can be explained in terms of changes of the neutral winds pattern between the two hemispheres (Oyama et al. 1996;, which causes large asymmetries in the plasma fountain over all sites for both daytime and nighttime hours . At least for daytime hours, the observations of values at the Northern middle latitudes greater than those at the Southern ones for all seasons are consistent with the results of Guo et al. (2011) obtained through the investigation of maps obtained by GPS ground-based receivers and COSMIC radio occultation measurements.

The Intrinsic Variability of the Slab Thickness
In our analysis we also highlighted the dispersion of τ values around the medians as represented by the boxes (IQR) and the whiskers in the boxplots shown in Sect. 5, a topic often overlooked by past studies. The investigation of τ dispersion around median values shed light on its intrinsic variability due to both physical phenomena and to the definition of τ itself as a function of both a punctual quantity (NmF2) and an integrated one (vTEC).
From our analysis it emerges that: • a greater dispersion is clearly observed at night and at both dawn and dusk hours. Nevertheless, in some summer cases, at low and middle latitudes, the dispersion characterizing daytime hours is larger, or at most comparable, than that observed during nighttime and solar terminator hours; • at high latitude, except for summer, a dispersion exceptionally high characterizes the nighttime and solar terminator hours. In particular, for HSA, a large dispersion is visible during daytime in the vernal equinox and winter; • at polar cap latitudes, except for some cases, no significant difference among nighttime, daytime, and solar terminator hours is observed.
The most relevant fact is the very high dispersion of τ values at nighttime and for solar terminator hours. The main reason of this behavior lies in the τ sensitivity to vTEC changes at small NmF2 values, like discussed in Sect. 2 (see Fig. 4). In fact, the nighttime and solar terminator hours are those for which the plasmasphere gives the most relevant contribution to the vTEC (Klimenko et al. 2015), and NmF2 is low. In such circumstances, even small vTEC changes cause large τ variations thus increasing its corresponding dispersion. Such conditions are even exacerbated at high latitudes, due to the lowest NmF2 values, and produce extremely high τ values as those recorded at Tromso. Moreover, part of the dispersion has to be ascribed to the intrinsic day-to-day variability exhibited by both vTEC and NmF2. As a consequence, stations located in regions characterized by large day-to-day ionospheric variations, such as the equatorial and low-latitude regions, as well as the auroral ones, are particularly affected by an increased dispersion of τ values around the median climatological trend. There is no doubt that this very high intrinsic variability of τ represents one of the toughest challenges in the τ modeling.

Current Gaps and Future Developments
In this work, an outstanding dataset of vTEC and NmF2 values has been considered to study the ionospheric equivalent slab thickness over the globe. The variability of τ in terms of local time, season, and solar activity has been studied with great accuracy at each selected location due to the availability of extremely long and complete time series, spanning almost two solar cycles. Differently, due to the lack of relevant instruments and due to the very strict requirements adopted for the definition of the dataset used, it is difficult to achieve the same level of specification with reference to the geographic coordinates. It would be therefore desirable to increase the knowledge of τ at low latitudes, for example in the African and Asian longitudinal sectors, with particular focus on the trough and the crests of the EIA, because of the relationship of τ with the ionosphere plasma scale height. Then, as it has been demonstrated by the very peculiar case of Thule, specific efforts should be devoted to study the behavior of τ at extremely high latitudes, where the slab thickness patterns are completely different from those at the other latitudes. The ever growing number and distribution of ground-based GNSS stations will allow covering most of the globe in the coming years; however, the lower availability of ionosonde stations compared to GNSS stations, with their uneven spatial distribution, represents a strong limitation for the improvement of the description of the τ spatial variations. The current distribution of ionosondes (see https://lgdc.uml. edu/common/DIDBFastStationList) is heavily unbalanced towards the Northern hemisphere middle latitudes, with a poor coverage of the low and high latitudes, and, most importantly, of the oceans where only a few island stations are installed. Then, in order to properly describe the τ spatial variations and for the implementation of a global τ model based on these observations, a higher number and a more uniform distribution of both ground-based GNSS and ionosonde stations is a requirement to achieve in the near future. A different approach relies on the use of near real-time NmF2 and vTEC global maps to obtain a global representation of τ (Froń et al. 2020;Galkin et al. 2022); however, the effectiveness of methodologies based on maps is still to be validated. Compared to the robust procedure based on co-located ground-based GNSS and ionosonde stations to obtain τ (as described in Sect. 5), methodologies based on global maps are possibly prone to misrepresentations due to the inherent smoothing of the mapping procedures, to the uneven distribution of stations used to produce corresponding maps, and to the data-assimilation schemes applied to produce maps (Pignalberi et al. 2021c). This is why a growing number of ionospheric facilities with a more uniform spatial distribution is a key point for a better comprehension of many τ variations.
As a further step, a deeper analysis should be carried out on the available data to better understand the τ dependence on solar activity since no clear schemes, common to all locations, have been identified so far. An evident problem in the study of solar activity variations on τ is the unavailability of time series covering enough solar cycles to properly catch this faint dependence, if any. As evidenced in Sect. 2 and by Pignalberi et al. (2021a), it is not possible to merge the oldest τ time series based on vTEC from Faraday rotation technique with newest ones based on vTEC from GNSS satellites; as a consequence, time series embracing at most three solar cycles can be obtained from both historical (Faraday rotation, from 60s to 90s) and modern observations (GNSS receivers, from 90s onwards). The shortness of available time series limits the application of spectral time series analyses to highlight solar activity variations around the eleven-year period. Different statistical approaches have to be explored in the future to solve the puzzle of the τ solar dependence. An attempt should also be made to establish a relationship between τ and geomagnetic activity, being aware of the intrinsic limitations posed by the scarce availability of case studies. Finally, an assessment about a possible influence on τ of the Ionospheric Winter Anomaly (Yasyukevich et al. 2018) or the Weddell Sea Anomaly (Chen et al. 2019) could also be envisaged.
In Sect. 4 it has been underlined the importance of the slab thickness parameter because of its relation to several physical processes in the ionosphere, to its capability to contribute to Space Weather monitoring and to ionospheric modeling, including the realization of fast data-ingestion procedures allowing, e.g., the retrieval of foF2 when vTEC is available. As far as future developments are concerned and considering that very few models of slab thickness currently exist, an attempt to build an improved global τ climatological model is then planned. Due to the geographic distribution of the available data and the remarkable difficulties to be overcome in order to implement a global model, a stepwise approach might be undertaken, with the initial implementation of local (single station) models, followed by the implementation of regional ones (based on multiple stations over a continental or subcontinental area). It is expected that diverse mathematical techniques could be used to accomplish the proposed task, especially in relation to the construction of a global model. Indeed, a classical approach based on Fourier series expansion for time evolution, and functions similar to Spherical Harmonics (Jones and Gallet 1962) to describe the spatial behavior of τ could be studied. Nevertheless, taking advantage of the newly developed technologies, the utilization of Machine Learning techniques (Yu and Ma 2021) should also be investigated, to go beyond the classical concept of median model and to likely introduce improved forecasting capabilities.
An intriguing and so far poorly explored field of research is the inclusion of a τ model into an electron density empirical model like NeQuick, even if this might appear as a contradiction because, as it has been pointed out in Sect. 4, empirical models like IRI and NeQuick are already able to provide estimates of τ as an indirect value obtained from the ratio of vTEC and NmF2. As a first extent, such a τ model could be used as a reference, or background solution, to provide a "warning" should the electron density model τ be far from the expected climatological value. As a further step, a procedure to modify NeQuick specific parameters allowing the model to match the expected τ value could be implemented. This task could be accomplished by taking advantage of the concepts used to develop the dataingestion techniques based on the use of multiple effective parameters (Nava et al. 2011). It is understood that the implementation of such procedure would require heavy modifications of the NeQuick source code, but it would allow correcting the model τ a posteriori, keeping the internal consistency of all the parameters values used in the model without modifying its analytical formulation.

Conclusive Remarks
In the present paper, first a revision of the literature about the ionospheric equivalent slab thickness has been done and then its climatological behavior has been investigated on the basis of thirty-two ionospheric stations distributed all over the world, for the last two solar cycles, which is the largest and most complete dataset compiled so far. The statistical treatment applied to both NmF2 and vTEC datasets makes the datasets of τ obtained for the various stations extremely reliable and valuable for the characterization of the τ behavior at very-large spatial and temporal scales. The climatological analysis highlighted very specific diurnal, seasonal, and solar activity variations depending on the location of the observing station, that were critically discussed in the light of the existing literature.
The extension and completeness of the dataset used in this investigation makes it appropriate for the development of a global empirical climatological model of τ which, besides representing the main climatological features, would overcome also the limitations characterizing the models developed in the past. For example, the Fox et al. (1991) regional model, being based on the observations of the Faraday rotation of signals from geostationary satellites, provides a relatively large underestimation of τ measurements, as recently highlighted by Pignalberi et al. (2021a). The Jakowski and Hoque (2021) global model, besides describing the main spatial, diurnal, and seasonal τ variations, can hardly represent all the variation and dispersion exhibited by τ , particularly during non-equilibrium conditions. Therefore, an empirical climatological global model of τ based on a very-large dataset derived by co-located ionosondes and ground-based GNSS receivers would be of utmost importance considering that the most important global empirical climatological models of the ionosphere, as NeQuick (Nava et al. 2008) and IRI (Bilitza et al. 2017), do not provide any direct τ modeling. Moreover, Space Weather applications relying on data-assimilation techniques of either vTEC or NmF2 ionospheric parameters would benefit of a reliable and robust climatological global model of τ .

Appendix: Supplementary Material Description
The statistics shown in Sect. 5 have been calculated for each ionospheric station listed in Table 1, and are available in the Supplementary Material.
Specifically, for each ionospheric station and level of solar activity the Supplementary Material includes: a) grids of fifteen-minutes binned counts of measurements on which the median value is calculated; b) grids of fifteen-minutes binned medians of measurements; c) grids of fifteen-minutes binned IQR values; d) boxplots of fifteen minutes binned measurements.
These plots are collected in 32 folders, one for each ionospheric station, with the corresponding station name. Funding Note Open access funding provided by Istituto Nazionale di Geofisica e Vulcanologia within the CRUI-CARE Agreement. This research is partially supported by the Italian MIUR-PRIN grant 2017APKP7T on Circumterrestrial Environment: Impact of Sun-Earth Interaction.

Author Contribution
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.