Assessment of Simulated Solar Irradiance on Days of High Intermittency Using WRF-Solar

: Improvements in the short-term predictability of irradiance in numerical weather prediction models can assist grid operators in managing intermittent solar-generated electricity. In this study, the performance of the Weather Research and Forecasting (WRF) model when simulating di ﬀ erent components of solar irradiance was tested under days of high intermittency at Mildura, a site located on the border of New South Wales and Victoria, Australia. Initially, four intermittent and clear case days were chosen, later extending to a full year study in 2005. A speciﬁc conﬁguration and augmentation of the WRF model (version 3.6.1) designed for solar energy applications (WRF-Solar) with an optimum physics ensemble derived from literature over Australia was used to simulate solar irradiance with four nested domains nudged to ERA-Interim boundary conditions at grid resolutions (45, 15, 5, and 1.7 km) centred over Mildura. The Bureau of Meteorology (BOM) station dataset available at minute timescales and hourly derived satellite irradiance products were used to validate the simulated products. Results showed that on days of high intermittency, simulated solar irradiance at ﬁner resolution was a ﬀ ected by errors in simulated humidity and winds (speed and direction) a ﬀ ecting clouds and circulation, but the latter improves at coarser resolutions; this is most likely from reduced displacement errors in clouds.


Introduction
Solar-generated electricity is expected to proliferate in the future due to declining costs of photovoltaics and its increased penetration into the grid [1,2]. The rapid expansion of the solar energy market supports ongoing efforts to reduce global emissions of carbon dioxide from fossil generated electricity reducing the profound threat caused by long-term changes in climate [3]. To meet the demand of solar-generated electricity, the supply chain in the solar industry must scale up accordingly. One of the major problems hindering the reliable supply of solar energy is intermittency. The temporal variability of the solar resource due to changes in cloud cover at several timescales immensely affects solar power production, which can be unpredictable due to the weather-climate interactions. Intermittency causes stress on utility operators managing load balancing, grid scheduling, and transmission operations. For example, opaque clouds over photovoltaic arrays can reduce solar-generated electricity by as much as 50%-80%, causing short-term voltage fluctuations. Voltage fluctuations increase maintenance costs by managing worn down line equipment on distribution feeders [4]. These fluctuations may require a complementary power source [5] for stable power production, which adds to the cost of operations. Thus, solar irradiance predictions are needed for effective grid management.
The structure and composition of clouds are stochastic in short-time (>10 min to 5 h) scales [6], imposing strong constraints on the modelling [7] and validation [8] of solar irradiance. Short-term Energies 2020, 13 irradiance prediction is conducted through advection of clouds from ground-based sky cameras or satellites using cloud motion vectors. This is vital for the detection of ramps in voltages and load following [9]. On the other hand, longer-term (>5 to 48 h) irradiance prediction is critical in operating and maintaining system reserves and requires numerical weather prediction models that account for clouds formation and dissolution, not just motion [9]. The financial value of improved irradiance prediction has been recently investigated in Australia for a concentrating solar thermal plant operating in the Australian National Electricity Market. For one year, the financial value increased by $1.04-$1.13 million and the forced outage rate reduced by 20%-21% for a 50 MW plant with 7.5 h of storage [10]. Taken together, irradiance prediction at all time scales is crucial for grid operators in the energy market. However, the Australian National Electricity Market (NEM) requires power generators to submit dispatch offers up to 40 h ahead every day [11]; thus, improved prediction of solar irradiance from weather models is essential for the Australian energy sector.
To predict clouds attenuating sunlight requires a sophisticated computer model with the accurate representation of the physics of clouds. The interplay between clouds, circulation, and climate sensitivity is key to understanding and anticipating changes in climate [12]. Current weather and climate models misrepresent low-level boundary layer clouds [13] and convection in the tropics [14], having the most significant impact on solar irradiance. Most problems in representing clouds in mesoscale models are due to shortcomings and deficiencies in representing convection initiation, evolution, and propagation, and the interaction between convection and its large-scale atmospheric environment and the underlying land surface [15,16]. These shortcomings in mesoscale models affect the timing and locations of clouds and radiation incidents on the ground. Furthermore, the cloud-radiation feedback due to changing microphysical properties of cloud particles further complicates the evolution of clouds and their interaction with radiation [17,18]. Similarly, the atmospheric state during the start of the simulation and the development of turbulent fluxes in the boundary layer also have enormous impacts on cloud formation and circulation in mesoscale models [19][20][21][22].
In Australia, prediction of irradiance using numerical weather prediction models has been tested in several studies [23][24][25][26][27] leading to similar conclusions: there exists a misrepresentation of clouds in the model either due to resolution or physics of clouds in the model. This issue has also been noticed in other parts of the globe [28,29], but recently, several physical and dynamical changes in representing clouds have been invented that improved solar irradiance forecasts using mesoscale models. These include calculations of direct and diffuse components from shortwave radiative transfer schemes [30]; increases in radiation calls and inclusion of fast radiative transfer calculations [31]; other particulates' direct effect in the shortwave spectrum [32]; cloud particles' interaction in microphysical schemes [33]; and sub-grid scale feedback produced by clouds in shortwave irradiance through shallow cumulus parametrization [34]. All these changes have been incorporated in a specific augmentation of the advanced research version of the WRF Model [35] designed for solar energy predictions known as WRF-Solar [36]. WRF-Solar has been extensively tested in the USA [36][37][38][39][40][41] and other countries, such as Spain [42], Singapore [43], Kuwait [44], and Saudi Arabia [45]. Most studies have reported significant improvements in solar irradiance predictions under different sky conditions with WRF-Solar in comparison to standard WRF simulations.
Specifically, our current understanding of the representation and movement of low-level clouds in models is limited, leading to errors in solar irradiance forecasts affecting solar power generation. However, recently several new physical and dynamical developments in clouds enhancing irradiance forecasts have been implemented in WRF-Solar that are yet to be explored in the Australian region, especially under days of high intermittency. This study aims to explore the potential of WRF-Solar for simulating irradiance on days of high intermittency at different horizontal resolutions. The paper is organized as follows: The data and model used in this study are described in Section 2. Results are shown in Section 3. Discussion and Conclusion are Sections 4 and 5, respectively.

