Fog Forecast Using WRF Model Output for Solar Energy Applications

Abstract

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 benefit 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 efficiency, 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 influence, 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 specific site under investigation. In this work, the thresholds are site-specifically 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 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 photovoltaic panels surface.

1. Introduction

Fog is defined as consisting of tiny droplets of water or ice crystals with a diameter between ~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 on life and human activities, particularly in meteorology, climatology, transports, energy, agriculture, economy and ecology [3,4,5,6]. The main consequence of fog is visibility reduction that can lead to financial damages, severe accidents and loss of lives [7,8,9]. Thus, fog prediction has always been an extremely important activity for land, sea and air transport sectors [10]. Moreover, the presence of fog tends to alter the radiative budget of the Earth–atmosphere system, hence it is important in meteorology because it can perturb the local temperature and humidity [11]. It impacts air quality because it develops when thermal inversions persist; in this situation, air pollutants are trapped in the boundary layer and not scattered in the atmosphere, especially in urban and industrialized areas [12]. Fog also scatters solar radiation, decreasing the power produced by photovoltaic panels. Fog scattering, in fact, reduces to zero the direct component of incoming solar radiation, which represents the main contribution to the electric power produced by photovoltaic panels [13,14]. Thus, fog forecast is useful to correct solar irradiance prediction and, accordingly, the power produced by solar systems. This certainly represents one of the frontiers of meteorology applied to renewable energy. In this sense, the fog forecast assumes a value that economically impacts the supply of solar energy on the market because the evaluation of the goodness of the forecast is related to an economic index, i.e., the lowest possible imbalance between the expected energy offered on the trading market and the measured final balance. Early attempts in forecasting fog are found in [15], where a method for predicting fog during night based on observations at 8:00 p.m. was presented. Cooling of moist air by radiative flux divergence, vertical mixing of heat and moisture, vegetation, horizontal and vertical wind, heat and moisture transport in soil, advection and topographic effect are all listed in [16] as important processes to consider in the fog formation, whilst longwave radiative cooling at fog top, gravitational droplet settling, fog microphysics and shortwave radiation are important processes driving fog duration.

In general, there are two broad categories of fog forecast: numerical (using computer simulations) and observational (diagnosing the likelihood of a future event based on current observations). An approach to implement fog forecast is to use statistics to define threshold values for key variables involved in the fog formation and to verify its predictability using observations and NWP model outputs. Among the statistical methods, different data mining, machine learning and multivariate techniques have been developed [17,18,19]. In [20], the proposed statistical method assigns a likelihood of fog occurrence from 0 to 1 based on comparing observations of key variables to predefined thresholds.

Another methodology is based on NWP fog simulation in terms of liquid water content or other thermodynamic quantities and processes (horizontal pressure gradient, advection, diffusion, etc.). These models are able to simulate radiation fog when configured with sophisticated options and high horizontal and vertical resolutions and can be implemented for fog forecast in one dimension (1D) similarly to other works [19,21,22,23,24,25,26,27]. The main limitation of these 1D methods is their poor representation of the large-scale situation. Thus, three-dimensional (3D) mesoscale models such as the Fifth-Generation NCAR/Penn State Mesoscale Model (MM5) and the Weather Research and Forecasting (WRF) model are also used [8,19,28,29,30,31]. Fog forecast with NWP models is still a challenge, with limited success in terms of accuracy and precision. This is mainly due to the intrinsically complex nature of the fog phenomenon and to the limited availability of observational and computational resources.

Some empirical methods, instead, are mainly based on thermal and hygrometric gradients between soil and medium-low troposphere [32,33,34,35,36,37,38]. However, these indices need to be adapted on the basis of the historical series and the local conditions, in order to include small scale factors like local advective transport, radiation, turbulent mixing, orography and to produce a multiapproach scheme for fog forecast [34,39].

