A discussion about the role of shortwave schemes on real WRF-ARW simulations. Two case studies: cloudless and cloudy sky

A wide range of approaches for radiative transfer computations leads to several parameterizations. Differences in these approximations bring about distinct results for the radiative ﬂuxes, even under the same atmospheric conditions. Since the transfer of solar and terrestrial radiation represents the primordial physical process that shapes the atmospheric circulation, these deviations must have an impact on the numerical weather prediction (NWP) model performance. In this paper, an analysis of the role of shortwave schemes on the Weather Research and Forecasting (WRF-ARW) model is presented. The study compares the effect of four parameterizations (Dudhia, New Goddard, CAM and RRTMG) in two cases: i) cloudless and ii) cloudy sky situations for a domain deﬁned over Catalonia (northeast of the Iberian Peninsula). We analyze the direct and the indirect feedback between the dynamical aspects and the physical parameterizations driven by changes in the radiative transfer equation computation. The cumulative effect of these variations are studied through three simulation windows: current day (0-23 h), day-ahead (24-47 h) and two days ahead (48-71 h). These analyses are focused on several NWP model ﬁelds. From the most directly related to shortwave schemes such as global horizontal irradiance or the heating rate proﬁle, to apparently secondary outcomes such as wind speed or cloud composition among others. The differences observed between model runs using different solar parameterizations increase with the simulation horizon, being more important in the cloudy scenario than in the cloudless sky.


