Decoupling the effects of clear atmosphere and clouds to simplify calculations of the broadband solar irradiance at ground level

. In the case of inﬁnite plane-parallel single- and double-layered cloud, the solar irradiance at ground level computed by a radiative transfer model can be approximated by the product of the irradiance under clear atmosphere and a modiﬁcation factor due to cloud properties and ground albedo only. Changes in clear-atmosphere properties have negligible effect on the latter so that both terms can be cal-culated independently. The error made in using this approximation depends mostly on the solar zenith angle, the ground albedo and the cloud optical depth. In most cases, the maximum errors (95th percentile) on global and direct surface irradiances are less than 15 W m − 2 and less than 2–5 % in relative value. These values are similar to those recommended by the World Meteorological Organization for high-quality measurements of the solar irradiance.


Introduction
Solar radiation drives weather and climate and takes part in the control of atmospheric chemistry. The surface solar irradiance (SSI) is defined as the power received from the sun on a horizontal surface at ground level. Of concern here is the SSI integrated over the whole spectrum, i.e. between 0.3 and 4 µm, called total or broadband SSI.
Numerical radiative transfer models (RTMs) simulate the propagation of radiation through the atmosphere and are used to calculate the SSI for given atmospheric and surface conditions. RTMs are demanding regarding computer time and in this respect are not appropriate in cases where operational computations of the SSI are performed such as at Deutscher Wetterdienst (Mueller et al., 2009), the Royal Netherlands Meteorological Institute (KNMI) (Deneke et al., 2008;Greuell et al., 2013), the MINES ParisTech  or prepared within the MACC European project (Granier et al., 2010). Several solutions have been proposed in order to speed up calculations of the SSI, such as abaci -also known as look-up tables (LUTs; Deneke et al., 2008;Huang et al., 2011;Mueller et al., 2009;Schulz et al., 2009).
The present work contributes to the research of fast calculations of the SSI under all sky conditions. It does not propose a new model but an approximation that can be adopted by models for calculations of the SSI. More precisely, it examines whether in the case of infinite plane-parallel singleand double-layered cloud, the SSI computed by an RTM can be approximated by the product of the SSI under clearsky and a modification factor due to cloud properties and ground albedo only. If this approximation were accurate enough, i.e. if the modification factor did not significantly change with clear-atmosphere properties, it would be possible to construct two independent models, possibly LUTbased models -for example, one for clear-sky conditions and the other for cloudy conditions. Recently, for example, Huang et al. (2011) used such an approximation with a very limited justification. This Technical Note aims at holding this Published by Copernicus Publications on behalf of the European Geosciences Union. justification by (1) exploring the influence of the properties of the clear atmosphere on the SSI in cloudy atmosphere, (2) proposing a general equation that decouples the effects of the clear atmosphere from those due to the clouds, and (3) computing the errors made with this approximation.

Objective
Let G denote the SSI for any sky. G is the sum of the beam component B of the SSI -also known as the direct component -and of the diffuse component D, both received on a horizontal surface. In the present article, following the RTM way of doing, B does not comprise the circumsolar radiation. Let G c , B c and D c denote the same quantities but for clear sky. The ratios K c and K cb are called clear-sky indices (Beyer et al., 1996): (1) K c is also called cloud modification factor in studies on UV or photosynthetically active radiation (Calbo et al., 2005;den Outer et al., 2010).
The indices K c and K cb concentrate on the cloud influence on the downwelling radiation and are expected to change with clear-atmosphere properties P c since the clouds and other atmospheric constituents are mixed up in the atmosphere. Equation (1) can be expanded: where θ S is the solar zenith angle, ρ g the ground albedo, and P c is a set of seven variables governing the optical state of the atmosphere in clear sky: (i) total column contents in ozone and (ii) water vapour; (iii) elevation of the ground above mean sea level; (iv) vertical profile of temperature, pressure, density, and volume mixing ratio for gases as a function of altitude; (v) aerosol optical depth at 550 nm; (vi) Ångström coefficient; and (vii) aerosol type. P cloud is a set of variables governing the optical state of the cloudy atmosphere: (i) cloud optical depth (τ c ), (ii) cloud phase, (iii) cloud liquid water content, (iv) droplet effective radius, and (v) the vertical position of the cloud. The objective of this article is to quantify the error made in decoupling the effects of the clear atmosphere from those due to the clouds in cloudy sky, i.e. if changes in P c are neglected in K c , respectively K cb in Eq. (2). This is equivalent to say that the first derivative ∂K c / ∂P c , respectively ∂K cb / ∂P c , is close to 0. In that case, Eq. (2) may be replaced by the following approximation: where P c0 is an arbitrarily chosen but typical set P c . The objective is now to quantify the error made when using Eq. (3) instead of Eq. (2).