Surface Observations
The ground based solar irradiance data measured at the Australian Bureau of Meteorology (BOM) station at Mildura site was used for model evaluation. Mildura is located on the border of NSW and Victoria (located at a latitude of −34.236 • and longitude of 142.087 • ) and has been identified as a potential solar farm location. The ground station takes measurements of global horizontal irradiance (GHI), direct normal irradiance (DNI), and diffuse horizontal irradiance (DHI) using a shaded pyranometer and shaded pyrgeometer. Apart from solar irradiance data, surface observations obtained from Automatic Weather Station at Mildura measuring relative humidity, sea-level pressure, temperature, windspeed, and wind-direction were also used for model evaluation. The quality-controlled minute measurements provided by the BOM are the mean of one-second measurements of the proceeding minute. The ground measurements' uncertainties given in BoM [46] were calculated using ISO guidelines.

Observed Cases
To demonstrate the ability of WRF-Solar in simulating irradiance under extreme cloudy and clear cases, we chose four highly intermittent and clear case days using the daily variability index (DVI) and the daily clearness index (DCI) for model evaluation based on observations of GHI at Mildura. DVI and DCI have been calculated using formulations of Huang et al. [47]: where GHI i is the measured GHI within time ∆t with n i=1 ∆t i = 1 day, and CSI i being the clear sky model derived GHI using [48]. Note, the clear days have been chosen as control days to compare simulated accuracy from highly intermittent days. In the following thresholds developed by Huang, Troccoli and Coppin [47], highly intermittent days were chosen with DVI > 30 and 0.95 > DCI > 0.69 and clear days with DVI < 1.2 and DCI >1, as described in Table 1.

Satellite Data
The Australian Bureau of Meteorology (BOM) produces hourly estimates of solar irradiance variables (DNI and GHI) at a resolution of 0.05° (≈5 km) over Australia. Since this study focused on observed cases of highly intermittent and clear days, we only used satellite data from selected dates in 2005, as shown in Figures 1 and 2. GHI estimates were derived from the Geostationary Meteorological Satellites GMS-4 and GMS-5, Geostationary Operational Environment Satellite (GOES-9), and Multi-functional Transport Satellites (MTSAT), which is a series of geostationary operated by the Japan Meteorological Agency (JMA). GHI estimates were computed from raw satellite reflectance based on the two-band physical model [49]

Satellite Data
The Australian Bureau of Meteorology (BOM) produces hourly estimates of solar irradiance variables (DNI and GHI) at a resolution of 0.05° (≈5 km) over Australia. Since this study focused on observed cases of highly intermittent and clear days, we only used satellite data from selected dates in 2005, as shown in Figures 1 and 2. GHI estimates were derived from the Geostationary Meteorological Satellites GMS-4 and GMS-5, Geostationary Operational Environment Satellite (GOES-9), and Multi-functional Transport Satellites (MTSAT), which is a series of geostationary operated by the Japan Meteorological Agency (JMA). GHI estimates were computed from raw satellite reflectance based on the two-band physical model [49] and corrected for biases. A modified form of the Ridley, et al. [50] model converted the bias-corrected GHI values into DNI. Earlier comparisons conducted with older datasets [51] showed significant improvements in the quality of