Introduction
Since the transfer of terrestrial and solar radiation represents the primordial physical process that shapes the atmospheric circulation, an understanding of numerical weather prediction (NWP) model performance also requires a detailed understanding of the interaction of these processes with the other mechanisms represented by the model.
For most of the weather applications, some modeling elements such as cloud physics, convection or initial and boundary conditions, among others, have a high impact, leaving solar radiation in second place. Furthermore, in a first-order approach (i.e. synoptic scale and mesoscale when short simulation horizons are considered), the diabatic contribution in the Euler equations can be neglected (Holton, 2004).
The development of new applications during the last decade is increasing the interest in shortwave scheme performance (Ruiz-Arias et al., 2013). Specially, in renewable solar energy to model the resource for forecasting as well as for prospecting applications.
A simple, fast and accurate computation of the radiative transfer equation (RTE) becomes an important challenge in NWP models. First, from a mathematical frame due to the complicated integro-differential form of the RTE (Chandrasekhar, 1960;Liou, 1980). Secondly, from a computational point of view, since the RTE in plane-parallel atmospheres involves integrals over the entire spectrum, A. Montornès et al.: A discussion about the role of the shortwave schemes on real WRF-ARW simulations with the solid angle and the vertical coordinate for each grid-point increasing the computational specifications (Lu et al., 2012). Finally, from a physical frame, because radiative parameters may take on a wide range of values increasing the complexity for a good fit (Joseph et al., 1976).
Typically, shortwave and longwave schemes in mesoscale models reduce the RTE to a set of ordinary differential equations in terms of the optical depth, the single-scattering albedo, the asymmetry factor and, in some cases, in terms of the second moment of the scattering phase function. These variables are parameterized as a function of the available fields in the model (e.g. air density, ice crystals or the water vapor mixing ratio among others) and by using look-up tables from experimental sources (e.g. HITRAN data set). The other required but not directly available from the mesoscale model, such as the ozone or the trace gases, is reduced to some pre-computed values defined inside the code or in auxiliary files. Finally, the continuum spectrum must be reduced to a few bands since line-by-line computations are not feasible with the current hardware resources (Lu et al., 2012). In recent years, some attempts to reduce the computational time have been performed by using Graphic Processing Unit, GPU-based architectures instead of the current Central Processing Unit, CPU (Michalakes and Vachharajan, 2008;Mielikainen et al., 2012). Although this hardware improvement will allow for a deep treatment of the radiative processes in the future, this type of computation is often hindered by the adaptation of mesoscale models to the new hardware architecture and coding (Ruestsch et al., 2010).
As a result, a wide range of approaches and strategies for radiative transfer computations are available, leading to several parameterizations. Each scheme induces different absorption, reflection and transmission outcomes producing departures on heating rate profile and surface irradiance values. Both outputs play an important role in many physical processes occurring in the Earth's atmosphere. Heat generated by solar radiation absorption is a first order contribution to the diabatic heating in the free atmosphere with a wide range of values due to the atmospheric composition. In the stratosphere, the main molecular absorber is ozone producing a heating rate ranging from 10 to 30 K day −1 , being an important mechanism in maintaining the thermal structure of that layer (Ramanathan and Dickinson, 1979). Contrarily, the troposphere shows different factors absorbing shortwave radiation as air molecules, hydrometeors and aerosols. In a cloudless and clear-sky situation (i.e. without aerosols), water vapor is the main absorber leading to a  Figure 3. The study covers different simulation horizons: 0-23 h (current day), 24-47 h (day-ahead) and 48-71 h (two days ahead). Shaded colors indicate the data used for the analysis of this paper.
heating rate that varies from 0.5 to 2 K day −1 . Cloud layers increase the solar radiation absorption depending strongly on the solar zenith angle, cloud composition (i.e. water or ice) and surface albedo (Chou et al., 1995). Aerosols have a high impact on radiative transfer showing a high variability in space, time and composition. The role of the aerosols has been analyzed for years due to the interest of its impact on the climate system. During the last decade, different authors, such as Ruiz-Arias et al. (2013) have analyzed this element due to the growth of the renewable energy applications in which the industry requires high accuracy.
In addition, solar surface irradiance warms the soil and is transformed into sensible and latent heat leading to many processes occurring in the Planetary Boundary Layer (PBL) as vertical transport of heat, mass and momentum as turbulence motion (Stull, 1988).
Hence, two identical simulations using different solar schemes undergo changes in the diabatic term of the thermodynamic energy equation (i.e. heating rate contribution) and in the surface energy balance (i.e. surface irradiance contribution) that are propagated to the other meteorological fields over time following high non-linear relationships given by the governing equations (Section 4).
For many years, the feedback between solar energy and the climate system has been under study by many authors a Budyko (1969) or Wild et al. (2005), among others. However, the effects on mesoscale simulations have not been treated as a relevant issue until now. This paper offers a discussion about the impact and the role of the shortwave parameterizations in mesoscale simulations comparing four available schemes in version 3.5 of the WRF-ARW model: Dudhia, New Goddard, NCAR Community Atmosphere Model (CAM) and Rapid Radiative Transfer Model for GCMs (RRTMG). As a high number of degrees of freedom are introduced by a study of this magnitude, we have reduced the analysis to two well-defined scenarios: i) simple cloudless sky and ii) complex cloudy sky. The proposed experiments consider three telescopic nests with two-way nesting centered over the central coast of Catalonia (northeast of the Iberian Peninsula).
The discussion of the results in Section 3 covers several fields, from the most direct variables to the solar parameterizations such as shortwave heating rate and global horizontal irradiance (GHI), to fields apparently as far removed as wind speed profile, for instance. Effects over time produced by these variations are analyzed throughout three simulation windows: current day (0-23 h), day-ahead (24-47 h) and two days ahead (48-71 h) with respect to the initialization time. Finally, in Section 4, we propose a conceptual frame explaining the relationships between the shortwave parameterization and the other processes represented by the NWP model.

