Fog Forecast Using WRF Model Output for Solar Energy Applications

: The occurrence of fog often causes errors in the prediction of the incident solar radiation and the power produced by photovoltaic cells. An accurate fog forecast would beneﬁt solar energy producers and grid operators, who could take coordinated actions to reduce the impact of discontinuity, the main drawback of renewable energy sources. Considering that information on discontinuity is crucial to optimize power production estimation and plant management e ﬃ ciency, in this work, a fog forecast method based on the output of the Weather Research and Forecasting (WRF) numerical model is presented. The areal extension and temporal duration of a fog event are not easy to predict. In fact, there are many physical processes and boundary conditions that cause fog development, such as the synoptic situation, air stability, wind speed, season, aerosol load, orographic inﬂuence, humidity and temperature. These make fog formation a complex and rather localized event. Thus, the results of a fog forecast method based on the output variables of the high spatial resolution WRF model strongly depend on the speciﬁc site under investigation. In this work, the thresholds are site-speciﬁcally designed so that the implemented method can be generalized to other sites after a preliminary meteorological and climatological study. The proposed method is able to predict fog in the 6–30 h interval after the model run start time; it has been evaluated against METeorological Aerodrome Report data relative to seven selected sites, obtaining an average accuracy of 0.96, probability of detection of 0.83, probability of false detection equal to 0.03 and probability of false alarm of 0.18. The output of the proposed fog forecast method can activate (or not) a speciﬁc fog postprocessing layer designed to correct the global horizontal irradiance forecasted by the WRF model in order to optimize the forecast of the irradiance reaching the photovoltaic panels surface. F.R.; resources, S.G., E.R. (Ermann Ripepi) and M.V.; software, S.T.N. and F.D.P.; supervision, F.R.; validation, D.C., D.G. and F.R.; visualization, S.T.N.; writing—original draft, S.T.N.; writing—review and editing, D.C., F.D.P., D.G., S.G., E.G., S.L., E.R. (Elisabetta Ricciardelli), E.R.


Introduction
Fog is defined as consisting of tiny droplets of water or ice crystals with a diameter betweeñ 5-30 µm suspended in the immediate vicinity of the Earth's surface, able to decrease the horizontal visibility down to less than 1 km [1]. Fog is a microscale phenomenon as it is directly influenced by local surface forcing and weather conditions, with a time scale of hours or less [2]. Fog has several implications on the atmosphere and the environment, leading to numerous direct and indirect effects

Materials and Methods
First, we introduce the METeorological Aerodrome Report, the surface SYNOPtic observations and the European Centre for Medium-Range Weather Forecasts reanalysis dataset. Successively, the specific configuration of the WRF NWP model used for the development of the fog forecast methodology is described.

METAR and SYNOP Dataset
METeorological Aerodrome Reports (METAR) data are collected in aeronautical weather stations. Data collection can be regular or special. Regular measurements are done at fixed temporal rate (10,20,30 or 60 min) while special measurements are taken when a significant variation of the weather between two regular reports occurs. For the fog presence verification purpose, we focus on the visibility METAR data. The Italian METAR network features both fully and partially automated stations [42]. This means that some METAR stations are supervised by an observer who manually reports the relevant information when necessary. Although the METAR can contain the specification about the direction in which the visibility measurements have been acquired, in our study, we considered only 360 • visibility data, with no indication on direction. Since vertical visibility is not always available, it was not considered in this study. Daylight hourly METAR reports are used as the source of Energies 2020, 13, 6140 4 of 28 predictand-i.e., the presence of fog. Fog is detected based on two conditions: (i) reported hourly visibility is less than 1 km (according to WMO fog definition) and (ii) the label "FG" alone (indicating fog condition) is reported in the METAR present weather field. In this way, the reduction of visibility has been attributed uniquely to the presence of fog and not to other weather conditions such as intense rain, thunderstorm, sand, or ash.
Moreover, METAR wind speeds recorded during the period 2011-2017 (in case fog condition was reported) were used in this study to define thresholds. METAR visibility data were used twofold: (i) for the climatological featuring of the main meteorological variables involved in the fog forecasting algorithm, and (ii) for the evaluation of the implemented method. To evaluate the proposed method, an independent METAR dataset was collected, covering a different time interval (January to April 2018) and relative to seven sites in Italy selected based on fog occurrences.
SYNOP bulletins are produced based on the same METAR measurements, but at higher signal resolution. SYNOP values of surface and dew point temperature, recorded with one tenth of degree resolution, were collected for the period 2011-2017 (in the case fog condition was reported) and used in this work to define thresholds. The SYNOP bulletins are issued worldwide at least on a six-hour basis; in Italy, the release frequency is three hours. We have therefore performed a temporal interpolation of these measurements to obtain the hourly frequency required in our study. METAR reports and SYNOP bulletins used in this work are property of the Italian Air Force Meteorological Service. Data were provided under the Educational/Research license.

