A simple parameterization of the short-wave aerosol optical properties for surface direct and diffuse irradiances assessment in a numerical weather model

Broadband short-wave (SW) surface direct and diffuse irradiances are not typically within the set of output variables produced by numerical weather prediction (NWP) models. However, they are frequently requested for solar energy applications. In order to compute them, a detailed representation of the aerosol optical properties is important. Nonetheless, NWP models typically oversimplify aerosol representation or even neglect their effect. In this work, a flexible method to account for the SW aerosol optical properties in the computation of broadband SW surface direct and diffuse irradiances is presented. It only requires aerosol optical depth at 0.55 μm and knowledge of the type of predominant aerosol. Other parameters needed to consider spectral aerosol extinction, namely, Angström exponent, aerosol single-scattering albedo and aerosol asymmetry factor, are parameterized. The parameterization has been tested using the Rapid Radiative Transfer Model for climate and weather models (RRTMG) SW scheme of the Weather Research and Forecasting (WRF) NWP model for data over the continental US. In principle, it can be adapted to any other SW radiative transfer band model. It has been verified against a control experiment and using data from five radiometric stations in the contiguous US. The control experiment consisted of a clear-sky evaluation of the RRTMG solar radiation estimates obtained in WRF when RRTMG is driven with ground-observed aerosol optical properties. Overall, the verification has shown satisfactory results for both broadband SW surface direct and diffuse irradiances. The parameterization has proven effective in significantly reducing the prediction error and constraining the seasonal bias in clearsky conditions to within the typical observational error expected in well maintained radiometers.