Model configuration
We use the WRF-ARW model for a dynamical downscaling consisting of three telescopic nests with two-way nesting consideration. All domains are centered in Barcelona and defined in a Lambert Conformal Conic geographical projection tangent to the standard latitude. The innermost domain has a horizontal resolution of 3 km with 100 points in the north-south and east-west directions. The relationship with the successive domains is 1:3 with respect to the previous one (i.e. 27, 9 and 3 km) with a buffer composed of ∼50 grid points in each direction ( Figure 1). In the vertical resolution, we consider 70 vertical levels fixed automatically with the top of the model (TOM) located at 10 hPa. We establish the default data sets provided by the model as static data (i.e. terrain, albedo, etc.).
Initialization and boundary conditions are set using the ERA-Interim reanalysis at 0.7 • x 0.7 • (Dee et al., 2011) whereas sea surface temperature SST are updated using the NOAA Optimum Interpolation (OI) SST V2 data set (Reynolds et al., 2002).
The WRF-ARW model includes seven shortwave schemes. However, for this paper only four parameterizations are considered: Dudhia (Dudhia, 1989), New Goddard (Chou and Suarez, 1999;Chou et al., 2001), CAM (Collins et al., 2004) and RRTMG (Iacono et al., 2008). The Goddard scheme (Chou and Suarez, 1994) is not examined since it is an old-fashioned version of the New Goddard with some relevant bugs and inconsistencies. The Geophysical Fluid Dynamics Laboratory, GFDL (Fels and Schwarzkopf, 1981) and Fu-Liou-Gu, FLG (Gu et al., 2011) are not analyzed to reduce the complexity of the discussion and because they do not add more knowledge beyond the particularities of different parameterizations. Nevertheless, the selection of these schemes is not arbitrary. The GFDL package was a solar scheme adapted from the ETA Model developed in the 1980s and, nowadays, it is deprecated in the WRF-ARW model because it has option 99 in the namelist.input file that is reserved for old parameterizations. In contrast, FLG is the newest scheme (since version 3.4) with few publications A. Montornès et al.: A discussion about the role of the shortwave schemes on real WRF-ARW simulations related to the WRF-ARW model.
Dudhia is the simplest shortwave scheme in the model. RTE is reduced to a broadband-integrated solar flux. Air scattering is tuned by an external parameter denoted as swrad scat defined in the namelist.input file of the model. Water vapor absorption is parameterized following an empirical expression from Lacis and Hansen (1974). Cloud albedo and absorption are computed using look-up tables from Stephens (1978) as a function of the cosine of the solar zenith angle and the liquid water path. Since ozone and trace gases are not considered, this scheme must be avoided in simulations with a TOM of less than 50 hPa. Nevertheless, it is used in the proposed experiments (i.e. TOM at 10 hPa) to discuss the ozone effect in the lower stratosphere.
New Goddard, CAM and RRTMG show common approaches in the radiative transfer solver. All use the two-stream approach reducing the RTE into two ordinary differential equations. In New Goddard and CAM, these equations are solved using the Delta-Eddington approximation while RRTMG uses the Practical Improved Flux Method, PIFM (Zdunkowski et al., 1980). Vertically, upward and downward fluxes are integrated using the Adding method (Liou, 1980). Differences in these schemes arise from the spectrum and optical property treatment as well as ozone, trace gases and aerosol composition.
New Goddard splits the solar spectrum into eleven spectral bands: seven in the ultraviolet (UV) band, one in the photosynthetic active region (PAR) band and three in the near-infrared (near-IR) band. The incorporation of non-gray gaseous absorption in multiple-scattering atmospheres is based on the K-distribution method (Liou, 1980) dividing each water vapor band into 10 sub-intervals. Cloud liquid, ice and rain particles are allowed to coexist at a layer parameterized in terms of the effective particle size. Clouds are clustered as low, middle and high, and the cloud overlapping is treated using the maximum-random approach (Chou and Suarez, 1999). The ozone profile is reduced to five data sets. Aerosols are not considered. The CAM scheme splits the shortwave spectrum into nineteen bands (seven in the UV band, three in the PAR band and nine in the near-IR band). Cloud scattering and albedo are parameterized following (Slingo, 1989) in which optical properties of the cloud droplets are represented in terms of the liquid water path and the effective radius. Several ozone profiles are available with a latitudinal resolution of 2.28 • for each month of the year. Background aerosol profiles are considered. Finally, RRTMG divides the spectrum into fourteen bands (three in the UV band, two in the PAR band and nine in the near-IR band). Water vapor absorption is based on the Correlated K-distribution, CKD method (Liou, 1980). This scheme interacts with cloud (liquid and ice) fractions A. Montornès et al.: A discussion about the role of the shortwave schemes on real WRF-ARW simulations using the Monte Carlo Independent Column Approximation, MCICA (Pincus et al., 2003). Ozone absorption is based on a single profile calibrated in mid-latitudes. The aerosol treatment can be coupled with WRF-CHEM (Grell et al., 2005). Since version 3.5, this scheme has been able to work with ozone and aerosol profiles available in the CAM scheme.
The latest version of the WRF-ARW model available when preparing this paper was the 3.6.1 (available since summer 2014). From version 3.6 a new module for the aerosol treatment using ancillary data for New Goddard and RRTMG was included. This kind of data can be provided using the namelist.input file to set the aerosol properties (assumed to be constant throughout the simulation) or using a gridded data set provided by the user as an auxiliary file (Ruiz-Arias et al., 2014).
Although aerosols have a non-negligible effect on the radiative transfer in the Earth's atmosphere, they are omitted from this paper for three reasons: i) they are an external factor to the NWP model that can increase the complexity in the interpretation of the results, ii) this kind of experiment is only available for two schemes and iii) these new options need a previous analysis in order to tune the aerosol representation for a good fit.  The set of simulations leading to the following results and conclusions have been performed using version 3.5 of the model to avoid the effect of new bugs or unknown issues in the newest versions.
In all cases, we call the radiative schemes (including the longwave) every 3 minutes. The other parameterizations are not relevant for the purposes of the current study and they are briefly described in Table 1. These schemes are chosen following the configuration detailed in Mercader et al. (2010) and updated to the version of the model used in this case.

