## Appendix B Detailed Description of Swarm Magnetic Data Analysis  
In this appendix, we describe the Swarm magnetic data analysis and processing in detail. The same procedure is applied to each component of the magnetic field in NEC (North-East-Centre) frame and to the scalar intensity (F) of the Earth Magnetic Field measured by three Swarm satellites (Alpha, Bravo, and Charlie). The raw data used for this analysis is the LR (Low resolution) data from ESA and we show an example of the step-by-step processing in Figure 11 of track 25 of Swarm Bravo acquired on 18 June 2021. This is the one shown in the main text in Figure 26. The figure shows the example of Y-East component of the magnetic field. The algorithm used to process the data is called MASS and was used for the first time to analyse the M7.8 Nepal 2015 earthquake by [PERSON] et al. [66]. Following that, it has been widely applied to study several earthquakes in the world with a magnitude between 5.5 and 8.3 [4,5,40,67,68,95]. The operation made by the algorithm on the magnetic signal \(B(t)\) consists of the following:Estimation of the numerical first order derivative of the original data \(B^{\prime}(t)\). This step is done by calculating the differences sample by sample (\(S\)) and dividing the result by the time \"\(t\)\" (i.e., the time difference between the two consecutive data): \[\frac{dB}{dt}\cong\frac{B_{S}-B_{S-1}}{t_{S}-t_{S-1}}=B^{\prime}(t)\]
2. In order to remove the further trend from the data, a cubic spline \(S(t)\) was fitted on \(B^{\prime}(t)\), setting the knock points every 20 samples: \[\textit{residual}\ (t)=B^{\prime}(t)-S(t)\]
3. The geomagnetic latitude of each sample was computed considering the position of the magnetic pole provided by the IGRF-13 geomagnetic model [96], and the samples outside the range [\(-\)50, \(+\)50] of geomagnetic latitude were rejected.
4. The Root Mean Square \"_RMS_\" of the whole track is calculated and compared with the root mean square \"_rms_\" of a moving window with a length for this work of 2.0\({}^{\circ}\) latitude. If \[rms>kt\cdot RMS\] the window is defined as anomalous. \"\(kt\)\" is a threshold to define the anomalies and, in this work, has been set to 2.5.
5. A track is classified as anomalous if it contains at least one anomalous window whose centre falls inside the research area (a circle of 420 km of radius), and geomagnetic conditions were quiet (\(\,\mid\) Dst\(\mid\,\leq\) 20 nT and ap \(\leq\) 10 nT).
6. The cumulative curve of all the identified anomalous tracks was computed for each component of the geomagnetic field, combining the three satellites together.  
In order to establish if the cumulative trend is anomalous, a deviation from a linear increase is searched. In fact, this algorithm identifies anomalous signals in the ionosphere during quiet geomagnetic time and a background of anomalies is also expected without the geophysical event as for example simulated in [67]. On the other hand, a deviation in terms of acceleration from the linear trend can indicate that a critical process, such as an earthquake or volcano eruption, is approaching [97; 98].  
18 June 2021 from the original signal (**A**) that has been derived (**B**), a cubic spline was further removed (**C**) and analysed the root mean square by a moving window (**D**).