Introduction
Broadband short-wave (SW) surface downward total solar irradiance (also known as global horizontal irradiance, GHI) is the sum of broadband SW surface downward direct normal 35 irradiance (DNI, received from the sun's direction) projected onto a horizontal plane and broadband SW surface downward diffuse irradiance (DIF, received from other directions).In general, DIF may also include reflected irradiance from surrounding areas, such as from mountain sides.Surface ir-40 radiance is also referred to as downward (or downwelling) SW flux in other disciplines.Both DNI and DIF are rarely included in predictions made with Numerical Weather Prediction (NWP) models.However, they are necessary in multiple applications such as those requiring a precise represen-45 tation of surface solar radiation or solar energy applications (Geiger, 1965;Hay, 1993;Whiteman, 2000;Gu et al., 2002;Oliphant et al., 2003;Pierce et al., 2005;Stoffel et al., 2010;Kleissl, 2013).
Surface irradiance is a key component in the representa-50 tion of the energy balance at surface in numerical land and weather models.In the surroundings of gentle terrain, and provided the atmospheric state is known, surface irradiance can be calculated at reasonable accuracy using simple models that assume isotropic conditions in both the sky and the 55 surface.However, under scattered clouds or steep terrain, the isotropy assumption fails.In such a case, a 3D solar radia-tion model would provide the best predictions (Cahalan et al., 2005;Iwabuchi, 2006;Pincus and Evans, 2009).Nonetheless, these models are so computationally intensive that, in practice, their use is restricted only to concrete applications such as validation studies (Mayer et al., 2010) or development of simplified parameterizations (Lee et al., 2011).However, if in particular both DNI and DIF are known, the uneven distribution of solar irradiance over complex terrain areas can be determined.The projection of direct irradiance on tilted surfaces is a straightforward geometrical problem.The exact computation of DIF over the surface would still be unfeasible but, in practice, isotropic or quasi-isotropic assumptions can be used at reasonable accuracy (Ruiz-Arias et al., 2010, 2011;Manners et al., 2012).Energy applications demand an accurate modelling of surface solar irradiances.Both GHI and DNI are acquiring greater importance in the energy sector as the construction rate of solar systems is booming.On the one hand, traditional flat-photovoltaic (PV) systems -the more mature and widely-utilized solar energy technology-are driven primarily by the incoming global irradiance onto the PV plane of array.As this plane very rarely coincides with the horizontal plane (the common irradiance output in most of the NWP models), a transposition model from the horizontal to the PV plane is also required; but accurate transposition models need DNI and DIF irradiances.On the other hand, solar concentrating technologies -both concentrating photovoltaic and solar-thermal plants-are driven primarily by DNI.These technologies increase the overall efficiency of the systems by concentrating DNI using an optical assembly of mirrors or lenses.Overall, solar energy systems require long time-series of GHI and DNI irradiances over wide areas for a proper evaluation of the solar potential or precise simulation of the energy produced by solar systems.Such simulations are most frequently undertaken using historic data.However, and very importantly, large solar power installations now require forecasts that enable an improved operation of the plants and maximize the integration rate of solar systems in the power grid without putting the power supply at risk.This is best done with NWP models for forecast horizons from about 4 to 6 hours onward (Diagne et al., 2013;Inman et al., 2013).
Among the downwelling solar fluxes that can be predicted at the surface, most of the NWP models only provide GHI.It is very likely that this has been motivated by the fact that DNI and DIF are challenging to calculate since they are more sensitive to input parameters, such as aerosols or clouds.Also, surface processes affected by solar radiation can be reasonably well represented with GHI alone, as long as the spatial resolution is more than a few km, which has been the typical case so far.Moreover, accurate calculation of DIF fluxes is computationally intensive compared with the simple methods that can be used to obtain GHI (e.g., Dudhia, 1989).However, the computational capabilities have grown enough to allow the use of more rigorous and precise methods to solve the atmospheric radiative trans-fer equation.Ruiz-Arias et al. (2013c) provide a comprehensive benchmarking study of some of the short-wave radiation schemes available in the Weather Research and Forecasting

115
(WRF) NWP model and their ability to predict GHI, DNI and DIF under clear-sky conditions in the contiguous US region.Albeit the three evaluated models yielded GHI estimates within the observational error range, not all of them showed good skills at predicting DNI and DIF.The best re-120 sults were achieved with the Rapid Radiative Transfer Model for climate and weather models (RRTMG; Iacono et al., 2008).In particular, for the period evaluated, the mean and root-mean square DNI errors when the RRTMG model was run without considering aerosol extinction (default setting 125 in WRF) were 66 W m −2 (7%) and 72 W m −2 (8%), respectively (percent magnitudes are relative to the mean observed value).In contrast, when RRTMG was run with instantaneous observations of aerosol optical properties (hereinafter, AOP), the mean and root-mean square errors decreased to 130 0 W m −2 (0%) and 9 W m −2 (1%), respectively.In the case of DIF, the mean and root-mean square errors when the model was not driven by AOP observations were -26 W m −2 (-34%) and 28 W m −2 (37%), respectively.When AOP observations were used, the mean and root-mean square errors substan-135 tially decreased to 2 W m −2 (3%) and 5 W m −2 (6%), respectively.

The need for a AOP parameterization
Many NWP models solve, or may solve, the solar radiative transfer in the atmosphere using a two-stream approach, 140 which allows for a fast and approximate solution by assuming azimuthal isotropy in radiant fluxes (Ritter and Geleyn, 1992;Edwards and Slingo, 1996;Chou et al., 1998;Iacono et al., 2008).Radiative transfer solvers in NWP models have been simplified by assuming an infinite and horizontally uni-145 form atmosphere and treating each model column independently.The major practical consequence of the two-stream approximation is a reduction in accuracy at large solar zenith angles.However, it is accurate enough under other conditions for most of the current applications.This allows for a 150 sufficiently detailed description of the solar direct and diffuse irradiances at a low-to-moderate spectral resolution.
In the absence of clouds, aerosols become the dominant driving factor for DNI and DIF and the greatest source of uncertainty.In particular, the impact of aerosols on DNI's mag-155 nitude is about 3 to 4 times larger than it is in GHI (Gueymard, 2012;Ruiz-Arias et al., 2013a).This results from the basic compensation effect by which an increase (decrease) of aerosol extinction results in a decrease (increase) of DNI and an increase (decrease) of DIF, in the general case.Thus, er-160 rors in DNI and DIF fluxes caused by a misrepresentation of the aerosol load partly cancel out in GHI, in the general case.In part, this explains why many NWP models have traditionally neglected the direct impact of aerosol in the assessment of GHI, or why it has been simply accounted for by using climatological values.However, this may result in DNI prediction errors up to 20% (Ruiz-Arias et al., 2013a,c).
Extinction by aerosols is described in radiative transfer problems in terms of three spectral quantities, namely, aerosol optical depth (AOD or τ ), single-scattering albedo (SSA or ω 0 ) and asymmetry factor (ASY or g).Aerosol optical depth is the integral of the extinction coefficient over an atmospheric path.The single-scattering albedo is the ratio of the scattering and extinction efficiencies.It represents the relative importance of the scattering events within the total extinction.Finally, the asymmetry factor is the first moment of the scattering phase function.It accounts for the preferred direction in which radiation is scattered (Liou, 2002).It is usual to model the spectral variability of AOD using the Ångström law τ (λ) = βλ −α , where λ is the wavelength in µm, β is the AOD measured at λ=1 µm and α is known as Ångström exponent (AE) ( Ångström, 1961).
The number and variety of region-wide aerosol datasets has steadily grown in recent years, from worldwide ground datasets such as the Aerosol Robotic Network (AERONET; Holben et al., 1998) to sensors aboard satellite platforms that regularly sweep the globe -the most well-known being the Moderate-resolution Imaging Spectro-radiometer (Remer et al., 2005)-.Both sources of data provide AOP observations that could be used in NWP models to evaluate DNI and DIF irradiances.Ground observations, essentially from AERONET, provide a reliable and comprehensive AOP description, at a number of wavelengths.However, the AERONET' spatial density is scarce and its nearreal-time availability is limited.Thus, in practice, its use in NWP model applications is constrained to a reduced number of cases.Satellite retrievals, on the other hand, provide broad spatial coverage, but the accuracy of their current estimates is often only reasonable for AOD at single primary wavelenghts, normally 0.55 µm (Remer et al., 2005).

200
In recent years, coupled Atmosphere-Chemistry Numerical Weather Prediction (ACNWP) models have received considerable attention and improvements.They now benefit from the growing number of available ground and remotely-sensed datasets.Moreover, ACNWP models now routinely offer 205 global forecasts of many molecular and particulate components of the atmosphere.Such is the case of the Monitoring Atmospheric Composition and Climate project (MACC, 2013) or the Goddard Earth Observing System model version 5 (GEOS-5, 2013).They evaluate AOP from prognoses of the chemical composition of the atmosphere and use them to calculate DNI and DIF irradiances.Nonetheless, in general, ACNWP models are computationally expensive and complex to run compared with the regular limited-area NWP models.Moreover, since they are initialized using mostly satellite observations, they suffer from similar biases regarding optical properties of aerosols.
For those applications that depend on NWP models and focus on DNI and DIF irradiances, it is convenient to set up a means to use AOP inputs from diverse sources.The 220 intended benefit of this approach is to use the best aerosol optical source for each potential application.In particular, for long-term evaluations of the regional surface solar radiation potential (commonly referred to as "solar resource assessment"), combined measurements of AOP from ground 225 sites, satellites and/or aerosol transport models could be used (Kinne et al., 2013;Ruiz-Arias et al., 2013b).On the other hand, when the application requires forecasts of surface solar radiation, the AOP predicted by global ACNWP models could be used.Nonetheless, since AOD is typically the only 230 AOP that is available, at least accurately, the rest of the required parameters, namely, SSA, ASY and AE, have to be specified/parameterized based on additional information, in most cases.
In radiances predicted in a one-year WRF simulation using the AOP parameterization is compared against independent surface solar irradiance ground observations in the contiguous US.
Section 3 describes the approach taken for the AOP pa-250 rameterization in the RRTMG SW model.Sections 4 and 5 present the results of a benchmarking study against a control experiment and the validation against ground observations, respectively.Finally, Section 6 highlights the most important conclusions of this work.

3 The AOP parameterization
The RRTMG SW radiative transfer model solves the multiple scattering problem using a two-stream algorithm (Oreopoulos and Barker, 1999) over 14 spectral bands spanning from 0.2 to 12.2 µm (Table 1).It accounts for extinction by 260 water vapor, carbon dioxide, ozone, methane, oxygen, nitrogen, aerosols, Rayleigh scattering and clouds.Under clear skies, the expected accuracy of RRTMG with respect to lineby-line calculations is about 4 W m −2 for direct fluxes and about 5 W m −2 for diffuse fluxes (Iacono et al., 2008).

265
In this study, the aerosol optical properties, which must be provided to the radiative transfer routine at every grid-cell of the domain being simulated and each RRTMG spectral band, are parameterized in terms of the vertically-integrated (total) AOD at 0.55 µm (τ 0.55 ) and a reference aerosol type, in a 270 similar way as it is done in many detailed radiative trans- fer models (Ricchiazzi et al., 1998;Gueymard, 2001;Berk et al., 2005;Mayer and Kylling, 2005).The reason why AOD is not parameterized here is twofold.First, optical depth is the specific property that dominates aerosol extinction in the shortwave, and thus it is important to make use of the best estimate available.Second, unlike other aerosol optical properties, both satellite retrievals and ACNWP models provide reasonable estimates of AOD for many current applications.
The reason to choose the value at 0.55 µm is to be consistent with the values usually provided by these data sources and the ground observations at AERONET.The parameterized reference aerosols are used to provide spectral values for SSA, ASY and AE, which are afterwards modulated in terms of the ambient relative humidity to account for the aerosol hygroscopicity.
Two different reference aerosol models from Shettle and Fenn (1979), namely rural and urban, have been included so far in the parameterization.They are representative of broad continental climate conditions.The rural aerosol model is intended for situations where the aerosol is not expected to be significantly affected by urban or industrial sources.Thus, it is expected to be the typical choice for nearly all simulations since large solar systems are typically absent from urban or industrial areas.The rural aerosol model is composed of a mixture of 70 percent of water soluble substance and 30 percent dust-like aerosols.In contrast, the urban aerosol is a mixture of rural aerosol (80 percent) and soot-like particles (20 percent).The two reference mixtures define the absorption, scattering and extinction coefficients, single-scattering albedo and asymmetry parameter for a number of wavelengths and relative humidities from 0% to 99%.The selection of these two reference aerosol models is motivated by their demonstrated ability to correctly represent clear-sky surface solar irradiances in other radiative transfer models (Ricchiazzi et al., 1998;Gueymard, 2001Gueymard, , 2008)).

