Coupling sky images with radiative transfer models : a new method to estimate cloud optical depth

A method for retrieving cloud optical depth (τc) using a UCSD developed ground-based sky imager (USI) is presented. The radiance red–blue ratio (RRBR) method is motivated from the analysis of simulated images of various τc produced by a radiative transfer model (RTM). From these images the basic parameters affecting the radiance and red–blue ratio (RBR) of a pixel are identified as the solar zenith angle (θ0), τc, solar pixel angle/scattering angle (θs), and pixel zenith angle/view angle (θz). The effects of these parameters are described and the functions for radiance, Iλ (τc,θ0,θs,θz), and RBR(τc,θ0,θs, θz) are retrieved from the RTM results. RBR, which is commonly used for cloud detection in sky images, provides non-unique solutions for τc, where RBR increases with τc up to about τc = 1 (depending on other parameters) and then decreases. Therefore, the RRBR algorithm uses the measured Imeas λ (θs,θz), in addition to RBRmeas (θs,θz), to obtain a unique solution for τc. The RRBR method is applied to images of liquid water clouds taken by a USI at the Oklahoma Atmospheric Radiation Measurement (ARM) program site over the course of 220 days and compared against measurements from a microwave radiometer (MWR) and output from the Min et al. (2003) method for overcast skies. τc values ranged from 0 to 80 with values over 80, being capped and registered as 80. A τc RMSE of 2.5 between the Min et al. (2003) method and the USI are observed. The MWR and USI have an RMSE of 2.2, which is well within the uncertainty of the MWR. The procedure developed here provides a foundation to test and develop other cloud detection algorithms.


Introduction
Advances in solar photovoltaic (PV) technologies have paved the way for higher PV penetration but as we rely more heavily on solar energy it becomes ever more crucial to understand the characteristics of this energy source.Unlike fossil fuels, solar energy has an inherent variability that causes difficulty when integrating solar energy into the grid.Improving forecasting of the available solar irradiance will support management of the electric grid and electricity markets and therefore ensure a more economical integration of solar power (Mathiesen et al., 2013).Currently several different methods are used to forecast at different spatial and temporal resolutions including numerical weather prediction (e.g., Lorenz et al., 2009;Mathiesen and Kleissl, 2011) and satellite imagebased forecasting (e.g., Hammer et al., 1999).For short-term forecasting, whole-sky imagery has been used (e.g., Urquhart et al., 2013).
Physics-based solar forecasting using whole-sky imagery requires geolocating clouds in the sky images, estimating their optical depth, motion, and dynamics (Chow et al., 2011).To estimate a cloud's optical depth (τ c ), the most advanced methods separate the image into clear sky, thin cloud, and thick cloud and assign a τ c to each of these groups.To distinguish thin and thick clouds, the red-blue ratio (RBR) (or a function of RBR) has been used as the default method (Koehler et al., 1991;Shields et al., 1993;Chow et al., 2011;Ghonima et al., 2012;Roy et al., 2001).It is defined as the ratio of the signal from the red channel to the signal from the blue channel.The RBR method takes advantage of Rayleigh scattering being greater in the blue wavelengths than the red wavelengths.When Rayleigh scattering is the predominant form of scattering, such as in clear skies, the RBR for a given view angle is smaller than under cloud scattering.RBR successfully differentiates clear sky from thin clouds and to a more limited extent thick clouds, but the RBR has not been applied to differentiate τ c .It is also difficult to apply the RBR method in the circumsolar region as thick dark clouds have lower RBRs than clear sky (Chow et al., 2011).In fact we will demonstrate in this paper through radiative transfer modeling (Sect.4) that RBR by itself is ineffective for differentiating τ c even for homogeneous cloud layers.
Differences in τ c can greatly affect the irradiance available for solar energy production.For solar power integration into the markets the impacts of forecast errors at individual plants can be considered to scale linearly with forecast error.Any kW of power not produced at solar power plant A will have to be provided by another dispatchable generator B and, depending on market rules, will result in economic losses for the solar power plant.For this analysis we consider the accuracy requirement of global horizontal irradiance (GHI) to be ±5 % of clear-sky GHI (Fig. 1a). Figure 1b and c demonstrate the corresponding absolute and relative error in τ c for a 5 % error of the clear-sky GHI.Relative τ c accuracy required for solar forecasting is large for thin clouds (τ c ∼ 1) and thick clouds (τ c > 30).A minimum occurs at τ c = 16, where a 21 % error in τ c is permissible for avoiding an underprediction of GHI by 5 % of the clear-sky value.PV power production is a close to linear function of global irradiance.Therefore Fig. 1a expresses effects of cloud optical depth on PV power production.Even at τ c = 100 PV power output is still at over 10 % of clear-sky conditions, which is still important for solar energy production.
Most current cloud detection methods are designed empirically using lookup tables and/or thresholds that are adjusted to "work" with a specific imager and cloud conditions (see Sect. 2).The present paper breaks new ground in that it attempts to improve our fundamental understanding of the impact of radiative transfer (RT) and τ c on the radiance and RBR of a given pixel in a sky image.To analyze this relation, the spherical harmonic discrete ordinate method (SHDOM) (Evans, 1998;Pincus and Evans, 2009) is used to produce synthetic overcast sky images (Sect.3) and analyze the determinants of sky imager (USI) radiances (Sect.4).The results reveal nonlinearities and non-monotonic behavior in radiances and RBR that explain many of the challenges previously observed with empirical cloud detection methods.The insights gained through RT are utilized to develop a τ c retrieval algorithm for sky imagery (Sect.5).The algorithm is compared to other methods in Sect.6 and a discussion and conclusions are presented in Sect.7.

