Sentinel-2 L1C radiometric validation using deep convective clouds observations

ABSTRACT In the frame of the European Space Agency (ESA) Scientific Exploitation of Operational Missions (SEOM) program new algorithms are developed to validate the Sentinel-2 level 1C (L1C) product radiometry, beyond the baseline algorithms used operationally in the frame of the Sentinel-2 Mission Performance Centre (MPC). In this context this paper presents the implementation of a Sentinel-2 radiometric validation approach based on deep convective cloud (DCC) observations. Due to their physical properties DCCs can be used to monitor the radiometric response degradation of the reflective solar bands of optical sensors. Their observation allows interband radiometry validation in the visible-near infrared (VIS-NIR) domain relatively to an a priori well calibrated reference band. We first present the selection of Sentinel-2 data acquired over DCC targets, as well as the tools and assumptions used for the modeling of the theoretical DCC radiometric response. The validation methodology is then thoroughly described and justified. It is based on the comparisons between the observed and the simulated top-of-atmosphere reflectance spectrum. Interband radiometric validation is performed through the statistical analysis of a large collection of individual observations. Results show the very good radiometric performance of Sentinel-2 with interband gains lower than 2%.


Introduction
Within the Copernicus program of the European Union the Sentinel-2 mission is dedicated to the acquisition of high-spatial resolution optical imagery from the visible to the shortwave infrared domains of the electromagnetic spectrum. From early 2017 it provides the continuous monitoring of terrestrial surfaces and coastal waters at a global scale completing a 5-days revisit frequency thanks to two identical Sentinel-2 platforms launched respectively on 23 June 2015 and 7 March 2017. Both platforms carry the same imaging payload, the Multi-Spectral Instrument (MSI), which acquires measurements in 13 channels at spatial resolutions of 10, 20, or 60 m . Table 1 recalls the MSI spectral bands with nominal central wavelength, associated bandwidth, and resolution.
Accurate calibration of the MSI is a prerequisite for the quality of the products delivered to the scientific community. It is firstly operated at level 1 for geometrical and radiometrical corrections  prior to the public release of L1C products. The methodology presented in this paper, which is based on the L1C products, therefore consists more in a validation of the MSI calibration than in an actual calibration.
In the frame of the ESA Scientific Exploitation of Operational Missions (SEOM) program (http://seom. esa.int/), new algorithms are developed to validate the Sentinel-2 L1C product radiometry (http://s2radval. acri.fr), beyond the baseline algorithms used operationally in the frame of the S2 Mission Performance Centre .
As aging differently impacts the spectral bands of an Earth Observation (EO) sensor, the interband calibration shall be regularly monitored. The methodologies usually used for this purpose are based on acquisitions over spectrally stable and bright natural targets. Assumption is made that one reference band can be considered well calibrated.
Acquisitions over two natural target types are usually adequate for the interband monitoring: sunglint over deep oceans and deep convective clouds (DCCs), both targets presenting a stable optical behavior from the visible (VIS) to the near-infrared (NIR).
Sentinel-2 is a Land mission with acquisitions mainly over land, coastal areas and main islands. Consequently, few acquisitions are performed over deep oceans and sunglint acquisitions are not regularly available, thus limiting the possibility to perform interband calibration monitoring from sunglint.
On the other hand, DCCs are more often observable from MSI. Such clouds occur in relationship with deep convection within the intertropical convergence zone and can be observed over water, as well as over land (e.g. Aumann & Ruzmaikin, 2013). DCCs have extremely large vertical extent stretching from the near surface up to the tropopause region (about 16-18 km in the Tropics), their tops are therefore very cold (about 200 K). Their horizontal spread depends on the size and level of maturity of the convective system which is generating them (e.g. Futyan & Del Genio, 2007) so that DCCs appear over horizontal scales from few to hundreds of kilometers.
They are made of water droplets on about their lower half and ice particles on their upper part. They are also very dense clouds and consequently exhibit very-high optical depths (order-of-magnitude 100), which minimizes the contributions of molecular scattering, aerosols, and ground to the TOA signal, at the exception of strong stratospheric aerosol loadings from volcano eruptions.
Their high brightness shows up with DCC TOA reflectance reaching unity in all VIS bands. Using radiative transfer simulations, Chen, Hu, Xu, and Zhang (2013) show that cloud albedo tends to stabilize as DCCs become more optically thick. What is true from simulations has been confirmed from observed targets for instance in Fougnie and Bach (2009) utilizing POLDER monthly-mean DCC reflectance. The latter is also mentioned to exhibit regional variability in Sohn, Choi, and Ryu (2015).
A very useful feature of DCC reflectance is its whiteness (i.e. spectral flatness) from VIS to NIR with shape and amplitude mainly depending on the cloud optical thickness (COT) and, to a lesser extent, on the cloud microphysics. Moreover, DCCs behave mostly like Lambertian reflectors, which avoid the use of a bidirectional reflectance distribution function for geometries close to nadir. However, fine radiative transfer at actual geometries of acquisition is necessary to detect subtle changes in the reflectance for more accuracy, which includes the assessment of bidirectional effects (e.g. Sohn, Ham, & Yang, 2009).
It is agreed that DCCs can then be used as spectral invariant targets to monitor the radiometric response degradation of reflective solar bands of EO sensors in the VIS-NIR domain (Fougnie & Bach, 2009;Wang & Cao, 2016). The accurate detection of such targets from remote sensing in combination with their accurate radiative transfer modeling allows a very reliable assessment of the interband calibration of EO sensors.
Next sections are then firstly dedicated to the extraction of DCC reflectance statistics from Sentinel-2 MSI radiometry, notably through the filtering in the water vapor absorption bands at 945 and 1375 nm (resp. B09 and B10). A second part is dedicated to a description of the assumptions used for radiative transfer modeling of the DCC TOA reflectance, notably through the decomposition of the signal into independent top-of-DCC reflectance and (mostly stratospheric) transmittance.
Finally, interband calibration coefficients are analyzed from S2A MSI data covering the period January-February 2016 along with a full uncertainty budget, providing a reliable insight on the quality of the MSI data calibration made prior to the L1C data release.

Sentinel-2 data handling
MSI Level-1C reflectance data is available from the Sentinels Scientific Data Hub (https://scihub.coperni cus.eu/), we concentrate on data from Sentinel-2A. Products are formatted per tile (ortho-images in UTM/WGS84 projection) of about 100 × 100 km 2 . The MSI swath, made of independent acquisitions from 12 detector modules, is 290 km wide and about three tiles so that the reprojection of the original MSI acquisitions only displays about 6 or 7 of the 12 detectors per single tile.
To deal with the different spatial resolutions of the MSI L1C products a scaling of all the reflectances, as well as of the other necessary parameters is done to 60 m, which is fine enough for our purpose. It provides square images of 1831 × 1831 pixels.
Bands of interest for the methodology presented in this article are all bands from 443 to 865 nm for proper calibration, as well as the 945 and 1375 nm water vapor absorption bands for the detection of DCCs. Indeed, strong water vapor absorption allows to detect very high altitude clouds as only above such clouds tiny water vapor amounts are within the optical path, absorption is consequently much weaker than above other clouds. The impact on TOA reflectance is more effective where water vapor absorption is the strongest, therefore giving more credence in the use of the 1375 nm band where the strongest absorption occurs.
Absorption from water vapor at 1375 nm is therefore the determining criterion for the selection of DCC pixels within Sentinel-2 images: we select tiles for further processing if their amount of pixels with TOA reflectance at 1375 nm above a threshold of 0.4 is higher than 10 %. Eventually, selection of the largest DCCs may be done by processing if their amount of  (m)  B01  443  20  60  B02  490  65  10  B03  560  35  10  B04  665  30  10  B05  705  15  20  B06  740  15  20  B07  783  20  20  B08  842  115  10  B08a  865  20  20  B09  945  20  60  B10  1375  30  60  B11  1610  90  20  B12  2190  180  20 pixels with TOA reflectance at 945 nm (resp. 1375 nm) above a threshold of 0.4 is higher than 70% (resp. 20%). When a tile is selected, further selections of the areas of the image consist in selecting pixels with TOA reflectance at 865 nm higher than 0.7 and TOA reflectance at 1375 nm higher than 0.4. Selection according to geometry to avoid strong off-zenith illuminations and/or strong off-nadir observations is performed though not necessary as deep convection is concomitant with high solar zenith angles (lower than 40°) and MSI acquisitions are very close to nadir (lower than 15°).
After this prior selection, DCC macro-structures (or "clusters") can be isolated. Neighboring DCC pixels belonging to a same cluster are labeled to eventually identify more than one macro-structure within the tile. The area of each DCC cluster is computed by conversion of the number of pixels which individual size corresponds to 60 × 60 m 2 . Clusters which area is larger than a circular area with typical radius 10 km are kept for further processing.
Due to radiometric noise the selection of pixels through the use of thresholds leads to eventual gaps within the finally selected areas, which should be avoided as one expects that the horizontal spread of DCCs rather cause continuous coverage. To avoid such effects a series of 30 pixel dilations (about 1.8 km) of the first selection mask is performed.
At such point, one can appreciate the variability of the TOA reflectance acquired by MSI over DCC targets in Figure 1 at B01 (443 nm, the other bands exhibiting the same dynamics are not displayed). The most striking aspect is the heterogeneity, as well as the dynamics of the reflectance with values sometimes as low as 0.5 and sometimes as high as 1.5 which traduces three-dimensional optical effects of DCCs. One-dimensional radiative transfer simulations, as used in this study, cannot capture such effects so that one must consider a certain amount of pixels over which a statistics of an average reflectance should be made to compare the measured reflectance to a modeled, theoretical, reflectance.
To do so, further subdivision of the selected pixels is made into blocks of 100 × 100 pixels (about 6 × 6 km 2 ) to provide a larger amount of DCC samples and a better account of the low frequency spatial variability (e.g. borders, convective core), resulting in a better retrieval of the COT per block of pixels.
The size of the blocks is defined to capture intra-DCC high frequency heterogeneity (i.e. sample regions with both high and low values relatively to the mean) more prone to be related to three-dimensional optical effects. Square blocks are determined nominally with 100 pixels width and height to the exception of margin blocks. Indeed, the reflectance images scaled to 60 m are 1831 pixels wide and long, which is not a multiple of 100, remaining margin blocks are consequently 31 pixels wide and/or long. Blocks of pixels are also attributed to any detector which footprint lies within the block as shown in Figure 1. Eventually, blocks are further subdivided when two different detectors share a common block. In all cases only pixels belonging to the original DCC cluster are considered within the blocks.
After final selection, statistics of the MSI radiometry and of the other necessary parameters (geometry of acquisition, ozone content, etc.) is performed for each block of pixels. For a given wavelength median MSI TOA reflectance is taken to represent the TOA reflectance which will be used as input of the methodology. In addition, a full variance-covariance matrix of all reflectances within each block of pixels will be exploited in the total uncertainty budget.

Radiative transfer modeling
The TOA reflectance acquired at sensor level over a DCC target can be modeled as the sum of the contributions from the atmosphere and the DCC (a "topof-DCC" (TODCC) reflectance), weighted by the gaseous transmission: Þ is the atmospheric signal (essentially Rayleigh scattering); • ρ DCC θ s ; θ v ; Δφ; λ ð Þis the reflectance at the top of the cloud resulting in the interactions of light within the cloud, including multiple-scattering effects; Þ is the gaseous transmission affected by concentrations in water vapor, NO 2 , O 2 , and O 3 , mostly made of contributions from the stratosphere above the cloud. Therefore only above-DCC transmission is computed. ρ t RTM denotes the TODCC reflectance simulated by the radiative transfer model, similarly a measured TODCC reflectance ρ t MSI describes the measured TOA reflectance corrected for gas absorption. Multiple scattering is addressed in the joint modeling of ρ DCC , ρ atm and their potential coupling in such TODCC reflectance.
The model is neglecting contributions from the ground and aerosols which have very little effect within and above DCCs (which is confirmed a posteriori). Only molecular scattering has a slight effect at the smallest DCC COTs and is taken into account. The computation domain is therefore split into two domains: one purely absorbing domain containing the trace gases mostly affecting from above the DCCs, such as ozone, oxygen, and water vapor, the other (within the cloud) purely diffusive in the visible containing ice particles.
The whiteness of DCC can only be seen on the TODCC reflectance. This feature allows to simplify the computation of ρ t DCC θ s ; θ v ; Δφ; λ ð Þonly at nominal central wavelength λ MSI rather than to compute it over the MSI filters that would be more time consuming. On the other hand, instrument filters must be considered for finest accuracy of the transmittances. The latter are therefore computed per nm between 400 and 1500 nm and then integrated over the MSI spectral response functions through

Transmittance computation
For our needs the standard tropical atmospheric profiles from McClatchey et al. (1972) are suitable for all trace gases except ozone. Indeed, stratospheric ozone has quite a strong absorbing effect in the Chappuis continuum between 500 and 700 nm and the variability of its concentration is influencing the transmission, the worst effect being potentially a bias in the modeled TOA reflectance. As most of the ozone is in the stratosphere McClatchey et al. (1972) profiles have been adapted to cope with stratospheric ozone variability by weighting the whole profile using total-column ozone concentrations (weight is unitary for the nominal profile) from a climatology. Also, as the standard profile is defined with a tropopause altitude of about 16 km, the weighted profiles are adjusted to the actual tropopause altitude in order to be conservative of the ozone content above the tropopause. The rationale for such adjustments is the use of the most appropriate ozone profiles above DCCs which tops are assumed to be at the tropopause level so that stratospheric ozone absorption is most accurately taken into account. Ozone total-column content climatological means and standard deviations are built upon monthly climatologies provided by ESA's Ozone Climate Change Initiative (CCI) project (www.esa-ozone-cci.org). Tropopause pressure daily climatologies are provided by National Centers for Environmental Modeling (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis (https://www.esrl.noaa.gov/psd/ data/gridded/data.ncep.reanalysis.tropopause.html). Tropopause altitude is built upon the McClatchey et al. (1972) temperature/pressure tropical profile.
Nadir-zenith two-way above-DCC transmittances T gas 0; 0; λ ð Þ have been tabulated for a series of ozone concentrations and tropopause altitudes to be exhaustively representative of the conditions in the Tropics (ozone concentrations between 200 and 350 DU every 5 DU, tropopause altitudes between 8 and 19 km every 1 km). In the retrieval, actual values are bilinearly interpolated using the ancillary data.
Nadir-zenith two-way transmittances T gas 0; 0; λ ð Þ are finally transformed into the geometry of acquisition using the total airmass m as: A sensitivity analysis provides evidence that such approximation (rather than a computation for series of angles) is valid up to about 0.02 % at all wavebands for DCC targets in the Tropics, which is sufficient in our study.
Sensitivity of nadir transmittance is shown in Figure 2 by varying either the set of ozone concentrations (left) or the set of tropopause altitude (right) over the MSI spectrum used in this study (up to B10, 1375 nm), nominal values being 250 DU and 16 km. Transmission is mainly modified in the Chappuis ozone absorption continuum if ozone concentration varies (left), which is impacting all bands in the VIS-NIR and especially B2 (490 nm), B3 (560 nm) and B4 (665 nm) as high as 4%.The tropopause altitude (and consequently the modeled DCC top altitude) variability (right) impacts the amount of O 2 and water vapor available above the cloud and consequently only affects transmission in the absorption bands but with very strong effect as high as 50%. Consistently with its functional design, the MSI spectral response functions avoid strong absorption features to the exception of water vapor absorption at 945 and 1375 nm.
Þis simulated for various COT and ice crystal effective diameters (De). Ranges are taken from the results of Sohn et al. (2015). The full dependencies then account for: • COT between 40 and 200 by step of 20; • De between 40 and 80 microns by step of 10 microns; • θ s between 0°and 60°by step of 2.5°; Figure 2. Above-DCC transmission for tropopause at 16 km and total-column ozone content from 200 to 350 DU (left) and for ozone content of 250 DU and tropopause altitude from 8 to 19 km (right). Wavebands between 430 and 1412 nm.
• θ v between 0°and 15°by step of 2.5°; • Δφ between 0°and 180°by step of 2.5°; • all nominal MSI wavelengths between 443 and 1375 nm. Figure 3 shows the variability of the simulated tropospheric reflectance for representative examples of sun and viewing conditions (Sun at 0°and 10°, varying viewing and relative azimuth angles). This shows that bidirectional effects have to be taken into account as they can affect the modeled reflectance up to 1%. These results are shown at higher angular resolution although interpolation of the 2.5°steps simulations to the actual geometry of acquisition is done in the methodology. Such interpolation provides enough precision, thanks to the Lambertian behavior of DCC targets, to the exception of scattering angles larger than about 173°, such exceptions are filtered out in the method. Figure 4 shows the behavior of the modeled reflectance with respect to COT or De change.
The dependency with COT (at fixed De, left) shows a strong monotonous response in all wavebands whereas the dependency with De (at fixed COT, right) is only strongly effective in the SWIR (1375 nm). Through the interband methodology there is no way to guide the retrieval by assessing precisely the ice crystal effective diameter; we only possess a rough range of the values. Interestingly, an increase of De provides a slight increase of the simulated spectral slope. We therefore prefer to use De as a fixed parameter for independent assessments of the gains because its uncertainty is more likely to provide bias in the final results.
The spectral dependency is more precisely shown in Figure 5 for the bands of interest up to 865 nm either from ρ TOA θ s ; θ v ; Δφ; λ MSI ð Þ (left) or ρ t RTM θ s ; θ v ; ð Δφ; λ MSI Þ (right). Simulations are displayed for 6 COT values within {40, 60, 80, 120, 160, 200} using different markers (lowest COT leads to lowest reflectance) and for 5 De values within {40, 50, 60, 70, 80} using different colors (deep blue = 40, deep red = 80). As expected, DCC whiteness is better captured on the ρ t RTM θ s ; θ v ; ð Δφ; λ MSI Þ spectrum, rather than on the TOA spectrum, up to 740 nm whatever the cloud model employed. Moreover, a slight contribution from Rayleigh scattering can be seen at the lowest COT through the inflexion of the ρ t RTM spectrum in the blue bands.

Principle of the method for individual interband calibration gains
The DCC interband calibration methodology relies in two steps assuming a reference band well calibrated. Nominally, the method has been developed using the  reference band λ ref ¼ 665 nm although we later provide evidence that results from any band in the VIS leads to similar results within 1/1000 to the exception of the water vapor absorption bands which must in any case be avoided as reference. The method is applied per detector and per DCC target, that is over a block of 100 × 100 pixels where only pixels belonging to a continuous DCC cluster and to a same detector are preselected. The whole block may not be filled by means of the subdivision into detectors but also by means of the initial coverage of the DCC target. In a first step, DCC optical depth is retrieved in the reference band λ ref . It is done by correcting the measured TOA reflectance ρ MSI λ ref ð Þ from gas absorption using the transmittance T gas θ s ; θ v ; λ ð Þ computed in the conditions of observations. This corrected reflectance (i.e. the measured TODCC reflectance in the reference band) ρ t Þ is then confronted to the monotonous response of the modeled reflectance ρ t RTM λ ref ð Þ with respect to COT. This is done by fitting a polynomial relationship P λ ref ρ!COT between log(COT) and ρ t RTM λ ref ð Þ. The choice of a thirdorder polynomial in log(COT) is resulting from a sensitivity analysis. The target COT is the result of the fit P λ ref ρ!COT applied to ρ t MSI λ ref ð Þ. In a second step, the retrieved COT is propagated through the LUTs to retrieve the modeled TODCC reflectance in the other bands ρ t RTM λ ð Þ. This is done per band by applying a similar thirdorder polynomial relationship P λ COT!ρ between log (COT) and ρ t RTM λ ð Þ. Then, the modeled TODCC reflectance is converted back to modeled TOA reflectance ρ RTM λ ð Þ by application of the transmittance ρ RTM λ ð Þ ¼ ρ t RTM θ s ; θ v ; Δφ; λ ð Þ :T gas θ s ; θ v ; λ ð Þ Finally, ratios between modeled and measured TOA reflectances provide individual interband calibration gains per DCC target. Writing g λ ð Þ i;j , where i and j stand for the target and the detector, we have: By definition, the methodology equals the ratio in the reference band so that g λ ref ð Þ i;j ¼ 1. Example of a DCC spectrum from an MSI target and the resulting modeled spectra is displayed in Figure 6. Error bars in the ρ MSI λ ð Þ spectrum are built upon the standard deviation within the block of pixels.

Uncertainty budget at individual gain level
At individual target level, gain uncertainties come from: • uncertainties in ρ MSI λ ð Þ i;j , inputs of the algorithm, by means of statistics of selected pixels; • raw L1C radiometric uncertainties at pixel level, these ones can be assessed through the Sentinel-2 radiometric uncertainty tool (S2-RUT, Gorroño et al., 2017) version 1.1; • uncertainties from the radiative transfer model arising from the choices made in the simulations (microphysical and macrophysical models, uncertainties in the ancillary inputs, choice of the solver. . .).
Uncertainties in ρ MSI λ ð Þ i;j arise from the methodology as ρ MSI λ ð Þ i;j is computed through a statistics over a block of pixels. We use the associated variance-covariance matrix cov ρ MSI λ ð Þ; ρ MSI λ ref ð Þ À Á computed over the collected pixels to propagate these uncertainties through the algorithm.
Raw L1C radiometric uncertainties are evaluated using the S2-RUT on two very different images containing DCCs (one with small DCC coverage, one with high coverage). S2-RUT uncertainties are stable and equal over both images where DCC pixels are selected, which is sufficient to define one set of absolute radiometric uncertainties for all processed images. Values obtained from the S2-RUT are summarized in Table 2. Additional L1C radiometric uncertainty given by the S2-RUT is calculated by the product between MSI reflectance ρ MSI λ ð Þ and the percent errors α RUT λ ð Þ given in Table 2. To be applied on the scale of a block of pixels one must consider that independent radiometric uncertainties are given per pixels through α RUT λ ð Þ Â ρ p MSI λ ð Þ where p stands for a given pixel. It follows by application of the standard error that the block of pixels variance from S2-RUT is: is therefore given by the square of the harmonic mean of ρ p MSI λ ð Þ weighted by α RUT λ ð Þ 2 N p . For a large amount of pixels (between 10 2 and 10 4 for a block of size 100) this additional radiometric uncertainty becomes very small, applying var RUT ρ MSI λ ð Þ À Á would however rewrite: i;j is here denoting the empirical variance of the L1C reflectance within a block of pixels. There is no extra term in the covariance cov ρ MSI λ ð Þ; ρ MSI λ ref ð Þ À Á as no covariance is produced by the S2-RUT.
Finally, uncertainties from the radiative transfer model are considered through the computation of final interband gains independently for different choices of microphysical and macrophysical properties. Uncertainties of ancillary data are propagated in the algorithm.
Considering the formulation of an individual gain at waveband λ and removing the i,j notation the gain variance can be approximated using a second order development. It follows: By means of a first-order approximation of the polynomial relationships P λ ref ρ!COT and P λ COT!ρ an analytical formulation of ρ RTM λ ð Þ can be found as a function of ρ MSI λ ref ð Þ.  Let ρ MSI;0 λ ref ð Þ, ρ t MSI;0 λ ref ð Þ, COT 0 , and ρ t RTM;0 λ ð Þ be respectively the MSI inputs (TOA or TODCC), the retrieved COT and the retrieved modeled TODCC reflectance computed through the retrieval process.
A first-order development of the polynomials allows expressing them as linear functions around the retrieved values provided radiometric uncertainty is small (which is justified a posteriori): where dP λ COT!ρ is the first derivative over log COT ð Þ. log COT 0 ð Þ and ρ t RTM;0 λ ð Þ, as retrieved values, are considered constant with a view to approximately express ρ t RTM λ ð Þ as a linear function of log COT ð Þ around the retrieved solution.
In the other hand, the first derivative dP λ ref Using and B λ ð Þ is a useless constant for each wavelength. It follows that: Transmittance terms are put outside of the covariance term as there is no covariance between the ancillary data nor the geometry of acquisition and the radiometry. However, in the variance term, the contribution of ancillary data and geometrical uncertainty on transmittance uncertainty appears in the form of independent terms. Finally, we distinguish three contributors: • radiometric uncertainty arising from the pixel selection and propagated through the retrieval. It can be fully expressed by means of the variance-covariance matrix of the MSI radiometry over a block of pixels: • ancillary data/geometrical uncertainty propagated through the transmittances • coupling between radiometric and ancillary data/geometrical uncertainties The covariance term highlights the spectral correlation between bands. By means of the spectral flatness of DCC targets this contribution is very strong and cancels out most of the contributions from independent variances as seen in Figure 7.
One has var g λ ð Þ ð Þ ¼ var rad g λ ð Þ ð Þþvar trans g λ ð Þ ð Þþ var coupling g λ ð Þ ð Þ. By means of the calibration with respect to one reference band, uncertainty in the gain in the reference band is null: var g λ ref ð Þ ð Þ¼0. Sensivity analysis (which result is displayed in Figure 8) provides evidence that the radiometric uncertainty within a block of pixels is the main factor contributing to the total uncertainty, therefore it is not necessary to further describe nor consider other contributions in the following.

Computation of global interband gains and associated total uncertainties
Weighted mean of selected and filtered individual gains over a population Ω provides an integrated view of interband calibration gains that can be analyzed per detector, region of interest, time period etc. It can be formulated for any subset of individual gains by where i stands for any individual target and/or detector, g i λ ð Þ being the associated gain at λ and w i λ ð Þ ¼ being the associated weight, here taken as the inverse of the total uncertainty of the individual gain derived from the preceeding section. Such formulation assumes that bias (i.e. offset) is zero, which supposes that proper system calibration (e.g. dark offset calibration) is already performed. The number N Ω is the total amount of individual blocks of pixels used in the statistics. For a DCC covering one full MSI image (i.e. 1831 × 1831 60 m-resolution pixels) N Ω roughly reaches a maximum 100 blocks per detector. The method can therefore be employed on a single image or, preferably to eventually catch geophysical variability, more images.
Since the reference band is supposed to be wellcalibrated we have g λ ref ð Þ Ω =1, interband gains g λ ð Þ obtained through this methodology are related to absolute calibration gains G λ ð Þ by means of the absolute calibration coefficient of the reference band G λ ref ð Þ by:  Consistently with being a radiometric validation methodology, we also have G λ ref ð Þ ¼ 1, meaning that G λ ð Þ ¼ g λ ð Þ and that even the quality of the absolute calibration is assessed although the methodology is relative at first.
Total uncertainties on hg λ ð Þi Ω can be expressed either by the standard deviation of the individual gains (i.e. the dispersion of the individual gains) or by the standard error (i.e. the precision of the mean value).
We choose the standard error to represent the precision on the mean gain values as we want to assess the reliability of the computed mean after application of a high amount of data. The error is then expected to decrease with increasing statistics. Total uncertainties are computed using individual gains uncertainties through An order-of-magnitude can be assessed using the typical uncertainty var From Figure 8 σ g λ ð Þ % 10 À2 so that σ hg λ ð Þi Ω easily gets below 10 À3 per detector as soon as few DCC images are collected. This means that the total uncertainties become very small, if not negligible.
However, these total uncertainties do not contain uncertainty from the modeling choices, such as the macrophysical description of the DCC composition and the ice crystals effective diameter De. Results provided by simulations from a variety of choices should then provide an insight on the modeling uncertainty in the obtained interband gains, resulting in eventual biases rather than in dispersion.

Results and validation
Individual gains are extracted over the Indonesian region delimited by longitudes 90°E-165°E and latitudes 15°S-15°N over the period January-February 2016. This region provides many DCC observations both over water and over land thanks to the sparse display of terrestrial surfaces. Independent computations are performed for nine DCC models: three macrophysical models times three microphysical models.
Three macrophysical models are defined to assess the sensitivity to the vertical proportions of ice and water within the cloud, both layers being assumed independent and homogeneously dividing optical thickness: • "1/2": one half of the cloud is made of water at the bottom, the other half of ice on top; • "1/3": one third of the cloud is made of water at the bottom, the rest of ice on top; • "2/3": two thirds of the cloud is made of water at the bottom, the rest of ice on top.
Three microphysical models are defined to assess the sensitivity to ice crystal effective diameter: 40, 60, or 80 microns. Gains are displayed along with their associated uncertainties (which are negligible for high statistics as explained earlier) in Figure 9 for all detectors combined up to 865 nm.
As these results use statistical weighting from all the detectors there occurs an oversampling of the middle detectors due to the selection of images with high percentage of high B10 reflectance. Indeed, images from side detectors are reprojected over incomplete tiles, normalization of the percentage with respect to the size of the image (instead of the actual number of valid pixels) under samples these incomplete tiles. This normalization will be corrected to allow equivalent DCC statistics for all detectors. The amount of blocks of pixels used in the computation of the gains at 443 nm is displayed in Table 3 relatively to the detector.
From Figures 10-12 results are computed for each detector separately at all wavebands between 443 and 1375 nm (except 665 nm which is the reference). At wavebands up to 865 nm, these results show that all L1C reflectances are within 2%, this is less than the 3% goal specification and better than the results presented in Gascon et al. (2017) from another DCC methodology. The variability of the gains compare well in both methods for the VIS bands, less for the NIR bands. Per detector interband gain variability is especially important at B01, which is also in line with the findings in Gascon et al. (2017) although the latter shows a clear drop of the gain between detectors 9 and 10 whereas it is captured in Figure 10 between detectors 8 and 9.
The influence of the macrophysical model is clearer in the VIS bands and nearly non effective in the NIR bands, the impact is smaller than the one of varying the crystals size. Indeed, the influence of De in the spectral slope provides notable changes in the gains slope from the VIS to the NIR.

Geographical stability
The stability of the results is assessed over Central Africa for the same time range. This region is also a very large area over which strong convection takes place mainly over land. Testing over such site is twofold: verify the hypothesis that the influence of the ground is negligible and assess eventual impact of different cloud microphysics (ice crystals size for instance) arising from different generation mechanisms. Shown in Figure 13 results are remarkably comparable to the ones over Indonesia (all detectors combined, as well as per detector, shown only at 443 nm but consistent at other bands).
A summary of the consistency of these results between regions can be given by comparing maximum differences between results over the two regions of interest and maximum differences between results using different macro-and micro-physical models for each region separately.
Differences are all very stable per detector so that simple means over the detectors can be used to display discrepancies per waveband. Tables 4 and 5 summarize the values. Inter-regional differences are nearly one order of magnitude smaller than intra-regional differences according to different models (microphysics and macrophysics) for all bands except 945 nm where differences are however still smaller.
We conclude that, given a sufficient timespan to collect DCC observations over each detector, results are stable with respect to the region of interest that gives credence in the applicability of our methodology over the entire tropical band as soon as DCCs are observed.

Stability with respect to the reference band
Following the above conclusions, interband gains can be computed over the period January-February 2016 combining data over Indonesia and over Central Africa. Using De = 60 microns and the "1/2" macrophysical cloud models interband gains are computed for each waveband used as a reference between B1 and B10. Combining all detectors results are displayed in Figure 14. Because gains are obtained relatively to one reference band each time, the result is scaled to the mean gain obtained between 443 and 865 nm by translation using the mean gain.
The results are remarkably stable (except at B09 and B10, which was expected and must be avoided anyway) even using B07 that is because we have selected a specific De. Note that the vertical scale is reduced (0.98-1.02) compared to results previously shown. This is in line with Gascon et al. (2017) who only performed the analysis on B04 and B07. Results per detector provide the same conclusions.

Applicability for monitoring purposes
The daily occurrence of deep convection allows observations of DCCs throughout the year. Consequently, to assess operational applicability of the methodology is to assess the convergence and/or the stability of the gains with respect to "final" values obtained from a large statistics. As shown earlier, this must be done through an analysis per detector.
Using the datasets from Africa and Indonesia Figure 15 firstly shows, for detector 7 where statistics is the strongest, the relative difference (in percent) between the average gains computed at available dates (expressed in days from first date) and the final gains averaged over the whole period. Results are provided for all bands but only from one DCC  model (De = 60 microns, half water-half ice). As explained earlier, middle detectors provide more DCC statistics by means of the selection procedure. This procedure can be corrected to allow equivalent DCC statistics for all detectors. Therefore, results from other detectors are expected to follow the one from detector 7.
Despite few outliers, independent gains are within ±0.5% of the final gains, which already gives an order of magnitude of the sensitivity of the methodology with respect to a variety of geophysical conditions. As for the time coverage, daily monitoring may not be achieved. On one hand only Africa and Indonesia are covered by these statistics, but on the other hand, as MSI is a land mission, one should expect about one half more of this statistics from extra coverage over South America. We expect then that interband monitoring is possible on a five-day basis with a precision of about 0.5%. Figure 16 then shows, for the same detector 7, the absolute relative difference (in percent) between the average gains computed from all dates up to the given date (again expressed relatively to the first date) and the final gains averaged over the whole period. Results are shown up to 60 days but the period includes few other dates, last points in the graph consequently do not show 0% relative difference.
These results provide the precision expected when gains are computed from a period where radiometry is stable. Corroborating random geophysical conditions, gains are randomly distributed over time and convergence may not be especially linear from first dates. But the precision gets better very fast: lower than 0.2% over a monthly period, lower than 0.05 % over a 1.5 month period for all bands, which is way below the uncertainties due to the DCC model employed.
Sudden changes in interband gains can be monitored through comparisons between individual gains and sliding monthly or bi-monthly means. By removing outliers for computation of the mean over a considered period, operational procedures may be done as soon as new DCC targets are ingested in the system, that is about every five days.

Applicability to other sensors
By similarity of the instruments between Sentinels 2A and 2B this methodology is fully transferable to Sentinel-2B and should be tested against data from this platform in the very near future. As well, only few adaptations are necessary to transfer the method to other sensors acquiring in the VIS-NIR range (e.g. Sentinel-3 OLCI) as soon as an efficient detection of DCCs is possible (for instance through the use of water vapor absorption channels).

Conclusion
We have implemented an interband calibration validation methodology based on the observation of DDCs by Sentinel-2. This methodology is based on the statistics of the Sentinel-2 L1C TOA reflectance over DCC targets used as inputs for radiative transfer simulations. The derivation of COT at a reference band provides its propagation in the simulations to derive the corresponding simulated reflectance spectrum in the VIS-NIR domain. Per MSI tile and for a given detector many individual interband radiometric gains are obtained by confronting the observed and simulated spectra. A synthetical analysis from many observed DCC targets over only a few weeks of MSI observations validates the achieved calibration of Sentinel-2A MSI at level 1C in the VIS-NIR with interband gains lower than 2%. Results are built upon a variety of cloud models which induces most of the variability in the obtained gains, all results agree with the meeting of the specifications. Such methodology is applicable operationally to  spot sudden changes from any detector with a precision of about 0.5% for a short-term assessment (about a week of data) to less than 0.05% for a longer-term assessment (about a month of data) providing a relative radiometric stability within the period considered. This methodology is fully transferable to Sentinel-2B and should be tested against data from this platform in the very near future. As well, only few adaptations are necessary to transfer the method to other sensors acquiring in the VIS-NIR range (e.g. Sentinel-3 OLCI) as soon as an efficient detection of DCCs is possible (for instance through the use of water vapor absorption channels).    Figure 14. Gains and uncertainties from all data over Central Africa and Indonesia for the period January-February 2016. All detectors combined. Results presented as a function of the reference wavelength used in the retrieval.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
This Uwork was supported by the European Space Agency under Grant [4000116454/16/I-Sbo].