Aerosol optical depth and Ångström exponent
AOD has to be specified at each RRTMG spectral band.In real applications, even in the best cases, AOD is only known (through measurement or modeling) at a small number of wavelengths, and the Ångström law is often used to describe its spectral variability.However, for some aerosol particle ensembles, such as the reference aerosol mixtures used here, this spectral variability is best described using a 2-band version of the Ångström law (Gueymard, 2001) as follows: where λ is the wavelength in µm and α i is the Ångström exponent for each band, defined as α i = α 1 for λ <0.55 µm, and α i = α 2 otherwise.The coefficients α i are obtained from the selected reference aerosol models by linearly fitting (in 320 log-log coordinates) the spectral extinction coefficients tabulated in Shettle and Fenn (1979) for each aerosol mixture and relative humidity.The corresponding values of α i are given in Table 2.For α 1 , the extinction coefficients at 0.337 µm, 0.55 µm and 0.649 µm were used.The values at 325 0.55 µm, 0.649 µm, 1.06 µm and 1.536 µm were used for α 2 .This modelling approach is found better than the regular Ångström law at resolving the distinct spectral contribution of the fine and coarse modes of the aerosol size distribution.
The fact that α 1 and α 2 show distinct values suggests this 330 approach is pertinent, at least for actual aerosol mixtures that behave like the selected aerosol models.The limit for the calculation of α 1 and α 2 (λ=0.55 µm) is similar to the limit of 0.6 µm suggested by Dubovik et al. (2002) to distinguish between the fine mode and the coarse mode in bimodal size 335 distributions.The decreasing α i values for increasing relative humidities indicate a particle size increase by water uptake and a relative extinction decrease at lower wavelengths.
It is worth mentioning that, unexpectedly, Ångström exponents for the rural aerosol model are greater than for the ur-340 ban aerosol model, indicating that overall the particles in the urban mixture have a larger size.This is very likely due to the assumption made by Shettle and Fenn (1979) that the soot-like particles in the urban mixture have the same size distribution as the water soluble and dust-like particles in the 345 rural aerosol mixture, despite the fact that soot particles are generally of smaller size.
In order to provide a representative value of AOD over each RRTMG spectral band, it is averaged for each band using Eq.(1).To take into account that the solar spectral irra-350 diance changes abruptly in the ultraviolet and visible regions and that some RRTMG infrared bands are wide, the extraterrestrial solar spectrum, E 0n (λ), as described by Gueymard (2004), is used as a weighting factor to compute the band- where j stands for each RRTMG spectral band, which extends over the range ∆λ j , and τ r (α ri ; λ) is the aerosol optical depth calculated with Eq. ( 1) for the relative humidity r.Factorizing τ 0.55 out of τ r (α ri ; λ), Eq. ( 2) can be re-written where ρ rj is the spectral scale factor with respect to τ 0.55 for the spectral band j and relative humidity r.It is given by Equation ( 4) was numerically evaluated for each RRTMG spectral band and relative humidity according to the α i coefficients in Table 2.The so-computed spectral scale factor values ρ rj were grouped in two look-up-tables (LUT) for the two aerosol types (Tables A1 and A2).For each RRTMG spectral band, the spectral scaling factors are interpolated using a 4-point Lagrange interpolation at the relative humidity values predicted by the NWP model.AOD is then calculated using Eq. ( 3) and the input τ 0.55 .Figure 1 illustrates the interpolation results for the rural aerosol mixture.It also compares the E 0n -weighted average as defined by Eq. ( 2) with a regular (un-weighted) average.The largest discrepancies appear in the ultraviolet, visible and near-infrared regions (RRTMG bands 8-12) as well as in the mid-infrared region (RRTMG band 14).The weighted average shifts the averaged AOD value towards wavelengths with higher extraterrestrial solar intensity.This results in an enhancement of aerosol extinction in the visible and infrared bands, and a decreased extinction in the ultraviolet region.Shettle and Fenn (1979) provide spectral values of SSA and ASY up to 40 µm starting at 0.2 µm for each aerosol mixture and relative humidity value.SSA is spectrally weighted for each band as follows:

