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 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⁻², a root-mean-square difference of 93 W m⁻², and a correlation coefficient 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.

1 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 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 multi-decadal 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 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; 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 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

$$\begin{array}{}\text{(1)}& n={\displaystyle \frac{X_{\mathrm{sat}}-X_{\min }}{X_{\max }-X_{\min }}}.\end{array}$$

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) (Herman et al., 2018; 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 clear-sky 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 archive-based approach. This exploratory direction aims at reproducing the satellite measurements of reflectances in both clear-sky 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.

2 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 solar-viewing 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.

2.1 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

$$\begin{array}{}\text{(2)}& \mathit{\rho }={\displaystyle \frac{\mathit{\pi }L}{E_{\mathrm{0}}\cos \left({\mathit{\theta }}_{\mathrm{s}}\right)}},\end{array}$$

with L the upwelling radiance at TOA for a given spectral channel, E₀ 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₀ 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

$$\begin{array}{}\text{(3)}& n={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{sat}}-{\mathit{\rho }}_{\mathrm{clear}}}{{\mathit{\rho }}_{\mathrm{ovc}}-{\mathit{\rho }}_{\mathrm{clear}}}},\end{array}$$

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 values of n usually follows

$$\begin{array}{}\text{(4)}& G=G_{\mathrm{c}}(\mathrm{1}-n),\end{array}$$

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

$$\begin{array}{}\text{(5)}& K_{\mathrm{c}}={\displaystyle \frac{G}{G_{\mathrm{c}}}},\end{array}$$

so we can rewrite Eq. (4) as

$$\begin{array}{}\text{(6)}& G=G_{\mathrm{c}}{K}_{\mathrm{c}}.\end{array}$$

In this paper, we keep the original and simple relation $K_{\mathrm{c}}=\mathrm{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) (reported in 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.

2.2 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) (Emde et al., 2016) 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₂, CO₂, and NO₂ 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.

Figure 2 (a, b) Distributions of simulated TOA reflectance spectra in overcast conditions ρ_ovc for the different viewing geometries in the look-up table and for a solar zenith angle of 30^∘, with a thick liquid cloud (COT =150). (a) CTH =15 km; cloud base height (CBH) =2 km. (b) CTH =0.5 km; CBH =0.2 km. (c) Error on ρ_ovc caused by a misattribution of cloud height to the “low thick cloud” category. Green, red, and blue arrows indicate spectral regions with main absorption features from O₃, O₂, and H₂O, respectively.

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.

2.3 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₂ (mainly at 690 nm, O₂-B band, and 762 nm, O₂-A band) and H₂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

$$\begin{array}{}\text{(7)}& {\mathit{\rho }}_{\mathrm{ovc}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\rho }}_{\mathrm{ovc},\mathrm{high}}+{\mathit{\rho }}_{\mathrm{ovc},\mathrm{low}}\right),\end{array}$$

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 look-up table of the MODIS collection 6 cloud product (Amarasinghe et al., 2017).

Gueymard (2018)Shettle (1990)

Table 1 Characteristics of the look-up table of cloudy TOA reflectances.

(Gasteiger et al., 2014)Kurucz (1992)Lefèvre et al. (2013)

Table 2 Characteristics of the simulated reflectances at the top of the atmosphere in clear-sky conditions.

2.4 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⁻² and a standard deviation error below 25 W m⁻² 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).

2.5 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_{\mathrm{sat}}=g(C-C_{\mathrm{0}})$, where C₀ 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.

Table 3 List of BSRN stations used for validation.

Figure 3 Colored lines: spectral response functions of different sensors in the spectral range considered by Heliosat-V. Gray lines: TOA reflectance spectra of typical scenes with a high- (dashed line) and low-altitude (solid line) thick cloud.

Figure 4 BSRN ground stations used for validation in this study, in the field of view of Meteosat Second Generation (0.6 µm channel).

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 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.

3 Results and discussions