ECMWF ERA5 HRES Reanalysis
The European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 HRES dataset [43] provides several meteorological variables at one-hour time resolution. It is a global climate reanalysis dataset with 0.25 by 0.25 degrees spatial resolution, covering the period 1950 to present. The reanalysis is a numerical description of the atmospheric conditions obtained combining models with a comprehensive set of observations. ERA5 HRES is produced using 4D-Var data assimilation in CY41R2 Version of the ECMWF's Integrated Forecast System (IFS) model, with 137 hybrid sigma/pressure vertical levels, with the top level at 0.01 hPa. Atmospheric data are available for these levels and they are also interpolated to 37 pressure, 16 potential temperature and 1 potential vorticity level(s). Among the hundreds of ERA5 HRES variables we have selected the ones closely related to fog: relative humidity, temperature and dew point at 2 m and wind speed at 10 m.

WRF Numerical Weather Prediction Model
The Solar Version [44] of the WRF numerical weather prediction model is currently operative at the Institute of Methodologies for Environmental Analysis of the National Research Council of Italy (IMAA-CNR). WRF is a limited area model developed by the National Center for Atmospheric Research (NCAR). WRF-Solar Version has been released from the Version 3.6 [45] and is designed to be used in the field of solar energy, thanks to the addition of specific tools. The main features of the WRF model are briefly described below.
The WRF is a nonhydrostatic NWP that solves and integrates the atmospheric dynamics equations; these are formulated using the terrain following vertical coordinate eta [45]. The numerical scheme used to integrate the low frequency modes (meteorological phenomena) is the third order Runge-Kutta method, while the forward/backward one is used for the high-frequency modes (acoustic and gravitational waves). As far as physics is concerned, WRF contains several options for each parameterized category, allowing them to be selected and combined according to the objective of the work. The parametrized physical categories are microphysics, convection, planetary boundary layer (PBL), surface model and radiation. Microphysics explicitly resolves the precipitation, vapor and cloud processes. The convective scheme reproduces vertical flows, and it is used only on low spatial resolution domains (>3-5 km) where convective eddies cannot be explicitly resolved. The PBL schemes Energies 2020, 13, 6140 5 of 28 reproduce the vertical and horizontal diffusion terms and consider the fluxes of latent heat and moisture calculating the vertical gradients. Finally, the radiative scheme provides atmospheric warming by considering the contributions of both shortwave and longwave radiation. In this study, we used WRF v3.8.1, which was released by the NCAR in August 2016, specifically configured for the purposes of this work. The parent domain covers part of the Mediterranean basin with a spatial resolution of 9 km, while the inner nested domain covers the whole Italian peninsula with a 3 km spatial resolution ( Figure 1).
Energies 2020, 13, x FOR PEER REVIEW 5 of 30 we used WRF v3.8.1, which was released by the NCAR in August 2016, specifically configured for the purposes of this work. The parent domain covers part of the Mediterranean basin with a spatial resolution of 9 km, while the inner nested domain covers the whole Italian peninsula with a 3 km spatial resolution (Figure 1). The large domain allows us to capture the synoptic structures entering the Mediterranean from the west side that often drive the meteorological situation on the Italian peninsula. The domains are represented with a Lambertian projection and employ land use and the DEM (digital elevation model) from the 2008 updated MODIS database with spatial resolution of 30 s of arc (about 900 m).
The output of the Global Forecast System (GFS) model with 0.25° horizontal resolution were first used to create the initial and border conditions for the WRF-Solar IMAA-CNR numerical model.
The WRF in the described configuration performs a daily processing and, starting from midnight, provides the weather forecast for the next 120 h with 1 h output writing frequency (tunable to higher frequency at the expense of increased computational time and storage capacity). For example, Figure 2 shows forecasts of the shortwave surface downward global horizontal irradiance (W/m 2 ), the 3 h cumulated total precipitation (mm/3 h) and the wind speed (knots) and direction, the 2 m temperature (°C) and the cloud cover. The large domain allows us to capture the synoptic structures entering the Mediterranean from the west side that often drive the meteorological situation on the Italian peninsula. The domains are represented with a Lambertian projection and employ land use and the DEM (digital elevation model) from the 2008 updated MODIS database with spatial resolution of 30 s of arc (about 900 m).
The output of the Global Forecast System (GFS) model with 0.25 • horizontal resolution were first used to create the initial and border conditions for the WRF-Solar IMAA-CNR numerical model.
The WRF in the described configuration performs a daily processing and, starting from midnight, provides the weather forecast for the next 120 h with 1 h output writing frequency (tunable to higher frequency at the expense of increased computational time and storage capacity). For example, Figure 2 shows forecasts of the shortwave surface downward global horizontal irradiance (W/m 2 ), the 3 h cumulated total precipitation (mm/3 h) and the wind speed (knots) and direction, the 2 m temperature ( • C) and the cloud cover.