Satellite Data
The Australian Bureau of Meteorology (BOM) produces hourly estimates of solar irradiance variables (DNI and GHI) at a resolution of 0.05 • (≈5 km) over Australia. Since this study focused on observed cases of highly intermittent and clear days, we only used satellite data from selected dates in 2005, as shown in Figures 1 and 2. GHI estimates were derived from the Geostationary Meteorological Satellites GMS-4 and GMS-5, Geostationary Operational Environment Satellite (GOES-9), and Multi-functional Transport Satellites (MTSAT), which is a series of geostationary operated by the Japan Meteorological Agency (JMA). GHI estimates were computed from raw satellite reflectance based on the two-band physical model [49] and corrected for biases. A modified form of the Ridley et al. [50] model converted the bias-corrected GHI values into DNI. Earlier comparisons conducted with older datasets [51] showed significant improvements in the quality of the derived data. For each day, the hourly measurements started at 18 UT on the proceeding day and ended at 11 UT of that day. The satellite observations were made once every hour, and the time of observation varies smoothly with latitude and differs between satellites. Further information about this dataset can be found in Prasad et al. [52] and Dehghan, Prasad, Sherwood and Kay [23].

Reanalysis Forcing Data
The initial and lateral boundary conditions data include ERA-Interim data from the European Center for Medium-Range Weather Forecasts (ECMWF). ERA-Interim is a global reanalysis data available at 0.71 • grid spacing on 60 vertical levels from the surface up to 0.1 hPa, at six-hourly time steps. The reanalysis dataset was downloaded from the NCAR Research Data Archive. The sea-surface temperature (SST) data was downloaded from archives of daily real time-global SST analysis data available at lower-resolution of 0.5 • from the National Oceanic and Atmospheric Administration (NOAA).

Model Setup
For this study, we used the WRF version 3.6.1 with the Advanced Research WRF (ARW) dynamical solver [53,54] designed for solar energy applications (WRF-Solar) [36]. The models' domain configuration is shown in Figure 3. The setup was nested with a 3:1 parent grid ratio and used during simulations starting at different resolutions described in Table 2.

Reanalysis Forcing Data
The initial and lateral boundary conditions data include ERA-Interim data from the European Center for Medium-Range Weather Forecasts (ECMWF). ERA-Interim is a global reanalysis data available at 0.71° grid spacing on 60 vertical levels from the surface up to 0.1 hPa, at six-hourly time steps. The reanalysis dataset was downloaded from the NCAR Research Data Archive. The seasurface temperature (SST) data was downloaded from archives of daily real time-global SST analysis data available at lower-resolution of 0.5° from the National Oceanic and Atmospheric Administration (NOAA).

Model Setup
For this study, we used the WRF version 3.6.1 with the Advanced Research WRF (ARW) dynamical solver [54,55] designed for solar energy applications (WRF-Solar) [36]. The models' domain configuration is shown in Figure 3. The setup was nested with a 3:1 parent grid ratio and used during simulations starting at different resolutions described in Table 2.  Previously, multiple studies have investigated different combinations of physics schemes to study solar radiation using WRF. In order to evaluate GHI values over a domain in the Iberian  Peninsula in Spain, Rincón et al. [55] used an optimal combination of physical schemes. Similarly, the evaluation of WRF solar irradiance forecasts conducted by Lara-Fanego et al. [56] used similar parameterization schemes as Rincón, Jorba, Baldasano and Monache [55] except for the Microphysical scheme, where the Thompson scheme was used. Otkin and Greenwald [57] have shown that the Thompson scheme simulates clouds more accurately; therefore, it was used by Mathiesen et al. [58] to develop a cloud-assimilating WRF-based model that improved the forecast of solar radiation in California. Also, Lopez-CotoI. et al. [59] used 72 combinations of WRF physical parameterizations from various schemes to assess solar forecasting errors using WRF in coastal California. Another study in Romania [60] used an optimal combination of schemes which produced the smallest rRMSE based on the GHI data obtained for five stations.
For this study, an optimal combination to simulate irradiance was adapted from Evans et al. [61], who showed that Mellor-Yamada-Janjic Planetary Boundary Layer (PBL) scheme and the Kain-Fritsch cumulus scheme worked well in simulating of a series of rainfall events near the south-east coast of Australia. Therefore, for this study, the model was run with two-way nesting using 51 terrain-following (ETA) vertical levels, of which ten are located in the lowest two km of the atmosphere, with a model top set at 50 hPa. The common physics options chosen included a WRF double-moment 5-class microphysics [62], Kain-Fritsch cumulus parametrization [63], Mellor-Yamada-Janjic PBL scheme [64], Dudhia shortwave radiation [65], rapid radiative transfer model longwave radiation [66], and the Noah land-surface model [67]. The model time step used for the outermost domain was 120 s, and the model output interval was chosen as 10 min for easier comparison with the satellite data. The radiation scheme was called every 10 min, and the shortwave radiation was interpolated based on the updated solar zenith angle between radiation calls [68]. The ERA-Interim initial and boundary conditions were updated every 6 h with default aerosol climatology [36].
Two separate configurations were designed for short-term (cases) and long-term (full year) simulations. The case experiments were run for 48 h with the first 12 h excluded from results due to model spin-up. Each simulation was initialised at 12:00 UTC the day before the case of interest. The full year experiment was only run with the first two domains due to increased computational costs.
The full year runs were updated with daily SST [69] and spectral nudging was used to constrain the simulations to the large-scale circulation for wavenumbers greater than 1000 km for wind, potential temperature, water vapor mixing ratio, and geopotential height above the boundary layer [70].

