Cloud optical properties retrieval and associated uncertainties using multi-angular and multi-spectral measurements of the airborne radiometer OSIRIS

In remote sensing applications, clouds are generally characterized by two properties: cloud optical thickness (COT) and effective radius of water/ice particles (Reff). Most of the current operational passive remote sensing algorithms use a mono-angular bi-spectral method to retrieve COT and Reff. They are based on pre-computed lookup 15 tables while assuming a homogeneous plane-parallel cloud layer and without considering measurement errors and the choice of ancillary data. We use the formalism of the optimal estimation method applied to near-infrared multi-angular measurements, to retrieve COT and Reff, and the corresponding uncertainties related to the measurement errors, the ancillary data, and the cloud model assumption. The used measurements were acquired by the airborne radiometer OSIRIS (Observing System Including PolaRization in the Solar Infrared Spectrum), developed by the Laboratoire 20 d'Optique Atmosphérique. It provides multi-angular measurements at tens of meters resolution, very suitable for refining our knowledge of cloud properties and their high spatial variability. OSIRIS is based on the POLDER concept as a prototype of the future 3MI space instrument planned to be launched on the EUMETSAT-ESA MetOp-SG platform in 2023. The used approach allows the exploitation of all the angular information available for each pixel to overcome the radiance angular effects. More consistent cloud properties with lower uncertainty compared to operational mono25 directional retrieval methods (MODIS-like methods) are then obtained. The framework of the optimal estimation method provides also the possibility to estimate uncertainties of different sources. Three types of errors were evaluated: (1) Errors related to measurement uncertainties, which reach 10% for high values of COT and Reff, (2) errors related to an incorrect estimation of the ancillary data that remain below 0.5%, (3) errors related to the simplified cloud physical model assuming homogeneous plane-parallel cloud and the independent pixel approximation. We show that not 30 considering the in-cloud heterogeneous vertical profiles and the 3D radiative transfer effects lead to uncertainties on COT and Reff exceeding 10%.