The implementation of a new fog forecast method represents the main objective of this work. Fog forecast is a rather complex task because of the extreme variability of the conditions that lead to its development, the relatively low knowledge about atmospheric phenomena and the coexistence and interaction of processes characterized by different time and spatial scales. Furthermore, fog is a local forced event that can rapidly grow or dissipate in conditions not uniquely defined. This is the main aspect we focus on in this research: the characterization of site-specific meteorological conditions that promote fog formation. The method to predict fog has been implemented using the output of a numerical weather prediction model. The selected domain is the Italian peninsula. In this territory, the Po Valley, the Apennines and Alpine valleys are often affected by fog during autumn and winter, especially in the morning hours, when a high-pressure stable condition and thermal inversion near the ground are present [2,40]. This kind of fog event has the highest impact on photovoltaic power production [13]. Solar energy is produced by converting the incoming solar radiation reaching the photovoltaic panel surfaces. Solar radiation enters the atmosphere and undergoes physical and photochemical interactions which determine a partial or total extinction. In particular, the atmospheric attenuation is due to the combined effects of backscattering and absorption by clouds and aerosol, the backscattering due to the air molecules and the absorption by the gases according to the Bouguer–Lambert–Beer law [14]. In absorption, a fraction of the energy that propagates through a layer of air is absorbed by the atmospheric constituents and can be re-emitted at a different wavelength. In addition to absorption, a fraction of the radiant energy that passes through the atmosphere is scattered. The scattering is mainly due to the impact with air molecules, dust and liquid water droplets which causes a portion of radiation to be reflected back in all directions (back scattering) or sent to the Earth in a diffuse form (forward scattering) [9]. Among the different types of clouds, we are interested in the specific effects of fog on the solar irradiance and, consequently, on photovoltaic power production. Fog affects solar radiation through Mie scattering and redistributing incident radiation in different directions according to the particle size. This means that in the fog condition, the diffuse component can account almost for all the incoming solar radiation. Hence, fog dramatically impacts the solar power production forecast when clear sky conditions are expected but a fog event occurs [13]. This happens especially in the first hours of daylight until solar heating triggers the fog dissolution process. The WRF model running at IMAA-CNR forecasts the global horizontal irradiance (GHI) and its direct and diffuse components. In particular, given the wide availability of measurements, we decided to focus on the GHI. Our idea is to adjust the WRF GHI output in those cases where the WRF model did not foresee fog while our multitest-based method for fog forecasting (MBFog) did it by activating a fog postprocessing layer. The fog postprocessing layer is designed to apply a fog attenuation transmittance coefficient to the GHI WRF outputs in order to optimize the forecast of the irradiance reaching the Earth’s surface. This transmittance is obtained using a parameterization of the fog extinction coefficient [41] and is a function of the solar elevation angle. In Section 2, the used dataset and the implemented methodology are introduced. Section 3 reports the main results and their discussion, while Section 4 reports the conclusions of the study.

2. 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.

2.1. 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 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.

2.2. 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.

2.3. 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 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).

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.

• Its configuration features the following parameterizations:
• the Thompson aerosol-aware microphysics [46];
• RRTMG (Rapid Radiative Transfer Model for GCMs) radiation scheme that explicitly calculates solar radiative components [44];
• the nonlocal Yonsei scheme for the planetary boundary layer [45];
• the convective scheme Kain–Fritsch [45], only for the low spatial resolution domain;
• the shallow cumulus convection active scheme [45].

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²), 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.

2.4. 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:

$$FSI=2\left(T_s-T_d\right)+2\left(T_s-T_{850}\right)+W{S}_{850}$$

(1)

where:

• $T_s$ is the surface temperature [K],
• $T_d$ is the surface dew point temperature [K],
• $T_{850}$ is the temperature at 850 hPa [K],
• $W{S}_{850}$ is the wind at 850 hPa [m/s].

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:

• test #1—FSI_test: fog stability index under a threshold (FSI_thresh);
• test #2—RH_test: surface relative humidity over a threshold (RH_thresh);
• test #3—Tdepr_test: difference between surface temperature and surface dew point under a threshold (Tdepr_thresh);
• test #4—WS_test: wind speed at 10 m between two thresholds (WS_min_thresh and WS_max_thresh);
• test #5—RHDiff_test: relative humidity difference between the first two vertical levels under a threshold (RHDiff_thresh).

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. List of sites used for WRF outputs performance evaluation and relative number of samples used for validation.