Evaluation Procedure
The model was primarily evaluated for solar irradiance variables (GHI, DNI, and DHI) against observations from the Mildura site at common times. Additionally, other surface variables measured from the site, such as relative humidity (RH), sea-level pressure (SLP), temperature (T), windspeed (WS), and wind-direction (WD), were also compared with the model to evaluate the simulated meteorological conditions at every model time step (120 s). The evaluation metrics used in this study include the mean bias error (MBE), root mean square error (RMSE) and the Pearson's correlation coefficient (r): where I i m and I i o are respectively modelled and observed values at the i th index from n number of samples with mean values I m and I o . Note, for WD, MBE and RMSE were calculated using directional differences [71,72] defined as: where positive (negative) differences relates to the modelled direction rotated clockwise (counterclockwise) with respect to the observed values.

Evaluating Solar Irradiance Variables
The performance of solar irradiance variables including GHI, DNI, and DHI for cases described in Table 1 are shown in Figure 4. The errors show two significant patterns. Firstly, errors (RMSE) in GHI mostly increase with resolution, especially on intermittent days (cases 1-4). This also holds true for DHI and DNI, but GHI is mostly underestimated by the model at coarser resolution, while DNI is overestimated at finer resolution. Secondly, clear days (cases 5-8) show far lower errors (MBE, RMSE) compared to the intermittent days with very high correlations (r > 0.95). Both GHI and DNI are slightly underestimated on clear days.
where and are respectively modelled and observed values at the index from number of samples with mean values and . Note, for WD, MBE and RMSE were calculated using directional differences [72,73] defined as: where positive (negative) differences relates to the modelled direction rotated clockwise (counterclockwise) with respect to the observed values.

Evaluating Solar Irradiance Variables
The performance of solar irradiance variables including GHI, DNI, and DHI for cases described in Table 1 are shown in Figure 4. The errors show two significant patterns. Firstly, errors (RMSE) in GHI mostly increase with resolution, especially on intermittent days (cases 1-4). This also holds true for DHI and DNI, but GHI is mostly underestimated by the model at coarser resolution, while DNI is overestimated at finer resolution. Secondly, clear days (cases 5-8) show far lower errors (MBE, RMSE) compared to the intermittent days with very high correlations (r > 0.95). Both GHI and DNI are slightly underestimated on clear days.
The errors on intermittent days likely result from the misrepresentation of clouds and circulation, while clear day related errors may relate to aerosol loading in the model. Note, compensating effects of errors in DNI and DHI produce lesser biases in GHI, especially on clear days (MBE < 15 Wm −2 ). DNI values on clear days are likely underestimated from lower aerosol loadings provided in the forcing data. This study focuses mostly on GHI; thus, to further explore GHI related errors, surface variables related to clouds and circulation need further investigation.  values on clear days are likely underestimated from lower aerosol loadings provided in the forcing data. This study focuses mostly on GHI; thus, to further explore GHI related errors, surface variables related to clouds and circulation need further investigation.

Evaluating Surface Variables
The surface interactions with the planetary boundary layer are crucial for the transfer of heat, moisture, and momentum in driving atmospheric circulation [73]. The interaction of clouds with radiation maintains atmospheric circulation and keeps the Earth's energy in balance. Therefore, clouds and circulation are critical in determining the amount of solar irradiance incident on a surface. Cloud transmission would reduce incoming solar irradiance, whereas circulation determines the location and height of clouds. Most of intermittent GHI observed at stations can be directly related to errors in the movement and the types of clouds overhead.
To understand the simulated GHI intermittency, surface variables influencing cloud formation and cloud movement (circulation) were evaluated at the site. Cloud formation is conducive to moisture (humidity) which also is affected by temperature. The performance of surface T and RH for cases described in Table 1 are shown in Figure 5.

