## 4 Methods and Data Analysis
### Data Analysis
Rinex-format daily GPS observation files, with 30s sampling rate, were used as input; a cut off angle of \(7^{\circ}\) was set for the processing. The ocean load was modelled using FES2004 coefficients [35]. GAMIT/GLOBK uses a typical two-step procedure to process the data, following the DGNSS differential technique. Firstly, the \"ionofree\" linear combinations (LC) are computed to account for the ionospheric delay, then, with IGS-precise satellite orbits, 10 ppm loose constraints coordinates of GPS stations are retrieved. In the second step, the Kalman filter procedure, implemented in the GLOBK program, is applied to retrieve the final coordinates in ITRF2014 [36]. Eight IGS stations (AJAC, MATE, NOT1, POTS, GRAS, GRAZ, WTRZ, ZIMM) were included in the analysis to achieve the frame referencing. Station-specific discontinuities due to updates of antenna calibration parameters or associated with changes in instrumentation (antenna and/or receiver) are accounted for in the first step (GAMIT) and those found unaccountable were corrected in the GLOBK time series. For the purpose of estimating tropospheric delay parameters during episodes of heavy rainfall, in addition to GPS data, meteorological observations (atmospheric pressure, temperature, relative humidity) acquired at the meteorological station collocatedwith the MAFE station were incorporated in the analysis. Moreover, according to [37], GNSS and meteorological observations from a remote station (distance > 200 km), MATE (Matera, Italy), were included in the processing to avoid correlation of the tropospheric delay parameters within the small NeVoCGPS network. Besides quasi-observations, the daily GAMIT solutions contain the parameters of the piecewise-linear model estimated from the data, hence, the zenith total delay (ZTD) is obtained by interpolation. PWV is extracted from the GAMIT solutions with \"sh_metutil\" [33], available in the script library of GAMIT software sub-package. This tool extracts the ZTD values from the daily solutions, applies corrections for the ZHD [23] to retrieve the ZWD values and converts them into PWV; the NS and EW horizontal gradients with associated errors are computed as well. In studying the thunderstorms, we have even used the freely available zenith total delay (ZTD) and PWV solutions produced by the Nevada Geodetic Laboratory (NGL), which employs the GIPSY/OASIS-II software in a PPP strategy [34]. The available solutions have a 5-min temporal resolution and the observation cut-off angle was 7\({}^{\circ}\). The satellite clocks, receiver clocks, station coordinates, integer ambiguity and total zenith delays were estimated. Finer details on the processing strategies and setting parameters can be found at NGL website ([[http://geodesy.unr.edu/magnet.php](http://geodesy.unr.edu/magnet.php)]([http://geodesy.unr.edu/magnet.php](http://geodesy.unr.edu/magnet.php)), accessed on 15 September 2021). As a further alternative method of analysis, the CSRS-PPP service was used. It relies on precise satellite orbit, clock and bias corrections, derived from a global network of receivers, to determine accurate user positions. For the purpose of this study, the latest version of the service (version 3) was used, which implements several ambiguity resolution (AR) algorithms (PPP-AR), allowing faster convergence with external ionospheric information and the processing of multi-GNSS observations. Currently, this new algorithm achieves millimetre-level accuracies for long observation sessions (daily or longer) in static mode, but an accuracy of a few centimetres can typically be attained in 1-2 h. Further details can be found in [32]. The CSRS-PPP solutions have a 30-s temporal resolution and the observation cut-off angle was 7.5\({}^{\circ}\).