Single-scattering albedo and asymmetry factor
where ωo,rj is the average SSA value for relative humidity r and spectral band j.The tabulated values of SSA for each relative humidity are interpolated to the wavelengths at which E on (λ) is known using cubic splines, which results in a series of values noted ωo,r (λ).Equation ( 5) assigns a 395 higher weight to the wavelengths at which extraterrestrial solar spectral irradiance and aerosol extinction are greater.The results of this calculation are grouped in one LUT for each aerosol mixture (Tables A3 and A4).These values are then interpolated for each spectral band and relative humidity us-400 ing a 4-point Lagrange interpolation.Following a similar approach, the spectrally-averaged asymmetry factor is calculated as: where ḡrj is the average ASY value for relative humid-405 ity r and spectral band j.The tabulated values of ASY for each relative humidity are interpolated to the wavelengths at which E on (λ) is known using cubic splines, which results in a series of values noted ĝr (λ).In this case, a higher weight is assigned to those wavelengths with greater E on (λ) and scat-410 tering coefficient.The resulting values of ḡrj are grouped in two LUT's corresponding to the two aerosol mixtures (Tables A5 and A6).These values are then interpolated for each spectral band and relative humidity using a 4-point Lagrange interpolation.

415
Figure 2 shows the parameterized SSA and ASY values for the two reference aerosol models and a relative humidity of 80%.The thin line is the resulting interpolation from the tabulated values (cross marks) in Shettle and Fenn (1979), both for SSA and ASY.The thick line is the resulting E 0n -420 weighted average for each model band after applying Eqs.(5 and 6).The shaded region represents the range of variability within each band due to relative humidity, from 0% to 99%.In general, the urban aerosol's SSA (Fig. 2c) has a smaller   value at all wavelengths and a higher sensitiviy to relative humidity than the rural aerosol model (Fig. 2a).Thus, the latter scatters more radiation but responds less to changes in humidity.Note that, for wavelengths above 4 µm, the band-averaged SSA remains close to the SSA value between 4 and 5 µm because the extraterrestrial solar intensity is very small 430 beyond 5 µm.Dubovik et al. (2002) presented an evaluation study of the aerosol optical properties observed during several years at various AERONET sites characterized by distinct aerosol types.At most of these sites, the measured SSA at midvisible wavelengths was about 0.95 (such as in Greenbelt, US; Crete-Paris, France; Bahrain; Solar Village, Saudi Arabia;...) which is a value roughly coincident with that provided by the rural aerosol mixture (see band 10 in Table A3).Only at those sites affected by high pollution (such as in Mexico City) or biomass burning (such as Zambia) the measured SSA at mid-visible wavelengths was smaller than 0.90.To describe such conditions, the urban aerosol mixture would be more appropriate (see band 10 in Table A4), even though it seems to be too much absorbing.
The asymmetry factor values are very similar for the two aerosol mixtures considered here (Figs.2b and 2d), with decreasing forward scattering in the ultraviolet and visible bands, as opposed to increasing values in the infrared up to 3 µm.Beyond 3 µm, it remains constant at about 0.75.