Evaluating Surface Variables
The surface interactions with the planetary boundary layer are crucial for the transfer of heat, moisture, and momentum in driving atmospheric circulation [74]. The interaction of clouds with radiation maintains atmospheric circulation and keeps the Earth's energy in balance. Therefore, clouds and circulation are critical in determining the amount of solar irradiance incident on a surface. Cloud transmission would reduce incoming solar irradiance, whereas circulation determines the location and height of clouds. Most of intermittent GHI observed at stations can be directly related to errors in the movement and the types of clouds overhead.
To understand the simulated GHI intermittency, surface variables influencing cloud formation and cloud movement (circulation) were evaluated at the site. Cloud formation is conducive to moisture (humidity) which also is affected by temperature. The performance of surface T and RH for cases described in Table 1 are shown in Figure 5. Both T and RH errors (RMSE) increase with the resolution on intermittent days (cases 1-4). Although the correlations of T and RH are high (r > 0.8) for most cases, T is mostly underestimated by the model by as much as 1.5 K and RH is overestimated by as much as 15%. A cooler and moist surface in the model will most likely influence the diurnal cycle of convection and the formation of clouds. This will influence GHI on overcast and intermittent days. For example, the intermittent case 2 shows the highest error (RMSE) in GHI coincident with errors observed with T and RH.
Similarly, cloud motion is determined by circulation, which can be characterized by large-scale motions represented by SLP and local advection based on WS and WD. The performances of these variables at the surface for cases described in Table 1 are shown in Figure 6. A strikingly clear pattern Both T and RH errors (RMSE) increase with the resolution on intermittent days (cases 1-4). Although the correlations of T and RH are high (r > 0.8) for most cases, T is mostly underestimated by the model by as much as 1.5 K and RH is overestimated by as much as 15%. A cooler and moist surface in the model will most likely influence the diurnal cycle of convection and the formation of clouds. This will influence GHI on overcast and intermittent days. For example, the intermittent case 2 shows the highest error (RMSE) in GHI coincident with errors observed with T and RH.
Similarly, cloud motion is determined by circulation, which can be characterized by large-scale motions represented by SLP and local advection based on WS and WD. The performances of these variables at the surface for cases described in Table 1 are shown in Figure 6. A strikingly clear pattern is observable with SLP for all the cases: an increase in RMSE with increase in the resolution, with modest correlations with observations (r > 0.6). Irrespective of the type of case investigated, the model consistently underestimates SLP by 0.5-2 hPa. Lower pressure in the model may influence the strength of winds produced, and thus the movement of clouds. This is also shown by evaluating WS and WD at the site.
Errors (RMSE) in WS also increase with the resolution, but the bias is significantly larger on intermittent days (cases 1-4) by as much as 2 ms −1 . Note, observed winds shows modest correlations with observations on most days (r > 0.6), except for a few days. The overestimation of WS may also relate to the underestimation of SLP, but it may also be locally driven by land-surface characteristics. Nonetheless, biases in WS added together with truth data can introduce displacement errors in observed and modelled clouds which can cause significant biases in GHI. Surprisingly, WD errors only slightly increased with resolution, but MBE is mostly negative (modelled direction is rotated counterclockwise relative to the observed directions). Note, errors in WD on intermittent days would also affect the displacement of clouds.
Energies 2020, 13, x 9 of 23 modest correlations with observations (r > 0.6). Irrespective of the type of case investigated, the model consistently underestimates SLP by 0.5-2 hPa. Lower pressure in the model may influence the strength of winds produced, and thus the movement of clouds. This is also shown by evaluating WS and WD at the site. Errors (RMSE) in WS also increase with the resolution, but the bias is significantly larger on intermittent days (cases 1-4) by as much as 2 ms −1 . Note, observed winds shows modest correlations with observations on most days (r > 0.6), except for a few days. The overestimation of WS may also relate to the underestimation of SLP, but it may also be locally driven by land-surface characteristics. Nonetheless, biases in WS added together with truth data can introduce displacement errors in observed and modelled clouds which can cause significant biases in GHI. Surprisingly, WD errors only slightly increased with resolution, but MBE is mostly negative (modelled direction is rotated counterclockwise relative to the observed directions). Note, errors in WD on intermittent days would also affect the displacement of clouds.

Bias Exploration
To further explore the biases observed in solar radiation variables, a few of the cases are explored in greater detail in terms of the representation of clouds and circulation, as described in Table 3. Table 3. Description of cases for bias exploration.

Bias Exploration
To further explore the biases observed in solar radiation variables, a few of the cases are explored in greater detail in terms of the representation of clouds and circulation, as described in Table 3. The representation of clouds and circulation described in Table 3 is also highlighted in Figure 7. The case on 2 December shows a good representation of clouds, although with biases in surface T and RH during the day (Figure 7-left). Similarly, WS and WD show good agreements during the day. It is not surprising to see that the intermittency in GHI from clouds is captured by the model at all resolutions due to the improved simulation of the type and location of clouds.
Energies 2020, 13, x 10 of 23 The representation of clouds and circulation described in Table 3 is also highlighted in Figure 7. The case on 2 December shows a good representation of clouds, although with biases in surface T and RH during the day (Figure 7-left). Similarly, WS and WD show good agreements during the day. It is not surprising to see that the intermittency in GHI from clouds is captured by the model at all resolutions due to the improved simulation of the type and location of clouds.  Table 3. This is also highlighted by comparing the spatial distribution of simulated and satellite-derived GHI at 5 km resolution shown in Figure 8. The model simulates broken clouds influencing GHI, similarly to the observed values at higher resolution. There are minor errors likely from misrepresentations of cloud types due to biases in T and RH. This is a classic case when the intermittency in the model slightly matches the observations due to better representation of broken clouds and winds.  Table 3. This is also highlighted by comparing the spatial distribution of simulated and satellite-derived GHI at 5 km resolution shown in Figure 8. The model simulates broken clouds influencing GHI, similarly Energies 2020, 13, 385 11 of 22 to the observed values at higher resolution. There are minor errors likely from misrepresentations of cloud types due to biases in T and RH. This is a classic case when the intermittency in the model slightly matches the observations due to better representation of broken clouds and winds. On the contrary, the case on 6 January is a classic case of failed representation of circulation which creates displacement errors in otherwise well simulated clouds, leading to significant biases in the simulated GHI. The model completely misses capturing any variability likely from displacement errors caused by WS (Figure 7-centre). The model simulated much stronger winds than those which were observed, while the clouds were likely better simulated due to reduced bias in T and RH. This is further elaborated in Figure 9, which shows clouds influencing GHI are displaced in the model due to the difference in the timing of frontal pass of clouds over the site in the model and observation. On the contrary, the case on 6 January is a classic case of failed representation of circulation which creates displacement errors in otherwise well simulated clouds, leading to significant biases in the simulated GHI. The model completely misses capturing any variability likely from displacement errors caused by WS (Figure 7-centre). The model simulated much stronger winds than those which were observed, while the clouds were likely better simulated due to reduced bias in T and RH. This is further elaborated in Figure 9, which shows clouds influencing GHI are displaced in the model due to the difference in the timing of frontal pass of clouds over the site in the model and observation.
Moreover, the case on 12 June is a classic case of proper representation of circulation on a mostly clear day observed at the site. The model reproduces most of the observations, including WS (Figure 7-left). This is further elaborated in Figure 10 with good spatial correlations amongst the model and observations. Evidently, a better representation of WS in the model avoids the influence of advection of unrealistic clouds much further away from the site. Note the noise observed with WD ( Figure 7-bottom left) is a plotting artefact due to values closer to 359 and 0 degrees.