Review of sky imager cloud detection methods and geometrical factors
Individual pixel cloud detection using the output image from sky imagers is based on either the radiance measurement or the ratio between radiance measurements for different wavelengths.Cloud detection algorithms using single channel radiance have found limited success (McGuffe et al., 1989;Kegelmeyer, 1994) due to the similarities in radiance values between clear skies and thick clouds in the visible spectrum.More success has been obtained when cloud detection uses the ratio between radiance measurements at different wavelength bands.One such algorithm is the RBR method, which uses the ratio of camera measurement in the red channel to the blue channel to classify a pixel as cloudy or clear.A fixed RBR threshold between clear sky and cloudy sky (Koehler et al., 1991) led to successful identification of opaque clouds but consistently failed to distinguish thin and clear skies.However, in a study of contrail clouds, Koehler et al. (1991) observed that the ratio of RBR to the clear-sky RBR was similar between contrail cases and motivated a method for identifying thin clouds.In other words, knowing the clear-sky background value aids in the detection of thin clouds.The main factors affecting the clear-sky RBR were found to be the solar zenith angle (θ 0 ), solar pixel angle/scattering angle (ϑ s ), pixel zenith angle/view angle (ϑ z , see Figs. 2, 3b, c for illustrations of these angles), and changes in aerosol properties.This lead to the development of clear-sky libraries (CSL) (Shields et al., 1993;Chow et al., 2011) to express clearsky RBR values under any condition.CSL are constructed by binning pixel values from clear-sky images into matrices as a function of θ 0 (Fig. 2), ϑ s , and ϑ z .From the CSL it is then possible to simulate a clear sky (Fig. 3a) for any given day, allowing the calculation of the ratio of measured RBR to clear-sky RBR.
The red-blue difference (RBD; Heinle et al., 2010) uses the same principles as the RBR for cloud detection but attempts to mitigate the strong directional variability in the RBR due to variability in the radiance, I λ , of the blue chan-Figure 2. Diagram of the UCSD sky imager (USI) and related solar and sky geometries.θ 0 is the solar zenith angle.The ϑ s is the angle subtended by the vector pointing at the sun and the vector pointing at the pixel in question.The ϑ z is the angle formed by the vector pointing at the pixel in question and zenith.It is important to note that ϑ s = θ 0 + ϑ z only holds in 2-D but not in 3-D because the incident and scattering directions may not be in the same azimuthal plane.An example of this is illustrated in Fig. 3.
where I r is the I λ in the red channel, and I b is the I λ in the blue channel.However, Ghonima et al. (2012) found minimal differences in performance between RBD and RBR retrieval with RBR outperforming RBD.Gauchet et al. (2012) used RBD combined with a different approach to account for the directional effects in cloud detection, in which they segmented images into five zones, solar disk, circumsolar disk, extended circumsolar disk, main zone, sky horizon, and oro-graphic horizon.The radiance and RBD thresholds to separate clear sky, bright cloud, and dark cloud varied by zone.These approaches have led to improved accuracy of cloud detection, yet limited progress has been made towards understanding the phenomena that influence the performance of these methods.Although a direct relationship with aerosol optical depth (τ a ) and RBR is observed for small τ a , (τ a < 0.3) (Ghonima et al., 2012) no direct relationship has been found between RBR, or other variables determined from sky imagers, and larger optical depths (τ > 0.3) such as those found typically in clouds.This has limited sky imager cloud detection to a binary classification in which the image is segmented into cloud or clear sky.The lack of research on τ c classification also stems from the fact that τ c is challenging to measure accurately and large spatio-temporal variability.Instead, a radiative transfer model is applied here to investigate the interrelationships between radiances, radiance ratios (RBR), and τ c and to devise a method to detect τ c .
3 Radiative transfer modeling of sky images and comparison to measurements

SHDOM model and input parameters
Radiance measurements can be obtained from a 1-D model for homogeneous clouds but in anticipation of future work a 3-D RT model was used.SHDOM is an explicit 3-D RT model that uses discrete ordinates to integrate the radiative transfer equation spatially, while spherical harmonics are used to save memory when computing scattering.SHDOM is more computationally efficient compared to Monte Carlo (MC) methods when solving the whole-sky radiance field.SHDOM is also found to be within 2-3 % (close to the noise level) of the MC models in the intercomparison of 3-D radiation codes (I3RC) (Marshak and Davis, 2005;Cahalan et al., 2005).Because of its computational efficiency and accuracy, SHDOM is selected for this analysis.SHDOM radiative transfer calculations are performed for 161 liquid water overcast skies with homogeneous τ c , ranging in τ c from 0 to 80 at solar zenith angles ranging from 21 to 70 • and for wavelengths corresponding to the peaks of the USI camera's red (620 nm), green (520 nm), and blue (450 nm) channels.
To calculate the single scattering properties of aerosols in the SHDOM simulation such as aerosol effective radius and refractive index, yearly average AERosol RObotic NETwork (AERONET) data from the Atmospheric Radiation Measurement (ARM) Southern Great Plains (SGP) site for the year of 2013 are used (Holben et al., 1998(Holben et al., , 2001) ) (Table 1).Background Rayleigh and aerosol optical depths are also obtained from yearly averages taken from the sun-tracking photometer at the ARM SGP site.Spectral surface reflectances of 0.043, 0.068, and 0.071 were used for the blue, green, and red channel simulations, respectively (Marchand et al., 2004).A cloud droplet effective radius of 8 µm (Min et al., 2003) is used to obtain the single scattering properties of the clouds in the SHDOM simulations.Given the desired τ c , cloud liquid water content (LWC) for input to SHDOM is computed as (Stephens, 1978) LWC where ρ l is the density of liquid water and z is the cloud geometric thickness.In this study, LWC is assumed constant between a cloud base of 1 km and cloud top of 2 km, giving a z of 1 km.
The SHDOM output radiance field is used to reproduce a sky image that would be obtained through a fisheye lens with an equisolid angle projection (Miyamoto, 1964): where f is the focal length, and r is the distance from the principal point in the image plane.