Method
The methodology used for assessing Eq.
(3) is of statistical nature. For a given condition related to the position of the sun, the ground albedo and the clouds (θ S , ρ g , P cloud ), several sets of clear-sky properties P c are randomly built. Each quadruple (θ S , ρ g , P cloud , P c ) is input to an RTM to compute G, B, D, K c and K cb . The variances of the K c and K cb series are then computed. The lower the variance, the lower the changes in K c or K cb with respect to the changes in P c and the more accurate the approximation given by the Eq.
(3). The errors made on the retrieved G and B when using Eq. (3) are quantified. The RTM libRadtran version 1.7 (Mayer and Kylling, 2005) is used with the DISORT (discrete ordinate technique) algorithm (Stamnes et al., 1988) to solve the radiative transfer equation. libRadtran needs input data of the atmosphere and surface properties. When not provided, data are replaced by standard assumptions. Atmosphere and clouds are assumed to be infinite plane-parallel. Table 1 reports the range of the 10 values taken respectively by θ S and ρ g . For computational reasons, θ S is set to 0.01 • , respectively 89 • in place of 0 • , respectively 90 • .
Cloud properties input to libRadtran are τ c , phase (water or ice clouds), heights of the base and top of cloud, the cloud liquid content and effective radius of the droplets. Default values in libRadtran for cloud liquid content and droplet effective radius are used: 1.0 g m −3 and 10 µm for water cloud, and 0.005 g m −3 and 20 µm for ice cloud. In a preliminary study (Oumbe, 2009, Fig. 4.6, p. 53), the influence of the changes in effective radius, from 3 to 50 µm, was found negligible for ice clouds. For water clouds, the smaller the radius, the greater the influence, though this influence is still negligible with respect to other variables.
The cloud properties are linked together. Table 2 presents the typical height of the base of cloud, geometrical thickness, and τ c for the different cloud types and is established after Liou (1976) and Rossow and Schiffer (1999).
A total of 10 values of τ c are selected in this study for water clouds and 10 others for ice clouds (Table 3, left column). Ranges of τ c are related to types of clouds to reproduce realistic conditions. Each τ c defines a series of seven couples cloud base height thickness for water clouds and three for ice clouds (Table 3).
According to Tselioudis et al. (1992), 58 % of the clouds are single layered and 28 % are double layered. The results presented hereafter are for single layer; the case of doublelayered clouds is briefly discussed at the end of Sect. 4 as results are similar in both cases. For a given cloud phase, there are 1000 (10×10×10) possible combinations of θ S , ρ g (Table 1) and τ c ( Table 3). The selection of a given τ c leads to the additional selection of a series of cloud base heights and thicknesses as shown in Table 3, i.e. seven for water clouds and three for ice clouds. At that stage, there are 7000 triplets (θ S , ρ g , P cloud ) for water clouds and 3000 for ice clouds.
Each triplet (θ S , ρ g , P cloud ) gives birth to 20 quadruples (θ S , ρ g , P cloud , P c ) by adding 20 P c randomly selected in Table 4. Similarly to what was done by Lefèvre et al. (2013) and Oumbe et al. (2011), the selection takes into account the modelled marginal distribution established from observation. More precisely, the uniform distribution is chosen as a model for marginal probability for all parameters except aerosol optical thickness, Ångström coefficient, and total column ozone. The chi-square law for aerosol optical thickness, the normal law for the Ångström coefficient, and the beta law for total column ozone have been selected. The selection of these parametric probability density functions and their corresponding parameters have been empirically determined from the analyses of the observations made in the AERONET network for aerosol properties and from meteorological satellite-based ozone products (cf. Table 4). For the sake of avoiding non-realistic cases, the allocation of the aerosol types is empirically linked to the ground albedo (Table 5).
Each combination (θ S , ρ g , P cloud , P c ) is input to libRadtran, yielding G, B, and D. In addition, G c and B c are obtained by libRadtran using θ S , ρ g , and P c as inputs. K c and K cb are obtained using Eq. (1). A series of 140 000 values for G, B, D, K c and K cb is thus obtained for water clouds and 60 000 for ice clouds. For each triplet (θ S , ρ g , P cloud ), the variances v(K c ) and v(K cb ) are computed over the 20 values K c and K cb . Since the clouds and other atmospheric constituents are mixed up in the atmosphere, K c , or K cb , is expected to change with varying P c . It is observed that v(K c ) and v(K cb ) are very small with respect to the squared mean values of K c and K cb for each triplet (θ S , ρ g , P cloud ), meaning that changes in K c and K cb with varying P c are small, thus supporting Eq. (3).
In order to illustrate this and to present this vast amount of results in a synthetic manner, it is firstly observed that these quantities v(K c ) and v(K cb ) do not vary noticeably with the cloud geometry for a given triplet (θ S , ρ g , τ c ): among the cloud properties for a given phase, the cloud optical depth τ c is the most prominent one. As a consequence, it is possible to illustrate the findings by averaging v(K c ) and v(K cb ) over the cloud geometry properties for each triplet (θ S , ρ g , τ c ). One obtains mean(v(K c )) and mean(v(K cb )). The positive root means of these averages are denoted RM(v(K c )) and RM(v(K cb )): RM(v(K c )) gives at a glance the influence of P c on K c for all cloud geometries. A small RM(v(K c )) means that the mean of v(K c ) is small. The variance v(K c ) and consequently RM(v(K c )) are linked to ∂K c / ∂P c . The lower RM(v(K c )) is, the lower the mean of v(K c ), the lower the change in K c with P c , and finally the lower the error made when using Eq. (3). The same reasoning holds for RM(v(K cb )). RM(v(K c )), respectively RM(v(K cb )), can be considered as a measure of the error made on K c , or K cb , when using Eq. (3). These quantities are also expressed relative to the mean K c and K cb for a given triplet, yielding relative values, noted rRM(v(K c )) and rRM(v(K cb )).