Evaluation of Long-Term Simulations
Running longer-term simulations, such as for a full year, with four nested domains ranging from 45 to 1.67 km, was computationally expensive and the results from short-term (cases) simulations ascertained that increased resolution does not necessarily improve errors from clouds and circulation at higher spatial scales. Therefore, longer-term simulations were conducted only with the first two

Evaluation of Long-Term Simulations
Running longer-term simulations, such as for a full year, with four nested domains ranging from 45 to 1.67 km, was computationally expensive and the results from short-term (cases) simulations ascertained that increased resolution does not necessarily improve errors from clouds and circulation at higher spatial scales. Therefore, longer-term simulations were conducted only with the first two domains (45 and 15 km) for statistical exploration of the biases observed in solar irradiance variables. Figures 11 and 12 show scatterplots of observed and modelled surface and solar irradiance variables at 45 and 15 km, respectively. This includes days with clear, intermittent and overcast conditions. Overall, the observed and modelled surface and solar irradiance variables align well with the line of agreement at both 45 km and 15 km. The largest spread is seen for surface variables RH and WD. WS measured at the site are clustered due to possible conversions from raw data, but still in reasonable agreement with modelled values. GHI also was in substantial agreement as compared to DNI and DHI.  The error and agreement metrics after evaluation of surface and solar irradiance variables for long-term simulations at d01 (45 km) and d02 (15 km) are reported in Table 4.   The error and agreement metrics after evaluation of surface and solar irradiance variables for long-term simulations at d01 (45 km) and d02 (15 km) are reported in Table 4.  The error and agreement metrics after evaluation of surface and solar irradiance variables for long-term simulations at d01 (45 km) and d02 (15 km) are reported in Table 4. At both resolutions, SLP and T are underestimated by the model, whereas RH and WS are overestimated. Note, WD has a directional offset which is negative (modelled direction is rotated counterclockwise relative to the observed directions). Errors in DNI and DHI compensate to produce better GHI at coarser resolution. Weak correlations in WS and WD (r < 0.4) ascertain circulation errors to be pivotal in comparing models with observations. Overall, errors (RMSE) increases with the resolution for most variables without much change in correlations. Simulations at the coarsest domain are quantitatively better when compared to observations; however, the model's ability to capture clear (DCI > 0.95), intermittent (0.5 < DCI < 0.95) and overcast (DCI < 0.5), cases daily and the associated DVI are demonstrated in Figure 13.  Out of the 343 days simulated, 175 days were identified as clear, 138 days were intermittent, and 30 days had overcast conditions using DCI calculated from observations. The number of days in agreement was determined to be within ±0.1 of the line of agreement. Thus, the model captured only 47 days as clear, five days as intermittent and three days as overcast equating to 27%, 4%, and 10% days being identified as clear, intermittent, and overcast, respectively. On the contrary, comparison of observed and modelled DVI to be in agreement within ±0.1 of the line of agreement showed the model captured intermittency of 52 clear days only without any success in capturing intermittency during DCI > 0.5 (intermittent and overcast days). Clearly, the model's ability to reproduce intermittency as cloudiness increases is severely limited. The model underestimates intermittency in observation by a mean factor of three for clear and overcast days, but that increases to 13 for intermittent days. A common problem identified on many days was the model's misrepresentation of cloudy days observed to be simulated as clear, or vice-versa. To further assess such misrepresentations, especially on days of high intermittency, Figure 14 illustrates the relationship of cloudiness increases is severely limited. The model underestimates intermittency in observation by a mean factor of three for clear and overcast days, but that increases to 13 for intermittent days. A common problem identified on many days was the model's misrepresentation of cloudy days observed to be simulated as clear, or vice-versa. To further assess such misrepresentations, especially on days of high intermittency, Figure 14 illustrates the relationship of errors (RMSE) from the surface and solar irradiance variables with respect to the varying degree of intermittency at the best-performing resolution (45 km). To the first order, errors in surface and solar irradiance variables are linearly related to DVI. RMSE in GHI, DNI and DHI increase strongly on highly intermittent days (DVI > 30) with errors over 200 Wm −2 , 400 Wm −2 , and 100 Wm −2 , respectively. Errors (RMSE) in solar irradiance variables are likely caused by cloud related errors on days of high intermittency, which can be as much as 20% in surface RH and 4 K in T. Similarly, circulation errors can be associated with surface SLP, WS, and WD by as much as 5 hPa, 3 ms −1 , and 60 degrees, respectively. This is further demonstrated in Figure  15 with effects on key variables in simulating GHI.  To the first order, errors in surface and solar irradiance variables are linearly related to DVI. RMSE in GHI, DNI and DHI increase strongly on highly intermittent days (DVI > 30) with errors over 200 Wm −2 , 400 Wm −2 , and 100 Wm −2 , respectively. Errors (RMSE) in solar irradiance variables are likely caused by cloud related errors on days of high intermittency, which can be as much as 20% in surface RH and 4 K in T. Similarly, circulation errors can be associated with surface SLP, WS, and WD by as much as 5 hPa, 3 ms −1 , and 60 degrees, respectively. This is further demonstrated in Figure 15 with effects on key variables in simulating GHI.
As expected, GHI errors are strongly correlated with DNI (r > 0.75) and DHI (r > 0.60) errors. Errors from other surface variables also show moderate correlations with errors in GHI, but only errors in RH, T, and WS showed significant correlations at a 95% confidence interval (p < 0.05). RMSE in GHI, DNI and DHI increase strongly on highly intermittent days (DVI > 30) with errors over 200 Wm −2 , 400 Wm −2 , and 100 Wm −2 , respectively. Errors (RMSE) in solar irradiance variables are likely caused by cloud related errors on days of high intermittency, which can be as much as 20% in surface RH and 4 K in T. Similarly, circulation errors can be associated with surface SLP, WS, and WD by as much as 5 hPa, 3 ms −1 , and 60 degrees, respectively. This is further demonstrated in Figure  15 with effects on key variables in simulating GHI. As expected, GHI errors are strongly correlated with DNI (r > 0.75) and DHI (r > 0.60) errors. Errors from other surface variables also show moderate correlations with errors in GHI, but only errors in RH, T, and WS showed significant correlations at a 95% confidence interval (p < 0.05).