Longitude (°)	Latitude (°)	IATA Code	City	Samples
8.72396 E	45.63 N	LIMC	Milano Malpensa	321
9.2626 E	45.46143 N	LIML	Milano Linate	318
11.29694 E	44.53083 N	LIPE	Bologna	321
11.6125 E	44.81556 N	LIPF	Ferrara	220
10.87228 E	45.38749 N	LIPX	Verona Villafranca	310
12.35194 E	45.50528 N	LIPZ	Venezia	313
12.72778 E	43.51694 N	LIVF	Frontone	86

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.

2.4.1. Evaluation of the IMAA-CNR WRF Variables Outputs

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. Statistical scores of WRF surface temperature outputs related to the seven selected sites measurements. (A negative BIAS means a WRF underestimation and vice versa).

		Surface Temperature
	N	BIAS (°C)	MAE (°C)	nRMSE	R	FB	FAC2
Milano Linate	12	0.21	0.88	0.08	0.99	−0.004	1
Milano Malpensa	15	0.29	0.74	0.08	0.99	−0.01	0.87
Verona	17	−0.15	1.39	0.11	0.98	0.03	0.94
Venezia	17	0.58	0.81	0.12	0.98	−0.02	1
Bologna	6	−0.30	0.7	0.1	0.98	0.01	1
Ferrara	6	−1.05	1.71	0.21	0.97	0.03	1
Frontone	4	1.17	1.99	0.29	0.89	0.07	1

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₂) and surface pressure (psfc) outputs according to the following formula [51]:

$$SurfaceDewPoint=\frac{B}{ln\left(\frac{A\ast eps}{q_2\ast psfc}\right)}$$

(2)

where $A=2.53\times {10}^8 \mathrm{kPa}$, $B=5.43\times {10}^3 \mathrm{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. Statistical scores of WRF surface dew point outputs related to the seven selected sites measurements. (A negative BIAS means a WRF underestimation and vice versa).

	Surface Dew Point
	N	BIAS (°C)	MAE (°C)	nRMSE	R	FB	FAC2
Milano Linate	15	0.98	2.03	0.63	0.91	−0.06	0.66
Milano Malpensa	12	1.49	1.68	1.6	0.93	−0.19	0.4
Verona	13	0.07	1.43	0.54	0.97	−0.01	0.93
Venezia	17	1.12	1.46	0.92	0.88	−0.09	0.78
Bologna	6	−0.07	0.79	0.35	0.96	0.007	0.83
Ferrara	6	0.63	1.62	0.17	0.88	0.17	0.75
Frontone	4	0.80	1.06	1.4	0.91	−0.07	0.72

).

The surface relative humidity is obtained through a formula that requires the surface mixing ratio and the surface pressure as input [36]:

$$R{H}_{surface}=100\frac{e\left(psfc-e_s\right)}{e_s\left(psfc-e\right)}$$

(3)

where:

• $e_s=e_0\ast exp\left(e_1\frac{t-c}{t-e_2}\right)$,
• $e=\frac{q_2\ast psfc}{eps+q_2}$,
• $e_0=6.112 \mathrm{hPa}$,
• $e_1=17.67$,
• $e_2=29.65$,
• $eps=0.622$
• $c=273.15$,
• $q_2$ is the surface mixing ratio
• $psfc$ 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. Statistical scores of WRF surface relative humidity outputs of the seven selected sites (a negative BIAS means a WRF underestimation and vice versa).

	Surface Relative Humidity
	N	BIAS (%)	MAE (%)	nRMSE	R	FB	FAC2