Cases studies and analysis
In trying to reduce the complexity of the problem, our analyses are simplified in two case studies: a cloudless sky scenario and a complex situation with cloud covering the innermost domain ( Figure 2). For each case, we consider three simulation windows: from 0 to 23 hours (current day), from 24 to 47 hours (day-ahead) and from 48 to 71 hours (two days ahead) as illustrated in Figure 3. For the current day simulations, the horizon from 0 h to 6 h is not considered in the discussion to avoid the spin-up of the model. All simulations start at 00 UTC (night in local time) avoiding the previous effects from the shortwave scheme in the current day.
For the cloudless case, we add the condition that the two previous days must show similar atmospheric situations to avoid any previous interaction with clouds in the day-ahead and two day-ahead simulations. The cloudless and cloudy scenarios are represented by 2011-10-12 and 2012-09-29, respectively ( Figure 2).
The analysis of the results are centered in a single node of the innermost domain located in Barcelona. For this discussion, we consider the model output with a 15 minute resolution.

Representation of the results
Based on this set of simulations, we will analyze the role of the shortwave schemes on real simulations and their impact on several fields in Section 3.
In general, the results are presented in terms of the departures from the mean value of all the simulations. Then, given one field f (e.g. temperature), the departure f for the scheme s at the time t is evaluated as where F is the average value of all schemes.
One of the aspects that we want to present in the next section is the variation of these departures as a function of the simulation horizon. For this reason, outcomes from different initializations are analyzed separately in Equation 1. In other words, if d is an index indicating the simulation horizon (i.e. current day, day-ahead or two day-ahead), Equation 1 can be rewritten as (2)

Results and discussion
The main research interests include a discussion about the impact of different solar transfer approaches on mesoscale simulations. The analysis presented is divided into two parts: i) an analysis of the model fields response (Section 3.1) and ii) study of the effects through the simulation horizon (Section 3.2).
The first part is focused on the current day outcomes because they simplify the discussion for two reasons. On the one hand, these simulations start at night (i.e. without any previous call of the solar codes) and hence, departures appear with sunrise. On the other hand, since the radiative processes have a cumulative impact (e.g. diabatic term on the energy equation), this simulation horizon offers a simpler analysis of the model response, being necessary for building a conceptual framework in order to explain the response of the model.
In the second part, we use day-ahead and two dayahead simulations to discuss the variations produced due to those cumulative effects through the execution horizon. In both parts, we will present a diagnosis of the response of different model fields from the most directly variables to the solar radiative transfer, such as heating rate or surface irradiance, to other meteorological fields, such as wind speed, temperature, surface fluxes or hydrometeors among others.