Discussion
Our results show that under all-sky conditions WRF-Solar simulated at 45 km resolution GHI, DNI and DHI produce errors (RMSE) of 134 Wm −2 , 248 Wm −2 , and 67 Wm −2 , respectively, which grew to 200 Wm −2 , 400 Wm −2 , and 100 Wm −2 under highly intermittent days (DVI > 30) at the Mildura site. Although WRF-Solar has not been extensively tested in Australia, similar error statistics were reported in other countries, especially under all-sky conditions. WRF-Solar GHI simulations from 2014 to 2016 in Singapore produced an RMSE of 242 Wm −2 [43]. Similarly, WRF-Solar simulations in 2017 in Kuwait produced RMSEs in GHI and DNI of 101Wm −2 and 137 Wm −2 , respectively [44]. Despite regional weather events dominating error statistics, results from WRF-Solar simulations in other countries show the consensus on misrepresentation of clouds and aerosols has not been fully resolved, but has improved from native WRF simulations [36][37][38]41,[43][44][45].
In Australia, native WRF simulations of DNI in 2009 at the Wagga-Wagga site (nearest to Mildura) showed MBE of ≈98 Wm −2 with days of extreme dust storm producing MBE of 388 Wm −2 [74]. In this study, DNI errors in MBE were about 24 Wm −2 . Direct comparisons to this study cannot be made, but the range of MBE errors produced at sites near Mildura is much higher as compared to the present study, suggesting the immense potential of WRF-Solar for simulating solar irradiance in regional Australia. Furthermore, our results are much improved over the errors reported by Dehghan, Prasad, Sherwood and Kay [23], who showed GHI errors (RMSE) at Mildura were of the order of 222 Wm −2 when simulated with the Air Pollution Model (TAPM). They also reported much improved performance in GHI from the coarsest domain of 45 km and attributed errors to misrepresentation of clouds in the model during deep low-pressure systems, the passage of cold fronts, easterly troughs, and cloud bands. High intermittent days explored in this study likely encompasses such weather events, but we also show that errors in reproducing circulation patterns may also exacerbate errors related to cloud misrepresentation.
Errors in cloud formation and timing (diurnal cycle) together with their horizontal and vertical positions need further exploration from the solar energy generation perspective. Issues of representing the diurnal cycle of clouds, the convection, and the circulation are prevalent in weather and climate models [16,19,22], but have gained less attraction with applications to solar irradiance forecasts [20]. This study demonstrates errors in simulated GHI to be significantly related to cloud formation (RH, T) and advection (WS) based on surface variables. The heat exchange between the surface and the planetary boundary layer drives circulation and convection influencing cloud formation and movement at various time-scales [73,75]. Errors in RH and T can make the atmospheric column drier or wetter influencing stability of the atmospheric column and the updraft of air parcels affecting cloud formation by limiting the height of supersaturation of parcels and the amount of convective available potential energy to drive convection [76]. In the tropics, satellite observations have shown that RH and thick clouds vary daily by as much as 20% and 10%, respectively, at any altitude over land [77]. Given RH errors in this study are as much as 20% which can be associated with errors in T [78], it is sufficient to say thick clouds will be influenced significantly during a daily cycle.
On the other hand, WS errors at the surface reflect the intensity of moisture mixing and atmospheric circulation, affecting the location and timing of the clouds [79]. Several studies have demonstrated the sensitivity of the diurnal cycle of clouds associated with the intensity of windspeeds driving thermal to mechanical induced lifting and triggering of convection [80][81][82]. Note errors in WS of 3 ms −1 can introduce displacement errors of 360 m in 120 s (model time step) at the surface. This will scale logarithmically with height and for typical low-level clouds (≈3 km) with high transmission produce displacement errors of around 800 m. Thus, average surface winds of 6 ± 3 ms −1 would geolocate low-level-clouds at most by 2.5 km or with a minimum of 0.8 km, which is likely to offset the location of clouds by at least one grid point at 5 and 1.67 km resolutions. Hence, errors in GHI at finer resolution increase with the strength of the mean winds advecting clouds. Spatial averaging would likely improve finer resolution results [26], but coarser resolution runs avoid such errors.
Therefore, misrepresentation of clouds and circulation in the model are severely undermining the reproduction of solar irradiance intermittency observed at ground stations. Efforts in simulating clouds at the correctly observed time and the location are sophisticated unless most of the physical processes are represented through better parametrizations in models. Possible improvements include the addition of sub-grid scale cloud feedbacks [37], accurate convective triggering [83], enhanced cloud microphysics [17,84,85], boundary layer inversions [86], and turbulence [87] in the model.

