## 2 Data and Methods
### Multi-Parameter Data and Processing  
The complete route of ionospheric disturbances is well explained in Figure 2. GNSS data from ground-based Continuously Operating Reference Stations (CORSs), GIMs from the Center for Orbit Determination in Europe (CODE), foF2, hnn2, and ionograms from Digital Ionosonde Database (DID)--Global Ionosphere Radio Observatory (GIRO), and Ne and Te derived from Swarm satellite measurements were used to investigate ionospheric disturbances in Afghanistan from a perspective of multi-parameter analysis. GNSS stations within the Dobrovolsky distance of the epicenter were considered [47]. Figure 1: Tectonic settings of Afghanistan earthquake cases along with ionosondes and GNSS ground stations. In the analysis of GPS data, a hybrid methodology involving machine learning, statistical analysis, and mathematical techniques was utilized. The 30 s interval TEC values were acquired from the available IGS databases of selected ground stations. GPS satellites transmit radio signals at two long-wavelength band frequencies: L1, which operates at 1.57542 GHz, and L2, which operates at 1.2276 GHz [37; 38]. Satellite and receiver biases were removed, and Vertical TEC (VTEC) was derived from absolute STEC, which is measured in TECU (total electron content unit = 10\({}^{16}\) e/m\({}^{2}\)). The STEC was calculated and converted into VTEC using GPS-TEC software v3.03. A hybrid machine learning algorithm, LSTM, was applied to train and predict the data over a 24 h period. By employing this machine learning-based envelope technique, a normal distribution of biased and unbiased VTEC values was opted for to establish a threshold of \(\pm\)1.34 u in order to achieve a higher level of confidence of almost 82% [48]. Any observed TEC values outside the range of \(\pm\)1.34 u signify the detection of abnormal signals [49]. Almost 10 GPS stations nearby the study area were processed, and those ground stations which were spatially closer to the epicenter, such as KIT3 (Uzbekistan), IISC (India), WUHN (China), LHAZ (China), CHUM (Kazakhstan), and POL2 (Kyrgyzstan), were considered for the actual anomaly estimation based on availability (see Table S2). This unique estimation method helped to accurately define abnormalities with the least data losses and RMS error rather than the historical methods of TEC anomaly estimations. GIMs, generated by CODE, provide daily updates with a resolution of 5\({}^{\circ}\) longitude by 2.5\({}^{\circ}\) latitude, offering valuable insights into ionospheric conditions. To clearly identify the ionospheric anomalies, the sliding interquartile range (SIQR) method was used for GIM-TEC processing. The SIQR method utilizes data from CODE to identify TEC disturbances associated with earthquakes, establishing a baseline reference from pre-earthquake or stable periods and calculating the interquartile range (IQR) to determine abnormal ionospheric parameter values [33]. The IQR approach identifies anomalies that may not have a significant impact on the quartiles, ensuring a thorough examination. Data from individual GNSS receivers and GIM were used to account for both local and large-scale fluctuations. This technology increases the dependability and accuracy of our anomaly detection system, offering important insights into ionospheric anomalies and their potential relationship with seismic activity [50]. For calculating the interquartile range (see Equations (1) and (2)), the  
Figure 2: The technical route of excluding, processing, and verifying the ionospheric anomalies prior to seismic activities. median of 2 consecutive days, i.e., Q\({}_{1}\), was subtracted from the 2-consecutive-day median of Q\({}_{3}\), and the upper and lower bounds were calculated by adding and subtracting the 1.4 IQR from Q\({}_{2}\). \[ \begin{split}\text{Q}_{1}=\frac{\text{x}_{7}+\text{x}_{8}}{2}, \text{Q}_{2}=\frac{\text{x}_{14}+\text{x}_{15}}{2},\text{Q}_{3}=\frac{\text{x}_ {21}+\text{x}_{22}}{2}\end{split}