USI hardware and calibration of the signal to radiance
On 14 March 2013 we deployed two USIs (serial numbers 1.7 and 1.8) at the ARM SGP site.The instrument domes were cleaned weekly.Daytime images from the USIs were collected continuously every 30 s for 220 days.Since USI 1.8 was located closer (at 200 m distance) to the instruments used for comparison, it is used for the analysis.The optical setup included a Sigma 4.5 mm fisheye lens, an infrared filter, and an Allied Vision GE2040 CCD camera (Fig. 2).
The fisheye lens creates an equisolid angle projection onto the CCD, resulting in an image where the solid angle subtended on each CCD cell (pixel) is approximately constant.Custom apertures were inserted into the lens of both USIs with diameters of 700 and 1000 µm for USI 1.7 and 1.8, respectively.A Bayer color filter on the CCD separates pixels into red, green, and blue pixels allowing for multispectral images.Three different images are taken at different exposure times and combined to create a high dynamic range image (Urquhart et al., 2015).The signal measured by each pixel is related to the amount of photons that are transmitted through the optics and converted to a voltage.The signal measured can therefore be calibrated to estimate the irradiance, E λ at a wavelength band, incident on a pixel.The radiance I meas λ observed by each pixel can then be calculated using where v is the camera measurement in counts at a given pixel, ) is a calibration factor between v and I meas λ , A in is the area of the pixel, is the solid angle, and λ is the wavelength band.Given the equisolid angle lens, A in , , and λ are assumed constant across the image sensor, resulting in a linear relationship (after correcting optical errors as shown in the Appendix and camera sensor nonlinearities which are negligible for our camera as seen in Fig. 4) between the camera signal v and the radiance I λ at the pixel's ϑ z as The calibration constant C 2 λ is obtained as the average (denoted as overbar in Eq. 6) of 131 overcast (cloud fraction (CF) is greater than 0.9) images on 98 different days.Overcast skies are preferred because the radiance is more homogeneous and since the method by Min et al. (2003) could be applied to obtain the τ c that is input to SHDOM.C 2 λ values are 1.16 × 10 −4 , 1.11 × 10 −4 , and 9.69 × 10 −5 W m −2 st −1 nm −1 for the red, green, and blue channels, respectively.Figure 4 demonstrates the three signal calibrations with a relative root mean square error (RMSE) of 0.155, 0.148, and 0.144 for the red, green, and blue channel respectively.This RMSE is within the range of the radiance variability expected in overcast clouds (Szczodrak et al., 2001).Field calibration to modeled SHDOM data was preferred here as lab calibrations of sky imagers are rarely available.Therefore the calibration method presented here is more widely applicable and provides a calibration that is consistent with the Min et al. (2003) method.Validation of the field calibration with independent lab calibration is left for future work.

Application to other imaging systems
Other ground-based sky imaging designs have also been developed (Seiz et al., 2007;Souza-Echer et al., 2006;Calbo and Sabburg, 2008;Cazorla et al., 2008;Heinle et al., 2010;Román et al., 2012;Gauchet et al., 2012) with the most dissimilar design consisting of a downward-pointing camera capturing the sky from a reflection off a spherical mirror (Pfister et al., 2003;Kassianov et al., 2005;Long et al., 2006;Mantelli et al., 2010;Martínez-Chico et al., 2011).Most ground imaging devices follow a relationship between the camera's signal and radiance similar to Eq. ( 5), differing only in the wavelength region λ, calibration factors C 1 λ , and C 2 λ and optical and sensor errors, with non-equisolid lens camera systems requiring to be specified per pixel.Therefore the method presented here can be adapted easily to images from other sky imaging systems.

Comparison of real and synthetic (SHDOM) images and stray light correction
Example measured images and their SHDOM equivalent images are illustrated in Fig. 5. Differences between the clearsky image (Fig. 5a) and the synthetic image (Fig. 5b) highlight the impacts of stray light as well as the vertical smear stripe caused by the CCD sensor (Fig. 5c).The stray light is particularly strong in the circumsolar region causing enhancement of the red radiance of up to 50 %.Stray light is caused by light from the direct beam being scattered through the optics (mainly the protective acrylic dome).This means that stray light is strongest for τ c = 0 and should decrease to 0 once clouds are thick enough to eliminate the direct beam, which occurs at roughly 5 < τ c < 12, depending on the solar zenith angle.Particular optical reflections are observed as circular patterns throughout the image that are aligned with the solar azimuth.While stray light patterns are often consistent for the same sun position, misalignments in the camera optics (e.g., during instrument maintenance) can lead to stray light changing under constant ϑ s and ϑ z , making it difficult to implement a general stray light correction, for example through a lookup table.Stray light leads to brighter pixel values than expected, which in turn can lead to misclassifications of clear sky as thin clouds (in the range 0 < τ c < 3).To mitigate some of the stray light effects the SHDOM results for clear sky (τ c = 0) are replaced by the measurements from the CSL for the rest of this analysis.Figure 5d-f demonstrate USI images and a synthetic image from SHDOM for τ c of 30.The cloud optical depth for input to SHDOM was determined from Min et al. (2003) measurements.The majority of the sky (ϑ z < 80 • ) red radiance differs by less than 5 %.At τ c = 30 direct normal irradiance (DNI) is absent and stray light can be neglected.