Conclusions
This paper assessed the performance of an optimum configuration and augmentation of the WRF model designed for solar energy applications (WRF-Solar) on days of high intermittency characterized by the daily variability index observed at Mildura, Australia. The performance was evaluated using incident solar irradiance variables (GHI, DNI, and DHI) with additional meteorological variables (RH, T, SLP, WS, and WD) measured at the surface.
Initially, results from four cases of highly intermittent and clear days mostly showed errors in simulating GHI increase with the resolution, especially on intermittent days. Similarly, errors in variables influencing cloud formation (RH and T) and circulation (SLP and WS) also increase with resolution. Both temporal (time-series) and spatial (satellite) exploration of GHI errors related to misrepresentation of clouds and circulation for specific cases were demonstrated using the best and worst-case scenarios. The worst scenario was associated with clouds located at different locations in the model and observations due to errors in cloud motion from differences in WS. On the other hand, the best scenario showed an improved simulation of cloud location as well as the cloud spatial structure, which coincided with agreements in cloud formation and circulation variables.
Owing to larger errors at finer resolutions, but for the robustness of results presented with only a few cases, more extended simulations conducted at coarser resolution confirm that misrepresentation of clouds and circulation in the model severely undermines the simulation of GHI. Errors in GHI showed significant correlations at 95% confidence interval (p < 0.05) with errors in RH, T, and WS at highly intermittent days. Overall, WRF-Solar performed better with errors in GHI (RMSE, 134 Wm −2 ) and DNI (MBE, 24 Wm −2 ) compared to standard WRF (DNI MBE, 98 Wm −2 ), and other models (GHI RMSE, 222 Wm −2 ) when compared to relevant studies conducted near the Mildura site. Despite improvements resulting from the optimum configuration and augmentation of standard WRF, Energies 2020, 13, 385 18 of 22 WRF-Solar still underestimates daily observed intermittency by an average factor of 13 related to the simulations of observed cloudy days as clear, or vice-versa.
Undoubtably capturing intermittency in models is a complex problem due to errors in clouds and circulation. To enhance future simulations of incident surface solar irradiance, the representation of clouds and circulation in models must be improved through physical parametrizations by capturing both the dynamic and thermodynamic formation of clouds with better representation of unresolved clouds and turbulence. The study is limited to only one site, but other sites are likely to be explored in future experiments. Similarly, the results presented in this study apply to weather conditions prevalent in the year 2005. Future simulations will explore more recent years to understand the impacts of current weather and climate extremes.