Influence of P c on clear-sky indices
Relative quantities rRM(v(K c )) and rRM(v(K cb )) depend on G and B. A large rRM(v(K c )) may not be important if G is very small. To better understand the results, Fig. 1 displays the averages of G and D for θ S = 40 • as a function of ρ g and τ c for water cloud. The beam irradiance B is not drawn as it does not depend on ρ g ; it decreases rapidly as τ c increases and the diffuse irradiance D tends towards G as  Fig. 1 shows that D increases with ρ g , and that both G and D decrease as τ c increases. Figure 2 exhibits rRM(v(K c )) for each couple (ρ g , τ c ) for θ S = 40 • for water cloud (left) and ice cloud (right). Each cell represents the changes in K c obtained for this θ S and this couple (ρ g , τ c ) when the geometrical parameters of the cloud and the other variables in P c vary. For both cloud phases, rRM(v(K c )) increases with τ c and ground albedo. As a whole, it is small. It is less than 2 % for the most frequent cases, i.e. τ c ≤ 20 and ρ g ≤ 0.7. It can be compared to the maximum relative errors (66 % uncertainty) recommended by the World Meteorological Organization (WMO, 2008) for measurements of hourly means of G or D which are 2 % for high quality, 8 % for good quality, and 20 % for moderate quality.
rRM(v(K c )) reaches a maximum of 8.5 % for τ c = 70 and ρ g = 0.9 for water cloud (7.5 % for τ c = 20 for ice cloud). Large ρ g and large τ c mean more reflected radiation by the ground and more backscattered radiation by clouds. This increases the path of the radiation in the atmosphere and, therefore, increases the influence of P c on K c . As G is small for (τ c = 70, ρ g = 0.9) (Fig. 1), a maximum of 8.5 % is not important in absolute value since it corresponds to approximately 30 W m −2 . This high relative deviation happens only for very high ρ g = 0.9. When ρ g ≤ 0.7, the corresponding error on G is less than 10 W m −2 .
The median and 5th (P5) and 95th (P95) percentiles of RM(v(K c )) for all corresponding couples (ρ g , τ c ) for a given θ S are computed and drawn in Fig. 3 for water cloud (left) and ice cloud (right) as a function of θ S . They are also expressed relative to the corresponding mean K c (Fig. 4) and are called relative median and relative P95. For both phases and for θ S from 0 • to 60 • , the relative median is less than 2 %, and the relative P95 ranges between 3.5 % and 5 %.
All three quantities increase sharply for θ S > 60 • . The relative median, respectively P95, reaches a maximum of approximately 8-9 %, respectively 11-12 % for θ S = 80 • . Then, a decrease is observed for θ S > 80 • . Further computations show that the increase in relative influence with large θ S is mostly due to the increase of the optical path in the atmosphere due to greater θ S and therefore a greater influence of P c and notably the aerosols.
Geosci. Model Dev., 7, 1661-1669, 2014 www.geosci-model-dev.net/7/1661/2014/  Overall, an increase in τ c or θ S increases the path of the sun rays in the atmosphere, and therefore the influence of changes in P c on K c increases along with τ c and θ S . This increase is compensated by a corresponding decrease in G. Since G c rarely reaches 120 W m −2 for θ S = 80 • , the error in G corresponding to P95 is less than 15 W m −2 . The diffuse irradiance D and therefore G are strongly influenced by ρ g . The influence of changes in P c on K c increases with ρ g . Deserts such as northern Africa and Arabia exhibit large ground albedo up to approximately 0.5 (Tsvetsinskaya et al., 2002;Wendler and Eaton, 1983); the error (P95) on G is of the order of 10 W m −2 . Fresh snow-covered or ice-covered areas may exhibit very large ρ g . For ρ g = 0.9, the error on G can be large for small θ S , i.e. 30 W m −2 . One has to be cautious in using Eq. (3) in such extreme cases.
Similar calculations are made for K cb . As expected with an RTM code, K cb changes neither with ground albedo nor with cloud phase; the cloud optical depth is the most prominent variable. Figure 5 exhibits the median and 5th and 95th percentiles of RM(v(K cb )) for all couples (ρ g , τ c ) as a function of θ S and their value relative to the corresponding mean K cb . The relative median, respectively relative P95, is less than 2 %, respectively 3 % for θ S ≤ 60 • . Then, it rises sharply. The relative median, respectively P95, reaches a maximum of approximately 8 %, respectively 17 % for θ S = 80 • . Then, a decrease is observed for θ S > 80 • . Large θ S values correspond actually to low irradiances. The clear-sky B c equals 53 W m −2 for θ S = 80 • , and therefore the corresponding median and P95 errors in B are approximately 4 W m −2 and 9 W m −2 . rRM(v(K cb )) has a tendency to increase as θ S increases. This increase is compensated by a corresponding decrease in B. The clear-sky irradiance B c rarely reaches 90 W m −2 for θ S = 80 • and the maximum error in B is less than 16 W m −2 .
If cases of large θ S and τ c for which the radiation is greatly attenuated are removed by considering only cases for which G > 100 W m −2 , the obtained rRM(v(K c )) and rRM(v(K cb )) are very small, even for large θ S . For θ S equal to 70 and 80 • , the medians are approximately 3 % of K c and K cb , and the P95 are 5 and 7 %, respectively.
It is concluded that -for all considered cloud properties and θ S , and for ρ g ≤ 0.7 -the influence of changes in P c on K c and K cb can be neglected. In these cases, Eq. (3) may be adopted with an error (P95) on G and B less than 15 W m −2 and most often less than 2-5 % in relative value. These results match the WMO requirements for high-quality measurements. However, in applications as discussed in the following section, there will be other sources of uncertainties, and the total uncertainty of any model using Eq. (3) will be greater and probably exceed these WMO requirements.
A similar analysis has been made for double-layered clouds with ice cloud topping water cloud. The water and ice cloud properties have been taken from Table 3, where only water clouds with a height top less than or equal to 5 km were considered since the minimum height of ice cloud base is 6 km. Accordingly, there were 5 (water cloud) × 3 (ice  cloud) cases. Results and conclusions are similar to those for single-layered clouds.

