## 2 Datasets and Methodology
### Model Setup and Lightning Data Assimilation Procedure  
The numerical model used in this study is the WRF model with advanced WRF dynamic (WRF-ARW), version 4.1.3 [72]. The simulations were conducted on one single domain, with 531 \(\times\) 531 grid points and 42 vertical levels, with a model top at 50 hPa. The model domain (Figure 1) covers the central Mediterranean and the whole Italian territory and has a horizontal grid spacing of 3 km. The physical schemes employed include the new Thompson microphysics scheme [73], the Yonsei University boundary layer scheme [74], the five-layer thermal diffusion for land surface processes scheme, the revised MM5 Monin-Obukhov scheme for surface layer physics [75], the Dudhia scheme [76], and the rapid radiative transfer model (RRTM) [77] as shortwave and longwave radiation schemes. No cumulus parameterization was activated. The 3 h 0.25\({}^{\circ}\) operational analysis/forecast cycle, issued at 12 UTC on the day before each day of forecast, from the integrated forecast system (IFS) global model of the European Centre for Medium-Range Weather Forecasts (ECMWF) was used as initial and lateral boundary conditions. To assess the impact of LDA on the precipitation forecast over the sea, we selected the period from 2 to 30 November 2019. This is because this period was characterized by abundant precipitation over Italy, especially in the first 20 days of the month. For all the simulations, we employed a 3 h rapid update cycle (RUC) data assimilation scheme. The RUC configuration was chosen because LDA impacts the weather forecast mainly in the short range [38; 39; 78]. We ran two configurations for each forecast: the first, named CTRL, without LDA, and the second, named LIGHT, with LDA. Each simulation spanned a 9 h period after its initialization time. The first 6 h of each simulation were used for LDA in the LIGHT simulations and to spin-up the model for both configurations, while the last 3 h were used for the forecast. Therefore, eight simulations were needed to cover the forecast for a full day. Lightning data used for assimilation were provided by the LIghnting detection NETwork (LINET) [79]). LINET includes more than 500 stations worldwide, more than 200 of which are in Europe. This network covers the Italian territory and the western Mediterranean Sea. LINET sensors detect very low frequency (VLF) and low-frequency (LF) waves emitted during the flash. LINET is a total lightning system since it is able to measure both intra-cloud (IC) and cloud-to-ground (CG) discharges. The data processing technique follows a 3D method (TOA, Time of Arrival, method, [80]), by which the height of IC lightning is also calculated. Position accuracy is about 75 m for CG lightning. All discharges recorded within a 1 s period and within a 10 km radius are considered as a single flash for data assimilation [81]. For LDA, we followed the method of [PERSON] et al. [29], based on nudging. The modeled water vapor is compared with that calculated by the following equation:  
\[q_{v}=Aq_{s}+Bq_{s}\text{tanh}(CX)\Big{(}1-\text{tanh}\Big{(}Dq_{g}^{\alpha} \Big{)}\Big{)} \tag{1}\]  
In this experiment, the coefficients were set to \(A\) = 0.95, \(B\) = 0.07, \(C\) = 0.25, \(D\) = 0.25, and \(\alpha\) = 2.2. The term \(q_{s}\) represents the saturation water vapor content at model atmospheric temperature, while \(q_{g}\) is the graupel mixing ratio (g kg\({}^{-1}\)). \(X\) is the number of total flashes per grid cell in each assimilation interval (10 min). The water vapor mixing ratio of Equation (1) increases with the flash rate (\(X\)) and decreases with the graupel mixing ratio (\(q_{g}\)). Simulated water vapor is nudged to the water vapor calculated by Equation (1)  
Figure 1: Weather research and forecasting (WRF) model domain. Yellow dots represent the small islands’ positions. when the latter is higher than the modeled one.