The MBFog Multitest Approach
In this study, a multitest-based method for fog forecasting (MBFog) is proposed. It was preferred to the single diagnostic approach because of limits that the latter shows in fog prediction [9,33]. Multitest-based approaches use two or more variables depending on the desired fog characterization. In particular, in [33], three tests were used to deal with different fog types: liquid water content (LWC), cloud base/top and surface relative humidity (RH)-wind rules. Since LWC is related to the horizontal visibility according to Kunkel equation [41], the LWC rule assumes fog condition, i.e., visibility of 1 km or less, when LWC is larger than 0.015 g/kg. Cloud base/top rule uses the nominal vertical features of fog (base of 50 m or less, top of 400 m or less) to identify a fog event. Finally, RH-wind rule identifies fog when two meteorological situations occur: a 2 m relative humidity of 90% or more and a 10 m wind speed of 2 m/s or less. In [39], the authors demonstrated that cloud base/top rule is good for a large-scale fog event like marine fog or coastal fog; for this reason, the RH-wind rule and a two-level approach using the temperature gradient rule were considered, in order to forecast shallow ground fog. The proposed method is a diagnostic forecast methodology based on the combination of different threshold tests applied to WRF model variables outputs, as shown in Figure 3. The main variables involved in the fog process, such as temperature, humidity, wind and dew point, are considered for threshold tests. The MBFog method is a combination of the tests implemented in [34,39] with the fog stability index [38,47] test. The fog stability index (FSI) is obtained according to the following formula:

The MBFog Multitest Approach
In this study, a multitest-based method for fog forecasting (MBFog) is proposed. It was preferred to the single diagnostic approach because of limits that the latter shows in fog prediction [9,33]. Multitest-based approaches use two or more variables depending on the desired fog characterization. In particular, in [33], three tests were used to deal with different fog types: liquid water content (LWC), cloud base/top and surface relative humidity (RH)-wind rules. Since LWC is related to the horizontal visibility according to Kunkel equation [41], the LWC rule assumes fog condition, i.e., visibility of 1 km or less, when LWC is larger than 0.015 g/kg. Cloud base/top rule uses the nominal vertical features of fog (base of 50 m or less, top of 400 m or less) to identify a fog event. Finally, RH-wind rule identifies fog when two meteorological situations occur: a 2 m relative humidity of 90% or more and a 10 m wind speed of 2 m/s or less. In [39], the authors demonstrated that cloud base/top rule is good for a large-scale fog event like marine fog or coastal fog; for this reason, the RH-wind rule and a two-level approach using the temperature gradient rule were considered, in order to forecast shallow ground fog. The proposed method is a diagnostic forecast methodology based on the combination of different threshold tests applied to WRF model variables outputs, as shown in Figure 3. The main variables involved in the fog process, such as temperature, humidity, wind and dew point, are considered for Energies 2020, 13, 6140 7 of 28 threshold tests. The MBFog method is a combination of the tests implemented in [34,39] with the fog stability index [38,47] test. The fog stability index (FSI) is obtained according to the following formula: (1) where: The first term, − , is the dew point depression and provides information regarding the availability of moisture content in the surface proximity. The second term, − , is the stability term and gives an estimation of the atmosphere steadiness. The FSI index can assume continue values between 0 and 100, an appropriate threshold identifies FSI values related to fog condition.
To summarize, the MBFog approach presented here considers the following criteria to derive the fog forecast: The last rule tends to identify the cases featuring negative relative humidity gradient between two layers next to the surface as "no fog". The main meteorological variables involved in fog onset Thresholds are fixed after site-specific climatological study (test #2, #3, #4), statistical study (test #1) or literature review (test #5).
The first term, T s − T d , is the dew point depression and provides information regarding the availability of moisture content in the surface proximity. The second term, T s − T 850 , is the stability term and gives an estimation of the atmosphere steadiness. The FSI index can assume continue values between 0 and 100, an appropriate threshold identifies FSI values related to fog condition.
To summarize, the MBFog approach presented here considers the following criteria to derive the fog forecast: The last rule tends to identify the cases featuring negative relative humidity gradient between two layers next to the surface as "no fog". The main meteorological variables involved in fog onset are temperature, relative humidity, dew point and wind speed. These quantities are derived from the output of a WRF operational chain implemented at IMAA-CNR to provide forecasts of all solar irradiance variables at high temporal and horizontal resolution for the benefit of solar energy industry [48]. The evaluation has been carried out against measurements from seven METAR sites in the Italian peninsula. The decision to focus on specific points of the area of interest followed by the localized nature of fog and its low probability of occurrence. This does not suggest general conditions to be valid into the whole domain. However, the method can be exported to other locations for which meteorological historical data are available. The METAR sites are listed in Table 1 and their location is shown in Figure 4. are temperature, relative humidity, dew point and wind speed. These quantities are derived from the output of a WRF operational chain implemented at IMAA-CNR to provide forecasts of all solar irradiance variables at high temporal and horizontal resolution for the benefit of solar energy industry [48]. The evaluation has been carried out against measurements from seven METAR sites in the Italian peninsula. The decision to focus on specific points of the area of interest followed by the localized nature of fog and its low probability of occurrence. This does not suggest general conditions to be valid into the whole domain. However, the method can be exported to other locations for which meteorological historical data are available. The METAR sites are listed in Table 1 and their location is shown in Figure 4.  The sites have been chosen for their relative high frequency of fog occurrence during the selected period. In particular, data collected between January and April 2018 relative to the first 30 h of the daily WRF run have been used for the evaluation: since the initial 6 h spin up time is discarded [49] and the WRF run starts at midnight each day, the first 30 h correspond to the forecast until the following day at noon. The number of samples for each METAR site are reported in Table 1 and refers to 35 WRF model runs for a total of 325 hourly data for the whole Italian peninsula. We evaluate the WRF outputs when clear sky conditions are identified, so to eliminate perturbations that could degrade the forecast. Considering this, a total of 73 samples are used for WRF output evaluation. Forecasts refer to different times during daylight. The sites have been chosen for their relative high frequency of fog occurrence during the selected period. In particular, data collected between January and April 2018 relative to the first 30 h of the daily WRF run have been used for the evaluation: since the initial 6 h spin up time is discarded [49] and the WRF run starts at midnight each day, the first 30 h correspond to the forecast until the following day at noon. The number of samples for each METAR site are reported in Table 1 and refers to 35 WRF model runs for a total of 325 hourly data for the whole Italian peninsula. We evaluate the WRF outputs when clear sky conditions are identified, so to eliminate perturbations that could degrade the forecast. Considering this, a total of 73 samples are used for WRF output evaluation. Forecasts refer to different times during daylight. The statistical analysis to evaluate WRF outputs of meteorological variables is based on the calculation of bias (BIAS), mean absolute error (MAE), normalized root-mean-square error (nRMSE), correlation coefficient (R), fractional bias (FB) and fraction of predictions within a factor of two of observations (FAC2). The statistics description and related formulas can be found in [50] and is reported in Appendix A. WRF outputs have been evaluated against METAR, SYNOP and ECMWF ERA5 datasets.
The surface temperature was investigated and the correlation between observed and predicted temperature ranges from 0.89 to 0.99, the nRMSE from 0.08 to 0.29 and the bias from −1.05 to 1.17 • C was analyzed (see Table 2 for the complete set of scores). The forecasted surface dew point temperature is not obtained directly as a WRF output, but it is calculated using surface mixing ratio (q 2 ) and surface pressure (psfc) outputs according to the following formula [51]: where A = 2.53 × 10 8 kPa, B = 5.43 × 10 3 K, eps = 0.622, q 2 is the surface mixing ratio expressed in kg/kg and psfc is the surface pressure expressed in hPa. Correlation coefficient varies in 0.88-0.97 range, the nRMSE between 0.17 and 1.6 and the bias from −0.07 to 1.49 • C (Table 3). The surface relative humidity is obtained through a formula that requires the surface mixing ratio and the surface pressure as input [36]: Energies 2020, 13, 6140 10 of 28 where: ps f c is the surface pressure in hPa.
Concerning the surface relative humidity, the forecasts are substantially in agreement with the measurements: Ferrara and Frontone sites are slightly underestimated whilst the remaining sites are overestimated (Table 4). The WRF model wind output is provided as Eastward (U) and Northward (V) components in the Arakawa C staggered grid [44] in which the U component refers to the center of the left grid face and the V to the lower grid faces. These components are relative to the model grid and not to the Earth coordinates. METAR data report wind speed and direction; hence, before comparison, U and V wind vectors have been converted in wind speed. With regards to wind speed, performances are worse than other variables, however reasonable if accounting for a systematic estimation error for all the selected sites (Table 5). Overall, the IMAA-CNR WRF model shows good performances for the selected sites; differences between observed and forecast values are accounted when tuning the thresholds of the MBFog method.

Definition of the Thresholds
Site-specific climatological and statistical analysis were carried out for each of the seven selected sites in order to derive appropriate thresholds. In particular, we performed a climatological characterization of surface temperature, surface relative humidity, 10 m wind speed and surface dew point during reported fog events, based on the historical series of METAR and SYNOP bulletins in the period January 2011-December 2017. Threshold values are computed via a statistical method aimed at maximizing the accuracy on the training dataset (climatological METAR and SYNOP dataset 2011-2017). Starting from this, we derived a corresponding analytical model given by: Threshold = MEAN ± STD + sign (BIAS)*MAE/2. This expression is only intended to provide a posteriori analytical reference for thresholds, explaining their order of magnitude in terms of three contributions related to the climatological dataset (i.e., MEAN and STD) and the NWP method (i.e., MAE) used. Thus, in this expression, RH_thresh and WS_min_thresh are set to the average minus the standard deviation values while TDepr_thresh and WS_max_thresh to the average plus the standard deviation values. These thresholds should be adapted considering the BIAS and the MAE of the WRF outputs presented above (see Section 2.4.1). All of them should be increased or reduced (depending on the BIAS sign) by a quantity equal to the half of the MAE. Regarding the relative humidity difference between the first two vertical levels (RHDiff) and the fog stability index, historical series are not available. Thus, in the case of RHDiff, we selected a threshold available and validated in literature, equal to −4.5% [39], while for the FSI we selected a subset of WRF outputs to derive an appropriate threshold. To this aim, we use the maximization criteria based on the area under ROC (receiving operating characteristic) curve and Youden's index [52]. Youden's index expresses the performances of a dichotomous diagnostic test (see Appendix B). We calculated this index using fixed values for relative humidity, surface temperature, surface dew point and wind speed thresholds but varying the FSI threshold value between its minimum and maximum value (0 and 100, assumed to be the highest and lowest probability of fog occurrence, respectively). The appropriate FSI_thresh was chosen as the FSI value corresponding to the maximum Youden's index. This permitted us to obtain the maximum area under ROC curve (AROC).

Results and Discussion
In this section, we report results of the site by site analysis. This is conducted to customize the fog forecasting method to the specific environment, in terms of values assumed by the variables of interest in the presence of fog and the relative derived thresholds.  Table 6, while Figure 5 shows the related histograms of measurements. It should be noted that METAR ambient and dew point temperatures are rounded towards the nearest integer, causing unrealistic gaps when T is close to Td and consequently in fog cases. This limitation has an impact also in relative humidity values that are obtained from the ambient and dew point temperatures [36]. Therefore, we used SYNOP data for the ambient temperature, dew point temperature and relative humidity climatological analysis in all the evaluated sites.
According to Youden's index, FSI_thresh equals 28, obtaining the ROC curve in Figure 5d, with AROC = 0.72. By applying these thresholds to the IMAA-CNR WRF variables outputs, the MBFog method has been evaluated against METAR, achieving the following performances: POD = 0.70, FAR = 0.12 and Accuracy = 0.94.   According to Youden's index, FSI_thresh equals 28, obtaining the ROC curve in Figure 5d, with AROC = 0.72. By applying these thresholds to the IMAA-CNR WRF variables outputs, the MBFog method has been evaluated against METAR, achieving the following performances: POD = 0.70, FAR = 0.12 and Accuracy = 0.94.   Table 7 while Figure 6 shows the related histogram of measurements. (RH), the difference between surface temperature and surface dew point (Diff_T_TDEW) and the wind speed at 10 m (WS10) have also been calculated. Values and adapted thresholds are reported in Table 7 while Figure 6 shows the related histogram of measurements.   According to Youden's index, FSI_thresh equals 26, obtaining the ROC curve in Figure 6d, with AROC = 0.71. By applying these thresholds to the IMAA-CNR WRF variables outputs, the MBFog methodology has been evaluated against METAR, achieving the following performances: POD = 0.88, FAR = 0.30 and Accuracy = 0.97. Verona is situated in the Po Valley next to Lake Garda. Its location is favorable for fog, and the near water basin represents a further source. Indeed, METAR during years 2011-2017 reports fog in 5.25% of measurements (6447 out of 118790). Among these, 50.58% are in winter, 46.01% in autumn, Energies 2020, 13, 6140 14 of 28 3.12% in spring and 0.29% in summer. Histograms and values for this site are reported in Figure 7 and Table 4. (45. 3875 N, 10.8723 E, 68 m a.s.l.) Verona is situated in the Po Valley next to Lake Garda. Its location is favorable for fog, and the near water basin represents a further source. Indeed, METAR during years 2011-2017 reports fog in 5.25% of measurements (6447 out of 118790). Among these, 50.58% are in winter, 46.01% in autumn, 3.12% in spring and 0.29% in summer. Histograms and values for this site are reported in Figure 7 and Table 4. In Figure 7, it is shown that the ROC curve obtained with an FSI threshold of 26 underlies an AROC of 0.72. With the thresholds reported in Table 8, MBFog forecast method achieves the following performances: POD = 0.80, FAR = 0.14 and Accuracy = 0.96. In Figure 7, it is shown that the ROC curve obtained with an FSI threshold of 26 underlies an AROC of 0.72. With the thresholds reported in Table 8, MBFog forecast method achieves the following performances: POD = 0.80, FAR = 0.14 and Accuracy = 0.96. The METAR 2011-2017 reports indicate the presence of fog in 4477 out of 120,877 cases (3.65% of the total observations). Of these, the 50.70% are in winter, 43.36% in autumn, 4.89% in spring and 2.14% in summer. Table 9 and Figure 8 report values and histograms for Venezia. The ROC curve is obtained using the FSI threshold equal to 27 and underlies an area of 0.71 as shown in Figure 8d. With the selected thresholds, the evaluation of the MBFog method in Venezia obtains the following performances: POD = 1, FAR = 0.14 and Accuracy = 0.99. Table 9. Mean, standard deviation and thresholds selected for relative humidity, surface temperature/dew point depression and wind speed measured in fog condition and FSI threshold calculated using Youden's index for Venezia. 3.4. Venezia (45.5053 N, 12.3519 E, 68 m a.s.l.) The METAR 2011-2017 reports indicate the presence of fog in 4477 out of 120,877 cases (3.65% of the total observations). Of these, the 50.70% are in winter, 43.36% in autumn, 4.89% in spring and 2.14% in summer. Table 9 and Figure 8 report values and histograms for Venezia. The ROC curve is obtained using the FSI threshold equal to 27 and underlies an area of 0.71 as shown in Figure 8d. With the selected thresholds, the evaluation of the MBFog method in Venezia obtains the following performances: POD = 1, FAR = 0.14 and Accuracy = 0.99. Table 9. Mean, standard deviation and thresholds selected for relative humidity, surface temperature/dew point depression and wind speed measured in fog condition and FSI threshold calculated using Youden's index for Venezia.

Bologna
Bologna is located in the Po Valley and, in the period between the 2011 and 2017, in 4366 cases out of 120,598 METAR reported fog (3.56%). Of these cases, the 53.39% were in winter, the 43.68% in autumn and the 2.93% in spring. Bologna's selected thresholds are reported in Table 10      index criteria, it has been selected the FSI threshold that permits to obtain the ROC curve of Figure  9d. This ROC subtends an AROC = 0.72. Calculated performances are: POD = 0.75, FAR = 0.40 and Accuracy = 0.98.  3.6. Ferrara (44.8156 N, 11.6125 E, 10 m a.s.l.) Ferrara, situated in the middle of Po Valley, is characterized by high frequency of fog occurrence, as reported in 2565 out of 32866 cases, i.e., 7.84% (the 41.92% during winter, the 53.7% in autumn, the 3.64% during spring and the 0.74% during summer) in the METAR 2011-2017 observations. Histograms and selected values for this site are reported in Figure 10 and Table 11. 3.6. Ferrara (44.8156 N, 11.6125 E, 10 m a.s.l.) Ferrara, situated in the middle of Po Valley, is characterized by high frequency of fog occurrence, as reported in 2565 out of 32866 cases, i.e., 7.84% (the 41.92% during winter, the 53.7% in autumn, the 3.64% during spring and the 0.74% during summer) in the METAR 2011-2017 observations. Histograms and selected values for this site are reported in Figure 10 and Table 11.
In Figure 10d, the ROC curve is shown, obtained with FSI threshold equal to 25 with an AROC of 0.74. Applying the thresholds reported in Table 12, the MBFog method achieves for Ferrara POD = 0.65, FAR = 0.08 and Accuracy = 0.90. (44.8156 N, 11.6125 E, 10 m a.s.l.) Ferrara, situated in the middle of Po Valley, is characterized by high frequency of fog occurrence, as reported in 2565 out of 32866 cases, i.e., 7.84% (the 41.92% during winter, the 53.7% in autumn, the 3.64% during spring and the 0.74% during summer) in the METAR 2011-2017 observations. Histograms and selected values for this site are reported in Figure 10 and Table 11.     Values and selected thresholds for this site are reported in Table 12, while in Figure 11 there are shown the related histogram of measurements in fog condition. With these thresholds, we obtained POD = 1, FAR = 0.08 and Accuracy = 0.97. = 0.65, FAR = 0.08 and Accuracy = 0.90. (43. 5169 N, 12.7277 E, 574 m a.s.l.) Frontone is located in the Apennines, in a small valley favorable to the development of radiation fog mainly during autumn and winter seasons. METAR 2011-2017 reported it in 2906 out of 60,125 cases (the 4.71%) distributed during the different seasons in this way: 39.82% in winter, 41.38% in autumn, 14.96% in spring and 3.84% in summer.

Frontone
Values and selected thresholds for this site are reported in Table 12, while in Figure 11 there are shown the related histogram of measurements in fog condition. With these thresholds, we obtained POD = 1, FAR = 0.08 and Accuracy = 0.97.  For the Frontone site, the ROC curve was obtained using an FSI threshold equal to 32 and underlies an area of 0.81 (see Figure 11d).
Note that all variables used in these comparisons (surface temperature, surface dew point temperature, wind speed at 10 m, temperature at 850 hPa, wind speed at 850 hPa and surface mixing ratio) are obtained from the same IMAA-CNR WRF runs used in MBFog performance evaluation, i.e., daytime hourly outputs of 35 runs selected in the period January-April 2018. In this framework, (see the Appendix C for all statistical scores), MBFog results in better performances with respect to the other methods mentioned above. This is mainly attributed to the fact that the thresholds are specifically tuned for each site. The evaluation has been carried out against METAR data on seven sites chosen for their frequent fog occurrence in the period January-April 2018. Table 13 contains the overall statistical scores for all the sites. Appendix B explains the definitions of the considered statistical indexes (see [50] for more details). In the evaluated sites, MBFog forecast method is able to predict fog in the 6-30 h after the WRF start time with an average accuracy of 0.96, an average probability of detection of 0.83 and an average false alarm ratio of 0.18. The average success ratio is 0.82; this means that the 82% of the fog forecasted have been actually observed. The goodness of the MBFog method has been evaluated with respect to the random chance by means of the odds ratio skill score obtaining average value of 0.994 (1 being the perfect score). MBFog suffers from a relatively high false alarm ratio (average value of 0.18); this means that the algorithm can produce a false underestimation of the irradiance reaching the photovoltaic panels when fog is expected while it does not occur.
Considering similar fog forecasting methods, in [40], an accuracy of 0.95 has been reported. In [34], the implemented multivariable-based diagnostic method scored an ETS = 0.334, while the single-rule methodology (LWC rule) scored an ETS = 0.063. Comparison against these methods [34,39] reveals that MBFog is aligned in terms of performances (average accuracy 0.96, average ETS = 0.409); however, it also reveals that this class of fog forecast methods experiences similar limitations, e.g., they predict fog formation fairly well, but they do not predict its duration as well.

First Results of the Application of the Fog Postprocessing Layer to the IMAA-CNR WRF Irradiance Forecast
In this subsection, we present a first evaluation of the application of the fog postprocessing layer to the WRF global horizontal irradiance forecast.
First of all, we selected those case in which IMAA-CNR WRF did not predicted fog while it occurs, after that we calculated statistical scores for 4 METAR site we have available the measured GHI values in the period between January and April 2018: Milano Linate, Milano Malpensa, Bologna e Ferrara. WRF GHI scores have been calculated in both the cases in which fog postprocessing layer is applied or not. In all the considered cases, we had an improvement of the statistical scores when using fog postprocessing layer, as reported in Table 14 below.

Conclusions
This work presents a multitest method to forecast the presence of fog using the output data of a numerical weather prediction model, namely WRF. The fog forecast method is based on several threshold tests applied to the main meteorological variables involved in fog processes obtained from the WRF NWP model at spatial resolution of 3 km and temporal resolution of 1 h over the Italian peninsula. The implemented method combines several deterministic tests in which meteorological variables involved in the fog process and the fog stability index (FSI) are evaluated with respect to empirical site-specific thresholds. The evaluation has been carried out against METAR data on seven sites chosen for their frequent fog occurrence in the period January-April 2018.
The MBFog method tends to maximize the accuracy and the probability of detection and at the same time to minimize the probability of false detection and the probability of false alarm with respect to single-test methods for the evaluated dataset. The strength of this method is its capability to adapt the thresholds to the specific site under investigation. This means that, if meteorological time series are available, the method can be adapted for other sites after a preliminary meteorological and climatological analysis. The output of the proposed fog forecast method can activate or not a specific fog postprocessing layer designed to correct the global horizontal irradiance forecasted by the WRF model in order to optimize the forecast of the irradiance reaching the PV panels surface. We can conclude that the proposed MBFog method is useful to forecast fog occurrence, and thus it can be used to improve the forecast of PV power production.

Acknowledgments:
We acknowledge the Italian Air Force Meteorological Service for the kind cooperation in providing us METAR reports and SYNOP bulletins used in this work. We are also grateful to Claudia Faccani and Gianluca Tisselli from the Italian air navigation service provider (ENAV) for their valuable support in the definition of the measurement dataset.

Conflicts of Interest:
The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript and in the decision to publish the results.

Appendix A. Standard Verification Methods for Continuous Variables Forecasts
In this appendix we describe the statistical indexes used for the validation of the WRF outputs forecasts: bias, mean absolute error (MAE), normalized root-mean-square error (nRMSE), correlation coefficient (R), fractional bias (FB), normalized mean square error (nMSE), fraction of predictions within a factor of two of observations (FAC2).

• Bias
Bias is defined as the sum of the difference between forecasted (F) and observed (O) values divided by the total number of samples.
Bias gives an indication of the forecast average error but does not measure the correspondence between forecasts and observations. Bias can assume values between −∞ and +∞, perfect score means a bias equal to zero.

Mean Absolute Error
Mean Absolute Error (MAE) is the ratio between the sum of the absolute value of the difference between forecasts (F) and observations (O) and the total number of samples.
It is used to address the average magnitude of the forecast errors but does not indicate the direction of them. It can be a value between zero and +∞ with perfect score zero.

•
Normalized Root-Mean-Square Error Normalized Root-Mean-Square Error (nRMSE) is calculated according to the following formula: It measures average error weighted according to the square of the error and is normalized by the mean observation value. In particular this index is influenced mainly by large errors encouraging conservative forecasts. It can assume values between 0 and 1 with perfect score 0.

• Correlation coefficient
Correlation coefficient is calculated using the following equation: It is a measure of the phase error. A good correlation coefficient means a scatter plot with values arranged around the diagonal but does not consider the bias and is sensitive to outliers. It can be a value ranging from −1 to 1 with perfect score equivalent to 1. Having a good correlation coefficient is a necessary but not sufficient condition for having a perfect forecast. •

Fractional Bias
Fractional Bias (FB) is defined as the sum of the difference between observed (O) and forecasted (F) values divided by the sum of the sum between observed (O) and forecasted (F), all divided by 2.
It is a measure of systematic errors and gives an indication on how the forecast underestimate or overestimate the measures. A good forecast means having a FB of 0. •

Normalized Mean Square Error
Normalized Mean Square Error (nMSE) is calculated according to the following formula: It quantifies random error beyond that systematic error. Perfect score is reached with value equal to zero. •

Fraction of predictions within a factor of two of observations
Fraction of predictions within a factor of two of observations (FAC2) is a robust measure because is not influenced by outliers and is obtained following this criterion:

Appendix B. Standard Verification Methods for Dichotomous Variables Forecasts
In this appendix are described the statistical indexes used for the validation of the multitest based approach fog forecast method: accuracy, bias score, probability of detection (POD), false alarm ration (FAR), probability of false detection (POFD), success ratio (SR), threat score (TS), equitable threat score (ETS), Hanssen and Kuipers discriminant (HKD), Heidke skill score (HSS) and Odds ratio skill score (ORSS).
In order to calculate categorical statistics, first, a contingency table must be defined showing the joint distribution, i.e., the frequency of positive and negative occurrences: where: • tp (true positive) or hits are those events that were forecasted and actually occurred; • tn (true negative) or correct negatives are the events that were not forecasted and did not occur; • fp (false positive) or false alarms are those events that were forecasted but not occur; • fn (false negative) or misses are the events occurred but not forecasted; and • n_total is the sum of true positives, true negatives, false positives and false negatives.
Keeping in mind that a perfect system would produce only true positives and true negatives and that in this work a positive event is represented by the correct detection/forecast of fog presence, below there are reported the description of the calculated categorical statistics and the criteria to their interpretation.

• Accuracy
Accuracy, or fraction correct, is the ratio between the sum of hits and correct negatives and the n_total. It can assume values ranging from 0 to 1 (best value). Accuracy gives the fraction of the overall correct forecasts. Accuracy = tp + tn n total (A8) • Bias score Bias score, also called frequency bias, measures the ratio of the frequency of forecasts to the frequency of observations. Can assume values over zero (perfect score is 1). If BIAS < 1 the system has a tendency of under forecast, the opposite if BIAS > 1.

Probability of detection (POD)
Probability of detection (or true positive rate or sensitivity) returns the frequency of observed event correctly forecasted. Values ranges from 0 to 1 with 1 perfect score. POD is sensitive to hits but not consider false alarms.
• True negative rate (Specificity) True negative rate (or specificity) is the ratio of the true negative by the sum of all negative observed events. It is useful in the determination of ROC (Receiving Operating Characteristic) curve and to calculate the Youden's index.
speci f icity = tn tn + f p (A11) • False alarm ratio (FAR) False alarm ratio gives an indication on the fraction of positive predicted events actually not occurred. It can assume values between 0 to 1 with 0 perfect score. It is sensitive to false alarms but not to misses.

FAR =
f p tp + f p (A12) • Probability of false detection (POFD) Probability of false detection (or false alarm rate) is the fraction correct negatives event not forecasted. Values ranges from 0 to 1. POFD should be 0 in a perfect system.

Success ratio (SR)
Success ratio measures the fraction of hits among all the positive forecasts (sum of hits and false alarms). SR is equal to 1 − FAR. Value range is between 0 and 1, with 1 best performance indicator.

SR =
tp tp + f p (A14) • Threat score (TS) Threat score (or critical success index-CSI-) is the measurement of the fraction of observed and forecasted events that were correctly predicted. Thus, it differs from the accuracy because it does not consider true negatives. It is sensitive to hits, A TS of 0 indicates that the system has no skill while 1 that system is perfect.
It is important to remark that all variables used in these comparisons (surface temperature, surface dew point temperature, wind speed at 10 m, temperature at 850 hPa, wind speed at 850 hPa, surface mixing ratio) are obtained from the same IMAA-CNR WRF runs used in MBFog performance evaluation, i.e., daytime hourly outputs of 35 runs selected in the period January-April 2018.
In the following tables all the statistical performances for each of the selected sites are reported. Overall, based on the IMAA-CNR WRF runs outputs, we can conclude that MBFog improves the performance compared to the other considered methods.