Vertical distribution
The vertical distribution of AOD is modelled after the spectral disaggregation has been completed.The latter is made following Eq.( 3) with spectral scale values ρ rj interpolated according to the relative humidity predicted by the NWP model, but only at the surface level.Then, the spectrally disaggregated τj values at the surface for each band are distributed vertically according to an exponential profile (Ruiz-Arias et al., 2013c) as follows: where z sf c and z toa are the altitudes at the surface and the top of the atmosphere, respectively.The scale height parameter Z h is set to 2.5 km (Gueymard and Thevenard, 2009).By following this procedure the vertically-integrated profile of AOD is consistent with the τ 0.55 value provided as input.
The vertical distribution of SSA and ASY is based only on the relative humidity profile in the NWP model.Therefore, the SSA and ASY vertical profiles resemble the model moisture profile.

Parameterization benchmarking
The consistency of the AOP parameterization at predicting clear-sky surface solar irradiance is first benchmarked with reference to a case study (hereinafter referred to as control experiment) in which the WRF's RRTMG model is driven using AOP and precipitable water data observed at five AERONET sites that have collocated surface solar irradiance observations.The control experiment represents a best-case estimate of the expected model performance at predicting clear-sky GHI, DNI and DIF irradiances.

Control experiment 480
In the control experiment, WRF is run using the RRTMG SW scheme.Clear-sky estimates of GHI, DNI and DIF are computed every 10 minutes for five completely cloudless days at five different locations in the contiguous US, namely, Bondville (IL), Sioux Falls (SD), Table Mountain 485 (CO), Desert Rock (NV) and Southern Great Plains (OK).All these stations are located far from urban or industrial centers, so that their predominant aerosol regime is likely to be well represented by the rural aerosol model.At all sites, concurrent observations of GHI, DNI and DIF, as well 490 as AOP and precipitable water observations from collocated or nearby AERONET sunphotometers, were available.Four of the experimental surface solar irradiance sites (Bondville, hereinafter referred to as BON; Sioux Falls, hereinafter referred to as SXF; Table Mountain, hereinafter referred to 495 as TBL; and Desert Rock, hereinafter referred to as DRA) belong to the Baseline Surface Radiation Network (BSRN; Ohmura et al., 1998) and to the Surface Radiation Network (SURFRAD; Augustine et al., 2005).The fifth site is at the Atmospheric Radiation Measurement Central Facility (here-500 inafter referred to as ARM), near Lamont, Oklahoma.The observed AOPs and precipitable water at the AERONET sites were ingested every 10 minutes at exactly the same time steps at which solar irradiance was computed in the model.The few traces of clouds generated in the WRF model during the 505 simulations were cleared up by setting the cloud mixing ratio to zero in order to ensure completely clear-sky conditions.Note that, since all AOPs were ingested from ground observations, there was no need to parameterize any aerosol property.Thus, the control experiment gives a fair estimate of 510 the RRTMG model performance at computing the clear-sky GHI, DNI and DIF irradiances using ideally accurate inputs.The control experiment is fully described in Ruiz-Arias et al. (2013c).

