An improved cloud index for estimating downwelling surface solar irradiance from various satellite imagers in the framework of a Heliosat-V method

We develop a new way to retrieve the cloud index from a large variety of satellite instruments sensitive to reflected solar radiation, embedded on geostationary as non geostationary platforms. The cloud index is a widely used proxy for the effective cloud transmissivity, also called clear-sky index. This study is in the framework of the development of the HeliosatV method for estimating downwelling solar irradiance at the surface of the Earth (DSSI) from satellite imagery. To reach its versatility, the method uses simulations from a fast radiative transfer model to estimate overcast (cloudy) and clear-sky 5 (cloud-free) satellite scenes of the Earth’s reflectances. Simulations consider the anisotropy of the reflectances caused by both surface and atmosphere, and are adapted to the spectral sensitivity of the sensor. The anisotropy of ground reflectances is described by a bidirectional reflectance distribution function model and external satellite-derived data. An implementation of the method is applied to the visible imagery from a Meteosat Second Generation satellite, for 11 locations where high quality in situ measurements of DSSI are available from the Baseline Surface Radiation Network. Results from our preliminary 10 implementation of Heliosat-V and ground-based measurements show a correlation coefficient reaching 0.948, for 15-minute means of DSSI, similar to operational and corrected satellite-based data products (0.950 for HelioClim3 version 5 and 0.937 for CAMS Radiation Service).


Introduction
Downwelling surface solar irradiance (DSSI) is one of the Essential Climate Variables defined by the Global Climate Observing 15 System (GCOS, 2016). It is the solar part of the downwelling irradiance at the surface of the Earth and on an horizontal unit surface. The solar irradiance is defined as the integration on the spectral interval 290-3000 nm, accordingly to WMO (2014).
DSSI considers the irradiance coming from all directions of the hemisphere above the surface: the irradiance coming from the direction of the Sun, plus a diffuse component due to scattering caused by the atmosphere (clouds, gases, aerosols) and reflection by the surface. 20 The knowledge of DSSI variations in space and time is of primary importance for various fields such as Earth sciences, renewable solar energy industries, agriculture, or some medical fields. To meet all these needs, an ideal information on DSSI would feature high spatio-temporal resolution, a coverage of the entire Earth's surface, and the longest time period possible. empirical proxy for effective cloud transmissivity. The latter, also named "clear-sky index" within the scientific community of solar energy, is defined as the ratio of the all-sky surface irradiance to the clear-sky surface irradiance (Long and Turner, 2008; 60 Beyer et al., 1996), i.e. the DSSI in cloud-free conditions. Figure 1. The calculation of a cloud index for a given location. The cloud index is the ratio between the distances "measurement to clear-sky" (red arrow) and "overcast-sky to clear-sky" (black arrow).
The cloud-index concept is based on the idea that the presence of a cloud brightens locally pixels of satellite shortwave imagery. In general, the value that quantifies reflectances of a given location observed from the top of the atmosphere (TOA), is comprised between a low and a high boundary values. The low boundary value X min is taken as the clear-sky case and the high one X max as the most cloudy case. The attenuation of downwelling surface solar irradiance by clouds is roughly given as a linear function of the difference between the measured value X sat and the clear-sky boundary, relatively to the cloudy caseclear case difference, as illustrated in Figure 1. The cloud index n is then given as: Differences between cloud-index methods mainly rely on : modifications of the relationship between the cloud index and the effective cloud transmissivity (Zarzalejo et al., 2009;70 Perez et al., 2002;Gupta et al., 2001;Rigollier and Wald, 1998); the choice of the variable used to calculate the cloud index, for example TOA albedo (Darnell et al., 1988), reflectance (Wang et al., 2014;Gupta et al., 2001;Möser and Raschke, 1984)), Lambert equivalent reflectivity (LER) (Herman et al., 2018;Dave, 1978) or raw satellite numerical counts (Pfeifroth et al., 2017;Perez et al., 2002); the way to retrieve the X max and X min for the chosen variable.
If very different approaches are used to estimate the upper boundary, the lower boundary is "archive-based", in most literature we reviewed: it is a minimum based on a time series of past satellite imagery. Such an approach is hardly appliable to non geostationary satellites due to variable viewing geometries and a low revisit time. As an example, Wang et al. (2014) use a climatology of surface albedo to derive DSSI from the Ozone Monitoring Instrument (OMI), embedded on the sun-synchronous satellite Aura.