Practical implications
A first practical advantage in adopting Eq. (3) instead Eq. (2) is that two independent models -one for modelling G c and B c , the other for modelling the effects of clouds -can be used. If the approach selected to assess the SSI is based on a LUT-based model, using Eq. (3) means that two LUT-based models for K c and K cb can be computed with only one typical set P c0 , therefore strongly reducing the number of runs of the RTM. One may select the following P c0 : -The middle latitude summer from the USA Air Force Geophysics Laboratory (AFGL) data sets is taken for the vertical profile of temperature, pressure, density, and volume mixing ratio for gases as a function of altitude.
-Aerosol properties are as follows: optical depth at 550 nm is set to 0.20, Ångström coefficient is set to 1.3, and type is continental average.
-Total column content in water vapour is set to 35 kg m −2 .
-Total column content in ozone is set to 300 Dobson unit.
-Elevation above sea level is 0 m.
It has been checked that the difference in K c and K cb using different typical sets P c0 was negligible, provided that the selected P c0 does not include extreme values.
As an example, this approach is that used in the MACC/MACC-II (Monitoring Atmosphere Composition and Climate) projects to develop the new Heliosat-4 method for a fast assessment of G, D and B (Qu, 2013;Qu et al., 2014).  cloudy conditions, the computation time for abaci combining the dimensions of McClear abaci and those for K c and K cb would have amounted to years. The immense gain in time justifies the slight loss in accuracy. Except θ S and ρ g , inputs to both models are independent. This is another practical advantage of Eq. (3) since it allows efficiently coping with the fact that P c and P cloud may not be available at the same spatial and temporal resolutions. This is exactly the case in the MACC/MACC-II projects. On the one hand, these projects are preparing the operational provision of global aerosol properties forecasts together with physically consistent total column content in water vapour and ozone (Kaiser et al., 2012;Peuch et al., 2009). These data are available every 3 h with a spatial resolution of approximately 100 km. They are inputs to the McClear model, yielding G c and B c . On the other hand, these projects are preparing the provision of P cloud at high temporal (15 min) and spatial (3 km at nadir) resolutions from an appropriate processing of images taken by the Meteosat Second Generation satellites. P cloud will be input to the K c and K cb models. Using Eq. (3) implies that the SSI may be computed at the best available time and space resolutions by resampling G c and B c , instead of resampling all variables contained in P c .