3.1 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

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathit{\delta }n& =\left({\displaystyle \frac{\partial n}{\partial {\mathit{\rho }}_{\mathrm{sat}}}}\right)\mathit{\delta }{\mathit{\rho }}_{\mathrm{sat}}+\left({\displaystyle \frac{\partial n}{\partial {\mathit{\rho }}_{\mathrm{clear}}}}\right)\mathit{\delta }{\mathit{\rho }}_{\mathrm{clear}}\\ & +\left({\displaystyle \frac{\partial n}{\partial {\mathit{\rho }}_{\mathrm{ovc}}}}\right)\mathit{\delta }{\mathit{\rho }}_{\mathrm{ovc}}.\end{array}\end{array}$$

This leads to

$$\begin{array}{}\text{(9)}& \mathit{\delta }n={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }}}\left(\mathit{\delta }{\mathit{\rho }}_{\mathrm{sat}}-\left(\mathrm{1}-n\right)\mathit{\delta }{\mathit{\rho }}_{\mathrm{clear}}-n\mathit{\delta }{\mathit{\rho }}_{\mathrm{ovc}}\right),\end{array}$$

where $\mathrm{\Delta }={\mathit{\rho }}_{\mathrm{ovc}}-{\mathit{\rho }}_{\mathrm{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.

3.1.1 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 EUMETSAT 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_{\mathrm{EUM}}-g_{\mathrm{D}\mathrm{2018}})/g_{\mathrm{D}\mathrm{2018}}$, 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.

3.1.2 Simulated reflectances at the top of the atmosphere in clear-sky conditions ρ_clear

As an intermediate assessment, simulated clear-sky reflectances at the top of the atmosphere (TOA) ρ_clear are compared to satellite measurements. Cloudy instants are manually filtered out of the satellite time series. Results comprising all the manually filtered clear-sky instants in 2011 for all 11 sites are shown in Fig. 6 as two-dimensional reflectance histograms.

Figure 5 Time series used for the 15 min mean statistics between satellite estimates and quality-checked BSRN measurements during the year 2011. In parentheses: percentage of data conserved.

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 6 Simulation of clear-sky reflectances at the TOA (ρ_clear) for MSG 0.6 µm (a) and 0.8 µm (b) spectral channels compared to actual satellite measurements. Represented data include simulations and measurements for all 11 locations, for the year 2011.

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.

Figure 7 Simulated and measured reflectances at the top of the atmosphere above the São Martinho da Serra (Brazil, a–c) and Camborne (United Kingdom, d–f) locations, for the MSG 0.6 µm channel and for the January, May, and September calendar months. Gray plus signs: MSG measurements (2011, Meteosat-9). Green asterisks: reflectances in overcast conditions ρ_ovc, derived from the liquid cloud look-up table. Blue asterisks: same from the ice-cloud look-up table. Purple asterisks: reflectances in clear-sky conditions ρ_clear, derived from radiative transfer simulations. Yellow and orange dots are, respectively, hourly percentiles 1 and 99 of MSG satellite measurements from years 2011 to 2019.

Figure 8 Comparison between simulations of clear-sky reflectances at the top of the atmosphere for the MSG 0.6 µm channel (ρ_clear, blue plus signs) and corresponding surface reflectances computed with the RossThick-LiSparse model applied to MODIS MCD43C1v6 BRDF parameters for the channel 646 nm (ρ_surface, red plus signs) for 5 d in June 2011. (a) Camborne station (CAM); (b) São Martinho da Serra station (SMS).

Figure 9 Blue plus signs: simulated reflectances at the top of the atmosphere in clear-sky conditions ρ_clear in January 2011 at Camborne station (CAM) and for MSG 0.6 µm. Red line: aerosol optical depth at 635 nm used for simulations.

Figure 10 Difference between simulated reflectances at the top of the atmosphere in overcast and clear-sky conditions $\mathrm{\Delta }={\mathit{\rho }}_{\mathrm{ovc}}-{\mathit{\rho }}_{\mathrm{clear}}$ for the São Martinho da Serra (Brazil, a–c) and Camborne (United Kingdom, d–f) locations and for the January, May, and September calendar months (a, d; b, e; c, f). In blue dots: MSG 0.6 µm channel; in red dots: MSG 0.8 µm channel.

3.1.3 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).