Test case 515
The simulations of the control experiment were repeated, this time using the AOP parameterization.That is, only the observed AOD at 0.55 µm at the AERONET sites and the type of aerosol were provided to WRF.The rest of the aerosol parameters, namely, AE, SSA and ASY, were parameterized 520 internally according to Sect. 3. Similarly to the control experiment, the model was driven by observations of precipitable water so that the real skill of the aerosol parameterization was better evaluated.Two different simulations, successively assuming rural and urban aerosols, were carried out at each 525 site.Note however that, following the discussion in Sect.3.2, the urban aerosol model is so absorbing that its use is not recommended over the US unless a highly absorbing aerosol is known to exist.An additional simulation for an ideal and perfectly clean atmosphere (i.e., no aerosols) was also conducted.
Figure 3 shows the relative errors of both the control experiment and the test cases as compared to the GHI, DNI and DIF ground observations at each site, and to the combination of all sites (referred to as case ALL in Fig. 3).If the parameterization were perfect, the grey blocks and the colour bars would match.Disagreements are caused by the prescription of the aerosol model, which sets synthetic values for the AOPs.
Figure 3a shows the relative errors in the case of DNI.The discrepancies between the control experiment and the test cases using the AOP parameterization are negligible (below 1% at all sites), regardless of the selection of aerosol mixture.This could be expected because, as far as aerosols are concerned, DNI is only impacted by optical depth, and the AOD at 0.55 µm is the same in both the control experiment and the test cases.The only distinction between the two experiments resides in the spectral distribution of AOD, which is observed in the control experiment, and rather inferred from the assumed aerosol type and the relative humidity predicted by the NWP model in the test cases.Nonetheless, since DNI is a broadband quantity, the impact of AE is reduced and so are the differences between the control experiment and the test cases.On the contrary, when an ideal aerosol-free atmosphere is assumed, the simulated DNI overestimates the observations beyond their expected uncertainty limit.
Figure 3b shows the relative errors in the case of DIF.Discrepancies between the control experiment and the test cases are greater than for DNI because DIF is also impacted by SSA and ASY, which now are parameterized.Specifically, for relative humidities below 90%, the urban aerosol is about 20% to 40% more absorbing than the rural aerosol.As a consequence, the urban aerosol model tends to systematically predict a 15-20% lower DIF than the rural aerosol model.Hence, unlike with DNI, the selection of the most appropriate aerosol type is important to accurately predict DIF.In particular, at four of the sites evaluated in this study, the rural aerosol model yields results that agree reasonably well with the control experiment.Conversely, at the TBL site, the urban aerosol model yields better results, even though this is normally an unpolluted site.The clear days that were originally selected for that site were associated, by chance, to low values of observed SSA.The anomalous preponderance of absorbing aerosols can be likely explained by wildfires, which originated from the Arapaho-Roosevelt National Forest (50 to 100 km away), and were particularly active during those days (Short, 2013).Nonetheless, this particular case serves to show that the urban mixture may be useful under specific circumstances, and not just to represent urban or industrial environments.Under ideally aerosol-free conditions, a systematic underestimation ensues, as could be expected, and reaches about 30%.
In the case of GHI (Fig. 3c), all the experiments provide estimates within the expected observational error range, even under an ideal aerosol-free atmosphere because, as already 585 commented, the large overestimation in DNI is partly compensated by the large underestimation in DIF.Overall, the rural aerosol fits the control experiment better.

Validation against ground observations
A major limitation of the benchmarking study described in 590 the former section comes from the fact that AOD, AE, SSA and ASY all need to be known simultaneously in the control experiment.However, measurement of SSA and ASY is limited by strong practical constraints (Dubovik et al., 2000) that drastically reduce their availability.Nonetheless, since 595 the only external and continuously variable input required by the AOP parameterization is the AOD at 0.55 µm, the validation period with the AOP parameterization can be extended as long as AOD and surface solar irradiance measurements are available.Thereby, two 1-year simulations, 600 with the same WRF model set-up described in Sect.4, have been conducted using the AOP parameterization, used alternatively with the rural and urban aerosol models.The observed AOD at 0.55 µm at the five test AERONET sites presented in Sect. 4 was ingested into WRF every 10 minutes 605 at exactly the same time steps at which GHI, DNI and DIF were computed.The predicted GHI, DNI and DIF are thus validated for the specific time steps having coincident AOD observations and under clear-sky conditions.The latter were identified based on the cloud screening method described in 610 Long and Ackerman (2000).
In addition, the simulation was repeated using WRF's Dudhia SW scheme as a skill reference for the case of GHI.The Dudhia SW scheme is the radiative transfer model most frequently selected by WRF users because of its large speed 615 gain compared to other selectable radiative transfer models in WRF, such as RRTMG.The Dudhia scheme (Dudhia, 1989) consists of a simple broadband parameterization of GHI (over the whole SW spectral range) that explicitly considers extinction by Rayleigh atmosphere and water va-620 por only.It does not account for multiple scattering effects.Extinction by ozone, aerosols, and other molecular absorbers is accounted for by using a bulk scattering parameter that was empirically fixed to represent average turbidity conditions (Zamora et al., 2003(Zamora et al., , 2005)).Further details may be found

