Satellite Retrieval of Downwelling Shortwave Surface Flux and Di ﬀ use Fraction under All Sky Conditions in the Framework of the LSA SAF Program (Part 1: Methodology)

: Several studies have shown that changes in incoming solar radiation and variations of the di ﬀ use fraction can signiﬁcantly modify the vegetation carbon uptake. Hence, monitoring the incoming solar radiation at large scale and with high temporal frequency is crucial for this reason along with many others. The European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Satellite Application Facility for Land Surface Analysis (LSA SAF) has operationally disseminated in near real time estimates of the downwelling shortwave radiation at the surface since 2005. This product is derived from observations provided by the SEVIRI instrument onboard the Meteosat Second Generation series of geostationary satellites, which covers Europe, Africa, the Middle East, and part of South America. However, near real time generation of the di ﬀ use fraction at the surface level has only recently been initiated. The main di ﬃ culty towards achieving this goal was the general lack of accurate information on the aerosol particles in the atmosphere. This limitation is less important nowadays thanks to the improvements in atmospheric numerical models. This study presents an upgrade of the LSA SAF operational retrieval method, which provides the simultaneous estimation of the incoming solar radiation and its di ﬀ use fraction from the satellite every 15 min. The upgrade includes a comprehensive representation of the inﬂuence of aerosols based on physical approximations of the radiative transfer within an atmosphere-surface associated medium. This article explains the retrieval method, discusses its limitations and di ﬀ erences with the previous method, and details the characteristics of the output products. A companion article will focus on the evaluation of the products against independent measurements of solar radiation. Finally, the access to the source code is provided through an open access platform in order to share the expertise on the satellite retrieval of this variable with the community.


Introduction
The downwelling surface shortwave flux (DSSF) refers to the radiative energy in the wavelength interval [0.3 µm, 4.0 µm] that reaches the Earth's surface per time unit and per surface area unit. DSSF essentially depends on the solar zenith angle, cloud coverage, aerosol load and type, and, to a lesser extent, gas absorption (water vapor, ozone content) and surface albedo. A study conducted by reference [1] showed the importance of clouds for the DSSF and the Earth's energy budget at the daily accumulated values of DSSF. Their distribution is performed in near real time with a timeliness of three hours, counted after the last satellite observation. The production is done every 30 min for the instantaneous MDSSF product, and daily for the DIDSSF product. Both products consider clear and cloudy skies, and provide total shortwave fluxes but without distinction of the diffuse fraction. Although these products have proven to be of high quality [31][32][33], they still show some limitations under clear-sky conditions, as they rely on the hypothesis of temporally and spatially constant aerosol properties [34]. The importance of aerosols on the DSSF estimation was underlined in numerous studies, especially in regions with high aerosol concentrations [35][36][37][38]. The main difficulty of using aerosols to derive DSSF was the general lack of accurate and dynamic information on the aerosol particles in the atmosphere. This limitation is nowadays less important thanks to the improvements in atmospheric numerical models combined with the use of an increase amount of satellite measurements [39][40][41]. Satellite products on aerosol properties (e.g., reference [42]) could represent alternatives. However, it was proved that their accuracy was still not enough to calculate accurately the radiative forcing due to aerosol at the surface [43], which is intimately related to DSSF.
This study describes the generation of a new LSA SAF product, referenced as LSA-207 and called the MSG-derived Downwelling Surface Shortwave Flux-Total and Diffuse (MDSSFTD). The two main advantages of this upgrade of the MDSSF operational product are (i) the improvement of the DSSF estimation (especially thanks to the consideration of the influence of aerosols on the atmospheric transmittance in case of clear-sky and cloudy-sky conditions), and (ii) the retrieval of the diffuse fraction at the surface level for all sky situations.
This article is organized as follows. Section 2 gives an overview of the methodology as well as the workflow for the estimation of the DSSF flux and its related fraction of diffuse radiation. Section 3 describes the input data that are needed for the calculation of these variables. Section 4 provides detailed descriptions of the algorithm under clear-sky and cloudy-sky conditions. Section 5 describes the characteristics of the product outputs and lists the known limitations before concluding in Section 6.

Downwelling Surface Shortwave Radiation
The downwelling surface shortwave radiation flux, E ↓ , is defined as the integral of the spectral irradiance E(λ) over the wavelength interval [λ 1 = 0.3 µm, λ 2 = 4 µm] DSSF is expressed in Watts per square meter (W/m 2 ). The spectral irradiance at wavelength λ is the hemispherical angular integral of the downwelling spectral radiance L(λ, θ s , ϕ s ) over the solar azimuth and zenith angles, noted respectively as ϕ s and θ s , weighted by the cosine of the solar zenith angle L(λ, θ s , ϕ s ) cos(θ s ) sin(θ s )dθ s dϕ s .
It includes contributions owing to the direct solar radiation (E dir ) attenuated by the atmosphere as well as diffuse radiation coming from other directions (E di f ): The DSSF can be expressed as where E 0 is the solar constant, θ s is the solar zenith angle, v(t) is the Sun-Earth distance factor, and T is the total effective transmittance of the entire atmosphere (including the cloud attenuation in the case of cloudy skies). The latter quantity can be decomposed into multiple attenuation processes of individual atmospheric components such as absorption by gases (T gases ), Rayleigh scattering (T Ray ), aerosols extinction (T aer ), and cloud extinction (T C ). The total transmittance is further detailed in Section 4.1.
The total DSSF, E ↓ , includes the contributions owing to the direct solar radiation attenuated by the atmosphere, as well as diffuse radiation. It can be approximated as where T dir and T di f are the effective transmittances of the atmosphere system (cloud-gas-aerosol) for the direct and diffuse radiation, respectively, and E dir and E di f are the corresponding direct and diffuse radiation that compose the total DSSF.