Response of the model fields (current day)
3.1.1 Effects on heating rate and surface global irradiance As discussed in Section 2.1, each parameterization assumes a wide range of approximations to model the solar transfer. Consequently, shortwave schemes lead to deviations in the radiative outputs even using the same atmospheric conditions. These changes are linked to: i) departures in the heating rate and ii) departures in GHI. Although recent WRF-ARW releases incorporate the direct and diffuse irradiance components, these fields have been avoided because they do not interact with other model elements.
The first kind are produced by changes in the absorbed energy. The magnitude of the deviations is different in the stratosphere (Figure 4) than in the troposphere (Figures 5 and  6). In the upper layer, the main absorber is the ozone. The differences in the modeled stratospheric heating rate using different schemes are observed in Figure 4. Without a specification of the ozone profile, Dudhia shows near-zero heating rate values. Consequently, the stratosphere experiences an important radiative cooling as noted in Figure 7. In contrast, the simulations using more sophisticated schemes are capable of maintaining the thermal structure of the stratosphere since they account for the ozone profile. These parameterizations undergo the highest heating rates at 10 hPa at noon (Figure 4). New Goddard shows the greatest heating, reaching a departure from the mean value of all simulations greater than 1.5 K day −1 . RRMTG and CAM experience similar heating rates with departures of around 0.5 and 1.5 K day −1 from the average, being more significant in the first scheme.
The troposphere requires a more complex analysis as a function of the cloud presence. In a cloudless and clearsky (i.e. without aerosols), the main shortwave absorber is the water vapor occurring in the near infrared bands. Since these solar bands have low available energy, the effect over the heating rate is less than 2 K day −1 with a similar distribution as the water vapor mixing ratio. As observed in Figure 5, the Dudhia scheme experiences the lowest heating rates with departures from the average that reach -1 K day −1 near the surface. In contrast, the highest heating rate values are undergone in CAM. The differences in the results between New Goddard, RRTMG and CAM are negligible with variations of less than 0.5 K day −1 . As a consequence, departures in the temperature profiles in the middle and in the upper troposphere are negligible (Figure 7).
In the cloudy scenario, two cloud regions can be observed ( Figure 8): i) discontinuous liquid clouds located below 600 hPa and ii) a thick ice cloud layer located around 400 hPa with the greatest thickness at 17:30 UTC. As discussed in the next paragraphs, hydrometeors show high sensitivity with the shortwave schemes that introduce more differences into the comparison of the heating rates ( Figure 6). For example, significant differences in the ice crystal mixing ratio are observed at 12:00 UTC comparing Dudhia and CAM simulations (Figure 8) with large departures in the radiative output ( Figure 6). For ice clouds, New Goddard shows the highest heating rates undergoing a difference with respect to the average of between 1 and 3 K day −1 with a maximum deviation reaching ∼8 K day −1 , 2 K day −1 in Dudhia, -3.5 K day −1 in CAM and 3 K day −1 in RRTMG when the cloud is thick. In contrast, for liquid clouds the highest heating rate is observed in RRTMG simulation.
Below clouds, the lowest heating rates are observed in Dudhia, which is the simplest parameterization (Section 2.1). In order to explain this behavior, we need to look in depth at the theoretical discussion about this scheme. As proposed by Dudhia (1989), this parameterization only considers the downward component of the solar flux from the top of the model to the surface. At each layer dz, the flux reaching the next layer is reduced assuming four contributions: clear-sky scattering dS cs , water vapor absorption dA wv and cloud absorption dA cld and reflection dR cld . These magnitudes are affected by the solar zenith angle µ 0 that increases the path length. Therefore, if µ 0 F is the monochromatic flux at top of the atmosphere (TOA), the downward component of solar flux F ↓ at the level z is evaluated as or normalizing with respect to µ 0 F , where lower case letters represent the normalized value of the respective capital letters. By construction, ds cs , da wv , da cld and dr cld are defined as positive values. Thus, they reduce the solar flux in Equation 4. Consequently, Dudhia does not have a representation of the multiscattering processes given by the source function in the RTE and the scattered radiation by clouds does not have a positive contribution to the downward flux below clouds.
Furthermore, the solar heating rate in this scheme can not be determined as the divergence of the flux because the upward solar flux is not computed. Therefore, the heating rate in a layer bounded by z 1 and z 2 is parameterized considering the contribution of water vapor and clouds as ∂T ∂t = z2 z1 ρc p (dA wv + dA cld + dR cld ) dz.
being ρ the dry air density of the layer and c p the specific heat. Below thick clouds, the solar beam is significantly reduced and water vapor absorption is the only contribution to Equation 5. Contrarily, the other schemes explicitly solve the upward and downward fluxes using the RTE and they can determine the divergence of the flux required to compute the physical definition of the heating rate, with F being the difference net flux (i.e. the difference between the downward and the upward components). The downward flux is composed of direct and diffuse components. Below thick clouds, the first one tends to be small and the main contribution is given by the diffuse radiation. This radiation is affected by several absorption processes, such as the water vapor and minor trace gases (i.e. the oxygen and the carbon dioxide) which increase the heating rate with respect to the Dudhia's case. The differences with the other schemes are around 0.5 and 1 K day −1 . New Goddard, CAM and RRTMG show the lowest deviations because they explicitly consider the diffuse flux contribution. High gradients are observed below and above clouds due to changes in the evolution and location of the hydrometeors between simulations.
The GHI (Figure 9) is strongly affected by cloud presence. In a cloudless sky, similar results are observed for all simulations. Dudhia leads to the most opaque atmosphere. Since low heating rates were experienced in Dudhia, this result may be a consequence of the tuning scattering parameter swrad scat, that attenuates the flux without any contribution to the heating rate. Conversely, the most transparent atmosphere is developed by the New Goddard scheme.
On the other hand, the results in the cloudy scenario require a more complex discussion. Dudhia and RRTMG experience the most opaque atmospheres. Results for the Dudhia simulation are linked to the non-consideration of the diffuse radiation. In contrast, New Goddard and CAM compute the greatest GHI values due to the variations in the cloud composition (different than in RRTMG) and due to the diffuse flux contribution.
The aforementioned discussion gains in strength by expressing the results in terms of the surface daily irradiation. In the cloudless scenario, all schemes show values from 15 MJ m −2 and 17 MJ m −2 (Table 2) with Dudhia being the case with the lowest energy amount and New Goddard the shortwave scheme that introduces the highest energy into the system. In the cloudy case, surface daily irradiation values range from 4.5 MJ m −2 to 7 MJ m −2 (Table 2). RRTMG and Dudhia are the simulations with the lowest energy input.