Impact of geometrical parameters and cloud optical depth on radiance and RBR
As described in Sect.2, individual channel radiance and RBR are the most fundamental parameters for cloud detection in sky images.To obtain τ c , the functions I λ (τ c ) and RBR(τ c ) must be parameterized.Furthermore, geometrical parameters (ϑ z , ϑ s ) and solar position (θ 0 ) have been found to affect RBR(τ c ) and I λ (τ c ) (Shields et al., 1993) such that we must obtain RBR(τ c , θ 0 , ϑ s , ϑ z ) and I λ (τ c , θ 0 , ϑ s , ϑ z ) to solve the inverse problem.SHDOM enables us to analyze each of these parameters individually.In this section we will demonstrate how each variable affects I λ and RBR.

Solar pixel angle
Figure 6a demonstrates for the red and blue channel that radiance decreases with increasing ϑ s for non-thick clouds (τ c < 30).For thin clouds (τ c = 1) the radiance peaks in the solar region as a result of the forward scattering peak of the cloud phase function.As τ c increases, radiance becomes constant with ϑ s .Figure 6a shows that the blue radiance is larger than the red radiance under clear skies, except for very small ϑ s , while for cloudy skies the two radiances are more similar.Therefore, most cloud detection methods assume that RBR is higher for clouds than for clear sky; however, Fig. 6b demonstrates that this is not always the case.At small ϑ s (ϑ s < 6 • ) the RBR of thick clouds is lower than that of clear sky.Moreover, for τ c = 1 thin clouds have a higher RBR at ϑ s < 30 • than thick clouds.For τ c < = 10 RBR increases as ϑ s approaches the solar region.At higher τ c (τ c > = 30) RBR becomes constant over ϑ s .Note that all of the statements in Sect. 4 strictly only apply for the θ 0 and ϑ z shown in the figure, but Figs.6-8 indicate that the conditions θ 0 = 60 • , ϑ s = 60 • , and ϑ z = 60 • are representative for a wide range of conditions.

Sky imager zenith angle
Near the horizon (large ϑ z ), diffuse irradiance is commonly observed to be enhanced.Horizon brightening is indeed observed for clear skies in Fig. 7a.As clouds become thicker the dependence of radiance on ϑ z is inverted and radiance decreases with increasing ϑ z .The radiative transfer transitions into the diffusion regime, where it only depends on ϑ z .In contrast, the RBR dependence has a similar shape independent of the τ c (Fig. 7b).Pixels near zenith have lower RBR than those near the horizon.

Solar zenith angle
The effects of θ 0 are intuitive and consistent with what is observed during a sunset and therefore not graphically presented; the red and blue radiance is observed to decrease with increasing θ 0 .The decrease in radiance is caused by the decrease in extraterrestrial horizontal flux as θ 0 increases.In contrast, RBR is found to increase with increasing θ 0 , reflecting the increase in air mass with increasing θ 0 .Increased air mass causes more blue light to be scattered back into space than red light.

Cloud optical depth
Figure 8a illustrates the ambiguity that arises when attempting to differentiate cloud optical depth with radiance.Radiance reaches a peak around τ c = 3.25 and almost for the entire range of τ c the solution is not unique; i.e., there are two τ c that lead to the same radiance.Figure 8b, however, demonstrates the ambiguity of τ c detection using only RBR.SHDOM simulations demonstrate that as τ c increases RBR increases until it reaches its maximum around τ c = 2 and then decreases until converging to a constant value for τ c > 20.This creates the following challenges: (i) RBR is insensitive to τ c for τ c > 20 and therefore thick clouds of different τ c cannot be distinguished and (ii) there is ambiguity because of the non-monotonic behavior.For example, clouds with a τ c of 1.5 have similar RBR values to clouds of τ c > 20.While (outside the solar region, see Fig. 6b) RBR is a useful differentiator between clouds and clear sky, more information is needed to differentiate between different τ c .
In addition, Fig. 8a highlights one of the main challenges of ground-based images compared to satellite-based cloud detection.In satellite-based τ c detection, the measured radiance can be used to calculate τ c (Nakajima and King, 1990) as the measured (upwelling) radiance monotonically increases with higher τ c .This same method cannot be used for ground-based imagery as radiance increases for thin clouds peaks and then begins to decrease.This means that two τ c can exist that produce the same radiance.It is again important to also note that the curves in Fig. 8 depend on ϑ s and ϑ z .For example, in the circumsolar (ϑ s < 30 • ) region red radiance peaks at τ c ∼ 0.75 while clear sky has a higher radiance and higher RBR than thick clouds (τ c > 5).Approaching ϑ z = 0 • and far from the sun (ϑ s > 60 • ), the red radiance peaks at τ c = 8.75.Near the horizon (ϑ z > 80 • ) the red radiance peaks again at lower τ c = 1.25.

Expressing cloud optical depth through geometrical and solar parameters
Since it is currently computationally infeasible to use SHDOM to solve τ c (I λ , RBR, θ 0 , ϑ s , ϑ z ) in real time (∼ 10 s) as required for sky imager solar forecasting, the homogeneous cases described in Sect.3.1 are used instead to create interpolants.As seen in Fig. 8a and b, τ c (I λ , θ 0 , ϑ s , ϑ z ) and τ c (RBR, θ 0 , ϑ s , ϑ z ) are multivalued functions.Therefore two separate interpolants are created for each function.τ c (I λ , θ 0 , ϑ s , ϑ z ) is split into τ c that is higher than the peak radiance and τ c that is lower than the peak radiance.τ c (RBR, θ 0 , ϑ s , ϑ z ) is similarly split into τ c that is higher than the peak RBR and τ c that is lower than the peak RBR.Section 5 will describe how these interpolants are used to find up to two τ c (I meas λ ) (one for the higher and one for the lower branch of τ c ) and how a unique τ c is obtained.
5 Radiance red-blue ratio (RRBR) method for cloud optical depth measurement