80
In this paper, we aim at finding an alternative to the need for archives of satellite imagery. It would then be easier to consider imagery from non-geostationary spaceborne platforms and produce a worldwide coverage. It would also let the results more reproducible, by making each instantaneous estimate independent of the preceding time series. Other drawbacks are frequent in archive-based methods such as systematic underestimations of the lower boundary X min (e.g. unusual dark shadows on the ground taken as the clear-sky reference), contamination of X min by clouds for cloudiest regions and of course, the difficulty to 85 find a trade-off between a large time window that ensures the observation of clear sky instants, and a small one that captures the temporal variability of X min .
But the development of an alternative to such archive-based approaches also means dealing with new issues: previous methods based on archives can be less dependent on absolute calibration of the original imagery (Pfeifroth et al., 2017;Perez et al., 2002) and consider implicitely the anisotropy of X min . The pixel-to-pixel estimation of X min is a surrogate for modeling the 90 influence of viewing geometry on measured reflectances, while the slot-by-slot temporal characterization of X min pictures the influence of varying solar-viewing geometry for the diurnal cycle of each pixel's reflectivity.
We aim at developing an alternative "stateless" method to extend the application field of the cloud-index approach to a wider variety of orbits and optical shortwave sensors. In order to limit the effects of molecular scattering, ozone absorption and 95 polarization present in the ultraviolet, and of the absorption of radiation by clouds in the near infrared, the method considers satellite imagery in the spectral range 400-1000 nm (λ < 1000 nm). This range is also wide enough to permit the use of imagery from many meteorological satellite imagers launched since the beginnings of spaceborne Earth observation (e.g. different generations of the Advanced Very High Resolution Radiometer and of GOES, Meteosat and Himawari radiometers).
Heliosat-V makes an extensive use of simulations from a radiative transfer model to estimate the upper boundary variables 100 X min and X max , which are here reflectances at the top of the atmosphere. This notably relies on the hypothesis that the absolute calibration of the satellite measurements is sufficiently well-known to ensure the quality of DSSI retrievals. It also means that the spectral sensitivity of the sensor and the anisotropy of reflectances caused by surface and atmospheric components have to be explicitely described to produce accurate estimates.
The impact of the anisotropy of surface reflectance has notably been shown for estimates of a cloud index derived from 105 measurements of the ultraviolet/visible Global Ozone Monitoring Experiment 2 (GOME-2) and OMI by Lorente et al. (2018).
The latter study also highlights the improvement of simulated shortwave clear-sky reflectances at the TOA, when using a model of bidirectional reflectance distribution function (BRDF) parameterized with data derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) spaceborne instruments. The method foreseen to compute the cloud index of Heliosat-V and eventually the DSSI is described in Section 2, along with the protocol of validation. Validation results are presented in Section 3 for simulated reflectances at the top of the atmosphere and for downwelling surface solar irradiance estimates. Section 4 is dedicated to the discussion of results, and Section 5 to the conclusion and perspectives. to the satellite radiometer viewing geometry and spectral sensitivity. Reflectances are defined by the relation: with L the upwelling radiance at TOA for a given spectral channel, E 0 the downwelling spectral solar irradiance at the top of the atmosphere on a perpendicular plan weighted by the spectral response function (SRF) of the channel, and θ s the solar zenith angle for a given location (i.e. latitude and longitude) and a given time. E 0 varies mainly with the Sun-Earth distance, computed here with the Solar Geometry 2 algorithm (Blanc and Wald, 2012). The cloud index is then defined as: where ρ sat is the reflectance measured by the radiometer for the given spectral channel, while ρ clear and ρ ovc are estimates of the reflectance that would be measured by the same sensor for, respectively, a clear-sky scene, and an overcast scene, i.e.
with an optically thick cloud covering the whole pixel considered. The notion of "optically thick cloud" will be described in detail in the subsection 2.3.
Because of its definition, the cloud index may also be calculated with radiances. We consider here reflectances in order to 130 visualize the anisotropic nature of different scenes. It has also the advantage to be a normalized quantity, so we can compare results for different radiometric channels and different SZAs.
The relationship between n and DSSI varies slightly from one method to another, in particular for highest and lowest values of n. The core of the relationship for intermediate values of n follows usually: where G is the all-sky DSSI and G c is the DSSI in clear-sky conditions and is provided by an external model. The external model used in this paper will be presented and discussed in section 2.4. The clear-sky index K c is largely used to simplify the reading and is defined as: so we can rewrite Equation (4) as: Several modifications of the relation K c (n) have been proposed, e.g. by Rigollier and Wald (1998)