Effects on other meteorological variables
Changes in GHI and heating rate profiles lead to indirect effects over the atmospheric available energy and, as a consequence, to the other weather fields.
The GHI is an input for the land surface model (LSM). Taking into account different grid-point information (e.g. soil layers, roughness, vegetation, or snow among others) this parameterization transforms the solar radiation into other kinds of energy such as latent LH, sensible SH and A. Montornès et al.: A discussion about the role of the shortwave schemes on real WRF-ARW simulations Figure 16. Cloud evolution during the day-ahead for the Dudhia and CAM simulations. ground GH heats.
The analysis of Figure 10 reveals a high relationship between the SH and the shortwave parameterization. In a clear sky, Dudhia (New Goddard), that shows the most (least) opaque atmosphere, the lowest (highest) values of SH are produced. At noon, differences between Dudhia and New Goddard reach 25 W m −2 . As in the previous section, RRTMG and CAM show intermediate cases with small departures from the mean value. On the other hand, in the cloudy scenario differences become greater. The largest deviations from the average are observed around 10:45 UTC in New Goddard and CAM with a positive maximum of greater than 90 W m −2 , while RRTMG reaches a negative value of near -100 W m −2 .
The LH flux undergoes negligible differences between simulations in the cloudless case. For the cloudy scenario, observed departures become significant (Figure 11) with the largest values achieved at noon (around 2 W m −2 ).
Additionally, departures in the shortwave heating rate introduce changes in the diabatic term of the energy equation. These variations are propagated with high nonlinear relationships by the model dynamics. Departures are observed over all predicted and diagnosed fields with higher values in the cloudy scenario. For example, the wind speed profile in the cloudy case ( Figure 12) undergoes departures that reach ±2 m s −1 or even greater in the current day. The planetary boundary layer (PBL) is driven, in part, by LH and SH. Therefore, changes in the surface energy balance lead to divergences in the PBL output (i.e. PBL height as well as surface fields). In the cloudless scenario, variations in PBL height (Figure 13), temperature at 2 m ( Figure 14) and wind speed at 10 m ( Figure 15) are negligible. The greatest departures in temperature are observed for Dudhia in which differences were between 0.25-0.5 K lower than for the average. For wind speed, the lowest values are experienced during the day, becoming greater than the average after sunset for the Dudhia simulation.
In the cloudy scenario, the lowest deviations in the temperature at 2 m are observed between New Goddard and CAM. The other simulations show differences of greater than 0.5 K with some peaks 1 K larger. The deviations in the wind speed at 10 m are low with values of less than 1 m s −1 . Changes experienced in LH and SH for the cloudy scenario induce large PBL height variations reaching departures of more than 250 m around 10:45 UTC.

