An alternative 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 of retrieving the cloud index from a large variety of satellite instruments sensitive to reﬂected solar radiation, embedded on geostationary and non-geostationary platforms. The cloud index is a widely used proxy for the effective cloud transmissivity, also called the “clear-sky index”. This study is in the framework of the development of the Heliosat-V 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 (cloud-free) satellite scenes of the Earth’s reﬂectances. Simulations consider the anisotropy of the reﬂectances caused by both surface and atmosphere and are adapted to the spectral sensitivity of the sensor. The anisotropy of ground reﬂectances is described by a bidirectional reﬂectance 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. For 15 min means of DSSI, results from our preliminary implementation of Heliosat-V and ground-based measurements show a bias of 20 W m − 2 , a root-mean-square difference of 93 W m − 2 , and a correlation coefﬁcient of 0.948. The statistics, except for the bias, are similar to operational and corrected satellite-based data products HelioClim3 version 5 and the CAMS Radiation Service.

Abstract. We develop a new way of retrieving the cloud index from a large variety of satellite instruments sensitive to reflected solar radiation, embedded on geostationary and non-geostationary platforms. The cloud index is a widely used proxy for the effective cloud transmissivity, also called the "clear-sky index". This study is in the framework of the development of the Heliosat-V 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 (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. For 15 min means of DSSI, results from our preliminary implementation of Heliosat-V and ground-based measurements show a bias of 20 W m −2 , a root-mean-square difference of 93 W m −2 , and a correlation coefficient of 0.948. The statistics, except for the bias, are similar to operational and corrected satellitebased data products HelioClim3 version 5 and the CAMS Radiation Service.

Introduction
Downwelling surface solar irradiance (DSSI) is one of the Essential Climate Variables defined by the Global Climate Observing 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, according to WMO (2018). DSSI considers the irradiance coming from all directions of the hemisphere above the surface: the irradiance coming from the direction of the Sun, usually referred to as beam horizontal irradiance, plus a diffuse component due to scattering caused by the atmosphere (clouds, gases, aerosols) and reflection by the surface, usually referred to as diffuse horizontal irradiance.
The knowledge of DSSI variations in space and time is of primary importance for various fields such as the Earth sciences, solar energy industries, agriculture, or some medical fields. To meet all these needs, ideal information on DSSI would feature high spatio-temporal resolution, coverage of the entire Earth surface, and the longest time period possible. Long time series of data are notably useful for identifying statistics of long-term inter-annual to multi-decadal variability and possible trends if bias and standard deviation of the error requirements are reached.
Different approaches already exist to produce such DSSI data. Sources of data mainly include ground pyranometric measurements (Driemel et al., 2018), numerical weather prediction modeling (Gelaro et al., 2017;Hersbach et al., 2020), and satellite-based remote sensing (Sengupta et al., 2021). Satellite-based methods are an efficient and accurate way Published by Copernicus Publications on behalf of the European Geosciences Union. B. Tournadre et al.: The Heliosat-V cloud index to produce kilometric and sub-hourly resolved multidecadal time series of DSSI. A more comprehensive review of pros and cons of different methods is notably described in Huang et al. (2019).
The imagery produced by satellite radiometers provides a unique perspective on DSSI. Upwelling radiances coming from each location on Earth are acquired several times per day by a wide set of satellite imagers. This can particularly be achieved thanks to imagers embedded on meteorological geostationary (GEO) and polar-orbiting satellites. Another approach has existed since 2015 thanks to the Deep Space Climate Observatory (DSCOVR) satellite mission: its Lissajous orbit around the L1 Lagrangian point between the Earth and the Sun makes it possible to picture the whole sunlit hemisphere of the planet, with a single satellite radiometer (Marshak et al., 2018;Hao et al., 2020).
Imagery of the Earth produced by satellite sensors has existed for about 6 decades and led early to the development of methods for estimating DSSI (Tarpley, 1979). Today, the information from multi-channel satellite measurements offers the possibility of deriving cloud physical properties and then computing cloud attenuation of the solar radiation with methods like the Fast All-sky Radiation Model For Solar Applications (FARMS) (Xie et al., 2016) and Heliosat-4 (Qu et al., 2017), Zhang et al. (2018), or Hao et al. (2019. Such methods are especially advantageous for highly reflective regions, where clouds are difficult to discriminate from the ground. Nevertheless, they require information on more than one spectral channel, limiting their versatility in the choice of satellite sensor. The use of radiative transfer models and look-up tables is also quite common in the field of satellite-based estimation of DSSI but usually requires pre-existing information on cloud properties or a cloud mask, e.g., ISCCP-F (Zhang, 2004), GEWEX-SRB (Pinker and Laszlo, 1992;Cox et al., 2017), CLARA (Mueller et al., 2009), Cloud_cci (Stengel et al., 2020;Stephens et al., 2001), and SICCS (Greuell et al., 2013).
Another group of methods, labeled "cloud index methods", is able to produce estimates of downwelling surface solar irradiance from the visible imagery of satellite radiometers without external knowledge of cloud physical and optical properties. This gives them potential to retrieve multidecadal time series, including from the imagery of the oldest two-channel sensors like the Meteosat Visible and Infrared Imager (MVIRI). Cloud index methods emerged quite early, notably with the seminal work of Möser and Raschke (1983) and the first Heliosat method (Cano et al., 1986;Cano, 1982). The cloud index quantity derives from the radiances measured by spaceborne sensors and relates them to the extinction of the DSSI caused by clouds. The greater the cloud index, the greater the extinction and the smaller the DSSI. More precisely, the cloud index can be used as an empirical proxy for effective cloud transmissivity. The latter, also named the "clear-sky index" within the scientific commu- 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). nity of solar energy, is defined as the ratio of the all-sky surface irradiance to the clear-sky surface irradiance (Long and Turner, 2008;Beyer et al., 1996), i.e., the DSSI in cloud-free conditions.
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 low and high boundary values. The low boundary value X min is taken as the clear-sky case and the high one X max as the cloudiest 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 relative to the difference between cloudy and clear-sky cases, as illustrated in Fig. 1. We name these variables X as they can be of a slightly different nature from one work to another (reflectance, albedo, radiance, digital count, etc.). 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;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) Dave, 1978), or raw satellite numerical counts (Müller et al., 2015;Perez et al., 2002;Cano et al., 1986), and the way of retrieving X max and X min for the chosen variable.
Very different approaches are used to estimate the upper boundary, but for the lower boundary, "archive-based" methods are used in most of the literature we reviewed: X min is a minimum based on a time series of past satellite imagery. This approach has some drawbacks. Firstly, it is hardly applicable 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. Secondly, systematic underestimations of the lower boundary X min are commonly detected, for example, due to dark shadowing caused by adjacent clouds on the surface and aerosol treatment (Mueller et al., 2015). On the other hand, contamination of X min by clouds in the cloudiest regions can lead to systematic overestimation of X min . Finally, ensuring the observation of clearsky instants by a sufficiently large time window and capturing the temporal variability of X min by a sufficiently small time window is a difficult trade-off that can lead to biases if not well respected.
In this paper, we propose a cloud index method based on radiative transfer modeling as an alternative to the archivebased approach. This exploratory direction aims at reproducing the satellite measurements of reflectances in both clearsky and overcast conditions based on description of surface, clear atmosphere, and cloud properties. Radiative transfer simulations are able to reproduce how TOA reflectances depend on viewing and solar geometries, also with their spectral distribution. In addition, it is possible to provide to the radiative transfer model input data that describe variations in space and time of clear atmosphere composition and of surface properties. Thus, our approach is useful for identifying and quantifying sources of errors in cloud index methods.
With a spectral and angular description, our method is also able 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 polarization present in the ultraviolet, and the absorption of radiation by clouds in the near infrared, the method focuses on satellite imagery in the spectral range 400-1000 nm. This range is wide enough to consider imagers on many meteorological satellites launched since the beginnings of spaceborne Earth observation.
The method foreseen to compute the cloud index of Heliosat-V and eventually the DSSI is described in Sect. 2, along with the protocol of validation. Validation results are presented and discussed in Sect. 3 for simulated reflectances at the top of the atmosphere and for downwelling surface solar irradiance estimates. Section 4 is dedicated to the conclusion and perspectives.