very important for constraining climate and meteorological models, improving the accuracy of climate forecast, and monitoring the cloud cover evolution. The instruments onboard Earth observation satellites allow continuous monitoring of the clouds and aerosols and retrieval of their properties from a regional to a global scale.
The cloud properties are retrieved using the information carried by the measurements of the reflected or transmitted 40 radiations by the clouds. Two main optical cloud properties are generally retrieved: the cloud optical thickness (COT) and the effective radius of the water/ice particles forming the cloud (Reff). These optical properties, along with the cloud altitude, allow to characterize the clouds at a global scale and help to determine the radiative impacts of clouds along with their cooling and warming effects (Platnick et al., 2003;Twomey, 1991;Lohmann and Feichter, 2005;Rivoire et al., 2020;Yang et al., 2010). Depending on the available information, various passive remote sensing methods are 45 operationally used in the retrievals of these optical properties. For instance, the infrared split-window technique (Giraud et al., 1997;Inoue, 1985;Parol et al., 1991) uses infrared measurements and is more suitable for optically thin ice clouds (Garnier et al., 2012). The bispectral method (Nakajima and King, 1990) which uses visible and shortwave infrared wavelength, is more suitable for optically thicker clouds. It is currently used in a lot of operational algorithms and in particular with the MODIS radiometer (Platnick et al., 2003). It is also possible to use a combination of multi-50 angular total and polarized measurements in the visible range, such as POLDER measurements (Deschamps et al., 1994), to retrieve COT and Reff (Bréon and Goloub, 1998;Buriez et al., 1997).
The above-mentioned methods are subject to several sources of error. Therefore, a moderate perturbation in the retrieved COT and Reff can, for example, cause variations of around 1 to 2 W/m 2 in the estimation of cloud radiative forcing (Oreopoulos and Platnick, 2008). The quantification of the retrieval uncertainties of these optical properties is 55 therefore critical. The sources of errors originating from the measurements can be quite well evaluated along the instrument calibration process and are often considered when developing a new algorithm (Sourdeval et al., 2015;Cooper et al., 2003) but the errors related to the choice of the cloud model to retrieve the parameters and the assumption made for the radiative transfer simulations should not be overlooked. Currently, computational constraints and lack of information in the measurements force the operational algorithms of cloud products (MODIS, POLDER, and others) 60 to retrieve the cloud optical properties with a simplified 1D-cloud model. In this model, clouds are considered flat between two spatially homogeneous planes in what is known as the plane-parallel and homogeneous (PPH) assumption (Cahalan et al., 1994). Another commonly used assumption is related to the infinite dimension of the PPH cloud and treats each pixel independently without considering the interactions that occur between neighboring homogeneous pixels, known as the independent pixel approximation (IPA, (Cahalan et al., 1994;Marshak et al., 1995b). The effect From medium to large-scale observations greater than 1 km (e.g. MODIS: 1×1 km 2 , POLDER: 6×7 km 2 ), the PPH approximation poorly represents the cloud variability. The subpixel horizontal heterogeneity and the nonlinear nature of the COT-radiance relationship create the PPH bias that leads to the underestimation of the retrieved COT (Cahalan 75 et al., 1994;Szczap et al., 2000;Cornet et al., 2018). The PPH bias increases with pixel size due to increase inhomogeneity. Using the bi-spectral method, the COT subpixel heterogeneity induces also an overestimation bias on the retrieved Reff (Zhang et al., 2012), while the heterogeneity effects appear limited with polarimetric observations (Alexandrov et al., 2012;Cornet et al., 2018). On the contrary, the microphysical subpixel heterogeneity leads to an underestimation of retrieved Reff (Marshak et al., 2006b). 80 At smaller scales, errors due to IPA become more dominant. At this scale, pixels can no longer be considered infinite and independent from their adjacent pixels. The radiations pass from one column to the others depending on the COT gradient. This leads to a decrease in the radiance of pixels with large optical thickness and an increase in the radiance of pixels with small optical thickness, which tends to smooth the radiative field and thus the field of retrieved COT (Marshak et al., 1995a, b). As a result, it can lead to a large underestimation of the retrieved optical thickness (Cornet 85 and Davies, 2008). Adding to these effects, for off-nadir observations, the tilted line of sight crosses different atmospheric columns with variable extinctions and optical properties which tend to additionally smooth the radiative field (Várnai and Davies, 1999;Kato and Marshak, 2009;Benner and Evans, 2001;Várnai and Marshak, 2003). In the case of fractional cloud fields not examined under nadir observations, the edges of the clouds cause an increase in radiances for high viewing angles, which by its turn, increases the value of the retrieved COT (Várnai and Marshak,90 2007), while overestimating the retrieved Reff (Platnick et al., 2003). They are often filtered out of cloud property retrievals especially under low sun angles (Takahashi et al., 2017;Zhang et al., 2019). The illumination and shadowing effects, on the contrary, lead to the roughening of the radiative field, their influence in over and under-estimating the cloud droplet size retrievals are documented in several papers (Zhang et al., 2012;Marshak et al., 2006a;Cornet et al., 2005). 95 The assumption of a vertical homogeneous profile inside the cloud is also questionable. The vertical distribution of the cloud droplets is important to provide an accurate description of the radiative transfer in the cloud (Chang, 2002) and obtain a more accurate description of the cloud microphysics such as the water content or the droplet number concentration. For simplicity reasons, classical algorithms assume a vertically homogeneous cloud model. However, several studies show a dependency of the effective radius retrievals on the SWIR band used. These differences are 100 explained by the non-uniformity of cloud vertical profiles and by the different sensitivities of spectral channels due to the absorption difference (Platnick, 2000;Zhang et al., 2012). Indeed, the absorption by water droplets being stronger at 3.7 µm, the radiation penetrates less deeply in the cloud than at 2.2 and 1.6 µm. Channel 3.7 is therefore expected to retrieve an effective radius that corresponds to a level in the cloud higher than that of channels 2.2 and 1.6. Considerable vertical variation along the cloud profiles is confirmed by many in-situ studies of droplet size profiles and water content 105 as summarized in (Miles et al., 2000). This vertical variation in liquid particle size is an important cloud parameter related to the processes of condensation, collision-coalescence, and the appearance of precipitation (Wood, 2005).
In the operational algorithms, the retrieval of COT and Reff is achieved through pre-computed Look-Up Tables (LUT). This method can be used to process large databases automatically. Its disadvantage is that a modification of the particle https://doi.org/10.5194/amt-2021-414 Preprint. Discussion started: 20 January 2022 c Author(s) 2022. CC BY 4.0 License. model or any other model parameter requires re-generating all these pre-calculated tables. But most importantly it lacks 110 also the ability to assess the uncertainties on the retrieved properties.
In this paper, we present a method to overcome these drawbacks and apply it to the measurements of the airborne radiometer OSIRIS (Observing System Including PolaRization in the Solar Infrared Spectrum, (Auriol et al., 2008) which was developed in the Laboratoire d'Optique Atmosphérique (Lille, France). OSIRIS is the airborne simulator of the 3MI (Multi-viewing Multi-channel Multi-polarization Imaging) radiometer, planned to be launched on MetOp-115 SG in 2023. It can measure the degree of linear polarization from 440 to 2200 nm and has been used onboard the French research aircraft, Falcon-20, during several airborne campaigns: CHARMEX/ADRIMED (Mallet et al., 2016), CALIOSIRIS and AEROCLO-sA (Formenti et al., 2019).
We couple the multi-angular multi-spectral measurements of OSIRIS with a statistical inversion method to obtain a flexible retrieval process of COT and Reff. The exploitation of the additional information on the cloud provided by 120 these versatile measurements implies the use of a more sophisticated inversion method compared to the LUT. The used algorithm is based on the optimal estimation method (Rodgers, 1976(Rodgers, , 2000 which has been widely used for applications in cloud remote sensing (Cooper et al., 2003;Poulsen et al., 2012;Sourdeval et al., 2013;Wang et al., 2016). In this method, the bayesian conditional probability together with a variational iteration method allows the convergence to the physical model that best fits the measurements. Therefore, it introduces the probability distribution 125 of solutions where the retrieved parameter being the most probable, with an ability to extract uncertainties on the retrieved parameters.
This article is organized as follows. Section 2 describes the basic characteristics of OSIRIS and some essential details of the campaign CALIOSIRIS-2. In section 3, a detailed description of the retrieval methodology is presented, including the mathematical framework needed to compute the uncertainties on the retrieved cloud properties. In section 130 4, a case study of a liquid cloud is presented and analyzed. We assessed the magnitude of different types of errors, such as the errors due to measurement noise, the errors linked to the fixed parameters in the simulations, and the errors related to the unrealistic homogeneous cloud assumption. The multi-angular retrievals and uncertainties are compared with the results obtained by the classical mono-angular retrieval algorithms (MODIS-like method) in section 5. Finally, section 6 gives a summary and some concluding remarks.

Instrumentation and airborne campaign
We use the new imaging radiometer OSIRIS. We will go through the main characteristics of the instrument and the airborne campaign CALIOSIRIS. More details about OSIRIS can be found in (Auriol et al., 2008).

OSIRIS
OSIRIS (Observing System Including PolaRisation in the Solar Infrared Spectrum, (Auriol et al., 2008)) is an extended 140 version of the POLDER radiometer (Deschamps et al., 1994) with multi-spectral and polarization capabilities extended to the near and short-wave infrared. This airborne instrument is a prototype of the future spacecraft 3MI (Marbach et https://doi.org/10.5194/amt-2021-414 Preprint. Discussion started: 20 January 2022 c Author(s) 2022. CC BY 4.0 License. al., 2015) planned to be launched on MetOp-SG in 2023. It consists of two optical sensors, each one with a twodimensional array of detectors; one for the visible and near-infrared wavelengths (from 440 to 940 nm) named VIS-NIR (Visible-Near Infrared) and the other one for the near and shortwave infrared wavelengths (from 940 to 2200 nm) 145 named SWIR (Shortwave Infrared). The VIS-NIR detector contains 1392×1040 pixels with a pixel size of 6.45×6.45 µm 2 while the SWIR contains 320×256 pixels with a pixel size of 30×30 µm 2 . Adding those characteristics to the wide field of view of both heads, at a typical aircraft height of 10 km, the spatial resolution at the ground is 18 m and 58 m respectively for the VIS-NIR and SWIR. This leads to a swath of about 25×19 km for the visible and 19×15 km for the SWIR. 150 OSIRIS has eight spectral bands in the VIS-NIR and six in the SWIR. Similar to the concept of POLDER, OSIRIS contains a motorized wheel rotating the filters in front of the detectors. The step by step motor allows only one filter to intercept the incoming radiation at a particular wavelength. The polarization measurements are conducted using a second rotating wheel of polarizers. Given the sensor exposure and transfer times, the duration of a full lap is about 7 seconds for the VIS-NIR and 4 seconds for the SWIR. Figure 1 shows the spectral response of each channel of OSIRIS. 155 The two channels (1240 and 2200 nm) used in this study are red-colored in the figure. OSIRIS is an imaging radiometer with a wide field of view. It has a sensor matrix that allows the acquisition of images 160 with different viewing angles. The same scene can thus be observed several times during successive acquisitions with variable geometries. The largest dimension of the sensor matrix is oriented along-track of the aircraft to increase the number of viewing angles for the same target. For example, when the airplane is flying at 10 km altitude with a speed of 200 to 250 m/s, the same target at the ground can be seen under 20 different angles for the VIS-NIR and 19 for the SWIR. Milieux, Observations Spatiales, Paris) and SAFIRE, the French Facility for Airborne Research. One objective of this campaign was the development of new cloud and aerosol properties retrieval algorithms in anticipation of the future 170 space mission of 3MI intending to improve our knowledge of clouds, aerosols, and cloud-aerosol interactions.
The data used in this work focuses on a cloudy case over ocean surface. A marine monolayer cloud that was observed on 24 October 2014 at 11:02 (local time). The aircraft flew at an altitude of 11 km above the Atlantic Ocean facing the French west coast (46.70°,-2.82°, red arrow in Figure 2a). The solar zenith angle was equal to 59°. In Figure 2 (b), the vertical profile measured by the LIDAR-LNG corresponds to the red rectangle in Figure 3b. The LIDAR-LNG detected 175 a monolayer cloud around 5.5 km. In the panels (c) and (d), we present colored compositions of total and polarized radiances obtained from three spectral bands of OSIRIS over this cloud scene. The white concentric contours represent the scattering iso-angles in a step of 10º.
The cloud strongly backscatters the total solar radiation at the three visible wavelengths, producing an intense white signal. On the polarized image (Figure 2d), we observe a strong directional signature of the signal, characteristic of 180 scattering by spherical droplets. The main structure is the peak of polarization around 140º (known as the cloud bow), which forms a white arc in Figure 2d. At larger scattering angles, we observe the supernumerary bows whose positions vary with the wavelength, alternating between the red, blue, and green channels. The measured polarized signal for scattering angles smaller than 130º is largely dominated by molecular scattering at 490 nm, hence the blue color. The ocean surface modulates the scattered total and polarized radiance and increased the signal in the specular direction as

Retrieval methodology
One of the most robust approaches in cloud property retrievals is the optimal estimation method (OEM). It is increasingly used in satellite measurement inversion (Cooper et al., 2003;Poulsen et al., 2012;Sourdeval et al., 2013;195 Wang et al., 2016). It provides a rigorous mathematical framework to estimate one or more parameters from different measurements. The OEM also characterizes the uncertainty on the retrieved parameters while taking into account the instrument error and the underlying physical model errors. A complete description of the optimal estimation method for atmospheric applications is given by Clive D. Rodgers (Rodgers, 2000). In this book, Rodgers described exhaustively the information content extraction from measurements, the optimization of the inverse problem, and the 200 solutions and error derivations. In the following, we will go through the basics of this method that define the core of our retrieval algorithm.

The formalism of the optimal estimation method
Considering a vector (of dimension ny) containing the measurements and a state vector (of dimension nx) 205 containing the unknown properties to be retrieved, these two vectors are connected by the forward model , which can model the complete physics of the measurements to an adequate accuracy. The errors associated with the measurement and the modelling are represented by the error vector . Eq. (1) states the relationship between these variables. Note that, hereafter, bold variables represent vectors or matrices.
The OEM aims to find the best representation of parameters that minimizes the difference between simulations ( ) 210 and observations while considering the linearity of the direct model near the solution. Therefore, the best estimate of the properties vector corresponds to the minimum of the so-called cost function ( ): The first term of ( ) represents the difference between the measurements and the forward model calculated for a given state vector , weighted by the covariance matrix associated with the measurement error and the forward model.
The second term represents the difference between the state vector and the a priori state vector a weighted by a 215 the covariance matrix associated with a . In line with the cost function, the optimal state estimation emerges from a balance between the information carried by the measurement about the state and what we already know about it before the measurement. In our case, we do not have a prior estimate of the state vector. The iterations are initiated by a first guess while applying a large a . The difference between the measurements and the forward model will be the decisive element in the minimization of the cost function. It will ensure that the estimated cloud properties have the optimal fit 220 with the observed system only.
The minimization is done through the Levenberg-Marquardt approach (Marquardt, 1963;Levenberg, 1944) based on the "Gauss-Newton" iterative method. Each iteration is calculated following the Eq. (3): where i is the state vector at the ℎ iteration, is the sensitivity (or Jacobian) matrix and is the covariance matrix 225 of the state vector defined in Eq. (4).
The parameter affects the size of the step at each iteration. If the cost function increases at an iterative step then is increased and a new smaller step ( +1 ) is calculated until the cost function decreases.
The iterative process stops when the simulation fits the measurement (Eq. (5)) or when the iteration converges (Eq. (6)). The left side of Eq. (5) represents the cost function without taking into account the a priori negligible contribution 230 and should be smaller than to stop the iterations. Eq. (6) deals with the iterative steps and will make sure that the iterations will stop when the difference between two successive steps weighed by is less than . In other words, when further changes in the state vector have small to zero changes in the minimization.

Basic setting of the retrieval algorithm
In order to apply this theoretical framework into our retrieval algorithm, we define next the basic elements stated in the previous subsection.
The state vector contains the properties to be retrieved. In our case, they are the cloud optical thickness (COT) and the effective radius of water droplets (Reff): The measurement vector contains the radiances (R) measured by OSIRIS at two wavelengths and for several view directions represented in Eq. (8) (8)by λ and respectively.
The forward model is based on the adding-doubling approach (De Haan et al., 1987;Van de Hulst, 1963) to solve the 245 radiative transfer equation and simulate the radiances measured by OSIRIS for the corresponding observation geometries and wavelengths. It is a major element of the retrieval and describes the radiation interaction with the cloud, the surface, and the atmosphere while fixing several parameters (e.g. wind speed, cloud altitude…). We assume a standard atmosphere with a mid-latitude summer McClatchey profile (McClatchey et al., 1972) for the computation of molecular scattering. As all the channels used in the retrieval are in atmospheric windows (as seen in Figure 1), the 250 atmospheric absorption is not accounted for. Our case study is purely above an ocean surface. The reflection by the surface can affect the measured radiances even in cloudy conditions and particularly for the optically thin clouds. The anisotropic surface reflectance of the ocean surface is characterized by a bidirectional polarization distribution function (BPDF). We used the well-known Cox and Munk model to compute the specular reflection modulated by ocean waves (Cox and Munk, 1954) with a fixed ocean wind speed based on measurements of the National Oceanic and Atmospheric 255 Administration (NOAA).
As in current operational algorithms, the cloud model used for the retrieval is a plane-parallel and homogeneous (PPH) cloud, which implies the independent column approximation (ICA). The case study is a liquid water cloud scene. Therefore, we used a log-normal distribution of liquid spherical particles (Hansen and Travis, 1974)  The Jacobian matrix includes the partial derivatives of the forward model to each element of the state vector (Eq. (9)). The columns of define then the sensitivity of the radiances (each with a specific wavelength -viewing angle 265 configuration) to COT or Reff. The rows of the Jacobian define the sensitivity of each radiance configuration to the two retrieved properties.

Error characterization
During the retrieval process, every element is associated with a random or systematic error embedded in the error 270 covariance matrix . The implantation of errors in the inverse problem adjust the solution from just a unique value of to a Gaussian probability distribution function (PDF) where is the expected value and is its covariance.
is calculated after a successful convergence using the Jacobian at the retrieved state and .
It leads to a representation of the uncertainty on a particular parameter defined as the square root of the 275 corresponding diagonal element σ = √ , where is the index of the parameter in the state vector (Eq. (10)).
We chose to express this uncertainty using the relative standard deviation (RSD) in % (Eq. (11)). The RSD will be used to characterize the quality of the retrieval. (11) regroups measurement errors ( ) and forward model errors. Indeed, the forward model uses ancillary information provided by a set of fixed parameters (listed in section 3.3.2). Errors related to an uncertain estimation 280 of these fixed parameters are represented by the covariance matrix described in the next section.

RSD = ( ) × 100
Besides the fixed parameters, the cloud model used in the radiative transfer computation adds a source of uncertainty.
The uncertainties on the retrieved parameters related to these approximations are regrouped in the covariance matrix described in the next section. is then addressed as a sum of these three components: Previous studies (Wang et al., 2016;Iwabuchi et al., 2016;Poulsen et al., 2012;Sourdeval et al., 2015)  the covariance matrices of errors from the measurement space into the retrieved state space (Rodgers, 2000). The gain matrix representing the sensitivity of the retrieval to the measurement is used: The total variance-covariance matrix of the retrieved state vector ( ) can then be decomposed into three contributions (Eq. (14)), with each term originating from its corresponding error covariance matrix.
Each term in this equation is developed and discussed in the following three subsections. 295

Uncertainties related to the measurements
Any type of measurement is subject to errors. It is necessary to apply calibration processes to study the relationship between the electrical signals measured by the detectors and the radiances and quantify its uncertainty. Calibration is done during laboratory experiments before the airborne campaign or the instrument launch into space (Hickey and Karoli, 1974) or in-situ if calibration sources are available onboard the sensor (Elsaesser and Kummerow, 2008). The 300 uncertainty of the measurements determined during calibration is usually random, uncorrelated between channels, and can be quantified as probability density function over the measurement space.
As errors between measurements are supposed to be independent, the covariance matrix of measurement noise ( mes ) is diagonal with dimensions equal to the measurements vector dimension ( × ). The diagonal elements 2 are the square of the standard deviation of measurement errors. In our retrievals, we calculated the covariance matrix based 305 on 5% of measurement errors that cover the measurement errors of OSIRIS: σ mes = R , × 5%.
The error covariance matrix for the retrieved parameters due to measurement errors is then expressed by mapping the covariance matrix mes from the measurements space to the state space by using the gain matrix : The uncertainty on a particular parameter originating from the measurements error is defined as the square root of the corresponding diagonal element corresponding to the standard deviation = √ . It is expressed using the RSD (mes) as in Eq. (11).

Uncertainties related to the fixed parameters
Every column in fp and σfp is dedicated to one fixed parameter. Therefore, we can separate the contributions of every element of the fixed parameters vector as follows: Each covariance matrix from the right side of Eq. (17)   The uncertainties related to the cloud vertical homogeneity and the cloud horizontal homogeneity are quantified separately. In the following, we present the elements of the forward model used to quantify the uncertainties of these 360 assumptions.

Non-uniform cloud vertical profile model
The vertically heterogeneous cloud model to assess the uncertainties of the assumed homogeneous cloud model is described by:

365
-an effective radius profile and possibly an effective variance profile but for simplification, we will consider that veff is constant over the entire vertical profile with a value of 0.02.
-an extinction coefficient (σ ext ) profile -a cloud geometrical thickness (CGT) characterized by the difference between the altitude of the cloud top ( top ) and the cloud base ( bot ). The values of CGT, top and bot are fixed based on the LIDAR measurements. 370 The effective radius and extinction coefficient profiles are computed using an analytical model already presented in (Merlin, 2016). It is based on adiabatic cloud profiles, which are described and used in several studies (Chang, 2002;Kokhanovsky and Rozanov, 2012). In the adiabatic scheme, the effective radius increases with altitude. However, several studies proved that a simple adiabatic profile is not sufficient to describe a realistic cloud profile (Platnick, 2000;Seethala and Horváth, 2010;Nakajima et al., 2010;Miller et al., 2016) because turbulent and evaporation 375 processes can reduce the size of droplets at the top of the cloud. The description of a more realsitic vertical cloud profile is obtained with two adiabatic profiles (Figure 3) that are joined at the altitude of maximum LWC ( max ): -The first profile from bot to max is considered adiabatic.
-The second profile from max to top follows an adiabatic LWC profile decreasing with altitude.

Figure 3: The heterogeneous vertical profile of effective radius (black line) and extinction coefficient (blue line) used to assess uncertainties due to the assumption used for the vertical profile. The equivalent homogeneous vertical profiles are shown in dashed lines.
Considering that LWC is equal to zero at the base and top of the cloud, and relying on the linear variation model of the 380 LWC with altitude (z) established in (Platnick, 2000), we can write that: The profiles of effective radius (Eq. (22)) and extinction coefficient (Eq. (23)) can then be computed by considering that the particle concentration is constant over the entire cloud which makes it possible to obtain analytical functions of LWC, R eff and σ ext . This unitless parameter varies from 0 to 1 representing the shape of the profile. The value 0 corresponds to max = top (adiabatic cloud) and the value 1 corresponds to max = bot (a reverse adiabatic cloud with a negative gradient 390 of water content). In the following results, a value of 0.15 is assigned which allows to have a profile close to the one studied in (Miller et al., 2016) from large-eddy simulations (LES) cloud scenes.
To assess the error due to the vertical heterogeneity of the cloud, we need to specify the maximum value of the extinction coefficient σ ext max and the effective radius R eff max of the vertically heterogeneous cloud, corresponding to the "equivalent" homogeneous clouds. We use Eq. (25) to assign σ ext max which leads to the same integrated exctinction 395 profile and Eq. (26) to assign R eff max to what ensures that the mean R eff of the heterogeneous vertical profile is equal to the R eff of the homogeneous cloud (R eff F ). σ ext max and R eff max are found analytically by integrating the profiles described in Eq. (22) and Eq. (23).
A vertically heteregenous cloud is computed for each pixel using the retrieved value based on the homogeneous assumption. The error-covariance matrix describing the error due to the simple homogenous cloud assumption (Eq. (20)) is calculated from the difference between radiances computed with homogeneous and heterogeneous vertical profiles, denoted and ′ respectively. 405 The 3D radiative transfer model The other assumption that might affect the retrieved cloud optical properties in the current operational algorithms is the horizontally plane-parallel and homogeneous (PPH) assumption for each observed pixel. It implies that each pixel is horizontally homogeneous and independent of the neighboring pixels (known as the independent pixel approximation 410 (IPA)). The homogeneous IPA assumption affects the cloud-top radiances and leads to differences between 1D and 3D radiances that are the result of several effects discussed in numerous publications and briefly summarized in section 1.
To assess the uncertainties in the retrievals arisen from this assumption, Eq. (20) is used. F is then the difference between the radiances computed with a 1D radiative transfer code (1D-RT), following the Adding-doubling scheme (Hansen and Travis, 1974), and the radiances computed with a 3D radiative transfer (3D-RT) code called 3DMCPOL 415 (Cornet et al., 2010). The 3D simulations used for each pixel, the COT and Reff retrieved using the PPH assumption.
Errors on cloud model assumptions are assessed independently so vertical homogeneous profiles are assumed. We also assume a flat cloud top, differences, and errors are thus minimized as cloud top variation may increase the differences between 3D and 1D radiances (Várnai and Davies, 1999;Várnai, 2000). Our strategy to assess the different types of uncertainty follows two steps. In a first step, we retrieve COT and R eff using a bispectral multi-angular method by considering the uncertainties related to the measurement errors alone. We use a weak absorbing channel centered at 1240 nm that is mainly sensitive to COT, and a channel partially absorbed by the water droplets centered at 2200 nm, and thus sensitive on Reff. In this case study, up to 13 viewing angles are available for each pixel. In the first step, only the measurement errors are accounted for and included in . This error 425 is straightforward and usually well characterized but does not change once the measurements are realized. Not considering the other errors at this stage allows benefiting from a faster retrieval algorithm without the calculation of and the heavy computation cost of heterogeneous cloud profiles and 3D-RT calculations. The second step consists of computing the error due to ertically and horizontally inhomogeneous cloud from the retrieved R eff and COT.
In Figure 4, COT (a) and Reff (b) retrieved from multi-angular SWIR radiances are presented. Spatial variations are 430 mainly due to variations in the observed cloud structures. The COT range is between 0.5 and 6 with a mean value of 2.1 while Reff varies between 2 and 24 µm around a mean value of 8.8 µm. As detailed in section 3.3, the final error is divided into three categories. Figure 5 shows the uncertainties originating from a 5% measurement error on the retrieved COT (RSD COT (mes)) and Reff (RSD Reff (mes)). RSD COT (mes) ranged from 0.5 to 5% with a mean value of 3.2% while RSD Reff (mes) ranged from 2 to 12% with a mean value equal 435 to 6.2%. Both cloud properties have relative uncertainties that increase mainly with the magnitude of the retrieved values. This is related to the quasi-linearity of the radiance as a function of COT and Reff in this cloud regime (small COT). When the radiance-COT (or Reff) relationship is quasi-linear, the sensitivity of the forward model to COT (or Reff) is high, which will consequently lead to parameters retrieved with high accuracy (low RSD). As COT (or Reff) increases, the gradient of the radiance-COT (or Reff) relationship decreases causing larger uncertainties.  The second type of uncertainty is related to the fixed parameters in the forward model. In Figure 6, we show the uncertainty on COT and Reff in % due to an incorrect estimation of each fixed parameter in the forward model. Panels The (c) and (d) panels in Figure 6 represent the uncertainties on the retrieved COT and Reff originating from the fixed effective variance of droplets. They are nearly null on COT with a mean value of 0.05% as the 15% uncertainty on the value of veff (0.02) does not modify the total radiances. On the other hand, RSD Reff (veff) reaches 0.5% with a mean 450 value of 0.25%. Indeed, veff modifies the width of the cloud droplet distribution and consequently slightly the absorption by cloud droplets, which gives information to retrieve Reff. For Reff higher than 15 µm, the relationship between SWIR radiances and Reff tends to flatten which make them less sensitive to veff and thus the uncertainties are smaller than 0.1%.
The panels (e) and (f) in Figure 6 shows that an error in the estimation of the ocean wind speed affects the retrieved 455 COT and Reff mainly for small COT. The water-air interface is reflecting mainly in the direction of specular reflection, but the ocean is not perfectly smooth; the surface is modulated by waves that enlarge the directions of reflection due to the specular reflection. These waves are mainly formed by the wind. The higher the surface wind speed is, the greater the amplitude of the waves leading to a larger reflection angle (wider sun-glint). The Sun-glint reflection is seen by OSIRIS only for very small values of COT and implies uncertainties on the retrieved parameters of about 0.5%. At 460 higher COT, the surface is non-apparent to OSIRIS measurements, and uncertainties are thus close to zero.
We note that all the uncertainties of the studied fixed parameters remain below 1%, which shows that retrieval of all the COT-Reff couples does not have a high dependence on the fixed parameters.

Figure 6: The uncertainties (RSD (%)) on COT (left column) and Reff (right column) originating from the ancillary parameter errors: altitude (a and b), the effective variance of water droplet distribution (c and d), and the surface wind speed (e and f).
The uncertainties due to the assumptions of the forward model are presented in Figure 7. The panels (a) and (b) 465 represent the uncertainties on COT and Reff respectively, originating from the vertically homogeneous assumption.
RSD COT (Fpv) ranges between 1 and 8% with a mean value of 4.9% while RSD Reff (Fpv) varies from 2 to 20% with a mean value of 13.3%. We note that when the cloud is optically thin, RSD COT (Fpv) and RSD Reff (Fpv) tend to be lower. When the extinction is small, the radiations penetrate deeper into the cloud and bring information on the whole cloud, similar to the one obtained with the homogenous vertical profile. The differences between radiances issued from  The uncertainties originating from the use of a 1D radiative transfer code instead of a more realistic 3D radiative 475 transfer are represented in Figure 7 (c) and (d) for COT and Reff respectively. RSD COT (F3D) ranges between 1 and 20% with a mean value of 4.35%, while RSD Reff (F3D) varies from 2 to 18% with a mean value of 9.25%. At high spatial resolution, these differences can be caused by several effects that can add or oppose each other. Considering the solar zenith incidence angle (59º), illumination and shadowing effects can be present depending on the viewing geometries and roughness of the radiative fields (Várnai, 2000). 480 However, in this work, we are dealing with multi-angular measurements where the same cloudy pixel is viewed under different viewing angles, which tends to mitigate the influence of illumination and shadowing effects. Therefore, contrary to MODIS mono-angular retrieval (Várnai and Marshak, 2002), 3D effects due to solar illumination do not appear in the retrieved cloud properties.
The plane-parallel and homogenous (PPH) approximation (Cahalan et al., 1994) strongly depends on the sensor spatial 485 resolution (Oreopoulos and Davies, 1998).  In Figure 9, the histograms of the relative difference between the radiances computed in 1D (R1D) and the radiances 495 computed in 3D (R3D) at 1240 nm for different bins of optical thickness are plotted. We can see the shift of the histograms from negative value for small optical thickness (R1D < R3D) towards positive difference for larger optical thickness (R1D > R3D) that is explained by the horizontal radiation transport between columns. Overall, we note that the uncertainties due to the forward model assumption are much more important than the one due to the fixed parameters. The forward model is not sensitive to small variations in the fixed parameters. However, while assessing uncertainties due to the vertical profile or radiative transfer assumptio, we change the parameters that our forward model is proven to depend on, thus minor changes in the integrated profile can lead to relatively large variations in the radiance fields, and consequently large uncertainties. 505

Advantage of using multi-angular versus mono angular information
The same strategy applied in section 4 is applied using a bispectral mono-angular method used for the MODIS instrument for example. Using the MODIS-like approach, the measurement vector for each pixel contains two monoangular total radiances, one at 1240 nm and the other at 2200 nm to retrieve COT and Reff. The retrieved COT over the whole field varies between 1 and 12 with a mean value equals to 3.44. The range of 515 retrieved Reff has a mean value of 15.65 µm affected by the high value of Reff retrieved around the scattering angles 130-140° where the sensitivity of 2200 nm radiances to the water droplet size is known to be small. As a matter of fact, (Cho et al., 2015) show indeed that in liquid marine cloud cases, the phase functions of different Reff converge to the same value for these scattering angle range leading to the failure of water droplets size retrieval from MODIS measurements. This reduced sensitivity explains also the high uncertainty on Reff due to measurement errors around 520 the cloud bow ( Figure 11). The reduced sensitivity on Reff, in this case, is not limited to the cloud bow directions and supernumerary bows but is also visible at small scattering angles (70-80°) that can be affected by specular reflection over the ocean.
525 Figure 11: Uncertainty on the effective radius originating from the measurement error, RSD Reff (mes) as a function of the scattering. The red line represents the mean RSD Reff (mes) = 12.55%.
Multi-angular retrieval presents the major advantage, that no aberrant values of Reff are retrieved near the scattering angles at 140° (comparing Figure 4b to Figure 10b). Clearly, the multi-angular measurements contain more information and allow to resolve the problem encountered with the MODIS-like method which is also clear in the reduction of 530 failed convergences from 5.9% to 2.4% (not shown here). In the overall scene, smaller Reff values (and consequently lower COT values) are obtained.
The relation between radiances and COT/Reff being monotonical, the mono-directional method always allows us to find a retrieved value as it is always possible to find a cloud model (a pair of COT and Reff) that matches a single measured radiance of a given target. However, this value can be more or less far from the real value and is not 535 necessarily an indication of an accurate retrieval, but only that a fit occurred. On the other hand, multi-angular retrieval makes the constraining of the forward model much more challenging. It uses jointly all the available measurements.
The retrieved state is then consistent at the best with all the measurements associated with different viewing angles. On this basis, it is not suitable to compare the cost function of a mono-directional method where one value of COT fits one radiance and the cost function of a multi-angular method where one value of COT corresponds to ~13 radiances. 540 To overcome this problem, we compare the accuracy of the retrievals by the response of the state vector to the errors, represented by the relative standard deviation (RSD). In Figure 12, we present the spatially average of the different types of errors, presented in section 4, for the mono-angular method (light green for COT and light blue for Reff) in comparison with the multi-angular method (dark green for COT and dark blue for Reff). We divide the source of errors in two panels, the left panel groups the lowest values of RSD and the right panel for the highest values of RSD. 545 Overall, Reff uncertainties are larger than the one on COT for any type of error. In the left panel of Figure 12, the three fixed model parameters errors related to an incorrect estimation of the fixed parameters of the model are weak compared to the others and remain below 0.3% for mono-angular retrievals. As explained in section 4, the fixed altitude does not contribute to the uncertainty on the two retrieved parameters. The average uncertainties originating from the fixed value of veff are about 0.05% for COT and slightly higher (0.15%) for Reff since veff affects the cloud bows that 550 are also sensitive to Reff. Concerning the surface wind speed, the uncertainties are around an average of 0.05%.
In the right panel, for mono-angular retrieval, the measurement error contributes to the uncertainty of about 8% on the retrieved COT and of about 12% on the retrieved Reff. The uncertainties are reduced by a factor two with multi-angular retrieval. The multi-angular approach leads indeed to more information available for each cloudy pixel and each additional information reduces the uncertainty on the retrieved parameters in the presence of the same 5% random 555 noise in the measurements.
The following two groups of bars correspond to the to the error introduced by the homogeneous assumption in the forward model. They are the main source of errors. For mono-angular retrieval, the assumption of a vertical homogenous profile contributes to an uncertainty of about 16% on COT and 54% on Reff. These uncertainties are reduced by a factor four in the case of multi-angular retrieval. As discussed previously, the principal effects of 1D 560 assumptions errors at the spatial resolution of OSIRIS come from the non-independence of the cloud columns that lead to smooth the 3D radiative fields. They lead to an uncertainty of 28% on COT and 45% on Reff when a mono-angular instrument is used.  The multi-angular approach provides additional information for each pixel and constrains the forward model to match 565 all the angular radiances at once. Indeed, the OSIRIS multi-angular characteristics have the advantage of decreasing the angular effects around the cloud bow directions by adding the contribution of other geometries. It avoids most of the failed convergences that occurred with the MODIS-like method and retrieved more homogeneous and coherent COT and Reff fields.

570
In this study, we present a method to retrieve retrieve two important microphysical and optical parameters of liquid clouds, COT and Reff and their uncertainties using NIR/SWIR multi-angular airborne measurements. The algorithm is based on the mathematical framework of the optimal estimation method (Rodgers, 2000) and focuses on assessing the different uncertainties on the retrieved properties originating from different sources of errors.
The studied case uses the measurements of the airborne radiometer OSIRIS obtained during the CALIOSIRIS 575 campaign. It consists of a monolayer water cloud located at 5 km altitude over the ocean with tilted solar incidence ( =59º).
In the first step of the retrieval, COT and Reff are retrieved by considering only the measurement error (without introducing any error linked to the forward model).The uncertainties originating from different sources of error are computed afterward by using the previously retreved COT and Reff, and are decomposed into 3 different sources related method and for comparison a mono-angular method, which is the usual approach in the operational algorithm called here the MODIS-like algorithm. 585 In the multi-angular retrieval, a 5% measurement error contributes to around 3% of uncertainty on the retrieved COT and 6% on the retrieved Reff. It increases with increasing values of COT and Reff for which the sensitivity of radiances starts to decrease. These uncertainties are doubled when considering a mono-angular retrieval.
The uncertainties related to the fixed parameters remain low with both mono or multi-angular retrieval. The largest one is due to the miss-knowledge of the effective variance of the droplet size and is respectively equal to 0.15% and 0.25% 590 for the mono and multi angular case.
This study clearly shows that the largest uncertainty is due to the homogeneous cloud assumption made in our forward model. First, the uncertainties related to the homogeneous vertical profile were quantified using a heterogeneous LWC profile with a triangle shape (known as quasi-adiabatic) composed of two adiabatic profiles. This more realistic profile takes into account the transition zone at the top of the cloud related to turbulent and evaporation processes. The scene 595 averaged values reach 5% and 13% for COT and Reff respectively in the multi-angular retrieval of our case study and go up to 13% and 54% for COT and Reff respectively when using mono-angular measurements. The largest uncertainties are obtained for the largest cloud optical thickness as the radiations sample only the higher layers of the cloud where the information is different between the homogeneous and heterogeneous vertical profiles.
The other source of uncertainty related to the simplified cloud physical model comes from the radiatively non-600 independence of cloudy columns that dominate at the high spatial resolution of OSIRIS. In the optically thin overcast cloud case, studied here, the scene average uncertainties originating from the 3D effects are 4% for COT and 9% for Reff when using multi-angular measurements, and 28% for COT and 45% for Reff when using mono-angular measurements. The non-independence of the cloud columns dominates and tends to smooth the 3D radiative field compared to radiances computed with the independent pixel approximation. Moreover, at the high solar incidence 605 (59º), illumination and shadowing effects can be present. However, the use of multi-angular measurements mitigated their effects. Consequently, the uncertainties are lower when using multi-angular measurements.
The method presented here addressed the exploitation of OSIRIS and can be easily adapted to the future 3MI imager.
The first step that consists of including the uncertainties related to the measurement errors is implementable in an operational algorithm. The second step that consists of computing the uncertainties resulting from the cloud model is 610 computationally expensive but can be applied to case studies to obtain a distribution of the error according to the cloud characteristic. In addition, the results obtained in this study show, not surprisingly regarding the numerous studies already published, that the vertical and horizontal homogeneity assumptions are major contributors to the retrieval uncertainties. One way to reduce it would be to define a more complex cloud model that can take into account the vertical and horizontal heterogeneity. This adds more complexity in the forward model as it would implies to retrieve 615 more sophisticated cloud parameters (e.g. extinction or effective size profile). It appears however possible given the important and complementary information provided by OSIRIS or 3MI measurements but more sophisticated retrieval methods need to be developed. Recent studies were proposed to retrieved vertical profile using cloud side information (Ewald et al., 2018;Alexandrov et al., 2020) or to realize multi-pixels retrieval to account for the non-independence of the cloudy pixels (Martin et al., 2014;Okamura et al., 2017;Levis et al., 2015).