RRBR Algorithm
We have shown that it is difficult to distinguish between different τ c by using RBR alone.As demonstrated in Fig. 8, radiance and RBR are non-monotonic functions of τ c with generally two τ c associated with the same radiance or RBR.However, for most cases, there is a unique τ c solution for a pair of RBRs and radiance.The RRBR method attempts to obtain this solution by first solving τ c (I λ , θ 0 , ϑ s , ϑ z ) and then substituting the (usually two) τ c solutions into RBR(τ c , θ 0 , ϑ s , ϑ z ) and identifying the correct τ c as the one with the smallest |RBR meas (θ 0 , ϑ s , ϑ z ) − RBR (τ c , θ 0 , ϑ s , ϑ z ) |, where RBR meas (θ 0 , ϑ s , ϑ z ) is the measured RBR.The algorithm for the RRBR method is depicted graphically in Fig. 9.I λ at a wavelength of 620 nm is used because its variations with τ c are larger than the other wavelengths.This larger dynamic range reduces the errors caused by instrument noise.The algorithm begins by comparing I meas 620 (θ 0 , ϑ s , ϑ z ) against max (I 620 (τ c , θ 0 , ϑ s , ϑ z )) (e.g., 0.19 W m −2 sr −1 nm −1 in Fig. 8a), where I meas 620 (θ 0 , ϑ s , ϑ z ) is the measured radiance in the camera's red channel.Heterogeneity in clouds can cause I meas 620 (θ 0 , ϑ s , ϑ z ) to be larger than max (I 620 (τ c , θ 0 , ϑ s , ϑ z )); in this case as the pixel conditions fall outside the range of the method the algorithm reverts back to τ c assignment solely based on RBR and τ c (RBR, θ 0 , ϑ s , ϑ z ) is used to find τ c .If there are two solutions, then the τ c associated with the maximum red radiance is used as there is no way to differentiate between the multiple solutions.Clouds brighter than the SHDOM radiance peak were found to only occur in 5.4 % of all pixels.
An example τ c estimate is presented in Fig. 10.The darker clouds (for example for the clouds between the sun and the horizon) are correctly identified as higher τ c even though the RBR is lower than for the thinner clouds.In the circumsolar region, the RBR is largest but the RRBR method correctly identifies a thinner cloud.The black points in Fig. 10b corresponds to undetermined τ c due to signal saturation.Since saturated pixel values near the sun exceed the dynamic range of the USI sensor, the RBR defaults to 1, the red radiance is unknown, and no τ c can be assigned.In practice one could interpolate across the saturated region, but we prefer showing the raw results in this paper.

Impact of 3-D effects
Although the RRBR method is developed from overcast scenarios, we also apply this method to broken cloud scenes.The largest 3-D effect is the geometric difference in a broken cloud's optical path (τ p ) compared to an overcast cloud's τ p (Hinkelman et al., 2007).Figure 11 illustrates the definition of τ p as the optical thickness along the path of the direct solar beam, while τ c is the optical thickness integrated along the vertical direction.For overcast clouds the τ p is simply related to the τ c as However, for partial cloud cover the optical path changes along the cloud and as a result it affects I λ , which in turn affects the RRBR retrieval.Ignoring horizontal photon transport, the RRBR's τ c is then a function of the τ p as in Eq. ( 7), which, unlike the actual τ c , changes across the square cloud.Figure 12 demonstrates how the RRBR method retrieves τ c for a 1 km × 1 km square cloud with a 0.2 km cloud geometric thickness.τ c is observed to increase in the same way that τ p increases.Therefore differences between the actual τ c and the RRBR's τ c will occur based on the geometry of the cloud.Again ignoring horizontal photon transport, in ideal cases, such as a cubic cloud, the region of uniform path length where the actual τ c and the RRBR's τ c are similar is limited to θ 0 45 • .For square clouds with a small vertical extent, the region of uniform path length is increased, while for a square cloud of large vertical extent the region of uniform path length is decreased compared to the cubic cloud.For a parallelogram cloud aligned with the solar beam, the local homogeneity is extended to include most of the cloud base.The specifics of defining when clouds can be consid-  ered locally homogeneous will be left for future work.When horizontal photon transport is included it would be expected to decrease the area of homogeneity.The second major 3-D effect is that heterogeneous clouds are brighter than homogeneous clouds under the same τ c .This is caused by increased upwelling solar irradiance from the unshaded part of the scene, illuminating the cloud from below.The reflected light from the cloud underside increases brightness.This is demonstrated in Fig. 13 where overcast and square clouds were compared for two different spectral surface reflectances (R) for τ c = 10.The results demonstrate that the cloud bottom radiance increases 5 % due to a spectral surface reflectance of 0.08 at a wavelength of 620 nm.Adjusting the RRBR method to account for these effects is left for future work.

Impact of aerosols
The chosen τ a is an additional source of error in the reference SHDOM simulations.Higher actual τ a values than those in the simulations may lead to τ a being classified as τ c , while smaller τ a lead to a reduced τ c estimate.This error is small since most τ c values are much larger than the variations in τ a in the USA.Furthermore, this error is not important for solar forecasting as -spectral effects aside -only the total atmospheric optical depth is of interest to estimate ground irradiance, not the partition between τ a and τ c .As demonstrated in Fig. 14, variations in aerosol optical depth from 0 to 0.2 lead to changes in I λ and RBR of less than 5 %.
The RRBR method was derived based on SHDOM results for a single layer liquid clouds but the model could be extended to ice clouds with additional SHDOM runs. Figure 14 demonstrates results from ice cloud simulations, with an effective radius of 100 µm.Ice clouds are not assessed in this paper as none of the methods used for comparison provide information for ice clouds.As for cloud scenes with multiple layers, the RRBR method represents the additive τ c of all cloud layers.