Diffuse Fraction
The diffuse fraction, f d , is the ratio between the diffuse incoming radiation and the total incoming radiation and is calculated as follows This ratio is calculated at the surface level and for the visible broadband domain ([λ 1 = 0.3 µm, λ 2 = 4 µm]).

Equivalent Aerosol Optical Depth at 550 nm
The equivalent Aerosol Optical Depth (AOD) at 550 nm is related to the extinction of the solar beam by aerosol particles in the wavelength 550 nm. It takes into account the extinctions from different aerosols components that may co-exist in the same aerosol mixture (sulfates, sea salt, black carbon, etc.). Therefore, the equivalent AOD is linked to the extinction of the total column of atmosphere in presence of an unclassified aerosol mixture at 550 nm.

Opacity Index
The opacity of the atmosphere is by construction the opposite of the clearness. The Opacity Index (OI) is defined as where Kt is the clearness index, in turn defined as the fraction of the solar radiation that is transmitted through the atmosphere to finally strike the surface of the Earth. The clearness index (Kt) is high under clear and sunny conditions, and low under cloudy conditions. The clearness index can be calculated as the incident radiation at the surface level divided by the extraterrestrial radiation at the top of the atmosphere Finally, it stands from Equation (4) that Kt is, in our case, nothing else than T, the total transmittance of the atmosphere. Analogously, OI, which characterizes the opacity of the total column of the atmosphere, could be also referred as the total extinction of the atmosphere.

Overview of the Retrieval Method
An overview of the retrieval method is presented in Figure 1 with a discrimination between the clear and cloudy cases. The list of the retrieved variables is in yellow color. The main two variables are the total DSSF and the diffuse fraction. However, two additional quantities characterizing the atmosphere are also estimated: an equivalent AOD at 550 nm and the Opacity Index. These four variables are all useful to understand the mutual dependencies between them and hence they will all be made available to the user and embedded in the final product (Section 5.1) together with the associated quality flags. Each output estimate is delivered on a pixel-wise basis that propagates the processing quality information of each module (clear-sky and cloudy-sky). The green color is used for the inputs, which will be further explained in Section 3.
Remote Sens. 2019, 11, x FOR PEER REVIEW 5 of 26 processing quality information of each module (clear-sky and cloudy-sky). The green color is used for the inputs, which will be further explained in Section 3. The method used for the retrieval of DSSF and is composed of three main steps (see the blue circles in Figure 1): first, the estimations of the individual atmospheric contributions and the total atmospheric transmittance (Section 4.1); second, the parameterization for clear-sky conditions (Section 4.2); and third, the parameterization for cloudy-sky conditions (Section 4.3). The approach chosen for the retrieval is based on explicit physical parameterizations of the contributions of each atmospheric component to the light attenuation. Moreover, the decoupling hypothesis between clearsky irradiance and cloud attenuation is discussed in reference [44]. Therefore, this first step of the MDSSFTD processing is to determine the transmittance of atmospheric gases, the Rayleigh scattering, the aerosol mixture, and clouds. These individual contributions are then used to estimate the total atmospheric transmittance. A Digital Elevation Model (DEM) is used (see Section 3.2.1) to perform an elevation correction. A discrimination with different parameterizations of the total atmospheric transmittance are therefore used for clear and cloudy pixels. Clear and cloudy pixels are determined through the cloud mask which, therefore, represents an important information. The cloud mask used by MDSSFTD is the satellite derived product described in Section 3.3.1. Despite this clear-sky/cloudysky discrimination, a particular attention is paid to ensure the spatial continuity of the radiation estimates at the borders between cloudy-and clear-sky conditions. The method used for the retrieval of DSSF and f d is composed of three main steps (see the blue circles in Figure 1): first, the estimations of the individual atmospheric contributions and the total atmospheric transmittance (Section 4.1); second, the parameterization for clear-sky conditions (Section 4.2); and third, the parameterization for cloudy-sky conditions (Section 4.3). The approach chosen for the retrieval is based on explicit physical parameterizations of the contributions of each atmospheric component to the light attenuation. Moreover, the decoupling hypothesis between clear-sky irradiance and cloud attenuation is discussed in reference [44]. Therefore, this first step of the MDSSFTD processing is to determine the transmittance of atmospheric gases, the Rayleigh scattering, the aerosol mixture, and clouds. These individual contributions are then used to estimate the total atmospheric transmittance. A Digital Elevation Model (DEM) is used (see Section 3.2.1) to perform an elevation correction. A discrimination with different parameterizations of the total atmospheric transmittance are therefore used for clear and cloudy pixels. Clear and cloudy pixels are determined through the cloud mask which, therefore, represents an important information. The cloud mask used by MDSSFTD is the satellite derived product described in Section 3.3.1. Despite this clear-sky/cloudy-sky discrimination, a particular attention is paid to ensure the spatial continuity of the radiation estimates at the borders between cloudy-and clear-sky conditions.

Satellite TOA Radiance
The LSA SAF operational system provides TOA radiances from SEVIRI in the channels centered at 0.6 µm, 0.8 µm, and 1.6 µm (see Figure 2). These data are only used in the cloudy-sky method for the estimation of the cloud albedo and cloud transmittance, as it will be explained in Section 4.1.5. Figure 2 shows the radiation reflected back to space (top-of-atmosphere radiances). The position of clouds for the 15 August 2017 at 12:30 is given in Figure 3. Cloudy areas usually correspond to the areas with highest TOA radiances. It illustrates that the brighter the clouds appear in the satellite images, the more radiation is reflected by them, and therefore the less radiation reaches the ground. The LSA SAF operational system provides TOA radiances from SEVIRI in the channels centered at 0.6 μm, 0.8 μm, and 1.6 μm (see Figure 2). These data are only used in the cloudy-sky method for the estimation of the cloud albedo and cloud transmittance, as it will be explained in Section 4.1.5. Figure 2 shows the radiation reflected back to space (top-of-atmosphere radiances). The position of clouds for the 15 August 2017 at 12:30 is given in Figure 3. Cloudy areas usually correspond to the areas with highest TOA radiances. It illustrates that the brighter the clouds appear in the satellite images, the more radiation is reflected by them, and therefore the less radiation reaches the ground.