Dynamical range performance
The performance of the AOP parameterization for each aerosol type has been analysed throughout the entire range of the observed AOPs at the five experimental sites.Figure 4(a-630 c) shows the relative frequency distribution of the observed AOD at 0.55 µm, the observed and parameterized SSA val-ues, and the observed and parameterized ASY values, respectively.Overall, the AOD values observed at the validation sites are relatively small, although the evaluation period spans an entire year and includes all available observations at the validation sites.The mean value is 0.06, the median is 0.05, and 95% of the values are smaller than 0.12.The mean observed SSA value is 0.92, with 95% of the values greater than 0.75.Very distinct estimations of SSA are made with the rural and urban mixtures.For the rural aerosol, 95% of the SSA values are between 0.92 and 0.94, with a mean value of 0.93.For the urban aerosol, 95% of the SSA values are smaller than 0.68, with a mean value of 0.62. Figure 4c shows the relative frequency distribution of the observed and simulated ASY values.Ninety-five percent of the observations span the range from 0.61 to 0.75, with a mean value of 0.67.The values simulated by the rural aerosol model have a mean of 0.66, and 90% of the data lies between 0.63 to 0.67.In the case of the urban aerosol model, 90% of the ASY values vary from 0.66 to less than 0.67, and their mean is also 0.66.
Since AE is not directly parameterized (note that it has been approximated by means of a two-band model), it has not been shown here for conciseness.However, its effective value can be estimated from the spectral distribution of AOD throughout the RRTMG bands.In so doing, it is found that 99% of the AE values for the rural aerosol model are contained in the range 1.19-1.22,and 99% of the AE values for the urban aerosol model are in the range 1.00-1.06.In contrast, the observations have a much wider range, 90% of them being between 0.72 and 2.59.This means that the effective AE values obtained by the parameterization span a much shorter range than their observed counterparts.This could be expected since only two reference mixtures are considered in the parameterization, as compared to the infinite number of possible mixtures that can exist in the real world.
Figure 4d-f shows the validation results for DNI.In each case, the relative error is within the expected DNI observational error.However, as it can be seen in Fig. 4d, for AOD above 0.05, there is a systematic bias of about 4 Wm −2 between the estimates obtained with the rural and urban aerosol mixtures.A side experiment (not shown here for the sake of conciseness) conducted with the SMARTS radiative transfer model (Gueymard, 2001) has revealed that this discrepancy is compatible with the different AE values obtained for each aerosol type.For AOD values below 0.05, the disagreement with the observations increases slightly.As shown in Ruiz-Arias et al. (2013c), this might be related to the observational uncertainty of the AOD observations from AERONET sites.
The larger uncertainty in AE at low AOD is also a probable contributor.As expected, DNI does not show any apparent trend with SSA and ASY (Fig. 4e-f).
Figure 4g-i shows the validation results for DIF.For all the test sites combined, the DIF estimates obtained with the rural aerosol model are within the expected range of the observational error.However, selecting the urban aerosol yields a negative bias that, interestingly, increases in magnitude for increasing AOD.The reason is that there exists a positive correlation between AOD and SSA in this experimental dataset 690 (not shown here) such that an increase of AOD results in an increase in SSA.In addition, as shown in Fig. 4h, a systematic underestimation of about 15% in the estimated DIF values results from selecting the urban aerosol model.No such bias exists with the rural aerosol model.No trend is observed 695 in the simulated DIF values with respect to ASY (Fig. 4i).
Figure 4j-l shows the validation results for GHI.In addition to GHI evaluated with the RRTMG model assuming either rural or urban aerosols, GHI calculated with the Dudhia SW scheme is also shown.The latter does not make use 700 of any aerosol optical variable as input.In any case, all the simulated values are within the range of the expected observational error.In particular, the RRTMG predictions of GHI have no discernible bias when obtained with the rural aerosol model.On the contrary, when the urban aerosol model is 705 assumed, the bias in DIF (Fig. 4g-i) propagates into GHI, but with a reduced relative impact (about 3%).The Dudhia scheme shows an increasing trend with respect to AOD, from an underestimation of about 5% (or, equivalently, 25 Wm −2 ) for very clean conditions, to unbiased estimates for AODs 710 about 0.12.This could be expected for this scheme because of its fixed aerosol scattering parameter.No trend is observed with respect to SSA and ASY.

Seasonality
One of the particular benefits of having a method that takes 715 aerosol extinction into account is to consider the impact of the seasonal variability of AOD in surface fluxes.Specifically, if AOD is not considered in the calculation of clear-sky surface irradiance, or if that is done using a fixed AOD value, a seasonal bias may appear in the computed irradiances at the 720 surface, which can become considerably large depending on the simulated region.Figure 5 shows the daily mean relative error of the computed DNI, DIF and GHI (simulated values minus observations, relative to the observations) in the 1-year simulations using the RRTMG model and either the rural or 725 urban aerosol model.The error series combine the predictions at the five validation sites.A 15-day moving average filter is used to underline the bias trend.For GHI, the calculated values with the Dudhia scheme are also shown.The expected observational error region for the surface solar irradi-730 ance observations, roughly estimated as ±5%, is highlighted in yellow.
Figure 5a and b show the DNI and DIF estimates, respectively.Overall, both the rural and urban aerosol mixtures produce unbiased DNI values during the entire simulated year.

735
The small difference between them is due to the different AE values that result from the selection of differing aerosol mixtures.Regarding DIF, the urban aerosol yields a sustained bias around -15%, with no seasonal trend, whereas the bias using the rural aerosol stays within the expected observa-tional error region, also without clear seasonal trend.This proves that the rural aerosol model fits the observations for the present experimental sites.
Figure 5c shows the results for GHI.The values computed with the RRTMG model for the rural aerosol are unbiased throughout the entire simulated year, whereas an assumed urban aerosol mixture would introduce a negative bias of about -2%.However, no seasonal trend is observed in either case.In contrast, the Dudhia scheme shows a clear seasonal trend in its bias: it underestimates by up to 5% in winter, apparently because of its reliance on mean annual conditions.This results in too much scattering in winter and too little in summer, to the point that it cannot reproduce the natural seasonal variability in GHI.