Other cloud optical depth measurements
The Min et al. (2003) method (Min and Harrison, 1996b;Min et al., 2003) is designed to estimate τ c for conditions with homogenous clouds using the measured atmospheric transmittance of global radiation (also referred to as clearness index).The atmospheric transmittance is obtained using a multifilter rotating shadowband radiometer (MFRSR) as where GHI 415 nm is the global horizontal irradiance, and GHI 415 nm 0 is the top of the atmosphere GHI, both at a wavelength of 415 nm.The MFRSR measurements at a wavelength of 415 nm is used in Eq. ( 8) to reduce effects of gaseous absorption.GHI 0 is adjusted from the true top-ofatmosphere GHI to remove τ a influences on T , by applying Langley regression calibrations from the DNI on clear skies to the GHI (Harrison and Michalsky, 1994;Min and Harrison, 1996b).A discrete ordinate radiative transfer model is applied to identify the τ c corresponding to the measured T (Min and Harrison, 1996a).By default a cloud effective radius (r e ) of 8 µm is assumed in the Min et al. (2003) method, but when liquid water path (LWP) values are available from a microwave radiometer (MWR), then r e is solved iteratively.r e is first solved for with Eq. ( 9 2003) method only works for liquid clouds, only data with estimated cloud base temperature above 5 • were used from the comparison dataset.The cloud base height is determined from ceilometer measurements and the cloud base temperature from the most recent sounding.Accurate τ c is obtained with this method for τ c > 10, but for τ c < 10 the Min method is no longer valid (Turner et al., 2004).

F. A. Mejia et al.: Coupling sky images with radiative transfer models
A true validation requires independence of datasets.Due to the need to calibrate USI radiances, results from the Min method were required to develop the RRBR method for overcast conditions.Therefore, at least for overcast conditions, the datasets are not independent and we call it a comparison.
τ c is also measured by an MWR.The MWR is a microwave receiver that detects the microwave emissions of the vapor and liquid water molecules.It measures cloud LWP in the zenith direction within a field of view of 6 • (Liljegren, 2000; Cadeddu et al., 2013).τ c can then be estimated as (Stephens, 1978) where LWC is the cloud liquid water content and ρ l is the density of liquid water.A r e of 8 µm is assumed, as in the Min et al. method.The MWR has an irregular timestep ranging from 20 to 40 s.The uncertainty in the LWP obtained from the MWR is ±0.03 mm (30 g m −2 , Morris, 2006), which corresponds to a τ c of ±5.6 with Eq. ( 8).

Comparison in overcast conditions with Min algorithm
Data from the Min et al. (2003) algorithm are compared to the average τ c from an entire USI image.Since the Min et al. (2003) method assumes overcast skies, only conditions with CF > 0.7 are used for this analysis yielding 5197 data points (about 43 h of data).The mean transmission of horizontally heterogeneous clouds is higher than the transmission of a uniform cloud with the same mean optical depth (Hinkelman et al., 2007).This is caused by the nonlinear relationship between τ c and radiance.
To adjust the heterogeneous USI τ c retrieval to be consistent with the Min et al. method, the USI τ c was converted to irradiance for each pixel using a lookup table, averaged over the entire image in irradiance space, and then converted back to τ c . Figure 15 compares results from both methods.An R 2 of 0.99 reflects the high correlation between the two methods.The relative RMSE decreases as τ c increases as demonstrated in Table 2, with thin clouds (τ c < 10) having an RMSE of 27.2 % and thick clouds (τ c > 30) having an RMSE of 5.8 % with the overall RMSE being 8.2 %.RMSE at τ c > 10 is well below the 21 % required for solar energy applications (Fig. 1) and validates the RRBR method for thick overcast clouds (τ c > 10), but for τ c < 10 the Min et al. (2003) method is no longer valid (Turner et al., 2004) and the relative RMSE increases drastically.These differences in RMSE with τ c highlight the difficulties in detecting thin clouds correctly.
Note that the zero mean bias error (MBE) between Min and RRBR for CF > 0.7 is partially a result of the cross-calibration in overcast conditions (because the crosscalibration data required CF > 0.95, the cross-calibration and validation data are different).However, the other errors be-  2003) method applied to MFRSR measurements for USI cloud fractions greater than 0.7.USI results are averaged in irradiance space over the hemisphere as shown in Eq. ( 9).tween Min and RRBR in Table 2 are non-zero (i) because some of the data used for the comparison were not used for radiometric calibration, (ii) due to differences in the wavelengths of measurement (620 nm for RRBR and 450 nm for USI), and (iii) because radiation differences as USI use diffuse radiance while Min uses global irradiance to derive τ c .Nevertheless the cross-calibration is expected to reduce the mean bias difference between Min and RRBR methods and the related error metrics are not representative of a truly independent dataset, especially for the CF > 0.7 scenario.