3.1.4 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.

3.2 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.

Figure 11 Two-dimensional histograms of satellite-based DSSI estimates from the Heliosat-V method versus ground-based BSRN measurements for the MSG 0.6 µm channel (a) and 0.8 µm channel (b).

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 McClear 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.

Table 4 Validation results for 15 min means of all-sky DSSI for the year 2011. Results based on the imagery of the Meteosat-9/SEVIRI 0.6 µm channel.

3.3 Comparison of satellite-based estimates of DSSI to the HelioClim3 and CAMS radiation products

The results of the method are also compared to satellite-based 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.

Table 5 Comparison between validation results of HSV and those of HC3v5 and CAMS-RAD, each one versus BSRN measurements. Statistics on 15 min means of DSSI for the stack of 11 stations and the year 2011. N=135107; BSRN mean = 424 W m⁻².

4 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.

Appendix A: Methods

A1 Setup of validation

Figure A1 First two rows: calibration gains provided by EUMETSAT (black stars) and by CERES-SYN1deg (Doelling et al., 2018) (red stars) for the 0.6 µm channel (first row) and 0.8 µm (second row) of the SEVIRI sensor aboard Meteosat-9 between 2011 and 2019. Third row: ID of the operational satellite at longitude 0^∘.

Appendix B: Results

B1 Simulated TOA clear-sky reflectances ρ_clear

Figure B1 Relative mean bias errors of simulated clear-sky reflectances at the top of the atmosphere ρ_clear (simulation measurement$/$measurement) for channels 0.6 µm (a) and 0.8 µm (b) and for each BSRN site.

Figure B2 Simulation of clear-sky reflectances at the top of the atmosphere (ρ_clear) for MSG 0.6 µm (a) and 0.8 µm (b) spectral channels compared to actual satellite measurements. The comparison is done for all 11 locations, for the year 2011. Only instants with BRDF data of the best quality are used (quality flag 0 of MCD43C1, “Best quality, 100 % with full inversion”).

Figure B3 Land cover types around measurement stations São Martinho da Serra (Brazil, upper panel) and Camborne (United Kingdom, lower panel) for 2011. In red: urban and built-up lands; in gray: croplands/natural vegetation mosaics; in light yellow: croplands; in dark yellow: savanna; in beige: grasslands; in blue: water bodies. Data from Terra + Aqua MODIS product MCD12Q1 version 6, following the International Geosphere-Biosphere Programme classification scheme (credit: NASA Worldview).

Figure B4 Comparison between simulated (red plus signs) and measured reflectances (blue plus signs) at the top of the atmosphere for 1 d in clear-sky conditions, for the 0.6 µm (first column) and 0.8 µm channels (third column). NDVI computed from satellite measurements is shown in the second column. Rows from top to bottom: locations of the São Martinho Da Serra, Camborne, Payerne, and Tamanrasset BSRN stations.

Figure B5 Difference between simulated reflectances at the top of the atmosphere in overcast and clear-sky conditions $\mathrm{\Delta }={\mathit{\rho }}_{\mathrm{ovc}}-{\mathit{\rho }}_{\mathrm{clear}}$ for the Tamanrasset (Algeria, a–c) and Sede Boker (Israel, d–f) locations and for the January, May, and September calendar months (a, d; b, e; c, f). In blue dots: MSG 0.6 µm channel; in red dots: MSG 0.8 µm channel.

Table B1 Validation results for 15 min means of the all-sky DSSI for the year 2011. Results based on the imagery of the Meteosat-9/SEVIRI 0.8 µm channel.

B2 Comparison of satellite-based estimates with ground-based measurements

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.