Cumulative effects (day-ahead and two days ahead)
In the stratosphere, the ozone layer representation has a relevant impact. The radiative cooling observed in Dudhia (Section 3.1.1) has a cumulative effect that becomes critical for simulation horizons of longer than 24 h (Figure 7). During the day-ahead, a cooling of around 5 K with respect to the current day is experienced at 10 hPa, whereas for the two days ahead it increases to 7 K.
In the troposphere, the analyses of the day-ahead and the two days ahead results reveal an amplification of the differences between simulations for all compared fields and scenarios. These variations are more important in the cloudy scenario than in the cloudless one.
In the cloudless case, differences between parameterizations are negligible. Generally, all schemes show similar values of surface daily irradiation with respect to the current day case ranging from 15 to 17 MJ m −2 ( Table 2). As observed in Figure 5 and Figure 9, heating rate and GHI experience low variations between schemes. In the result analyzed in this study, a slight tendency is observed for all shortwave schemes to become more opaque reducing the GHI (Table 2). Analyzed surface fields (Figures 14 and 15) undergo small changes between New Goddard, RRTMG and CAM simulations. Conversely, Dudhia shows an increment of the differences throughout the simulation window. The most significant departures are observed in the wind speed at 10 m reaching variations of near 1 m s −1 at some moments of the day-ahead and even a different evolution in some hours of the two days ahead.
For the cloudy scenario, the most significant differences between simulations are in the hydrometeors and, as a consequence, in the radiative output. Following the previous example detailed in Figure 8 for the current day simulation, Figure 16 shows the same comparison for the day-ahead. While in the current day, variations were observed in magnitude, in the day-ahead deviations were experienced in magnitude as well as in distribution. The largest differences were observed in the two days ahead comparison but are not presented in this paper.
Variations in the radiative results lead to large departures to the other fields as we will discuss in detail in Section 4. These deviations are greater in the two days ahead simulation window than in the day-ahead. In the day-ahead, all parameterizations tend to produce lower surface daily irradiation values (Table 2) within the range from 2.5 to 4.4 MJ m −2 . In contrast, higher surface daily irradiation values are observed in the two days ahead window ranging from 7 to 13.2 MJ m −2 due to a lower cloud cover over the analyzed grid point during the entire day.
Changes in the GHI cause departures in SH ( Figure 10) with departures reaching ±40 W m −2 in the day-ahead, drifting to more than ±100 W m −2 in the two days ahead window. These changes produce irregular deviations in PBL height (Figure 13), which are more significant in RRTMG simulation during the two days ahead.
For the surface fields, similar patterns are observed with an increment of the departures in the two days ahead. For temperature at 2 m ( Figure 14) these changes lead to deviations of 1 K in both windows, with New Goddard being the coldest while CAM and RRTMG are the warmest ones. Wind speed at 10 m ( Figure 15) differences are around ±2 m s −1 or less with peaks reaching +6 m s −1 .
High departures in wind speed profiles are observed in the cloudy scenario with deviations reaching ±10 m s −1 (Figure 17) in some cases.