Heterogeneous and homogenous cloud conditions with the microwave radiometer
The RRBR method is compared to τ c estimates from the MWR using the 12 422 pixels in each USI image with ϑ z < 6 • .Figure 16 shows the comparison of the two methods.All conditions result in a RMSE of 3.6 or 19.0 % and R 2 of 0.98, again well within the minimum error requirement of 21 %.Since overcast data were already compared in Sect.6.2, we now focus on cloud fractions of less than 0.7.The RMSE is 2.23 for the heterogeneous cases, which is well within the uncertainty of the MWR measurements of ±5.6 but which corresponds to a relative RMSE of 85.0 % that exceeds the objectives set at the beginning.Just like in overcast conditions (Sect.6.2), RMSE is highest for thin clouds (τ c < 5) at 71.7 %, decreases at medium cloud thickness (10 < τ c < 30) to 31.9 %, and increases once again for thick clouds (τ c > 30) to an RMSE of 56.8 %.The heterogeneous cases are associated with a higher relative RMSE of 85.0 % compared to 8.2 % reported in Sect.6.2 for the homogeneous Min et al. (2003) method.The lower correlation of 0.66 between the two methods is probably related to (i) the uncertainty of the MWR, (ii) random errors in τ c re-  A MBE of −14.4 % is observed demonstrating a tendency for the RRBR method to underpredict τ c .This can further be analyzed when MBE is split into τ c categories.Thick clouds (τ c > 30) have the highest MBE of −52.8 % compared to thin clouds (τ c < 10) that have an MBE of −11.3 %.As de-scribed in Sect.5.2 heterogeneous clouds are brighter than homogeneous clouds because of the reflected light from the ground surface, leading to higher radiance measurements and lower τ c .Another factor that causes higher USI radiance measurements in heterogeneous clouds are cloud sides.Since clouds sides are no longer obscured such as those in overcast clouds, an increase in cloud illumination relative to overcast clouds increases radiance.For clouds that are thicker than the red radiance peak (τ c = 7.25) this increased radiance along the sun-facing edge of the cloud results in an underprediction of τ c .The fact that MBE becomes more negative with increasing τ c could be a result of neglecting 3-D effects in the RRBR method.

Discussion and conclusions
This paper presents an analysis of the atmospheric radiative transfer effects on sky imager cloud detection and retrieval of τ c using synthetic images produced by the SHDOM.Synthetic images demonstrate that θ 0 , τ c , ϑ s , and ϑ z all significantly and often nonlinearly and non-monotonically affect radiance I λ and RBR of sky image pixels.For thin clouds, I λ (ϑ s ) increases rapidly as it approaches the sun, as a result of the strong forward peak in the cloud phase function.In contrast, for thick clouds I λ (ϑ s ) is found to be near constant with solar pixel angle for τ c > 30.ϑ z has two main effects: horizon brightening for thin clouds and horizon darkening for thick clouds.Thick clouds fall in the diffusion regime, where I λ decreases with ϑ z but is independent of other parameters.
At a ϑ s of 45 • I λ (τ c ) is demonstrated to increase with increasing τ c for thin clouds.It reaches a peak at a τ c between 0 and 5, depending on ϑ s and ϑ z .At τ c greater than 5, I λ (τ c ) decreases with increasing τ c .Similar characteristics are ob-served for the RBR although it does not decrease as much after reaching its maximum, making it an effective tool for distinguishing between clear sky and thick clouds.However, neither I λ (τ c ) nor RBR(τ c ) are monotonic, leading to the difficulties in cloud detection and τ c characterization with one parameter.The RRBR method combines the RBR and I λ to overcome the non-monotonic nature of each individual parameter.
Summary statistics of the different comparisons are presented in Table 2.For overcast skies the RRBR yields τ c that is consistent with the Min et al. (2003) method.For heterogeneous cloud fields (CF < 0.7), comparisons with MWR measurements of LWC at zenith demonstrated that the RRBR method provides τ c estimates with typical R 2 of 0.68 and RMSE of 2.2, which is well within the uncertainty of the MWR instrument (±5.6), but more work needs to be done to validate that heterogeneous clouds are within the 21 % uncertainty required for solar applications.As demonstrated by the relative RMSE in Table 2, the RRBR method provides accurate τ c for overcast thick clouds.The relative RMSE is larger for τ c < 10 for all comparison datasets and future development requires a new comparison method for thin clouds.These results validate the RRBR method for overcast clouds but consistent underpredictions of heterogeneous cloud optical depth require improvement in the method.
Characterizing the cloud heterogeneity effects may improve the RRBR method.As the RRBR method is based on interpolants developed from simulations of homogeneous overcast skies, cloud heterogeneity violates the assumptions and is likely the leading source of errors.Errors due to cloud heterogeneity have been analyzed mainly in the context of satellite remote sensing.Varnai (1998) and Chambers et al. (1997) observed that the cloud spatial reflectance variation is smoother than variations in τ c .They hypothesized that optically thicker clouds would scatter more light to their thinner neighboring clouds, causing the thinner clouds to appear brighter and thicker (looking from space), while the thinner clouds would scatter less light to the thicker clouds, making them appear darker and thinner than expected for a homogeneous cloud scene.A similar but opposite effect is observed in ground-based imagery, where thicker clouds shade their neighboring thinner clouds, making them appear darker and thicker, but this effect is moderated by the location of the sun relative to the clouds.Figure 9 also shows that sunfacing cloud edges scatter more light, increasing I λ and leading to thinner τ c estimates than the cloud edges on the opposite side.Cloud edges facing away from the sun will be shaded by the rest of the cloud and will be estimated as being thicker.These 3-D effects introduce noise in RRBR estimations of τ c .Although the comparison methods presented here are able to highlight some errors in the RRBR method no method was accurate enough to provide information about thin clouds (τ c < 10) and future development requires a new comparison method for thin clouds.
For solar applications, the pixel-by-pixel τ c allows an estimate of the contribution of each individual pixel to the DNI, GHI, and diffuse horizontal irradiance at the surface.Solar forecasting will benefit in two ways.(i) Static fields of τ c can be used to spatially modulate surface irradiance fields beyond the standard binary (clear -cloudy) solar irradiance, either in 1-D similar to Eq. ( 9) or using the homogeneous results from SHDOM simulations from a lookup table.(ii) Temporal evolution of τ c can be used to extrapolate future cloud τ c resulting from cloud dynamics.At present, τ c is assumed to be steady as the cloud field is advected.An analysis of τ c evolution especially for individual cumulus or stratocumulus clouds through a time series of the τ c could improve solar forecast accuracy and extend forecast skill to longer forecast horizon.For example, stratocumulus clouds that tend to form overcast layers in coastal Southern California thin during the mornings and then break up rapidly over land.The RRBR can detect the cloud thinning rate and forecast breakup time.

