Over the last decade, thanks to the availability of historical satellite observations that have begun to be significantly large and thanks to the exponential growth of artificial intelligence techniques, many advances have been made in the detection of geophysical parameters such as seismic-related anomalies. In this study, the variations of the ionospheric Total Electron Content (TEC), one of the main parameters historically proposed as a seismic-connected indicator, are analyzed. To make a statistically robust analysis of the complex phenomena involved, we propose a completely innovative machine-learning approach developed in the R programming language. Through this approach, an optimal setting of the multitude of methodological inputs currently proposed for the detection of ionospheric anomalies is performed. The setting is optimized by analyzing, for the first time, multi-year-mostly twenty-year-time series of TEC satellite data measured by global navigation satellite systems (GNSS) over the Italian region, matched with the corresponding multi-year time series of seismic events. Seismic events including all the countries of the Mediterranean area, up to Turkey, are involved in the analysis. Tens of thousands of possible combinations of input methodological parameters are simulated and classified according to pre-established criteria. Several inputs examined return clear results. These results combined with each other highlight the presence of anomalous seismic-related sequences that have an extremely low probability of having been detected randomly (up to 2 out of 1 million). The anomalies identified represent the most anomalous behaviors of the TEC recorded during the entire period under investigation (e.g., 20 years). Some of the main conclusions are that, at mid-latitudes, ? the detection of seismic-TEC anomalies can be more efficient looking for punctual rather than persistent phenomena; ? the optimal thresholds for the identification of co-seismic anomalies can assume different values depending on type of anomaly (positive or negative) and type of observation; ? single GNSS receiver data can be useful for capturing local earthquake-ionospheric effects and Global Ionospheric Maps (GIM) data can be functional in detecting large-scale earthquake-ionospheric effects; ? earthquakes deeper than 50 km are less likely to affect the ionosphere.