The clear-sky reflectances ρ clear
We use a radiative transfer model to estimate what a spaceborne optical imaging system would measure in clear-sky conditions, for a given radiometric channel. Using simulations in cloud indices has previously been done notably to retrieve effective cloud fractions from the OMI instrument (Lorente et al., 2018;Veefkind et al., 2016;Stammes et al., 2008). We apply the same approach to satellite radiometers. "ammonium" and "nitrate" added to the aerosol categories of the service. It is worth mentioning that we do not use it, even such change would not affect the method itself.
An algorithm developed by Lefèvre et al. (2013) translates MACC partial aerosol optical depths information into aerosol mixtures designed for the Optical Properties of Aerosols and Clouds (OPAC) software package (Hess et al., 1998). These mix-170 tures are associated to aerosol properties: scattering and absorbing coefficients, single scattering albedo, asymmetry parameter and the Angström coefficient. The total AOD at 550 nm is then calculated as the sum of partial AOD at 550 nm provided by CAMS. As libRadtran needs a total AOD input for the simulated wavelength, the OPAC Angström coefficient of the given mixture is used to estimate the AOD at the required wavelength.
The reflective properties of land surfaces are described with the RossThick-LiSparse (Ross-Li) model of bidirectional re-175 flectance distribution function (Roujean et al., 1992;Lucht et al., 2000). It is then possible to consider the variations of the surface reflectance depending on viewing and solar zenith angles and on the azimuthal difference of both geometries ∆φ. The Ross-Li model decomposes the BRDF of a surface into a sum of three components: an isotropic contribution, independent of viewing and solar geometries; a volumic contribution, following a mathematical model of an idealized canopy ; and a geometric contribution, considering the shadows induced by the roughness of the surface.

180
Algorithms have been developed to estimate the parameters f iso , f vol and f geo that weight respectively each of the contributors to the surface reflectance for all lands. This has notably been made with the imagery produced by the Moderate Resolution Imaging Spectroradiometer (MODIS) embedded on Terra and Aqua satellites (Wanner et al., 1997;Lyapustin et al., 2018).
We test here simulations with data from the Algorithm for Modeling MODIS Bidirectional Reflectance Anisotropies of the Land Surface (AMBRALS) (Wanner et al., 1997) with its derived product MCD43C1 v6 (Schaaf et al., 2002).

185
This product provides f iso , f vol and f geo parameters with a 0.05°resolution (about 6 km at the equator), a daily sampling rate, 16-day average and for seven spectral channels, including 4 channels in the 400-1000 nm spectral interval considered for the Heliosat-V method (Fig. 3). Owing to libRadtran documentation (Mayer et al., 2017), the values of each parameter are assigned to the central wavelength of its channel and a linear spectral interpolation is applied for the radiative transfer calculations. For wavelengths shorter than the 0.47 µm MODIS channel, values are considered spectrally constant. For wavelengths longer than 190 0.85 µm, the interpolation is made between parameters at MODIS channels 0.86 µm and 1.24 µm.