Data availability
The raw image data used as input can be made available upon request to the corresponding author.In fact, he would be happy to provide any data from the ARM field trial from 11 March to 4 November 2013 from either or both camera systems that were deployed.(Each day, a single camera's imagery is about 3 GB and the request volume would have to be reasonable, e.g. for the 10 days cited.)However, it is expected that other authors interested in this paper would be users of their own sky imaging system and will likely want to perform the calibration on their own system.
The sky imagery data used in this work were collected as part of the UCSD Sky Imager Cloud Position Study (https: //www.arm.gov/campaigns/sgp2013sicps)(Urquhart, 2013).Data from that study are available upon request to the corresponding author.

Figure 1 .
Figure 1.(a) GHI divided by clear-sky GHI as a function of τ c for homogenous clouds as derived from SHDOM.The black line represents the results while the blue and red are 5 % offsets in clearsky GHI.(b) Error bounds of ±5 % clear-sky GHI in (a) converted to absolute intervals for τ c .For example, for GHI to stay within 5 % clear-sky GHI at τ c = 30, τ c cannot be more than 7.7 below 30 and not more than 15.6 above 30.(c) Same as (b) but the y axis is divided by τ c .

Figure 3 .
Figure 3. (a) Clear USI image created from a clear-sky library for a solar zenith angle of 60 • on 26 March 2013, 15:00:00 UTC.For the image in (a) solar pixel angles and pixel zenith angles are shown in (b) and (c), respectively.The red lines in (b) and (c) highlight the pixels with ϑ s = 60 • and ϑ z = 60 • , respectively, which are often used in the following chapters to illustrate relationships with cloud optical depth and radiance.North is located on the bottom right corner of the image.

Figure 6 .
Figure 6.(a) SHDOM red and blue channel radiance over various sun pixel angles (ϑ s ) at ϑ z = 60 • , and θ 0 = 60 • (pixels used for Fig. 6 are highlighted as a red line in Fig. 3c).Results are shown for different cloud optical depths from clear (τ c = 0) to thick clouds.(b) RBR as a function of ϑ s at constant ϑ z = 60 • and θ 0 = 60 • .

Figure 7 .
Figure 7. (a) Red and blue channel radiances and (b) RBR over various ϑ z at constant ϑ s = 60 • and θ 0 = 60 • .Pixels used for Fig. 7 are highlighted as a red line in Fig. 3b.

Figure 10 .
Figure 10.(a) USI image for 25 March 2013, 22:10:00 UTC.(b) τ c retrieval from the RRBR method.Pixels inside the black ring are the pixels used for averaging and comparison with the MWR (Sect.6.4).(c) RBR meas (ϑ s , ϑ z ).(d) I meas 620 (θ 0 , ϑ s , ϑ z ).For this scene, the MWR measured a τ c of 0.56 and the USI measured a τ c of 0.20, the highest τ c readings within 10 min of this image are 19.4 and 15.3 for the MWR and USI respectively.

Figure 11 .
Figure 11.Illustration demonstrating differences between RRBR measured cloud optical depth, cloud optical depth, and cloud optical path.

Figure 12 .
Figure 12.(a) SHDOM simulated sky image of a 1 km × 1 km square cloud with cloud geometric thickness of 0.2 km at a θ 0 = 60 • and τ c = 10 (b) and RRBR τ c retrieval.The RRBR τ c is 8, 0.2 km from the cloud edge facing the sun.
) using LWP from the MWR and the Min et al. (2003) τ c .Once obtained r e is used as an input in the discrete ordinate model, which provides a different τ c , which leads to a different r e ; this process is repeated until the changes in τ c are below a threshold value.Min et al. concluded that the uncertainty in the inferred cloud properties caused by r e was less than 5 %.Since the Min et al. method uses GHI measurements to estimate τ c , the τ c is representative of the sky hemisphere.At the ARM site the Min et al. (2003) τ c is sampled and reported every 20 s.Since the Min et al. (

Figure 15 .
Figure 15.Comparison of RRBR τ c retrievals from the sky imager vs. the Min et al. (2003) method applied to MFRSR measurements for USI cloud fractions greater than 0.7.USI results are averaged in irradiance space over the hemisphere as shown in Eq. (9).

Figure 16 .
Figure 16.Comparison of USI RRBR vs. MWR measurements of cloud optical depth for CF < 0.7 in black and CF > 0.7 in red.

Table 1 .
Atmospheric radiative properties for the ARM site used as input to SHDOM.τ a and Rayleigh optical depth are averages for the year 2013 from AERONET data.

Table 2 .
Statistics of RRBR comparison against the Min et al. method in overcast skies (cloud fraction > 0.7) and microwave radiometer (MWR) measurements.RMSE[-]is the absolute root mean square error, RMSE [%] is the relative root mean square error, MAE [%] is the relative mean average error, and MBE [%] is the relative mean bias error.