Discussion and conclusions
A parameterization of the aerosol optical properties for the prediction of all three components of short-wave surface solar irradiance in NWP models has been proposed.It has been implemented and verified in the RRTMG SW scheme of the WRF NWP model.The verification has been conducted at five radiometric stations in the contiguous US with nearby or collocated AERONET sites.A previous experiment, using observed aerosol optical properties, has been used here to serve as control experiment.The latter consists of a best-case clear-sky evaluations obtained with two of the WRF shortwave solar radiation schemes when forced with observed aerosol optical properties taken at test AERONET sites.Thus no aerosol optical property is parameterized in the control experiment.In contrast, the aerosol optical parameterization only uses observations of AOD at 0.55 µm, whereas AE, SSA and ASY are parameterized based on the (assumed) predominant type of aerosol and relative humidity.
The approach to parameterize the aerosol optical properties is versatile since the only mandatory parameter is AOD at 0.55 µm, which can be provided either as a fixed value or as a time and space varying field.The remaining aerosol optical parameters, namely, AE, SSA and ASY are parameterized from the known properties of two different bimodal aerosol mixtures, namely rural and urban, dominated by the accumulation mode.The urban aerosol model is basically a more absorbing version of the rural aerosol.However, as for AOD at 0.55 µm, AE, SSA and ASY can also be provided to WRF (starting with version 3.6) as either a fixed value or a time and space varying field.This allows the undertaking of sensitivity studies or the use of external data sources.
The proposed parameterization has been evaluated over a period of one year, which is considerably longer and significant than the 5-day control experiment.Overall, the verification has shown satisfactory results.Regardless of the reference aerosol that is invoked, virtually no difference is found in DNI when evaluated using the AOP parameterization or the control case.In the case of DIF, the selection of the most appropriate reference aerosol is important because it conditions SSA and ASY, which in turn affect DIF.At four of the five experimental sites, the rural aerosol model results in very 795 good agreement with the control experiment.At the remaining site, the SSA values observed at the AERONET station were anomalously low during the five test days.This explains why the urban aerosol model is better there under such exceptional conditions.Its use can be effective to consider the 800 effect of highly absorbing aerosols caused by pollution or smoke from wildfires.Based on the 1-year simulation conducted at the five test sites, it is found that the use of the AOP parameterization to simulate the time variations of AOD at 550 nm contributes to effectively removing seasonal biases 805 in the predicted DNI, DIF and GHI.For the latter irradiance component, the performance of the AOP parameterization has been illustrated by comparing its results against the Dudhia short-wave scheme, which considers aerosol extinction only by assuming a bulk extinction parameter.

810
Arguably, a major limitation of the proposed AOP parameterization might be in its requirement to select one of the only two aerosol models currently offered.This method still makes sense for limited-area models, such as WRF, under the assumption that any significant changes regarding the aerosol 815 mixture occur at spatial scales larger than the domain being simulated with the NWP model.
At this point, the two aerosol mixtures do not allow the simulation of aerosol situations having a dominant coarse mode, such as what is typical with sea salt or desert dust.The 820 inclusion of these other types of aerosol mixture is already in progress.When completed, this will extend the applicability of the method, particularly to arid regions, where higher AODs are the norm, and where considerable solar development is already taking place.Finally, it is worth mentioning 825 that the modelling of aerosol extinction for the prediction of GHI, DNI and DIF based on reference aerosol mixtures, as here presented, is an approach to the general problem of surface irradiance prediction at high spatio-temporal resolution.An alternative approach to obtain irradiance predictions at 830 larger scales and coarser resolutions would be to use climatological aerosol optical properties, if their temporal variability over the region of interest is known to be small.
Fig.1.AOD spectral scale factor interpolated using 4-point Lagrange interpolation for relative humidities from 0% to 99% for each RRTMG spectral band and the rural aerosol model.For the sake of comparison, the results using Eon-weighted and un-weighted spectral scale factors are shown.

Fig. 2 .
Fig.2.Parameterized SSA and ASY parameters for the rural and urban aerosol mixtures and a relative humidity of 80% (thick line).TheShettle and Fenn (1979) spectral values are shown with cross marks.They are interpolated using cubic splines (thin line).The grey region encompasses the variability range of the parameters with different values of relative humidity.

Fig. 4 .
Fig. 4.Error analysis with respect to the variability range of AOD, SSA and ASY parameters for GHI, DNI and DIF that results from the 1-year WRF simulations.(a-c) shows the relative frequency distribution of the observed AOD at 0.55 µm, the observed and parameterized SSA values, and the observed and parameterized ASY values, respectively.(d-l) shows the observed and simulated DNI, DIF and GHI values (upper half of the panels) as well as their relative errors (lower half of the panels) as a function of the observed AOD at 0.55 µm, SSA and ASY values.The expected observational error region for the surface solar irradiance observations, roughly estimated as ±5%, is highlighted in yellow.
this work, a parameterization approach for the AOPs

Table 1 .
Spectral bands distribution in RRTMG.From top to bottom rows, λ's (in nm) are band mean, band minimum and band maximum values, respectively.Note the band numbering does not follow increasing or decreasing wavelength values.The band naming convention follows the RRTMG's definition.

Table 2 .
Ångström exponents for the two spectral bands of the modified Ångström's law, the aerosol mixtures and relative humidity values.Ångström exponents are computed as described in Sect.3.1

Table A1 .
AOD spectral scale factor ρrj for the rural aerosol mixture.

Table A2 .
AOD spectral scale factor ρrj for the urban aerosol mixture.

Table A3 .
Single-scattering albedo for the rural aerosol mixture.