Milano Linate	12	3.87	9.16	0.21	0.86	−0.01	1
Milano Malpensa	15	3.28	4.63	0.09	0.98	−0.01	1
Verona	13	3.11	7.87	0.46	0.91	0.001	1
Venezia	17	3.22	8.03	0.16	0.79	−0.01	1
Bologna	6	0.36	1.14	0.02	0.99	−0.001	1
Ferrara	6	6.5	6.5	0.13	0.92	0.02	1
Frontone	4	3.61	8.55	0.10	0.89	0.01	1

).

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. Statistical scores of WRF wind speed at 10 m outputs of the seven selected sites. (a negative BIAS means a WRF underestimation and vice versa).

	Wind Speed at 10 m
	N	BIAS (m/s)	MAE (m/s)	nRMSE	R	FB	FAC2
Milano Linate	12	−0.13	0.93	1.15	0.26	0.17	0.67
Milano Malpensa	15	−1.12	1.67	0.82	0.45	0.10	0.67
Verona	13	−0.42	1.14	0.62	0.63	0.05	0.69
Venezia	17	0.39	1.16	0.72	0.70	0.18	0.41
Bologna	6	−0.44	0.93	0.48	0.29	0.05	0.66
Ferrara	6	1.73	1.81	1.19	0.38	−0.16	0.67
Frontone	4	−0.19	1.25	1.04	0.64	0.04	0.72

).

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.

2.4.2. 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).

3. 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.

3.1. Milano Linate (45.46143 N, 9.2626 E, 103 m a.s.l.)

Milano Linate is situated in the Po Valley and is characterized by a relative high probability of fog formation. This is confirmed by METAR and SYNOP 2011–2017 observations: during the whole period, fog has been reported in 5311 of 122075 measurements (4.33% of total observations number) distributed unequally in all seasons: 51.16% in autumn, 46.96% in winter, 1.73% in spring and 0.15% in summer. Values and selected thresholds are reported in

Table 6. 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 Milano Linate.

	Mean	Standard Deviation	Threshold
RH	99.30%	3.40%	Min: 98%
Diff_T_TDEW	0.10 °C	0.21 °C	Max: 1.2 °C
WS10	1.43 m/s	0.75 m/s	Min: 0 m/s Max: 1.3 m/s
FSI	-	-	28

, 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.

3.2. Milano Malpensa (45.62 N, 8.7231 E, 211 m a.s.l.)