TOA Incoming Sun Radiation
The solar constant used in this algorithm, for both clear-sky and cloudy-sky modules, covers the shortwave band between 0.250 μm and 4 μm, and equals 1367 W/m 2 , as given by the World Meteorological Organization [45]. However, this value will probably need to be adjusted in the next years. The spaceborne radiometers are becoming more and more accurate, and recent measurements seams to yield 1362 Wm −2 with an uncertainty of order of 2 Wm −2 [46]. The Sun-Earth distance factor v(t) can be estimated as with the day angle Г (in radian) defined by reference [47] as Г = 2 ( − 1) 365 ⁄ , being the day of the year.

Digital Elevation Models
The MDSSFTD algorithm for clear-sky retrievals relies on two Digital Elevation Model (DEM) files, that are used to correct the aerosol layer height and distribution. The first DEM, noted DEMMSG, is derived from the USGS GTOPO30 DEM by bilinear interpolation to the SEVIRI grid resolution. The

TOA Incoming Sun Radiation
The solar constant used in this algorithm, for both clear-sky and cloudy-sky modules, covers the shortwave band between 0.250 µm and 4 µm, and equals 1367 W/m 2 , as given by the World Meteorological Organization [45]. However, this value will probably need to be adjusted in the next years. The spaceborne radiometers are becoming more and more accurate, and recent measurements seams to yield 1362 Wm −2 with an uncertainty of order of 2 Wm −2 [46]. The Sun-Earth distance factor v(t) can be estimated as with the day angle Γ (in radian) defined by reference [47] as Γ = 2Π(k i − 1)/365, k i being the day of the year.