The overcast-sky reflectances ρ ovc
Cloud-index methods in the literature use various ways to estimate the TOA reflectances in overcast conditions ρ ovc (Perez et al., 2002;Lefèvre et al., 2007;Pfeifroth et al., 2017). One way to approximate it without the use of archives of satellite imagery has been proposed within the Heliosat-2 framework (Lefèvre et al., 2007). An empirical relation based on the work of 195 Taylor and Stowe (1984) was developed, considering a dependency of ρ ovc with the single solar zenith angle θ s .
However, spectral radiative transfer simulations of ρ ovc show that there is also a significant dependency between the TOA cloudy reflectances and other variables. In Figure 2, we represent 2-dimension histograms of TOA reflectances calculated from such simulations, with a solar zenith angle set to 30°as an example. For most wavelengths, a significant spread of 200 the distribution is observed (Fig. 2, two upper rows), corresponding only to different viewing geometries defined by a linear meshgrid in cosine of viewing zentih angle (θ v ) and difference ∆φ of solar and viewing azimuth angles (Table 1).
In this paper, we assume a cloud optical thickness (COT) of 150 to define optically thick clouds and overcast conditions.
This assumption relies on COT statistics from retrievals by the International Satellite Cloud Climatology Project (ISCCP) and surface measurements of irradiance shown by Trishchenko et al. (2001). The simulations for a low thick cloud (cloud top  reflectances with low clouds can be much lower than for high clouds, for a given cloud optical thickness. But outside these specific spectral regions, the height of clouds will not affect significantly the results of the method.

210
An alternative way is therefore to produce look-up tables (LUT) from radiative transfer simulations, an approach notably applied in the framework of the Heliomont cloud-index method to the Spinning Enhanced Visible and Infra-Red Imager (SE-VIRI) High Resolution Visible channel (HRV) (Stöckli, 2014). It is then possible to take into account the viewing geometry and also the spectral variability of ρ ovc . Assumptions have to be made on the properties of the optically thick clouds as the Heliosat-V method is designed to work by using only one spectral channel in the range 400-1000 nm: cloud top height, phase 215 of cloud, cloud optical thickness, cloud droplet radius or ice crystal shape and size.
Here we construct a liquid cloud LUT of ρ ovc , setting different cloud and atmosphere properties, geometry and spectral grids, as described in Table 1. The optical properties of the clouds come from the precalculated Mie tables provided by the libRadtran software package.
As no information is provided on the actual cloud vertical structure, ρ ovc are calculated as : where ρ ovc,high and ρ ovc,low are respectively derived from the high and low liquid cloud LUTs, interpolated on the viewing and solar geometries of the satellite time series and adapted to the spectral response function of radiometric channels.
An ice cloud LUT is also produced, to study the sensitivity of surface irradiance estimates to the assumed cloud phase. The ice cloud characteristics follow the parameterization by Yang et al. (2013). We use the "aggregate of 8 columns" ice crystal 225 habit and the "severe" degree of roughness, which are notably used for the description of ice clouds in the look-up table of the MODIS collection 6 cloud product (Amarasinghe et al., 2017).

The clear-sky model of surface irradiance G c
The clear-sky surface irradiance is given by the the version 3 of the McClear model (Gschwind et al., 2019

Set-up of validation
To study the validity of the method, we compare DSSI estimates from MSG satellite measurements with pyranometric DSSI data retrieved from measurement stations part of the reference Baseline Surface Radiation Network (BSRN) (Ohmura et al., 1998;Driemel et al., 2018). Considered stations are listed in Table 2 and displayed in the MSG field of view in Figure 4. Only Water vapour total column 20 kg m −2 Aerosols default aerosol described in Shettle (1990) Temperature and pressure profiles AFGL midlatitude summer  Figure 5 shows the time series when data are considered valid, for each station.   The method has been tested on images from the SEVIRI sensor, aboard the Meteosat-9 meteorological geostationary satellite belonging to the family of Meteosat Second Generation (MSG). We consider measurements in the solar channels 0.6 µm and 0.8 µm channels, for the year 2011 and for 11 locations in the field of view of the satellite, and corresponding to locations of pyranometric in situ sensors from the BSRN network. 245 We use the calibration coefficients provided by EUMETSAT that operates MSG. This is worth noting as some calibration methods recommend to use significantly different gain factors to compute radiances from raw numerical counts (e.g. Doelling et al. (2018)). The use of optimal calibration is out of the scope of our work. Still, we compared gains coefficients proposed by EUMETSAT g EUM with those provided by Doelling et al. (2018) g D2018 for the measurements produced by the Meteosat-9 0.6 and 0.8 µm channels in 2011. They show a mean disagreement, calculated as (g EUM − g D2018 )/g D2018 , of about -9 % for 250 0.6 µm and -8 % for 0.8 µm during this period (also illustrated on Fig. A1). Such errors will affect with the same magnitude the agreement between numerical simulations and measurements of clear-sky TOA reflectances, underlining the importance of absolute calibration for the Heliosat-V method.
To compare the results of our method to operational satellite-based products of surface irradiance, we use data from Helio-Clim3 version 5 (HC3v5) and CAMS Radiation (CAMS-RAD) DSSI databases. Both are derived from the imagery of the SE-

255
VIRI sensor and are produced by a Heliosat method: a modified version of Heliosat-2 for HC3v5 (Qu et al., 2014) and Heliosat-4 for CAMS-RAD. Boths products and their descriptions are provided by the SoDa service (http://www.soda-pro.com/).
All validation results thereafter are produced for solar zenith angles lower than 80°.
3 Results in the near infrared 0.8 µm channel are significantly higher, so is the absolute value of STD. Both higher STD and bias for the 0.8 µm will cause a lower precision in the calculation of the cloud index than for 0.6 µm. Correlation coefficients are significantly high, both larger than 0.9 but the correlation is much better for 0.6 µm with 0.974: the variability of ρ clear is 270 significantly better represented (almost 95 % of the total variance) than for 0.8 µm (82 % of the total variance). Figure 7 shows that estimates are able to reproduce partly the diurnal variability observed in clear sky conditions. When studying station by station, highest mean biases reach +0.035 for 0.6 µm (Payerne) and -0.07 for 0.8 µm (Camborne) (see also https://doi.org/10.5194/amt-2020-480 Preprint. Discussion started: 12 February 2021 c Author(s) 2021. CC BY 4.0 License. Fig. B1). It is worth noting that we use MCD43C1v6 BRDF data regardless of their quality flags. We observe although that keeping only the highest quality data slightly improves statistics.

Simulated reflectances at the top of the atmosphere in overcast conditions ρ ovc
The validity of ρ ovc is more difficult to test than that of ρ clear by comparing with satellite measurements as the occurence of optically thick clouds can be rare depending on the location, the season and the hour of the day. We therefore use 9 years of Meteosat measurements, between 2011 and 2019 to extract the 1% most reflective scenes for each station, month and hour of the day (Fig. 7). The first row of Figure 6 shows a good agreement between most reflective satellite scenes of the São Martinho 280 da Serra pixel and ρ ovc . On the other hand, some stations show regularly values of measured reflectances beyond the ρ ovc simulated boundary, as in the exemple of Camborne (Fig. 7, second row). Figure 7 illustrates also how ρ ovc depends on the liquid or ice phase of the cloud, due to their different scattering phase functions.

Comparison of satellite-based estimates with ground-based measurements
Validation results are shown in Table 3, for 15-min averaged DSSI estimates. Satellite-based estimates are obtained with MSG 285 0.6 µm imagery. Results for MSG 0.8 µm imagery show generally lower quality in terms of correlation and STD (Fig. B2).
14 https://doi.org/10.5194/amt-2020-480 Preprint. Discussion started: 12 February 2021 c Author(s) 2021. CC BY 4.0 License. We tested the sensitivity of the DSSI estimates to the cloud phase by using in one case the reference look-up table, featuring a liquid cloud, and for the test case, an ice cloud as described in Section 2.3. Results show only minor differences, pointing out a limited influence of the cloud phase on DSSI estimates (Fig. B3).

290
The results of the method are also compared to satellite-based DSSI products HelioClim3 version 5 (HC3v5) and CAMS Radiation Service (CAMS-RAD) on Table 4. Results for the new HSV method show statistics similar to HC3v5 and CAMS-  Better results from the channel 0.6 µm could be attributed to a smaller influence of the cloud top height, compared to the 0.8 µm channel which is affected by water vapour absorption (Fig. 3). Also, the choice of a spectral linear interpolation 300 between MODIS channels to simulate surface reflectances in SEVIRI channels is supposed to contribute significantly to biases observed in ρ clear simulations, in particular for the 0.8 µm channel with vegetated surfaces due to the red edge spectral pattern.
Such biases could affect significantly DSSI estimates.
The surface reflectivity is lower for shorter wavelengths in general. Selecting a channel for which the surface reflectivity is low will favor a high contrast between clear-sky and overcast scenes, and improve the precision in the computation of the cloud 305 index.
Our ability to reproduce reflectances at the top of the atmosphere in overcast conditions depends also on our knowledge of cloud properties, including their scattering phase function, tridimensional structure and top height.
But the introduction of radiative transfer simulations in the computation of the cloud index also enhances the importance of an accurate calibration of satellite radiance measurements. Using the gain coefficients developed by Doelling et al. (2018) 310 for CERES-SYN1deg instead of EUMETSAT operational coefficients is sufficient to remove most of the mean bias observed between simulations and measurements of reflectances at the top of the atmosphere in clear-sky conditions, for the channel 0.6 µm. Besides, this increases the mean bias for the 0.8 µm. It is possibly due to a weaker compensation of the errors on the calibration and caused by the linear spectral interpolation applied between MODIS channels to reproduce reflexion properties of vegetated surfaces.

315
The simple relationship between the cloud index and the clear-sky index used here explains the significant amount of negative values of DSSI estimates. The improvement of this relation will be the object of a future study.
Finally, the quality of the results depends also on the quality of the clear-sky surface irradiance model : the example of the McClear model shows typical biases of 3 % for the studied stations, when compared to BSRN irradiance data (Fig. C1). The improvement towards a least biased estimation of the downwelling surface solar irradiance based on a cloud index will require 320 better estimates of the attenuation of the solar radiation by the clear atmosphere.

Conclusion and perspectives
A method to compute the cloud index is described in the framework of the development of the future Heliosat-V method for estimating downwelling surface solar irradiance from satellite imagery. The cloud index is built to deal with a single radiometric channel in the spectral range 400-1000 nm. It also does not need archives of data to quantify the cloud effective transmissivity.

325
This approach has advantages. First, the concept of the Heliosat-V cloud index enables the use of imagery from geostationary and non-geostationary platforms, an asset to reach an extended spatial coverage. Moreover, the approach has the potential to deal with long time series of imagery from radiometers characterized by different spectral sensitivities and viewing geometries.
Validation results using SEVIRI imagery show that DSSI can be estimated by a cloud index method that does not rely on archives of imagery, with a quality similar to operational satellite-based data products like CAMS Radiation Service and 330 HelioClim3, in terms of RMSD and correlation. This is an encouraging step toward the application of a Heliosat method to non geostationary satellite sensors. However, we note that there are differentiated errors depending on the spectral channel considered. This could be attenuated notably by an external knowledge on cloud top height and by improving the spectral interpolation of reflexion properties of vegetated surfaces.
To clarify the potential of the method for long time series of imagery, we will need to explore how sensitive to the quality 335 of input data the results are. The knowledge on atmospheric composition in absorbing and scattering species and on surface reflectivity properties is notably lower for past periods like 1980's than for today. Also, the absolute calibration of satellite imagery can be more uncertain, without on-orbit calibrated instruments. Many inputs of the method have very different degrees of quality, depending on the period considered: the composition of the clear-sky atmosphere (aerosols and gases), surface properties, external clear-sky irradiance model. Further work is still to be done on multidecadal time series to study how the 340 quality of such ancillary data affect the estimates of DSSI.
Also, producing global maps of DSSI requires to deal with non geostationary satellite imagery. First tests of the method have been made on the imagery of the Earth Polychromatic Imaging Camera (EPIC) embedded on the DSCOVR platform. They show encouraging results that will be extended and detailed in a future publication.
Global coverage of DSSI information obviously requires also to deal with ocean surfaces and snow covered regions, and this 345 will need to be treated in the future.
Code availability. Excerpts of code are available at https://cloud.mines-paristech.fr/index.php/s/HAWmw7Fs927EtME Data availability. DSSI results derived from the implementation of Heliosat-V for validation on all 11 stations are available at https://cloud. mines-paristech.fr/index.php/s/HAWmw7Fs927EtME, along with simulated and MSG measured reflectances, cloud and clear-sky indices and clear-sky irradiance estimates from the McClear model. The manually filtered clear-sky instants are also provided for all 11 locations.