The METAR 2011–2017 observations at Milano Malpensa reported fog in 2797 out of 121987 cases (2.29% of total observations) occurring in autumn (42.97%), winter (56.67%), spring (0.25%) and summer (0.15%), respectively. The mean and the standard deviation of the surface relative humidity (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. Mean, standard deviation and thresholds selected for relative humidity, surface temperature/dew point depression and wind speed measured in fog condition and fog stability index (FSI) threshold calculated using Youden’s index for Milano Malpensa.

	Mean	Standard Deviation	Threshold
RH	98.6%	4.5%	Min: 99%
Diff_T_TDEW	0.22 °C	0.37 °C	Max: 1.6 °C
WS10	1.07 m/s	0.68 m/s	Min: 0 m/s Max: 1.3 m/s
FSI	-	-	26

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.

3.3. Verona (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. 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 Verona.

	Mean	Standard Deviation	Threshold
RH	95.54%	3.58%	Min: 96%
Diff_T_TDEW	0.66 °C	0.57 °C	Max: 1.9 °C
WS10	1.01 m/s	0.84 m/s	Min: 0 m/s Max: 1.3 m/s
FSI	-	-	26

, MBFog forecast method achieves the following performances: POD = 0.80, FAR = 0.14 and Accuracy = 0.96.

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. 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.

	Mean	Standard Deviation	Threshold
RH	96.15%	9.24%	Min: 91%
Diff_T_TDEW	0.65 °C	0.87 °C	Max: 2.3 °C
WS10	1.44 m/s	0.9 m/s	Min: 1.1 m/s Max: 2.9 m/s
FSI	-	-	27

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.

3.5. Bologna (44.5308 N, 11.2969 E, 42 m a.s.l.)

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. 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 Bologna.

	Mean	Standard Deviation	Threshold
RH	96.86%	2.75%	Min: 95%
Diff_T_TDEW	0.45 °C	0.4 °C	Max: 0.5 °C
WS10	1.82 m/s	1.05 m/s	Min: 0.3 m/s Max: 2.4 m/s
FSI	-	-	24

while in Figure 9 are shown the related histograms of measurements in fog condition. Following the Youden’s 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. 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 Ferrara.

	Mean	Standard Deviation	Threshold
RH	92.98%	4.75%	Min: 91%
Diff_T_TDEW	1.07 °C	0.75 °C	Max: 1 °C
WS10	0.49 m/s	0.71 m/s	Min: 0.7 m/s Max: 2.1 m/s
FSI	-	-	25

.

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. 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 Frontone.

	Mean	Standard Deviation	Threshold
RH	94.9%	3.71%	Min: 95%
Diff_T_TDEW	0.78 °C	0.6 °C	Max: 2.4 °C
WS10	1.78 m/s	1.39 m/s	Min: 0 m/s Max: 2.5 m/s
FSI	-	-	32

, the MBFog method achieves for Ferrara POD = 0.65, FAR = 0.08 and Accuracy = 0.90.

3.7. Frontone (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.

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).

3.8. Discussion of MBFog Evaluation Results

MBFog performances were evaluated with respect to the following methods:

(a)

Fog Stability Index under a fixed threshold (31), as reported in [38] (i.e., FSI test);

(b)

RH larger than a fixed threshold (90%) and Wind Speed smaller than a fixed threshold (2 ms⁻¹), as reported in [34] (i.e., RH/Wind test);

(c)

Visibility obtained from Liquid Water Content (LWC) using Kunkel’s law [41] less than a fixed threshold (1 km), as reported in [37,41] (i.e., Visibility test).

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. Overall statistic performances of multitest fog forecast method for the selected sites.

	Milano Linate	Milano Malpensa	Verona	Venezia	Bologna	Ferrara	Frontone
tp	14	7	12	11	3	22	12
tn	119	109	104	122	127	104	16
fp	2	3	2	2	2	2	1
fn	6	1	3	0	1	12	0
n_total	141	120	121	136	133	140	29
Accuracy	0.94	0.97	0.96	0.99	0.98	0.90	0.97
BIAS	0.80	1.25	0.93	1.17	1.25	0.71	1.08
POD	0.70	0.88	0.80	1	0.75	0.65	1
FAR	0.12	0.30	0.14	0.14	0.40	0.08	0.08
POFD	0.02	0.03	0.02	0.02	0.02	0.02	0.06
SR	0.88	0.70	0.86	0.86	0.60	0.92	0.92
TS	0.64	0.64	0.71	0.86	0.50	0.61	0.92
ETS	0.59	0.61	0.67	0.84	0.49	0.54	0.87
HKD	0.68	0.85	0.78	0.98	0.73	0.63	0.94
HSS	0.75	0.76	0.80	0.91	0.66	0.70	0.93
ORSS	0.99	0.99	0.99	1	0.99	0.98	1

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.

3.9. 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. Comparison of GHI WRF output performances with or without the application of fog postprocessing layer (ppl).

	Milano Linate		Milano Malpensa		Bologna		Ferrara
	No fog ppl	With fog ppl	No fog ppl	With fog ppl	No fog ppl	With fog ppl	No fog ppl	With fog ppl
n_total	23	23	15	15	14	14	41	41
BIAS (W·m−2)	89.44	23.49	80.98	49.42	88.62	26.02	82.73	17.60
MAE (W·m−2)	89.96	31.11	80.98	49.42	88.62	27.12	102.12	44.03
RMSE (W·m−2)	122.68	39.97	118.58	76.4	125.47	36.36	139.28	59.89
R	0.83	0.96	0.80	0.86	0.79	0.98	0.63	0.80
FB	−0.15	−0.05	−0.25	−0.20	−0.14	−0.05	−0.14	−0.04
nMSE	0.45	0.06	1.22	0.63	0.46	0.05	0.70	0.19
FAC2	0.52	0.78	0.13	0.27	0.36	0.93	0.56	0.73

below.

4. 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.