Digital Elevation Models
The MDSSFTD algorithm for clear-sky retrievals relies on two Digital Elevation Model (DEM) files, that are used to correct the aerosol layer height and distribution. The first DEM, noted DEM MSG , is derived from the USGS GTOPO30 DEM by bilinear interpolation to the SEVIRI grid resolution.
The GTOPO30 DEM exploited is originally available at a spatial resolution of 30 arc-seconds (i.e., at approximately 1 km resolution, https://lta.cr.usgs.gov/GTOPO30). This DEM estimates the height at each SEVIRI pixel. The second DEM, noted DEM CAMS , is the bilinear interpolation to the SEVIRI grid of the coarser DEM used to model the atmospheric fields that are used as input to the MDSSFTD algorithm (see Section 3.3). Figure 3a shows DEM MSG over the MSG disk.

Land Surface Albedo
The surface albedo used as input for the MDSSFTD algorithm corresponds to the MSG Daily surface Albedo product (MDAL) distributed daily by the LSA SAF in near real time [10,48,49]. The total Broad-Band shortwave Bi-Hemispherical albedo (BB-BH dataset) is chosen as the most appropriate variable to estimate the albedo of the surrounding environment. This 'environmental' surface albedo considers the contribution of the area surrounding the processed pixel to the multiple scattering of radiation between the surface and the lower boundary of the atmosphere (indeed, BH offers an integration of directional effects over the target). The MDAL product covers the full MSG disk and is provided at the SEVIRI native spatial resolution. Due to near real time constraints, the albedo of the previous day is used as input of the algorithm. This hypothesis of a stationary albedo is plausible in most of situations. Figure 3b shows the MDAL BB-BH surface albedo for the 14th of August, 2017.

Cloudiness
Two distinct methods are used to estimate the DSSF, one for clear-sky conditions and one for cloudy-sky conditions (see Section 2.2). The cloud mask used in MDSSFTD is therefore a crucial variable as it allows identifying the cloudiness for all pixels with a high confidence. The absence or presence of clouds for each SEVIRI pixel is derived from the cloud mask that is generated and distributed by the SAF NWC (http://www.nwcsaf.org, reference [50]) to support nowcasting applications. This cloud mask is generated using the operational MSG cloud detection method initially developed by reference [51] and made available every 15 min over the full MSG disk. Figure 3c shows an example of this cloud mask for the 15 August 2017 at 12:30 UTC.

H 2 O and O 3 Columnar Content Modeling
The total columnar water vapor (u H 2 O ) and total columnar ozone (u O 3 ) content exploited in this study are obtained from the ECMWF Integrated Forecast System (IFS), made available globally with a 0.125 • × 0.125 • spatial resolution. The operational LSA SAF chain makes use of forecasts up to 36 h ahead available at hourly time-steps. The latest water vapor and ozone forecasts available are then interpolated temporally to a 15-min estimate, and spatially on the MSG disk for ingestion by the MDSSFTD algorithm. Figure 4 gives an overview of u H 2 O and u O 3 fields.

Aerosol Optical Depth and Speciation
The Copernicus Atmosphere Monitoring Service (CAMS, https://atmosphere.copernicus.eu) distributes the aerosol optical depth (AOD) at 550 nm forecasted by the European Centre for Medium-Range Weather Forecasts (ECMWF) for eleven natural and anthropogenic components. Different components are: Sulfate (SU), Sea-Salt (SS), Black Carbon (BC) for 3 classes of particle size, Organic Matter (OM) for 3 classes of particle size, and Dust (DU) for 3 classes of particle size. These aerosols optical depths are available in near-real time and at a spatial resolution of approximately 0.4 • × 0.4 • (~40 km × 40 km). The delay of a few years between the forecast data and the reanalyses of CAMS is not acceptable in the LSA SAF due to the NRT constraints. Also, CAMS forecasts of AOD are used as input to the MDSSFTD algorithm to constrain the aerosol content and speciation. These aerosol data are interpolated spatially to the MSG/SEVIRI grid. The latest forecast available before the processed timeslot is used (but no longer than 3 h before).

Aerosol Optical Depth and Speciation
The Copernicus Atmosphere Monitoring Service (CAMS, https://atmosphere.copernicus.eu) distributes the aerosol optical depth (AOD) at 550 nm forecasted by the European Centre for Medium-Range Weather Forecasts (ECMWF) for eleven natural and anthropogenic components. Different components are: Sulfate (SU), Sea-Salt (SS), Black Carbon (BC) for 3 classes of particle size, Organic Matter (OM) for 3 classes of particle size, and Dust (DU) for 3 classes of particle size. These aerosols optical depths are available in near-real time and at a spatial resolution of approximately 0.4° × 0.4° (~40 km × 40 km). The delay of a few years between the forecast data and the reanalyses of CAMS is not acceptable in the LSA SAF due to the NRT constraints. Also, CAMS forecasts of AOD are used as input to the MDSSFTD algorithm to constrain the aerosol content and speciation. These aerosol data are interpolated spatially to the MSG/SEVIRI grid. The latest forecast available before the processed timeslot is used (but no longer than 3 h before).  In order to compute the radiative quantities related to each type of aerosol, we use the Global Aerosol Data Set (GADS) [52], which makes available the optical properties of a number of aerosol types. However, these types do not completely match the CAMS aerosol classification and, therefore, a correspondence between CAMS and GADS aerosol components was implemented, as proposed by reference [53]. At the end, the aerosol layer considered in MDSSFTD is a dynamic mixture of the In order to compute the radiative quantities related to each type of aerosol, we use the Global Aerosol Data Set (GADS) [52], which makes available the optical properties of a number of aerosol types. However, these types do not completely match the CAMS aerosol classification and, therefore, a correspondence between CAMS and GADS aerosol components was implemented, as proposed by reference [53]. At the end, the aerosol layer considered in MDSSFTD is a dynamic mixture of the following aerosol components: the insoluble particles (modeled by the GADS component "INSO"); the water soluble particles (GADS "WASO" component); the black carbon particles (GADS "SOOT" component); the fine and coarse sea-salt particles, referred to as SSALL (combination of GADS classes SSAM and SSCM); and the fine, medium-sized and coarse mineral particles, referred to as MIALL (combination of GADS classes MINM, MIAM and MICM). The properties of the aerosol components used as input of MDSSFTD algorithm are shown in Table 1.  Table) In order to estimate aerosols individual transmittances, the algorithm uses the SIRAMix method (Surface Incident Radiation estimation using Aerosol Mixtures) developed by reference [53]. SIRAMix is based on a physical parameterization of the direct and diffuse DSSF that is fed by a pre-computed Look Up Table (LUT) of aerosol optical properties. The LUT was generated using the radiative transfer model libRadtran [54]. The LUT file is less than 300 kB and allows fast identification of the aerosol optical properties with an optimized computation cost that is well suited for operational purposes. The LUT provides the following variables for each of the 5 GADS aerosol components (INSO, WASO, SOOT, MIALL, SSALL), indexed i in the following and listed in Table 1: (i) the shortwave spherical albedo of aerosols, A i aer (δ i 0 , u H 2 O ) and (ii) the direct and diffuse individual transmittances, T i aer,dir θ s , The computation of these quantities was made with libRadtran for a predefined range of values for the sun zenith angle (θ s ), varying between 0 • and 85 • ; for the total column content of atmospheric water vapor (u H 2 O ) (with the corresponding modification of the optical properties of hygroscopic aerosols) varying between 0 and 5 g/cm 2 ; and the aerosol optical depths at 550 nm for each aerosol component i (δ i 0 ) with values taken between 0 and 4. Interpolation techniques are used to identify the values of these variables for a given atmosphere/geometry combination. Table 2 summarizes the input data sets that are used for the clear-sky and cloud-sky cases. The only difference is the exclusive use of the SEVIRI TOA radiances for the cloudy-sky case.

Gas Transmittances
The transmittances from gases are used by both methods (clear and cloudy cases) to account for the absorption of gas molecules. They are parametrized empirically according to reference [55] as follows where the sub index gas stands for any of the seven gases mentioned in Table 3; the coefficients a, b, c, and d depend on the gas (see Table 3); the pressure-corrected optical air mass (m ) is described in the following; and u gas corresponds to the total column content for a given gas in units of atm-cm. Here, only the total column water vapor and ozone content are considered as changing variables. These values are derived from the ECMWF forecasts presented in Section 3.3.2, while other u gas values are considered static. Fixed concentrations are given in Table 3. For gases, there is no distinction between direct ('dir') and diffuse ('dif ')-same value of the transmittance is used in both configurations of the radiation.  (11) where p is the atmospheric pressure at the surface altitude (in Pa); p 0 is the mean atmospheric pressure at sea level (p 0 = 101,325 Pa), and H 0 is the altitude above sea level (in kilometers); and where m is the optical air mass at standard pressure p 0 , and is formulated according to reference [59] as: This equation was found accurate for any air mass up to θ s < 85 • with an error of less than 0.5%.

Rayleigh Direct and Diffuse Transmittances
The total (direct plus diffuse) transmittance due to Rayleigh scattering is formulated according to reference [60] as: where m is the pressure-corrected air mass as defined above in Equations (11) and (12). The diffuse transmittance due to Rayleigh scattering is formulated according to reference [61] as: where the 0.5 factor stands for the forward scattering fraction [62], i.e., the Rayleigh scattering is assumed to be isotropic and only half of the scattered radiation goes downward.

Aerosols Transmittances and Spherical Albedos
The characteristics of each individual aerosol component were presented in Section 3.3.3. The combination of several components results in the aerosol mixtures that are observed in reality. To derive the transmittance and the spherical albedo of the whole aerosol mixture, the individual spectral albedos, and the individual transmittances of each of the aerosol components are determined following reference [53]. These quantities have been pre-computed for different combinations of aerosol load, water vapor and solar zenith angle and used in the SIRAMix as Look Up Table (LUT) for each of the GADS components (see Section 3.3.4). The AOD values from each CAMS aerosol component are converted into AOD values corresponding to the GADS aerosol components. This is done based on the following equations [53]: The optical depths δ i 0 of the individual aerosol components i are first height-adjusted before being used for the computation of total aerosol transmittance. Indeed, the aerosol layer spreads vertically from the ground (H 0 ) to the top of the aerosol layer (H TOL ). The CAMS AOD is provided with a ground height (H CAMS ) defined at its own grid resolution (0.4 degrees), which can differ from the real ground height H 0 . Since aerosols are not homogeneously distributed on the vertical, the CAMS AOD converted to GADS types, determined as described in the equation above, may not be adequate to use directly in the computation. Appropriate correction is therefore executed as where i_GADS stands for the GADS components (INSO, WASO, SOOT, SSALL, MIALL) and Z is the scale height in kilometers. More details on this height correction are given in reference [53]. Values of Z and H TOL are provided for each GADS component in Table 4. The equivalent spectral AOD at 550 nm δ EQU550 0 is defined as: where i index stands for the 5 GADS aerosol components (INSO, WASO, SOOT, SSALL and MIALL), δ i,height_cor 0 the equivalent and height-corrected aerosol optical depth for each aerosol component. This equivalent spectral AOD at 550 nm (Equation (17)) is one of the output variables provided in the MDSSFTD product.
Finally, these spectral AOD at 550 nm are converted to total shortwave (broadband) AOD noted ∆ i 0 , as the DSSF is defined based on an integration over a given spectral range. The broadband conversion of the spectral aerosol optical depth of the individual aerosol components is defined as: where α and β are the regression coefficients determined by radiative transfer simulations for varying aerosol optical depths, as detailed in Table 5. The broadband spherical albedos and the transmittances can be obtained by extracting the LUT values that correspond to the broadband AOD values for each of the GADS component.
The total aerosol transmittances and total aerosol spherical albedos for the aerosol mixture are determined as the sum of the individual quantities of each component weighted by the corresponding total shortwave AOD ∆ i 0 , and normalized by the equivalent broadband AOD (∆ O ), as follows [63]: where ∆ o is the equivalent broadband AOD, defined as the summation of the individual broadband aerosol optical depths ∆ i 0 for the 5 components of the aerosol mixture presented above (i.e., INSO, WASO, SOOT, SSALL, and MIALL) where i index stands for the 5 aerosol components, δ i 0 the equivalent and height-corrected aerosol optical depths for each aerosol component, and α i and β i the AOD spectral to broadband conversion coefficients provided in Table 5.

Total Cloud-Free Atmospheric Transmittance
The total cloud-free atmospheric transmittance is equal to the clear-sky transmittance. In cloudy situations, the same value of cloud-free atmospheric transmittance is used for the transmittance happening below the cloud layer. This parameter helps providing an estimation of the total cloudy-case transmittance (Section 4.1.5). Both cloud and cloud-free transmittances and albedos are key parameters for the estimation of the incoming solar radiation at the surface level.
The cloud-free transmittance, T C− f reeA is delivered by using the output of SIRAMix module which calculates the clear-sky total DSSF (E clear tot , described in reference [53] and briefly discussed in Section 4.2.2): This cloud-free atmospheric transmittance T C− f reeA is an estimate of the combined effective transmittances of gas absorption, Rayleigh scattering, and aerosol extinction. The parameter T C− f reeA is combined with the transmittance of clouds in order to estimate the total DSSF for cloudy-sky case.
T C− f reeA includes single and multiple scattering contributions between the surface and the albedo of the cloudless atmosphere. One important variable for the cloudy-sky case, which is indirectly estimated here from Equation (21), is the single scattering clear-sky transmittance (T ss C− f reeA ). The two terms in Equation (21) express the product between the single scattering clear-sky transmittance (T ss C− f reeA ) and a multiple scattering factor (in the right hand side of the equation) between the clear-sky atmosphere albedo (A C− f reeA ) and the surface albedo (A S ). The albedo of the cloud-free atmosphere is estimated as follows where A aer is the total broadband aerosol albedo derived from Equation (19), and A Ray is the Rayleigh albedo approximated by reference [64] as a constant value (equal to 0.0685).

Total Cloudy-Sky Atmospheric Transmittance
For cloudy pixels the DSSF is estimated based on a simplified physical description of the radiation interactions in the atmosphere-surface system according to reference [65] and reference [66]. Here it is assumed that a homogeneous cloud layer covers the whole image pixel. The effective transmittance factor T C−A for cloudy atmosphere is therefore given by where T ss C− f reeA is the cloud-free single scattering atmospheric transmittance, T C is the cloud transmittance, A S is the surface albedo, T bc is the atmospheric transmittance between the surface and the cloud, and A C is the cloud albedo of the surrounding pixel approximated by the cloud albedo. The uncertainty due to this approximation can become important in the case of a broken-cloud sky -especially in presence of a dense cloud layer surrounding a partially cloudy contaminated area. The denominator of Equation (23) quantifies the multiple scattering between the surface and the bottom of the cloud layer.
The cloud transmittance T C and the cloud albedo A C are key parameters for the estimation of the incoming solar radiation at the surface level. They are highly variable with respect to time and space as they depend on the often-rapid evolution of clouds. The instantaneous values of these two cloud variables are determined from the satellite measurements and with the help of a physical model. The estimation is performed by inverting a system of two equations that approximate the radiative transfer happening within the atmosphere-surface medium including aerosols, gases, and clouds. The approach is an improved version of the one proposed by reference [10]. Because these two cloud variables have broadband characteristics, the satellite spectral reflectances in the spectral channels of the satellite instrument (here the MSG/SEVIRI channels at 0.6 µm, 0.8 µm, and 1.6 µm) are first transformed to broad-band top-of-atmosphere albedo A TOA . We apply the spectral conversion relations proposed by reference [67] based on CERES (Clouds and the Earth's Radiant Energy System) TOA broadband data and the angular reflectance model of reference [68].
The first equation of the system that is inverted to derive T C and A C expresses the albedo at the top of the atmosphere A TOA as the sum of the different contributions that are illustrated in Figure 6. The TOA albedo depends on (i) radiation coming from Rayleigh scattering of gases above the cloud layer (A R ), (ii) radiation reflected by the cloud layer that is attenuated by gases above (A c T SunCloudSat ), (iii) radiation reflected by the surface that is attenuated by the atmosphere and the cloud layer (A s T SunSur f aceSat T 2 c ), and (iv) radiation reflected by the aerosol layer that is attenuated by the gases and the cloud layer (A aer T SunCloudSat T 2 c ): where T bc is the transmittance below cloud, T SunSur f aceSat is the transmittance between the sun, the surface, and the satellite, and T SunCloudSat is the transmittance between the sun, the cloud, and the satellite. It is important to note that the third term of the equation above accounts for multiple scattering between the surface and the bottom of the cloud layer. Similarly, the fourth term accounts for multiple scattering between the aerosol layer and the cloud layer above it. The initial parametrizations of bc and (hereafter called ′ bc and ′ ) proposed by references [64,65], and also used by reference [10], are only considered as gas transmittances between the cloud and the surface. We improved this parametrization by also considering the influence of aerosols particles by doing = ′ × The initial parametrizations of T bc and T SunSur f aceSat (hereafter called T bc and T SunSur f aceSat ) proposed by references [64,65], and also used by reference [10], are only considered as gas transmittances between the cloud and the surface. We improved this parametrization by also considering the influence of aerosols particles by doing T SunSur f aceSat = T SunSur f aceSat × T 2 aer and T bc = T SunSur f aceSat /T SunCloudSat = T bc × T 2 aer . In comparison to previous formulation by reference [10], the fourth term in Equation (24) was added for the cloud-aerosol radiative coupling. The new formulation also accounts for the multiple scattering between the top of aerosol layer and the bottom of the cloud layer by using the aerosol albedo. We remember that the MDSSFTD algorithm assumes that aerosols only exist below the cloud layer.
The second equation of the inverted system gives the cloud transmittance T C expressed in terms of the cloud albedo A C and the cloud absorption a C as: The cloud absorption is modeled as a linear function of the cloud albedo by introducing the "cloud absorption factor" α. The currently employed numerical value of α = 0.11 was not derived from first principles, but has been adjusted by matching the final flux estimates with the help of a validation database [69]. This parameter therefore mainly serves for "absorbing" the methodological approximations and uncertainties, rather than for quantifying the physical cloud properties. Note that the cloud transmittance T C could be improved by using a radiative transfer model and taking into account the varying cloud optical thickness as a function of cloud type. Finally, combining the expressions Equation (24) and Equation (25) allows us to calculate the two unknowns T C and A C from the "observable" A TOA by solving a quadratic equation.
The total aerosol transmittance T aer is derived from the total DSSF formula (see Equation (4)) applied to the clear-sky case, as follows where the individual transmittances of gases and Rayleigh scattering are defined in Sections 4.1.1 and 4.1.2. Finally, the total cloudy-sky atmospheric transmittance is obtained by combining Equation (23) and Equation (21): where T C− f reeA is the cloud-free total atmospheric transmittance. The previous formulation assures the continuity at the transition between the cloud and cloud-free situations. In the case of no cloud (when T C = 1 and A C = 0), T C−A becomes equal to T C− f reeA , which corresponds to the cloud-free transmittance of the atmosphere.

Description
In the absence of clouds, the estimation of DSSF is carried out using the method developed by reference [53] referred to as SIRAMix in order to estimate the aerosols individual transmittances. SIRAMix is based on a physical parametrization method coupled with a pre-computed LUT of aerosol properties (see Section 3.3.4) that comprises mainly the direct and diffuse aerosol transmittances of different aerosol components and their corresponding albedos. The input data for aerosol and gaseous components are described in Section 3.3. The originality of SIRAMix is that the effective transmittance of the atmosphere (see Equation (4)) is calculated based upon a combination of the individual transmittances of each gaseous and each aerosol components (i.e., sulfate, soot, sea salt, organic carbon, and mineral dust) with a distinction between direct and diffuse components. The total atmospheric attenuations for direct and diffuse radiation are decomposed into multiple attenuation processes of individual atmospheric direct and diffuse components in order to consider the presence of atmospheric gases (T gases ), Rayleigh scattering (T Ray ), and aerosol extinction (T aer ). This approach is well adapted for the estimation of the diffuse DSSF. The calculation of the total cloud-free transmittance T C− f reeA is detailed in Section 4.1.4.

Total DSSF for Clear-Sky Case
The total DSSF, E clear tot , reaching the surface under clear-sky conditions is derived with the SIRAMix method as explained above (Section 4.2.1). Its estimation is decomposed in two components: one diffuse, E clear di f , and one direct, E clear dir , and technically computed as: The diffuse DSSF is also decomposed into two components again: a single scattering term and a multiple scattering term. The estimation of these terms mainly depends on the diffuse transmittance of aerosols (T aer,di f ; Equation (19)); the diffuse transmittance due to the Rayleigh scattering (T Ray,di f ; Equation (14)); the shortwave 'environmental' albedo of the surrounding surface (see Section 3.2.2); the gas transmittances (T gas ; see Section 4.1.1); and the albedo of the atmosphere under clear-sky condition (see Section 4.1.4). A more detailed description of the calculation of E clear tot and of its two sub-components (E clear di f , E clear dir ) is given in reference [53]. The most important feature to retain is that the approach is based on the individual transmittance contributions of the different atmospheric components as described above (Section 4.1). In addition, satellite radiances are not directly used here to estimate the incoming solar radiation, but only the sun radiation at the TOA (Section 3.1.2). However, the satellite derived surface albedo is used to retrieve the multi-scattered radiation from the atmosphere to the surface, which could become an important contribution over bright surfaces.

Equivalent AOD, Diffuse Fraction and Opacity Index for Clear-Sky Case
The diffuse fraction at the surface level in clear-sky conditions ( f clear D ) is derived as follows: This diffuse fraction at the surface level is one of the datasets provided in the MDSSFTD product. The clearness index (Kt) is approximated as the amount of the total clear-sky DSSF to the amount of solar flux that is observed at the top of the atmosphere The opacity index (OI), which is also provided as output of the MDSSFTD product, is derived from Equation (7).
The calculation of the equivalent AOD at 550 nm is detailed in Section 4.1.3 (Equation (17)).

Description
In cloudy conditions, the total DSSF is strongly anti-correlated with the observed reflectance at the TOA. The method for the total DSSF retrieval under cloudy conditions first determines the total shortwave TOA albedo (cloud albedo) and cloud transmittance from the TOA reflectances observed by SEVIRI (Section 3.1.1). Second, the same representation of gases and aerosol contributions than in the clear-sky case is used. Cloud-related terms are determined by inverting a simplified radiative transfer-based model (i.e., the system of two equations described in Section 4.1.5). Finally, the effective transmittance of the cloudy atmosphere, T C−A (Section 4.1.5), is used in Equation (4).
In a different way than the clear-sky method which calculates the diffuse DSSF as a prior step to estimate the diffuse fraction, the cloudy-sky method determines f d based on the computation of the clearness index (Kt). This index was defined by reference [70] and is equal to the overall transmittance T of the atmosphere (see Equation (4)).

Total DSSF in Cloudy Conditions
The total DSSF, E cloud tot , in a cloudy situation is determined injecting the total effective transmittance T C−A from Equation (27) in the following equation where E 0 is the solar constant; θ s the solar zenith angle, v(t) is the Sun-Earth distance factor (see Section 3.1.2).

Equivalent AOD, Diffuse Fraction and Opacity Index for Cloudy Case
The diffuse fraction at the surface level in cloudy-sky conditions ( f cloud D ) is estimated from the total DSSF using the clearness index (Kt) as defined in Equation (8). Variable f cloud D is estimated as proposed by reference [70]: To ensure realistic values, the diffuse fraction is set to a value of 1 when the cloud transmittance is zero (i.e., in the presence of very bright clouds). When the cloud transmittance is one (e.g., clear-sky situation incorrectly flagged as cloudy by the input cloud mask), then the clear-sky method is used. The case when K t ≥ 0.78 is quite uncommon and only occurs in case of very thin clouds and pure atmosphere with low quantities of atmospheric particles and gas (or false cloud detections). This empirical method was chosen as it was validated over a number of sites in Europe. Reference [71] reviewed a number of direct/diffuse irradiance separation methods. Many of them were complex methods with several predictors which proved good efficiency. However, these methods were developed for high-frequency irradiance data, or they rely on future state estimators so that they are not implementable in near-real time nor with 15-min instantaneous satellite measurements like SEVIRI. The methods developed by reference [26] or reference [72] also seem powerful with a low degree of complexity. The impact of such methods is planned to be further analyzed for the case of our MSG data in the future in the framework of the evolution of the MDSSFTD product.
The estimation of the opacity index (OI) of the total column of the atmosphere remains unchanged compared to the clear-sky case (Equation (7)).
The calculation of the equivalent AOD at 550 nm is strictly identical under clear-sky and cloudy-sky conditions and detailed in Section 4.1.3.

Comparison with the Previous LSA SAF DSSF Product Version
The MDSSFTD method differs from the MDSSF method (from reference [10] and currently used for the operational production of DSSF) in several aspects such as the parametrization of gaseous transmittances and the aerosol contribution. The only case when the two algorithms would give As it can be seen in Figure 7, the total DSSF decreases over the areas that are cloudy, whereas the diffuse fraction increases over these regions. However, diffuse fraction can also increase with the AOD (see in North Africa over the Sahara Desert due to dust and in Southern Africa over Angola and Congo due to biomass-burning aerosols).

Comparison with the Previous LSA SAF DSSF Product Version
The MDSSFTD method differs from the MDSSF method (from reference [10] and currently used for the operational production of DSSF) in several aspects such as the parametrization of gaseous transmittances and the aerosol contribution. The only case when the two algorithms would give similar results would be under clear-sky conditions when MDSSFTD aerosol inputs are equivalent to a horizontal visibility of 20 km corresponding to continental aerosols, which is the constant aerosol conditions set in the MDSSF algorithm. In the cloudy-sky case, the main difference lies in the impact of the aerosol layer below the cloud layer that is now considered in the MDSSFTD product. Figure 8 presents the difference of the total solar flux measurements from the two LSA SAF satellite-derived product versions (ref. LSA-201 and LSA-207 for MDSSF and MDSSFTD, respectively) for a particular day and time slot (15 August 2017 at 12:30 UTC). The spatial differences between the two estimates of the DSSF are many. For example, the area with the maximal difference between the two products is located in Central Africa (see red-coloured values on Figure 8), with values above 100 W/m 2 over Angola/Republic of Congo. This is due to the important amount of aerosols in the atmosphere over this region and for this period of the year (see equivalent AOD in Figure 7d), resulting from biomass burning (see Figure 5). The strong influence of aerosols is considered for the estimation of MDSSFTD total flux, while it is not the case for MDSSF. Moreover, Figure 7c shows that this region is cloudy, which contributes to increase the differences. Consequently, this is an interesting case that illustrates the major importance of considering the attenuation of the radiation due to aerosols not only in clear-sky conditions. A more detailed comparison between the operational MDSSF product and the upcoming MDSSFTD will be given in the companion article [73].  Figure 7d), resulting from biomass burning (see Figure 5). The strong influence of aerosols is considered for the estimation of MDSSFTD total flux, while it is not the case for MDSSF. Moreover, Figure 7c shows that this region is cloudy, which contributes to increase the differences. Consequently, this is an interesting case that illustrates the major importance of considering the attenuation of the radiation due to aerosols not only in clearsky conditions. A more detailed comparison between the operational MDSSF product and the upcoming MDSSFTD will be given in the companion article [73].

Known Issues and Limitations
The MDSSFTD product for clear-sky conditions has been updated compared to the previous operational method (MDSSF). The main evolution is the consideration of the presence of dynamic aerosol conditions and the delivery of diffuse fraction estimates. The first limitation is that the current version of the algorithm depends on the quality of the aerosol inputs. The uncertainty associated to the AOD coming from CAMS is the major source of error for the estimate of DSSF and for the diffuse fraction under clear sky conditions. Reference [53] showed that clear-sky retrievals of DSSF are accurate enough for CAMS reanalysis dataset but of much lower quality when using the forecast products as input to the processing.
The accuracy of the MDSSFTD product also depends on the uncertainty introduced by the dependencies on the surface albedo and on the approximations for the calculation of TOA albedo

Known Issues and Limitations
The MDSSFTD product for clear-sky conditions has been updated compared to the previous operational method (MDSSF). The main evolution is the consideration of the presence of dynamic aerosol conditions and the delivery of diffuse fraction estimates. The first limitation is that the current version of the algorithm depends on the quality of the aerosol inputs. The uncertainty associated to the AOD coming from CAMS is the major source of error for the estimate of DSSF and for the diffuse fraction under clear sky conditions. Reference [53] showed that clear-sky retrievals of DSSF are accurate enough for CAMS reanalysis dataset but of much lower quality when using the forecast products as input to the processing.
The accuracy of the MDSSFTD product also depends on the uncertainty introduced by the dependencies on the surface albedo and on the approximations for the calculation of TOA albedo from the SEVIRI spectral radiances. For example, in the presence of snow, the satellite values of surface albedo are often less accurate. This will affect the accuracy of the estimated DSSF, as surface albedo helps to estimate the multiple scattering of radiation between the surface and the atmosphere. Furthermore, multiple scattering can be very important in the presence of snow, which corresponds to a high value of albedo.
The cloud mask also represents a highly sensitive auxiliary input data. Missing cloud detections can significantly impact the product quality. Also, we know that the presence of broken clouds (i.e., clouds not covering a whole SEVIRI pixel) are not properly considered by the model (nor the cloud mask), which assumes a homogeneous cloud layer over the whole pixel [74]. MDSSFTD results on SEVIRI pixels switching in time from clear to cloudy conditions or reciprocally, from cloudy to clear, are very likely to correspond to conditions of broken clouds, and therefore shall be considered with caution by the users.
A limitation could be also found due to the underestimation of the multi-scattering processes in the case of a recent snow-fall that could therefore not be considered yet by using as input the yesterday's LSA SAF satellite surface albedos.
Finally, a restriction is applied in the look-up-table for aerosol properties to sun zenith angles smaller than 85 degrees. This is required as the radiative transfer model, and the MDSSFTD algorithm itself, are not well suited for extreme zenith angles. In conclusion, this restriction disables the data availability under twilight conditions (θ s > 85 degrees). For the MSG/SEVIRI platform, this occurs in wintertime over high latitude regions.

Access to the Code Sources and Data Policy
The MDSSFTD code is made available in open source under the CeCILL-C license. Users shall therefore follow the principles of this license and, in particular, the rules for the exploitation of the code. In addition, users are kindly requested to duly acknowledge authors of the code: "The MDSSFTD code (Carrer et al., 2019a, b) was made available to the community by CNRM/Météo-France thanks to the support of EUMETSAT." The code is accessible at the following location: https://opensource.umr-cnrm.fr/projects/mdssftd.

Conclusions
This article presents a method for the estimation of the total DSSF and the corresponding diffuse fraction for both clear-sky and cloudy-sky situations observed by the SEVIRI instrument aboard the MSG satellite. The method first determines the individual contributions of clouds, gases and aerosols. Second, two specific retrieval modules based on physical parameterizations are used to estimate the DSSF and the diffuse fraction for clear-sky and cloudy-sky situations. The code corresponding to the scientific algorithm is made available to the community. The algorithm is today implemented in the LSA SAF chain and the associated MDSSFTD product is available to the community since July 2019; back-processing may also be requested. The LSA SAF dissemination methods are described in reference [49]. The validation of the method presented in this article is reported in the companion paper [73]. The evaluation was conducted by comparing the instantaneously retrieved total DSSF and diffuse fraction (fd) against ground measurements and other satellite-based products. For all sky conditions, total DSSF obtained an average MBE of 3.6 W m −2 and rMBE of 0.3% when compared to in situ measurements, whereas MBE and rMBE for the diffuse fraction were found to be, respectively, −0.04% and −17.7%. Furthermore, the MDSSFTD product showed a good agreement when compared to other radiation products such as the one from CAMS and the original MDSSF product from the LSA SAF. Finally, it is important to notice that although the LSA SAF was initially designed to serve the needs of the meteorological community, there is no doubt that this MDSSFTD product can address a much broader set of applications, which includes users from the agricultural and forestry sectors, as well as from the renewable energy industry. Funding: This research was partially funded by EUMETSAT through the LSA SAF project and by the authors respective affiliations.