Conclusions
This Technical Note analyses the influence of the prominent atmospheric parameters on the SSI, with the objective of finding a practical way to speed up the calculations with an RTM. The presented results have been obtained by the RTM libRadtran. It has been checked that the results and conclusions do not depend on this model by obtaining similar results with the streamer RTM (Key and Schweiger, 1998).
It was found that -for all considered cloud properties, solar zenith angles and ground albedos -the influence of changes in clear-atmosphere properties on K c and K cb is generally less than 2-5 %, provided that the ground albedo is less than 0.7. This variation is similar to the typical uncertainty associated with the most accurate pyranometers. In these cases, Eq. (3) may be adopted with an error (P95) on G and B less than 15 W m −2 .
The longer the path of sun rays in the atmosphere is, the greater this variation and the greater the influence of clearatmosphere properties on the clear-sky indices. The mean error made when using Eq. (3)  and direct irradiances expressed as the 95th percentile (P95) is less than 15 W m −2 . The P95 can be greater than 15 W m −2 when the ground albedo is greater than 0.7. In that case, one should be cautious in using Eq. (3). Such high albedos are rarely found; they may happen in case of fresh snow. Like in other RTMs the beam irradiances are modelled by libRadtran as if the sun were a point source. On the contrary, pyrheliometers measure the radiation coming from the sun direction with a half-aperture angle equal to 2.5 • according to WMO standards. The diffuse irradiance in this angular region is called the circumsolar irradiance (CSI). If it were to be compared to measurements, irradiances estimated in this work have to be corrected by adding CSI to B, and by removing CSI from D. In clear sky, the CSI correction to B is approximately 1 % of B (Gueymard, 1995;Oumbe et al., 2012). Under cloudy skies, and especially thin clouds, the CSI can be greater than 50 % of B. A CSI correction needs to be applied only in cloudy skies. Therefore the CSI can be taken into account a posteriori by correcting K c and K cb obtained by Eq. (4) with a specific model.
The presented work has demonstrated that computations of the SSI can be made by considering independently the clear-sky conditions and the cloudy conditions as shown in Eq. (3). A first practical advantage is that two independent models may be developed and used: one for clear-sky conditions and the other for cloudy conditions with their own set of inputs. Another practical advantage is that it allows efficiently coping with cloud and clear-sky variables available at different spatial and temporal resolutions.
These results are important in the view of an operational system as it permits separating the whole processing into two distinct and independent models, whose input variable types and resolutions may be different. The benefit of this separation is not limited to LUT-based models. For example, one may combine LUT-based models for K c and K cb with an analytical model predicting G c and B c such as the ESRA model (Rigollier et al., 2000) or the SOLIS model (Mueller et al., 2004). When both models are LUT-based, using Eq. (3) means two ensembles of abaci: one for clear-sky and the other for cloudy skies. In doing so, the number of entries for each ensemble is reduced leading to the reduction of (i) the size of the abaci, (ii) the number of combination between parameters, and (iii) the total number of interpolations between nodes, thus increasing the speed in computation.