Feedback between solar transfer schemes and other NWP model components
NWP models are composed of several elements each with their scientific considerations with the resulting model being as the sum of these components and their interactions (Dudhia, 2014). The description of the results in the current day (Sections 3.1 and 3.2) allows us to propose a conceptual frame to figure out the interaction of the shortwave parameterization with the dynamical and numerical aspects and the other physical parameterizations.
As detailed in the model technical guide (Skamarock et al., 2008), the primitive or Euler equations governing the WRF-ARW model are given by ∂ t θ + (∇Vθ) = F θ (10) Expressions 7, 8 and 9 are the momentum equations describing the fluid movement, Equation 10 is the energy conservation, equations 11 and 13 describe the mass conservation for dry air and water vapor, and expression 12 is the geopotential equation.
A further description of these equations and variables is detailed in the technical guide. However, for the following discussion, the interpretation of the left-hand and righthand parts of the equations becomes relevant. The left-hand terms are related to the dynamical and numerical aspects of the model.
In contrast, the right-hand terms correspond to the unresolved physical processes that need to be parameterized. F U , F V and F W include the friction with the surface as well as the horizontal and the vertical transport of momentum by turbulent eddies. In mesoscale simulations, these mechanisms are modeled by the atmospheric surface layer and PBL schemes. F θ is the diabatic heating term in the energy conservation (Equation 10) taking into account the condensation (microphysics scheme), the convection (cumulus scheme), the radiative heating rate (shortwave and longwave schemes) and the heat transport in the boundary layer (PBL scheme). Finally, F Qm corresponds to the moisture variation due to the precipitation rate, either stratiform (microphysics scheme) or convective (cumulus scheme), due to the evaporation at surface (i.e. LSM) and due to the phase variation during the fall of the hydrometeors (microphysics scheme).
Focusing the discussion on the shortwave scheme, Figure 18 is presented. In the solar transfer parameterizations, optical properties are computed as a function of the air density ρ given by the primitive equations and as a function of the water mixing ratio Qv, liquid clouds Qc, ice clouds Qi and other hydrometeors such as rain Qr, snow Qs and graupel Qg given by the microphysics scheme. These fields are used to compute the radiative variables required in the RTE. As a result, upward R ↑ s and downward R ↓ s fluxes are computed at each layer. The divergence of the flux determines the heating rate (i.e. F θ ) as where g is the gravity, c p is the heat capacity at constant pressure and T and p are the temperature and pressure, respectively. As observed in Figures 4, 5 and 6, differences in modeling the RTE produce important variations on the heating rate profile. These departures are given by different ozone representations in the stratosphere (Montornès et al., 2015) and by different atmospheric compositions in the troposphere (i.e. water vapor and clouds, since aerosols are not considered). The vertical integration of the fluxes throughout the entire column leads to the irradiance reaching the ground (Figure 9). Solar energy at the surface interacts in the LSM with surface properties such as the vegetation or the soil use among others. As a result, the shortwave radiation is transformed into LH (Figure 11), SH ( Figure 10) and GH that are important to force the diurnal cycle.
Both LH and SH drive the atmosphere surface layer and the PBL schemes fixing the surface temperature (Figure 14), wind speed ( Figure 15) and water vapor mixing ratio.
Surface temperature defines the radiative emission, R ↑ l in the longwave spectrum following Plank's Law for black-body radiation. This value is the main input for the longwave scheme that works in a similar way to the shortwave schemes described previously. The divergence of the longwave flux is transformed into heat producing variations in the diabatic term of the energy equation F θ . Longwave downward flux R ↓ l reaching the surface interacts with the LSM producing changes in the surface energy budget.
Hydrometeors (Figures 8 and 16) are computed in the microphysics scheme mainly as a function of the air temperature and Q v provided by the Euler equations.
To conclude, given a set of fixed atmospheric conditions, a variation in the solar radiative approximations leads to changes in the heating rate profile and surface energy budget at each grid-point. The first one adds/removes energy in the Euler equations with an effect on the next time step. Hence, that energy is propagated throughout prognostic and diagnostic fields with high non-linear relationships leading to variations in the other physical schemes such as microphysics. The second one leads to deviations in the physical processes that occur at the surface. In a parallel way, these processes interact with the Euler equations as previously described. Consequently, in the next radiative call, the atmospheric conditions for the shortwave parameterizations are different. Thus the variations in the results through the simulation horizon are broadened. These variations are more important when clouds are considered on the simulation domain (e.g. Figures 12 and 17). On the one hand, due to their impact on radiative absorption and surface irradiance. On the other hand, because of the wide range of approaches and methods for the cloud treatment in the solar schemes.

Conclusions
The impact of four solar transfer parameterizations (Dudhia, New Goddard, CAM and RRTMG) in real WRF-ARW simulations has been analyzed with two study cases considering the impact with and without clouds over the simulation domain. Moreover, the effect through the simulation horizon has been studied in three model windows (0-23 h, 24-47 h and 48-71 h).
There is a wide range of approaches for solar transfer computation that leads to different results on heating rate and GHI (i.e. available energy in the simulation) even under the same initial and boundary conditions. These variables contribute to the dynamical aspects of the model through the diabatic term in the thermodynamic energy equation and to the other physical parameterizations. In consequence, variations in RTE computation introduce changes in the simulation results.
In the stratosphere, ozone profiles available in the shortwave schemes are critical in maintaining the thermal structure of the layer, as was illustrated in Figure 7.
The impact of the solar schemes at the troposphere shows differences between the cloudless situation and the cloudy scenario. In a situation without clouds and aerosols, the most important absorber is the water vapor, the absorption of which occurs mainly in the near-IR bands with low available energy. Therefore, differences in the RTE treatment introduce small variations in the other parts of the mesoscale model arising as a negligible factor.
In contrast, the interaction between clouds and solar radiation has a relevant role in the simulations. Modifications in modeling these interactions lead to changes in the heating rate with an effect on cloud composition for the successive radiative calls. A departure in the cloud composition produces more departures in the heating rate and the GHI. As a result, these deviations are propagated to the other parts of the model producing changes in the results.
In conclusion, this work clearly shows how shortwave radiation schemes play a significant role in the mesoscale NWP simulations that must not be omitted. A deeper analysis covering all the grid-points may be valuable as a future work.