Methods
Previous methods based on archives can avoid dependency on absolute calibration of the original imagery (Müller et al., 2015;Perez et al., 2002) and consider implicitly the anisotropy of X min . The pixel-to-pixel estimation of X min is a surrogate for modeling the influence of viewing geometry on measured reflectances, while the slot-by-slot temporal characterization of X min captures the influence of varying solarviewing geometry for the diurnal cycle of each pixel's reflectivity. The development of an alternative to archive-based approaches means dealing with new issues: a challenge is to reproduce explicitly and accurately the TOA reflectances. For this, input data and models used need to satisfy the requirements for accurate DSSI estimations, as will be discussed hereafter.

The cloud index n
As stated in the introduction, Heliosat-V (HSV) is a method approximating the attenuation of DSSI radiation by clouds with a cloud index, n. Here, the cloud index components are reflectances considered at the TOA and corresponding 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 plane weighted by the spectral response function 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 Sect. 2.3. Because of its definition, the cloud index may also be calculated with radiances. We consider here reflectances in order to visualize the anisotropic nature of different scenes. It also has the advantage of being a normalized quantity, so we can compare results for different radiometric channels and different solar zenith angles.
The relationship between n and DSSI varies slightly from one method to another, in particular for the highest and lowest values of n. The core of the relationship for intermediate 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 Sect. 2.4. The clear-sky index K c is largely used to simplify the reading and is defined as so we can rewrite Eq. (4) as In this paper, we keep the original and simple relation K c = 1 − n introduced by Darnell et al. (1988). Its improvement is beyond the scope of this work but has been explored by various studies, e.g., by Rigollier and Wald (1998) Rigollier et al., 2004), Gupta et al. (2001), Perez et al. (2002), and Zarzalejo et al. (2009), notably to better characterize cloudy situations with n ≈ 1. In the following subsections, we describe the method used to compute ρ clear , ρ ovc , and G c .

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 (Lorente et al., 2018;Veefkind et al., 2016;Stammes et al., 2008). We apply the same approach to satellite radiometers.
Radiative transfer simulations are able to estimate reflected radiation at the top of the atmosphere (TOA) considering the non-Lambertian nature of the atmosphere and of Earth's surfaces. For the implementation of the method applied here, we use the model uvspec within the software package libRadtran (version 2.0.2)  and the one-dimensional solver DISORT (Buras et al., 2011). We chose to use 32 streams for DISORT as a good compromise between time computation and a good angular representativeness of simulated radiances. For the spectral description, radiative transfer simulations are made following the so-called REPTRAN spectral approximation (Gasteiger et al., 2014). This parameterization enables the production of fast computations of radiative transfer adapted to the spectral sensitivity of satellite radiometric channels.
The spectral description of downwelling solar irradiance at the top of the atmosphere is provided by data from Kurucz (1992) for simulating ρ clear . The composition of the atmosphere is provided by time series of total atmospheric columns of ozone and water vapor and partial aerosol optical depths (AODs) from the Monitoring Atmospheric Composition and Climate (MACC) reanalysis (Inness et al., 2013) distributed by the ECMWF. Data from MACC are extracted from the McClear service (http://www.soda-pro.com/ web-services/radiation/cams-mcclear, last access: 16 January 2020). MACC values are originally given on a 3 h time step and with a spatial resolution of about 80 km (Inness et al., 2013;Lefèvre et al., 2013). The McClear service applies to MACC data a bilinear spatial interpolation onto the considered location and a linear interpolation in time to a 1 min time step (Lefèvre et al., 2013). The atmospheric abundance profiles of O 2 , CO 2 , and NO 2 are kept to the fixed values of the Air Force Geophysics Laboratory (AFGL) midlatitude summer profile (Anderson et al., 1986) along with the temperature, pressure, and air density profiles.
Partial AODs from MACC are provided at the wavelength 550 nm for five types of aerosols (black carbon, dust, sea salt, organic matter, sulfate). Even though two supplementary classes "ammonium" and "nitrate" are now included in the Copernicus Atmospheric Monitoring Service (CAMS) reanalysis, these do not impact the method proposed here and were, thus, not considered.
An algorithm developed by Lefèvre et al. (2013) translates MACC partial aerosol optical depth information into aerosol mixtures designed for the Optical Properties of Aerosols and Clouds (OPAC) software package (Hess et al., 1998). These mixtures are associated with 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.
An important component to simulate ρ clear is the reflection properties of surfaces. The impact of the anisotropy of surface reflectance has notably been shown for estimates of a cloud index derived from measurements of the ultraviolet/visible Global Ozone Monitoring Experiment 2 (GOME-2) and OMI by Lorente et al. (2018). This study also highlights the improvement of simulated shortwave clear-sky reflectances at the TOA when using a model of the bidirectional reflectance distribution function (BRDF) parameterized with data derived from the Moderate Resolution Imaging Spectroradiometer (MODIS) spaceborne instruments.
Here, we describe reflective properties of land surfaces with the RossThick-LiSparse (Ross-Li) model of the BRDF (Roujean et al., 1992;Lucht et al., 2000). It is then possible to consider the variations of the surface reflectance depending on the 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. 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 done with the imagery produced by MODIS embedded on the 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). 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, a 16 d average, and seven spectral channels, including four 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 0.85 µm, the interpolation is made between parameters at MODIS channels 0.86 and 1.24 µm. A summary of inputs used to produce simulations of clear-sky reflectances is provided in Table 2.

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;Mueller and Träger-Chatterjee, 2014). 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) with an empirical relation based on the work of Taylor and Stowe (1984). It considers 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 Fig. 2, we represent two-dimensional 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 the distribution is observed ( Fig. 2a and b), corresponding only to different viewing geometries defined by a linear mesh grid in the cosine of the viewing zenith angle (θ v ) and difference φ of the 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 Trishchenko et al. (2001). The simulations for a low thick cloud (cloud top height (CTH) at 500 m) and a high thick cloud (CTH at 15 km) show in general good agreement (Fig. 2c), except in absorbing bands of O 2 (mainly at 690 nm, O 2 -B band, and 762 nm, O 2 -A band) and H 2 O (mainly at 725, 820, and 950 nm) and for short wavelengths where scattering becomes increasingly significant (e.g., Jin et al., 2011). For these wavelengths, the TOA reflectances with low clouds can be much lower than for high clouds for a given cloud optical thickness. However, outside these specific spectral regions, the height of clouds will not affect significantly the results of the method.
An alternative way is therefore to produce look-up tables (LUTs) from radiative transfer simulations, an approach notably applied in the framework of the HelioMont cloud index method (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 about 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 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, geometries, 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 ρ ovc = 1 2 ρ ovc,high + ρ ovc,low , where ρ ovc,high and ρ ovc,low are reflectances, 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 eight columns" ice crystal habit and the "severe" degree of roughness, which are notably used for the description of ice clouds in the lookup 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 version 3 of the McClear model (Gschwind et al., 2019). The McClear model is a fast and accurate model that provides clear-sky estimation of DSSI with an absolute bias below 21 W m −2 and a standard deviation error below 25 W m −2 for six stations part of the reference Baseline Surface Radiation Network (BSRN) (Ohmura et al., 1998;Driemel et al., 2018), namely, Brasilia, Carpentras, Palaiseau, Payerne, Sede Boker, and Tamanrasset. The McClear model was fed with the partial optical depths at 550 nm for black carbon, dust, sea salt, organic matter, and sulfate from MACC reanalysis. It is also fed by water vapor atmospheric total columns and the ozone total columns provided by ECMWF. Data were downloaded from the McClear web service (http://www.soda-pro.com/ web-services/radiation/cams-mcclear, last access: 16 January 2020).

Setup and datasets for validation
The method has been tested on images from the Spinning Enhanced Visible and Infra-Red Imager (SEVIRI) 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 and 0.8 µm for the year 2011 and for 11 locations in the field of view of the satellite, corresponding to locations of pyranometric in situ sensors from the BSRN. We use the calibration gains provided by EUMETSAT that operates MSG. For sensors with a linear count response like MSG/SEVIRI (Doelling et al., 2018), the radiance L sat is related to digital count C via L sat = g(C − C 0 ), where C 0 is the so-called space count.
To study the validity of the method, we compare DSSI estimates from MSG satellite measurements with pyranometric DSSI data retrieved from BSRN measurement stations. Considered stations are listed in Table 3 and displayed in the MSG field of view in Fig. 4. Only the highest-quality BSRN measurements of surface irradiance are used, having passed a quality check (Lefèvre et al., 2013). Figure 5 shows the time series when data are considered valid, for each station.
We also compare the results of our method to operational satellite-based products of surface irradiance. For this, we use data from the HelioClim3 version 5 (HC3v5) and CAMS Radiation (CAMS-RAD) DSSI databases. Both are  Lefèvre et al. (2013) and OPAC Temperature and pressure profiles AFGL midlatitude summer derived from the imagery of the SEVIRI 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. Both products and their descriptions are provided by the SoDa service (http://www.soda-pro.com/, last access: 9 September 2019).
As this work is exploratory on a new method, we limit ourselves to conservative situations with solar zenith angles lower than 80 • , covering most cases. For higher angles, some effects not considered by the method can occur, including shadowing and high parallax effects.

Validity of cloud index components
The validity of cloud index components, ρ sat , ρ clear , and ρ ovc , defines the uncertainty of n. From Eq. (3), the uncertainty in the cloud index can be written as δn = ∂n ∂ρ sat δρ sat + ∂n ∂ρ clear δρ clear + ∂n ∂ρ ovc δρ ovc .
where = ρ ovc − ρ clear . It appears that the "clear-sky error" (1−n) δρ clear will be more significant in clear-sky conditions (i.e., n is close to 0), and the "overcast-sky error" nδρ ovc will be more important in overcast conditions (i.e., n is close to 1). Besides, the error in the cloud index will be inversely proportional to , the difference between overcast and clear-sky TOA reflectances. Because of this relationship between the errors in cloud index and reflectances, the discussions in this section are focused on absolute values of reflectance errors.

Measured reflectances at the top of the atmosphere
A potentially important source of the measurement error δρ sat comes from the calibration gain. The operational calibration gains that we use in this paper have a claimed uncertainty of around 4 % (EUMETSAT, 2019). On the other hand, Hewison et al. (2020) assert that the alternative method by Doelling et al. (2018), used for GSICS-corrected computation of calibration gain, limits its bias to below 1 %. The use of optimal calibration is beyond the scope of our work. Still, we compared gain coefficients proposed by EU-METSAT g EUM to 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 relative disagreement, calculated as (g EUM − g D2018 )/g D2018 , of about −9 % for 0.6 µm and −8 % for 0.8 µm during this period (also illustrated in Fig. A1). We expect that these errors will affect with the same magnitude the agreement between numerical simulations and measurements of clear-sky TOA reflectances. This underlines that an accurate source of absolute calibration is important for the Heliosat-V method.  Fig. 6 as two-dimensional reflectance histograms. For the 0.6 and 0.8 µm channels, correlation coefficients are both higher than 0.9, but the correlation is much better for 0.6 µm, with a value of 0.974. This means that the variability of ρ clear is significantly better represented, with almost 95 % of the total variance, for 0.6 µm than for 0.8 µm, with 82 % of the total variance. The root-mean-square difference (RMSD) between the simulated reflectance and measured reflectance in the clear-sky conditions is 0.03 (15 %) for the 0.6 µm channel and 0.04 (12 %) for the 0.8 µm channel. The bias is 0.02 (10 %) and −0.02 (−7 %) for the 0.6 and 0.8 µm channels, respectively, contributing a big part to the RMSD.
The standard deviation of the difference (SD) is 0.02 for the 0.6 µm channel and 0.04 for the 0.8 µm channel. Both higher bias and SD for the 0.8 µm will contribute to lowering the precision in the calculation of the cloud index based on this channel compared to 0.6 µm. When studying station by station, the highest absolute standard deviation of the difference between simulations and measurements is reached for Sede Boker with 0.03, while the lowest is reached for Tamanrasset with 0.008. For bias, the most significant mean values reach +0.035 for 0.6 µm (Payerne) and −0.07 for 0.8 µm (Camborne) (see also Fig. B1).
Using the gain coefficients developed by Doelling et al. (2018) for CERES-SYN1deg instead of EUMETSAT operational coefficients is sufficient to remove most of the mean bias observed between simulations and measurements of ρ clear for the channel 0.6 µm. Besides, it increases the mean bias for the 0.8 µm channel.
It is worth noting that we use MCD43C1v6 BRDF data regardless of their quality flags. We observe though that keeping only the highest-quality data improves significantly statistics (Fig. B2). Also, the choice of a spectral linear interpolation 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 (low reflectivity below around 700 nm, high reflectivity above around 750 nm). Another part of the bias, difficult to quantify, is linked to the accuracy of the calibration of satellite measurements. Figure 7 shows the diurnal variations of measured and simulated reflectances for the SMS and CAM stations. Both SMS and CAM are surrounded mainly by various types of vegetation and some urban area in the case of CAM (Fig. B3). We observe that simulations are able to reproduce partly the diurnal variability observed in clear-sky conditions (also refer to Fig. B4 for channels 0.6 and 0.8 µm under different surface conditions). In Fig. 8, we compare ρ clear values to the surface reflectance ρ surface , computed with the RossThick-LiSparse model applied to BRDF parameters derived from the MODIS 646 nm channel and using the viewing and solar geometries considered. Firstly, we see that ρ clear values are significantly higher than ρ surface , with a different diurnal pattern. This shows the importance of considering the atmosphere anisotropic reflectance to reproduce TOA reflectances. We can also see the contribution from the surface anisotropy in the ρ clear simulations. This appears in particular close to the backscattering direction, where surface reflectance is enhanced: around noon in Camborne and the morning in São Martinho da Serra.
For CAM, some higher values of ρ clear are observed in January. This can be attributed to high aerosol optical depth during this period, as illustrated in Fig. 9. It shows that ρ clear is not only sensitive to time variations of surface properties, but also to atmospheric composition changes.

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 comparison to satellite measurements as the occurrence 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 (orange dots in Fig. 7). In the first row of Fig. 7, we can see that some patterns are similar in simulated ρ ovc and the 99th percentile of measurements over the São Martinho da Serra pixel: in the forward-scattering conditions (evening on the western edge of the Meteosat disc), both agree on increased values of ρ ovc . On the other hand, some stations show regular values of measured reflectances beyond the ρ ovc simulated boundary, as in the example of Camborne (Fig. 7, second row). Figure 7 also illustrates how ρ ovc depends on the liquid or ice phase of the cloud, due to their different scattering phase functions. Our ability to reproduce reflectances at the top of the atmosphere in overcast conditions depends therefore on our knowledge of cloud properties, including their scattering phase function and top height.
Other effects like the tridimensional structure of clouds probably explain part of the discrepancies between measurements and plane-parallel simulations in overcast conditions (Horvath and Davies, 2004).

Difference between simulated overcast and clear-sky reflectances
The difference between overcast and clear-sky reflectances is bigger when the overcast reflectance is relatively low and clear-sky reflectance is relatively high. High values of mean a good quality of cloud index estimation (cf. Eq. 9). We study the dependencies of with the simulated reflectances to identify conditions that will cause high uncertainties in the computation of the cloud index. In general, we observe that the computed value of is higher for the 0.6 µm channel than for 0.8 µm as a combination of surface, cloud, and clear atmosphere spectral signatures. This is illustrated in Fig. 10 for stations SMS and CAM. We observe however for the desert stations TAM and SBO that both channels present similar values of (Fig. B5). also depends on the viewing and solar geometries because of ρ ovc and ρ clear different angular signatures. It leads, for example, for the SMS station and channel 0.8 µm to low values of in January mornings and high values of in the evening, which can be explained by the strong forward scattering of clouds occurring in these conditions.

Comparison of satellite-based estimates of DSSI to ground-based measurements
Validation results are shown in Table 4 for 15 min-averaged DSSI estimates. Satellite-based estimates are obtained with MSG 0.6 µm imagery. Results for MSG 0.8 µm imagery show generally lower quality in terms of correlation and SD, as shown in Fig. 11. The simple relationship between the cloud index and the clear-sky index used here explains the significant number of negative values of DSSI estimates. The improvement of this relation will be the object of a future study.
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 Sect. 2.3. Results show only minor differences, pointing out a limited influence of the cloud phase on DSSI estimates (Fig. B6).
Finally, the quality of the results also depends on the quality of the clear-sky surface irradiance model. Gschwind et al. (2019) report, for example, relative mean biases of the Mc-Clear model from −3.6 % (Barrow, Alaska, USA) to +3.2 % (Payerne, Switzerland) when compared to BSRN irradiance measurements. The improvement towards a least-biased estimation of the downwelling surface solar irradiance based on a cloud index will require better estimates of the attenuation of the solar radiation by the clear atmosphere.

Comparison of satellite-based estimates of DSSI to the HelioClim3 and CAMS radiation products
The results of the method are also compared to satellitebased DSSI products HelioClim3 version 5 (HC3v5) and CAMS Radiation Service (CAMS-RAD) in Table 5. Results for the new HSV method show statistics similar to HC3v5 and CAMS-RAD, for both estimates based on the 0.6 and 0.8 µm channels, in terms of correlation and of SD. One may note very low values of bias for operational products. This is expected because CAMS-RAD and HC3v5 estimates are calibrated with DSSI measurements from a similar set of BSRN stations.
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 vapor absorption (Fig. 3). Biases discussed for the computation of clear-sky and overcast TOA reflectances could also affect significantly DSSI estimates.

Conclusion and perspectives
Heliosat-V is a cloud-index method for estimating downwelling surface solar irradiance from satellite imagery. In the framework of its development, we proposed an alternative way of retrieving the components of the cloud index, this index being used to quantify the attenuation of DSSI by clouds. The method takes advantage of radiative transfer modeling to provide versatility to the concept of the cloud index. Thus, it does not need archives of data to quantify the cloud effective transmissivity. It is applicable for optical sensors on geostationary and non-geostationary orbits; flexible for future improvements to describe surface, clear atmosphere, and clouds; and investigates physical solutions for limitations observed in previous cloud index methods.
It is built to deal with a single radiometric channel in the spectral range 400-1000 nm. 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 the CAMS Radiation Service and 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 external knowledge of cloud top height and by improving the spectral interpolation of reflection properties of vegetated surfaces.
To clarify the potential of the method for long time series of imagery, we will need to explore how sensitive the results are to the quality of input data. The knowledge of atmospheric composition in absorbing and scattering species and on surface reflectivity properties is notably lower for past periods like the 1980s 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, and external clear-sky irradiance model. Further work is still to be done on multidecadal time series to study how the quality of such ancillary data affect the estimates of DSSI.  Also, producing global maps of DSSI requires users to deal with non-geostationary satellite imagery. The first tests of the method have been done 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 also requires users to deal with ocean surfaces and snow-covered regions, and this will need to be treated in the future.       Figure B6. Impact of the cloud phase on DSSI estimates. Two-dimensional histogram of satellite-based DSSI estimates from the Heliosat-V method versus ground-based BSRN measurements for the MSG 0.6 µm channel. The liquid cloud look-up table of overcast-sky TOA reflectances is replaced for the ice cloud LUT.
Data availability. DSSI results derived from the implementation of Heliosat-V for validation on all 11 stations are publicly available at https://doi.org/10.23646/tg31-1452 (Tournadre and Gschwind, 2022), 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.
Author contributions. Conceptualization was done by BT, YMSD, and PB. Investigation, validation, and writing of the original draft were done by BT. Visualization was done by BT and PB. BT and BG did the software, with contributions from PB. Supervision was by PB and BG. Methodology was done by BT, PB, and YMSD. All the authors brought contributions to the writing process. XC and RAES brought contributions to the writing process, as did other coauthors.
Competing interests. The Transvalor company that partially funds this work distributes HelioClim3 and CAMS Radiation Service.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Acknowledgements. We would like to thank the operators of the BSRN stations for producing and providing valuable validation measurements. We are grateful to the libRadtran team for their open radiative transfer model and interactions. The MODIS MCD43C1 version 6 data product was retrieved from the online NASA Earthdata Search, courtesy of the NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC) and USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, https://search.earthdata.nasa.gov/ (last access: 8 July 2019). We also thank EUMETSAT and the ECMWF for providing, respectively, Meteosat-9 and MACC reanalysis aerosol, water vapor, and ozone data. We used imagery from the NASA Worldview application (https://worldview.earthdata.nasa.gov, last access: 24 February 2022), part of the NASA Earth Observing System Data and Information System (EOSDIS). We thank Lionel Ménard for his help in the distribution of data linked to the paper (see the Data availability section). We thank the three anonymous reviewers for their helpful comments and suggestions. We finally thank Olivier Boucher, Alba Lorente, Review statement. This paper was edited by Piet Stammes and reviewed by three anonymous referees.