Combining visible and infrared radiometry and lidar data to test simulations in clear and ice cloud conditions

. Measurements taken during the 2003 Paciﬁc THORPEX Observing System Test (P-TOST) by the MODIS Airborne Simulator (MAS), the Scanning High-resolution Interferometer Sounder (S-HIS) and the Cloud Physics Lidar (CPL) are compared to simulations performed with a line-by-line and multiple scattering modeling methodology (LBLMS). Formerly used for infrared hyper-spectral data analysis, LBLMS has been extended to the visible and near infrared with the inclusion of surface bi-directional reﬂectance properties. A number of scenes are evaluated: two clear scenes, one with nadir geometry and one cross-track encompassing sun glint, and three cloudy scenes, all with nadir geometry.CPLdata is used to estimate the particulate optical depth at 532 nm for the clear and cloudy scenes and cloud upper and lower boundaries. Cloud optical depth is retrieved from S-HIS infrared window radiances, and it agrees with CPL values, to within natural variability. MAS data are simulated convolving high resolution radiances. The paper discusses the results of the comparisons for the clear and cloudy cases. LBLMS clear simulations agree with MAS data to within 20% in the shortwave (SW) and near infrared (NIR) spectrum


Introduction
A recent review of the light scattering properties of cirrus (Baran, 2009) points out that it is more desirable to construct cirrus ice crystal models that predict the light scattering properties of non-spherical ice crystals that can be applied at any wavelength rather than at particular wavelengths. It also points out that the choice of ice crystal model, beside its importance for climate modeling, is also important for the space-based remote sensing of cirrus properties, since inappropriate choice of the scattering phase function may lead to errors in retrieved optical depth of several factors. The review contains a wealth of references pertaining to this problem.
The issue of the quality of cloud products that are now routinely produced from satellite data is addressed in Ham et al. (2009) that examines the quality of the Moderate Resolution Imaging Spectroradiometer, MODIS, retrieved cloud products. These are used as input to a radiative transfer model to compute multiply scattered radiances at a number of MODIS channels and comparing the measurements to the simulations. The main findings are that radiances for shortwave bands between 0.466 and 0.857 µm appear to be quite accurate, while simulated radiances for the 1.24, 1.63 and 3.78 µm bands do not well agree with measurements. Large differences between simulations and data are also found in the infrared window bands (such those centered at 8.56, 11.0 and 12.0 µm).
In Yang et al. (2007) the differences of the bulk optical properties of ice clouds retrieved in MODIS collection 4 and 5 are investigated and it is shown that collection 5 optical thickness over ocean are a factor 1.9 higher than the collection 4 counterpart. Moreover it is stated that the differences can lead to either an enhancement or a reduction of the warming effect of ice clouds, depending on the specific ice cloud of interest.
In the Zhang et al. (2009) paper, the main concern is the influences of different ice particle micro-physical and optical models on the resulting optical thickness retrievals from satellite measurements of solar reflection. They find that the ice cloud optical thickness retrieved from POLDER is substantially smaller than that from MODIS, and the difference is attributed primarily to the difference of asymmetry factor used in the two retrievals. They conclude that ice cloud optical thickness retrievals based on satellite measurements of solar reflection are highly sensitive to the choice of the ice particle model assumed in the retrieval.
The present study was initiated to evaluate the quality of a forward modeling methodology, called Line-By-Line Multiple Scattering (LBLMS), that is an extension to the shortwave of a state-of-the-art methodology already used to simulate high resolution spectral data in the infrared spectral range, from 3000 to 50 cm −1 (Rizzi et al., 2001;Amorati and Rizzi, 2002;Maestri and Rizzi, 2003;Tjemkes et al., 2003;Rizzi and Maestri, 2003;Maestri et al., 2005).
Two diverse data sets are used. The first is a field study during the 2003 Pacific THORPEX (The Observing System Research and Predictability Experiment) Observing System Test (P-TOST, http://angler.larc.nasa.gov/thorpex/). During P-TOST the spectro-radiometric data are combined with lidar products that describe the particulate (aerosol and clouds) extinction profiles. Since our long-term interests are on the retrieval of cloud variables, the experimental cases selected for the core study involve radiances measured in presence of a cloud layer over the marine surface. However the methodology and results for two clear scenes over the sea are also included, since nadir and cross track simulations of a clear scene provide evidence of correct modeling of the surface and of the aerosol layers, especially in the shortwave (SW) part of the spectrum. This is especially important when dealing with thin clouds in order to avoid that incorrect simula-tions of surface or aerosols properties affect the cloud properties retrieval process.
The second data set is from the Moderate Resolution Imaging Spectroradiometer (MODIS) and corresponding retrieved cloud products. This data set is used to clarify the unexpected results obtained from the P-TOST cloud cases, and is not intended to present results of statistical significance.
The description of the experiment and of all the case studies is given in Sect. 2, together with the details of the modeling methodology. The results for the clear P-TOST case are presented in Sect. 3, the three P-TOST cloudy cases are presented in Sect. 4, and the MODIS computations are discussed in Sect. 5. An overall discussion of the results with conclusions follow in the last section.

Description of the experiment and instruments
The P-TOST data-set was measured on 22 and 23 February 2003 during a flight over the Pacific Ocean SE of the Hawaiian Islands, when the high altitude NASA aircraft ER-2 carried the three instruments of interest for our analysis: S-HIS, MAS and CPL.
The Scanning High-resolution Interferometer Sounder (S-HIS) (Revercomb et al., 1998) is a Fourier-transform spectrometer with laser-controlled sampling, operating in the thermal spectrum between 3.3 µm and 17.2 µm (3000-680 cm −1 ); it utilizes a 45 • scene mirror that rotates through a measurement sequence consisting of views of the earth and two calibration sources, one at ambient and another up to 60 K above ambient. During each scan, 11 cross-track Field of Views (FOVs) are sampled (±35 • total view angle), with a nadir spatial resolution of 2 km from the nominal 20 km ER-2 flight level and a spectral resolution of 0.5 cm −1 .
The MODIS Airborne Simulator (MAS) (King et al., 1996) was built as support to the development of the MODIS satellite instrument: it is a 50 channel scanning spectrometer that covers the spectral range from visible to thermal infrared and acquires 50×50 m (at nadir) pixel data across a 37 km swath (±43 • view angle) from the nominal ER-2 flight level. At the nominal ER-2 ground speed of 210 m/s, MAS FOVs show an along-track superposition of the scan lines of about 33% of the pixel width.
The Cloud Physics Lidar (CPL) (McGill et al., 2002) provides cloud and aerosol backscatter profile at 30 m vertical and 200 m horizontal resolution at 1064 nm, 532 nm, 355 nm. During the P-TOST experiment the CPL provided optical depths at 532 nm up to a saturation value of about 3.
In addition to the ER-2, the NOAA G-4 research aircraft flew carrying the Airborne Vertical Atmospheric Profiling System (AVAPS) to measure the temperature and humidity profiles from cruise level (12 km) to the ground.
The ER-2 Mission Report specifies that the aircraft has flown from 21:40 UTC on 22 February to 04:30 UTC on 23 February (11:40 to 18:30 local time) at a cruise altitude of 20 km and the NOAA G-4 took off 20 min after the ER-2. An overlook of the entire mission's route can be seen in the GOES image in Fig. 1. At the surface, a high pressure system extends north of the Hawaii islands, smoothly decreasing towards the equator. The analysis of the geopotential at higher levels (not shown) reveals an undulation of the tropical jet stream over the eastern Pacific. A ridge extends to the West of Hawaii whereas to the East a trough with an associated slowing down of the jet and upper-level divergence extends towards the equator. The inspection of the vertical velocities at 200 hPa (not shown) shows upwelling on the east side of the trough associated to the extended high cloud cover that can be observed in Fig. 1.
The ER-2 flew from Hickam AFB to get on WNW-ESE oriented track line (20 • N, 153 • W to 16 • N, 144 • W) designed to enter and transect subtropical jet zone running up from tropics between 150 • and 140 • W longitude. ER-2 did 2 back and forth runs (4 total segments) of this line with the G-4 doing profiling on a similar but shorter leg, releasing 11 drop-sondes. G-4 flew back after two ER-2 legs. The Western third of this line contained clear to partly cloudy (low cloudiness) skies. Middle and eastern ends were overcast with thick cirrus associated with the subtropical flow. Transect crossed core of subtropical jet positioned at about 145 • W.
The P-TOST campaign was designed specifically as a fine tuning opportunity of the MODIS and the Atmospheric Infrared Sounder (AIRS) science product algorithms. The decision to use the P-TOST data was taken shortly after the experiment was completed, and the choice was based on the following elements: 1) high spatial resolution radiometric data spanning a spectral interval from shortwave to infrared; 2) high spectral resolution infrared interferometric data colocated with 1); 3) a clear sky and an ice cloud scenery within the same mission; 4) detailed measurement of the atmospheric profile of temperature and humidity and 5) the availability of an airborne LIDAR (on same platform) to provide aerosol products for the clear case and some cloud products for the cloudy one. The campaign would have greatly benefited from a satellite overpass, but the closest Terra MODIS granule is at 22:25 of 22 February and covers the uppereastern part of the GOES image in Fig. 1, outside the experimental area. The closest Aqua MODIS granule is at 19:20 of 22 February and covers the eastern edge of the same GOES image, again outside the experimental area.
The P-TOST campaign provides, in the clear case, an interesting combination of wind pattern and sun glint to test the surface reflection properties of the modelling suite. Additional in-situ sampling of the microphysical features of the cloud layer would have been key observations but, to the best recollection of the authors, there were no dataset available at that time with the characteristics described above under 1-5), combined with in-situ microphysical measurements.

Modeling methodology
Although S-HIS and MAS flew together on the ER-2, the data-recording time of the two instruments had different reference, resulting in different time overpass over the same scene. Moreover, it had been noted (R. Holz, University of Wisconsin, Madison, personal communication, 2004) that S-HIS data had wrong geographic data positioning due to use of the inertial navigation system of the ER-2 as reference, which in the analyzed mission did not work properly. Therefore a scheme was developed to collocate MAS and S-HIS data using MAS channel 45 (a window centered at 11 µm). This channel performs a relatively more stable measurement of the scene radiance with respect to window channels located in NIR range, which are affected, moreover, by solar contribution. Two sets of virtual measurements are generated, having the spatial resolution of S-HIS and the spectral resolution of MAS: MAS pixels are averaged over the S-HIS footprint closest to nadir and the S-HIS data is convolved over the MAS spectral response function to produce the equivalent MAS spectral bands (Moeller et al., 2003). The minimisation of the mean square differences between the convolved S-HIS signal and the averaged MAS signal is performed over a 750 km long ER-2 track (corresponding to flight between 01:06 UTC and 02:00 UTC on 22 February) in which clear sky, broken clouds and overcast situations are present. The mean temporal displacement between MAS and S-HIS data is found to be 41±2 s. The uncertainty of 2 s produces a possible spatial displacement of about 400 m at the nominal ground speed of the ER-2, i.e. it is smaller than the linear dimension of the S-HIS nadir footprint.
The radiative transfer equation is solved for a planeparallel geometry and the atmosphere is divided into a number of layers (82), each assumed to have homogeneous scattering properties (that change from layer to layer). The simulations are done with a suite of codes collectively called Line-By-Line Multiple Scattering (LBLMS). In the present work the line-by-line computations of layer spectral optical depths are done using the Line-by-Line Radiative Transfer Model (LbLRTM) (Clough and Iacono, 1995). LbLRTM can solve the clear sky radiative transfer equation, but in our study is used to generate layer monochromatic optical depths (OD). These are interpolated at 0.002 cm −1 and then convolved to compute spectrally averaged OD.
Tests were done to find the appropriate spectral resolution for the specific sensors that are modeled (S-HIS and MAS in this work), so that the difference between radiances convolved using the averaged optical depths or directly the LbLRTM monochromatic radiances were below a given threshold. Results are dependent on the spectral widths over which the OD average is performed and also on the atmospheric layering. Three different spectral resolutions were used, as a compromise between a reasonable computing time and accuracy: 0.01 cm −1 from 580 to 3000 cm −1 (3.33-17.2 µm, indicated as IR in this paper), 0.05 cm −1 from 3000 to 7000 cm −1 (NIR, 1.43-3.33 µm), and 0.5 cm −1 from 7000 to 22 000 cm −1 (SW, 0.45-1.43 µm). The HITRAN 2004 (Rothman et al., 2005) spectroscopic database and the MT-CKD 1.3 water vapour continuum absorption model (Clough et al., 2005) are used.
The integration of the radiative transfer equation, including multiple scattering, is based on the code RT3 (Evans and Stephens, 1991). Layer spectral absorption optical depth is the sum of molecular and particle absorption and spectral total scattering optical depth is the sum of particle and Rayleigh scattering with the total phase function being a weighted mean of the two components.
The emissivity of the ocean surface in the IR is computed using the method described in Masuda et al. (1988). The ocean surface reflection properties take into account the contribution from the wind-roughened surface (Cox and Munk, 1954), the white caps and the up welling radiation scattered back from below the surface as function of the oceanic pigment concentration. The surface model implemented in RT3 is analogous to the one adopted in 6S (Vermote et al., 2006), with the only exception being the wind direction dependency that in RT3, due to limitations inherent to its computational structure, is not allowed. An in depth description of the methods adopted can be found in Bozzo (2009).
Post-processing of high resolution radiances to produce un-apodised spectra at S-HIS resolution is done as described in Rizzi et al. (2001). Simulated MAS radiances are obtained by convolving the spectral values at the indicated resolution with the MAS response functions (NASA, cited 2008).
Aerosol single scattering optical depth and phase function are computed using an in-house code, but the main physics contained (refractive index of components, aerosol as a collection of homogeneous spherical particles, size distributions, definition of external mixtures and growth coefficients) is analogous to the package OPAC (Hess et al., 1998) and so are the computed properties. The Mie scattering portion is handled by routines distributed together with the initial version of the RT3 code by Frank Evans.
The retrieval methodology (RT-RET) is described (and applied to Arctic and Mid latitude clouds respectively) in Maestri and Holz (2009) and Maestri et al. (2010). RT-RET uses LbLRTM and the same doubling and adding algorithm (RT3) described earlier. In this paper it is used to retrieve cloud optical depths and effective dimension of the cloud particle size distribution (PSD) from S-HIS radiances. The definition of effective dimension D e of a PSD made of nonspherical shaped particles is the same that was introduced by Foot (1988): where P (D) and V (D) are respectively the projected area and volume of a particle with maximum dimension D. The ice density ρ i is assumed constant with the particle dimension. Effective radius is defined as R e = D e /2. The single scattering properties for different ice habits in the short-wave (0.25-4.5 µm) and long-wave (4.5-100 µm) spectral ranges are taken from Yang and Liou (1998) and Yang et al. (2005) and will be referred to in the following by the acronyms SSP-SW and SSP-IR.
Previous work (e.g. Wendisch et al., 2005;Wyser, 1999) has shown large differences in simulations based on PSD composed of single shapes and there are no reasons to prefer one shape over another. Our choice is to use a mixture (called MIXML) of different habits, bullet rosettes, aggregates, droxtals and a small percentage of solid columns as described in Bozzo et al. (2008) and Fig. 1 of the quoted paper. MIXML is based on Lawson et al. (2006) measurements, developed in the context of mid latitude ice clouds. The major differences between mid-latitude and tropical cirrus clouds are found for cirri formed near strong tropical convective events, at the top of large anvils. Such ice clouds present usually higher fractions of larger particles than mid latitude-synoptic ice clouds, associated to strong updrafts . In case of synoptically generated midlatitude cirrus, the large ice crystals tend to subside quickly due to the weak updrafts. In our case, although the cirrus is associated with the subtropical jet-stream, it does not show the characteristics related to the tropical convective structures. MIXML has been tested in the analysis of infrared interferometric data collected at Mid Latitudes and in cloudy conditions during the italian phase of the EAQUATE experiment .
The ice particle size distribution (PSD) adopted for the forward and inverse computations is a Gamma distribution: where D is the particle maximum dimension, N 0 is the intercept value of the distribution, λ the slope (with unit of an inverse dimension) and µ is the dispersion (or shape) parameter. The PSD is extended to particles smaller than the one available in the quoted databases by assuming that spherical particles are present with diameter smaller than 2 µm, whose properties are computed using same Mie code as for the aerosol properties. Particular attention has been given to the treatment of the phase function to account for the sharp diffraction peak exhibited by the ice crystals. Since the scattering phase function is handled with an expansion in Legendre polynomials, such peaked functions would require thousands of components, which in turn would imply an extremely timeconsuming solution. Instead, the procedure proposed by Potter (1970) has been followed and the phase function peak is modified by using a spectrally varying algorithm which is tailored to minimise the number of Legendre terms required for an accurate reconstruction. The contribution of the deltatransmission peak associated with the forward scattering at a 0 • (Takano and Liou, 1988) is also accounted for. The number of Legendre terms used in the computation is linked to the number of angles employed in the zenith discretization in each hemisphere (N u ). After considerable tests we have used, for the ice cloud cases, N u = 32 (and 125 Legendre terms) in all spectral regions except from 600 to 4850 cm −1 where N u is set to 60 (237 Legendre terms). Since the CPU time required for LBLMS computation is proportional to the cube of N u , the accurate computation of high resolution radiances has required a massive computational effort that was made possible only after code parallelisation. As an example in Fig. 2 the original phase function, after truncation using Potter's method, is shown together with the one reconstructed with a number of Legendre coefficient as explained above. The figure shows that the fractional difference is an oscillating function that is generally below 10% except at some specific angles where it can reach values between 20% and 30%. The geometry of the nadir measurement in the cloudy case is such that the (single) scattering angle is 123 • . pixel long area, which covers a section 2 • wide (±1 • across the nadir line) and roughly 15 km long, along the blue line in Fig. 3. For the cross-track case MAS measurements are averaged over an area as wide as the full MAS swath and 120 pixel long (the red box on Fig. 3), which roughly corresponds to a 4 km long and 40 km wide region.

Measurement conditions for the clear cases
In both case studies the averaged MAS data is compared to the LBLMS simulation. The solar zenith angle is approximately 31 • and solar azimuth 204 • ; sensor azimuth angle swings between 204 and 24 • , hence aligned with the surface incidence plane of solar radiation.
Three drop-sondes were launched in coincidence with this track: one at 23:07 UTC, the second at 23:12 UTC and the third at 23:17 UTC. All drop-sondes measured a 10 m wind speed between 6 m/s and 7 m/s with azimuthal direction around 30 • , hence slightly skewed with respect to the instrument-sun plane: wind speed has been set in the simulation at 6.5 m/s. The vertical profile of temperature and humidity used for the simulation is a composite of the data from the 23:17 UTC drop-sonde and a standard tropical profile (Anderson et al., 1986) to fill the 10 km gap between the G-4 and the ER-2 cruise altitude. The CO 2 mixing ratio profile is modified to measured global mean value for the period of the campaign. Absorption from chloro-fluoro-carbon macro molecules (CFC) is also accounted for.
The CPL detected the presence of an aerosol layer distributed in the boundary layer from 1000 m down to about 500 m. Below this level there is no information on aerosol optical depth. Since the source of aerosol is the oceanic surface an exponential extrapolation of the measured optical depth profile down to the marine surface is assumed. The integrated optical depth at 532 nm obtained with this pro- cedure is 0.07. Because of the extremely low aerosol load measured by CPL, such an extrapolation must be considered only a minor weakness in the definition of our study case. As already mentioned there is no MODIS data over the experimental area and no aerosol product different from the CPL one is thus available. The aerosol layers (from surface to 1000 m height) are simulated with the optical characteristics of a maritime tropical aerosol model (same mixture defined in Hess et al. (1998), grown in an ambient with 80% relative humidity, in accordance with the drop-sondes measurements).
The pigment concentration of the oceanic water derived from the global products of the SeaWiFS satellite (Johnson et al., 1998), averaged over 8-days around the 22 February 2003 in the Pacific Ocean SE of the Hawaiian Islands, has a mean value of 0.07 mg/m 3 .

Results for the clear nadir case
Results for the nadir case are shown in the three panels of Fig. 4 in terms of radiance (all radiance values are given with unit W/(m 2 µm sr)) for the SW and NIR and brightness temperature (K) for the IR range. Figure 5 shows the relative errors in radiance between LBLMS and MAS data for the SW and NIR range and the absolute error in brightness temperature (K) for the IR range. The error bars added to the measured values represent one standard deviation of MAS data and denote the overall scene variability. In the MAS file no uncertainty is provided for the radiance data. The nominal errors associated with the MAS channels are listed in Table 2 of King et al. (1996). For a standard scene, the declared equivalent noise in the IR channels is smaller than 0.5 K, while in the solar range is of order of 0.2 W/(m 2 µm sr) between 1 and 0.5 µm and 0.02 W/(m 2 µm sr) between 1.6 and 2.4 µm. In this case the error bars are barely noticeable due to the stable signal provided by the clear sky scene. The channel centered at 1848 cm −1 did not work properly during the selected mission and has been ignored in the comparison.
The difference between the simulations and the MAS data is almost constant in the SW range between 0.45 and 1 µm with a value of about 1.5 W/(m 2 µm sr), hence the relative error increases with increasing wavelength since the radiance level is decreasing: it remains below 10% for wavelength shorter than 0.7 µm and it grows up to 35% approaching the channels at 0.95 µm. In the NIR range the relative error is around 10-15% with the exception of the strong H 2 O band, where it reaches values of about 90%.
In the IR, from 500 cm −1 to 2800 cm −1 , the absolute errors are between 0.5 and 2 K except in the 3 channels located in the strong CO 2 absorption bands (one at 700 cm −1 and the other two at 2250 and 2450 cm −1 ). The discrepancy in the opaque channels due to CO 2 absorption is certainly linked to the assumed atmospheric temperature profile . The surface skin temperature is set equal to the last temperature measured by the drop-sonde and this affects the simulation in the window regions. These differences could be easily eliminated to a large extent by improving the assumed atmospheric/surface temperature profile but in this case study we were not particularly concerned with the small deviations in the infrared range and we did not optimise those computations.
The largest relative differences in the SW and NIR are found in channels with important water vapour absorption. We have performed simulations also of the S-HIS data (not shown) and the results show a negative bias (the simulation being colder that the S-HIS data) in the infrared vibrorotational band of water vapour. Therefore the water vapour profile, and in particular the profile assumed between the G-4 and the ER-2 flight levels, is likely more humid than the true profile, which also explains the negative bias in the SW and NIR channels.
In conclusion the main causes of the observed discrepancies in the SW and NIR could be the modeling of the scattering by the oceanic surface, the assumed aerosol optical depth vertical profile, the aerosol mixture adopted and its growth properties. Each of these will be discusses in what follows.
To understand the importance of each process, the four panels of Fig. 6 show the upwelling radiance at MAS height simulated for four channels when a) the source function is the surface with a molecular atmosphere that only absorbs radiance (line labeled SUR in all panels), b) when molecular scattering is added as source of radiance (SUR+RAY) and finally c) when aerosol is also present (FULL). It is evident that in our case study the surface and molecular scattering are the dominant radiation sources.
Currently our modeling of the oceanic surface assumes that the sun glint pattern is independent on the wind direction. Comparisons with the 6S model indicate that this hypothesis gives an error in the range 10-20% of the signal for low solar zenith angles (Bozzo, 2009). In particular, in case of a zenith angle of 31 • the up-welling radiation is enhanced if the wind is blowing along the incident and reflected solar beam plane, whereas it is reduced in case of a wind blowing orthogonal to the sun's reflection plane. The wind's azimuthal angle retrieved from the drop-sondes measurements is around 60 • , hence 36 • from the sun's reflection plane, which lies along the 24 • -204 • azimuthal line.
The OPAC standard Maritime tropical aerosol model is an external mixture of three components, one of which, the sea salt coarse mode (SSCM) accounts for 1.4% of the mass of the mixture. The extinction optical depth has a maximum around 5 µm decreasing toward shorter wavelengths, markedly different from the sea salt accumulation mode component (SSAM) which has a maximum around 1 micron. Therefore one could devise a slightly different mixture, by increasing the mass of the SSCM and reducing the mass of SSAM, and the new mixture would reduce the underestimation observed in the NIR. Also data on condensational growth by aerosol components is sparse, and it is well known that assuming aerosols are spherically homogeneous particles is only a rough approximation. However the relevance of the aerosol contribution is so small in this case study that it is not necessary to dwell on the various sensitivity tests performed.
In conclusion the overall agreement of LBLMS simulations to MAS data appears to be good (mostly within 20% in the SW and NIR and 2 K in the IR), and certainly consistent with the approximations used. The small underestimation of LBLMS simulations with respect to MAS measurements in the SW can be related to the simplification of a wind direction independent bidirectional reflectance (BDRF). In the next section we examine the same scenario, but from different viewing angles.

Results for the clear cross-track case
LBLMS is tested with a full MAS cross-track scan swinging in the analyzed scene from the maximum to the minimum of the glint reflection region, along the reflection plane of the incoming solar radiation. Up-welling radiation at 20 km is simulated at 32 zenith angles for each hemisphere using 64 terms of the azimuthal-mode Fourier expansion. Pigment concentration, aerosol profile, model and relative humidity are same for the nadir case. Figure 7 shows the reflectance at 20 km for 4 MAS channels defined in the previous section. The BDRF for various wavelengths is computed from averaged MAS data (full lines) inside the red box in Fig. 3, and from LBLMS simulations (dashed lines).
The peak in reflectance is reached for all the 4 curves at around 36 • , slightly shifted toward the horizon from the Fresnel reflection point that is supposed to be at the incidence angle of 31 • . The glint pattern appears to be less steep and slightly more skewed towards the horizon in the MAS observations: in fact some difference in the glint pattern is expected due to the assumption of a wind direction-independent BDRF in the simulations.
The relative difference are between −10% to +10% over the whole scan for the SW channels, with a larger underestimation for the NIR channels, and are largest in the direction opposite to the reflection point, at minimum MAS radiance values.

Measurement conditions for the three cloudy cases
As already mentioned the ER-2 did 4 segments and the middle and eastern end of the flight leg were characterized by an extended cirrus associated with the subtropical jet-stream. The transect crossed the jet at about 145 • W and 19-20 • N. The ER-2 made two back and forth overpasses over the high, thick tropical cirrus.
As pointed out in the discussion of the clear cases, some assumptions adopted in the ocean's reflectivity model implemented in LBLMS could lead to spurious effects in the interpretation of the upwelling radiance, especially at low sun zenith angles. For this reason all measurements between the ER-2 take off time and 24:00 UTC (14:00 local time) were disregarded. The chosen transect is located between 01:26-01:34 UTC and is characterized by a solar zenith angle of 57 • .
Thermal infrared and visible MAS imagery show (see Fig. 8) that the cloud is fairly homogeneous only in the optically thickest part and quite variable elsewhere. The thinnest part, located at the edge of the very extended cloud layer, is rather inhomogeneous with many open gaps over the underlying ocean. Some small cumulus clouds can be spotted below the cirrus layer, although not directly beneath the flight line.
Being the flight of the NOAA G-4 much shorter than the ER-2's mission, the drop-sonde spatially closest to the data analyzed was launched at 23:27 UTC, hence 2 h before the chosen ER-2 sector, but from satellite imagery there seem to be no important cloud development in these two hours. In the humidity profile, the cloud layer position is characterized by a net increase of relative humidity between 200 hPa and 300 hPa up to a value of 60% in the middle of the layer for the section with the optically thickest cloud. As in the clear sky case, a sub-tropical climatological standard profile is used to fill the gap between the two aircrafts. As an estimate of the surface wind speed we use the same wind speed adopted in the clear sky scene (6.5 m/s), though the relative importance of the surface wind in absence of a sun glint is very small. The pigment concentration and the aerosol optical depth, profile and mixture are the same adopted in the clear sky case.
For layers with optical depth at 532 nm less than 3, CPL data is used to characterize the internal cloud structure and to define top and bottom cloud levels. Figure 9 shows the extinction cross section at 532 nm at nadir for the flight stretch chosen for the comparison. The cloud becomes thicker and the top height increases while flying from the edge of the cloud layer to the inner part.
The CPL optical depth at 532 nm is shown in Fig. 10 and it is seen that after 01:34 UTC the CPL signal is saturated, hence the cloud bottom information is not reliable. Three sectors are selected from the whole track: they are representative of a very thin (between red lines in Figure 10), a medium thin (green lines) and a moderately thick (blue lines) cirrus layer. The retrieval methodology RT-RET is used to determine cloud optical thickness (OT) and effective dimension D e from S-HIS data, collocated with MAS data, averaged over each of the three sectors selected. A-priori parameters (to be defined before the actual retrieval can be done) are: the µ parameter of the size distribution, the ice crystal shape and the bottom and top heights of the cloud layer. The retrieval, moreover, assumes that the IWC is distributed homogeneously with height.
In next section the retrievals and the forward computations for a standard case (called V0) defined with a specific choice of these parameters, are shown. Since no measured data is available from the field campaign about the micro-physical composition of the cirrus layer, we have used MIXML, for reasons already discussed. The size distribution is described by a Gamma with a value of µ = 0. Cloud top and cloud bottom are the average in each sector of the values determined from CPL measurements collocated with S-HIS FOVs. Other choices are possible: in regions where the CPL signal is not saturated, the extinction profile can be used to infer the IWC vertical distribution (assuming a vertically invariant constant composition and size distribution). A comparisons of results obtained with an IWC that varies with height versus the homogeneous case is shown in Sect. 4.3. In Sect. 4.4 results are presented that assume different values of the µ parameter of the PSD (cases V1 and V2). Table 1 summarises the characteristics of each cloud scene for the standard case (V0) and the other PSDs study cases discussed in Sect. 4.4. The effective diameter retrieved by RT-RET is in the range 64 to 80 µm, therefore the bulk of the mixture is a combination of the optical properties of irregular aggregate of hexagonal columns and 3-dimensional bullet rosette. The retrieved OT in the IR range is used to Table 1. Characteristics of the 3 cloudy sectors used for the comparisons. CLT and CLB are cloud top and bottom height obtained rom CPL; RT-RET OT (IR) and D e , hence IWC, are retrieved by RT-RET averaged in each sector; CPL OT is measured by CPL at 532 nm and CPL OT 1σ is one standard deviation of measured CPL OT in each sector.

Sector 1 (red)
Sector 2 (green) Sector 3 (blue)  V0  V1  V2  V0  V1  V2  V0  V1  V2 CLT (km) 11.6 11.6 11.6 11.9 11.9 11.9 12.1 12.1 12.1 CLB (km) 9.7 9.7 9.7 9.8 9.8 9.8 10.0 10.0 10.0 IWC (g/m 3 )*10 −3 8.9 7.3 8.9 11.7 10. determine the total ice mass for the whole depth of the cloud (IWP) and then used to compute the OT at 532 nm using the same PSD, type of mixture and optical properties database: such OT values agree within the natural variability of the OT values derived from CPL, also shown in Table 1 in the three sectors. Figure 11 shows the brightness temperature in MAS channels obtained by averaging MAS pixels over the S-HIS footprints in each sector, and the corresponding S-HIS data (averaged in each sector) convolved over the MAS spectral response function. The agreement between the instruments is quite good throughout most of the IR spectrum.

Results for the standard cloudy cases
LBLMS is used to simulate the upwelling radiance from the three sectors highlighted in Fig. 10 over the whole spectral range covered by MAS, using the same configuration as in the clear sky case. Figures 12, 13 and 14 show the comparison between MAS measurements and LBLMS simulations in the three cloud sectors. Figure 15 shows the relative difference between LBLMS simulations and MAS measurements.
Since the procedure RT-RET is based on an iterative least square fit of the S-HIS radiance between 820 and 980 cm −1 the agreement between the observations and the simulations in that range is obviously quite good. Agreement between data and simulations (accounting for one standard deviation of the MAS data) is in fact found between 1200 and 650 cm −1 , in all the sectors. The difference between LBLMS and MAS lies between −2 K and +2 K, with an overestimation for the thin cloud and a slight underestimation for the other two cases, but well inside the signal variability over the scene, represented by one standard deviation of the MAS data around the mean value. All the simulations show a narrow overestimation region at about 1050 cm −1 due to incomplete information about the O 3 atmospheric profile. For wavenumber greater than 1800 cm −1 the signal is the sum of the radiation emitted at terrestrial temperatures and of the  Fig. 11. MAS data averaged over S-HIS nadir FOV and S-HIS data convolved with MAS relative spectral response. The red solid line and yellow dashed refer to Sector 1 (the thin cloud layer); the two green shades are for Sector 2 (medium cloud layer) and the blue lines are for the moderately thick layer of Sector 3. The vertical bars are the one standard deviation of the MAS data around the average values and denote the natural variability as well as radiometric error. reflected solar radiation and the influence of the latter is evident in the spectral ranges where the brightness temperature is higher than that observed in the IR window. LBLMS underestimates MAS in the red and green sectors between 2000 and 2200 cm −1 , while overestimates in the blue sector between 2300 and 2700 cm −1 .
The SW upwelling radiance computed by LBLMS for the three cloudy sectors (Fig. 12, and top panel of Fig. 15) shows a strong underestimation with respect to MAS. The difference between LBLMS and MAS for the thick cloud is between 20% to 45% and the difference for the thin cloud is between 25% to 55% relative to MAS. A large underestimation is also seen in the NIR range ( Fig. 13 and middle panel of Fig. 15), in all sectors at all wavelengths, ranging from 50% to 60%.

Homogeneous versus inhomogeneous cases
The results for the standard (homogeneous) case are consistent with RT-RET that retrieves cloud parameters assuming that the cloud is vertically homogeneous. In this section the CPL derived cloud extinction profile is used to modulate the amount of ice mass within the cloud layer, in each sector. We have used the same D e constant with height, as in the standard case.
The cloud properties are then used to calculate the total ice amount so that the OT, computed at 532 nm, matches the CPL one (given in Table 1). The average CPL extinction profiles for the three chosen sectors are shown in Fig. 16. The extinction profiles are interpolated at the model levels (computations in sectors 1/2/3 are done with 11/14/14 cloud layers, of thickness 0.2 km, out of 81 layers to describe the atmosphere) and the IWC profile is computed so that it is proportional to the average extinction profile. The final IWC profiles are shown in Fig. 17.
It is well known, from first principles, that radiance in the infrared is very dependent on the vertical distribution of mass, since the latter changes the source function profile which, being mostly due to emission, is a function also of temperature. Our interest is however centered in the shortwave and near infrared. Fractional differences between the inhomogeneous and the homogeneous cases, for the three sectors, are shown for the shortwave and near infrared MAS channels in Fig. 18. Largest fractional differences are seen in sector 1 in the whole spectral range. The fractional differences increase with increasing wavelength in the shortwave and stay nearly constant in the near infrared, except between  . 17. Average IWC profile, for the LBLMS layering, for Sector 1 (red), 2 (green) and 3 (blue) derived as explained in the text. The vertical dashed lines indicate the profiles adopted for the homogeneous computations.
1.8 and 1.9 µm, where molecular absorption is relevant. For the case under study, the inhomogeneity reduces the scattered radiance to the sensor with respect to the standard case. These results show quite clearly that the difference in simulated radiance between the homogeneous and the inhomogeneous cases are, in any case, much smaller than the difference between measured and simulated radiances described in Sect. 4.2.

Results with different PSDs
Small particles are strong scatterers of solar radiation, though they bring small contribution to the total PSD's mass and one could imagine that one of the causes of the radiance underestimation is due to an improper treatment of the PSD in the small particle range. There is in fact much discussion on the role of small ice particles (whose effective diameter is smaller than about 50 to 60 µm) and on the physical mechanisms that produce and eventually maintain these particles inside cirrus clouds. Jensen et al. (2009) report on new aircraft measurements in anvil cirrus sampled during the Tropical Composition, Cloud, and Climate Coupling (TC4) campaign with the 2-Dimensional Stereo (2D-S) probe, which can detect particles as small as 10 µm. They suggest that micro-physical measurements in tropical cirrus clouds obtained with the CAS (Cloud Aerosol Spectrometer) should be considered suspect when large crystals are present, and that measurement made with instruments at the wing-hatch location are presumably also affected by shattering artifacts. They however point out that their findings are relevant for relatively low and warm tropical cirrus and do not imply that small crystals do not play a significant role in the radiative properties of other types of cirrus, such as anvils generated by continental convection, mid-latitude cirrus, or even anvil cirrus in other tropical regions. They also point out that small particles could persist in uppermost tropical tropopause which is often saturated with respect to ice.
The sensitivity of the results to the assumed set of PSDs is here investigated. The main parameter defining each set of PSDs is the dispersion µ. For each value of µ (i.e. each PSD set) the slope λ is varied so that 26 PSDs are defined by their effective dimension D e . Typical values of the dispersion range between −2 and 16, while typical slopes range from 0.1 to 10 −4 cm −1 (with D expressed in cm). These values are derived from data discussed in Baum et al. (2005) and Heymsfield et al. (2004), where measured PSDs of ice clouds from six field experiments are fitted to theoretical modified gamma type PSDs. The intercept N 0 determines the total number of particles in the PSD, constraining the IWC (and thus is a parameter derived from the retrieved OT and D e ).  Fig. 19. Comparison between radiance simulated for cases V0, V1, V2 for Sector 1. Top and middle panels: fractional radiance difference with respect to V0; lower panel: brightness temperature difference with respect to V0.
Two new sets of size distributions (each consisting of 26 PSDs for each set) were defined with the values µ = 7.0 (called V1) and µ = −0.1 (called V2), and for each set RT-RET is used to retrieve the corresponding OT and D e . Table 1 contains the retrieved cloud parameters for the three sectors and the two PSD sets (V1 and V2), besides the results for the standard case (V0).
The three retrieved PSDs are used to simulate the spectral radiance at high resolution in each of the three sectors, and the S-HIS and MAS convolved radiances. The difference among the results for cases V1 and V2 are plotted in Fig. 19, Fig. 20 and Fig. 21 respectively for Sectors 1 to 3. Fractional radiance difference (L(V x) − L(V 0))/L(V 0) is plotted for the SW and NIR ranges, while brightness temperature difference is plotted for the IR range. The maximum difference in the SW range is for case V1 in Sector 3 and is of about 0.02 (or 2%). In the NIR range maximum difference is around 2 microns and is less than 10%. In the IR range the absolute differences are below 0.3 K except at 2650 cm −1 in Sector 3 (less than 0.6).
Therefore, although the three sets of PSDs used are quite different, yet the difference in simulated radiances for the three PSDs are much smaller than the difference between measured and simulated radiances in the standard (V0) case.

A MODIS case study: a cross check for the LBLMS results in the SW and NIR
The MODerate resolution Imaging Spectroradiometer The scene over the Indian Ocean, observed by MODIS on Terra the 19 January 2009 (granule MOD021KM.A2009019.0320.005.2009019145552), is characterised by a large cirrus located SW of Australia, extending between latitude 40 • and 45 • (see Fig. 22).
The MODIS cloud mask ensures the ice phase of the cloud. Only pixels with OT much larger than 1 are considered, in order to minimise the influence of the surface properties and of the atmospheric profile below the cloud layer, given the difficulty to obtain an accurate description of their properties. Using MODIS cloud products three small areas are chosen, within the red circle, with different cloud OT; within each area several pixels with same retrieved OT and effective radius R e are averaged.
The relevant parameters of the three areas and the retrieved properties from MODIS, used for the LBLMS computations, are summarised in Table 2, where the effective diameter is defined as D e = 2R e and decreases as the cloud OT increases.  file, and it is understood to be comprehensive of all sources of errors in the retrieval of all parameters that are required to the determination of the OT: cloud top level, cloud phase and effective radius.
The database of ice crystals optical properties are SSP-SW and SSP-IR, same used for the simulations of the P-TOST case. The standard mid latitude summer profile (Anderson et al., 1986) is used as input to LBLMS, with a surface wind of 5 m/s and an ocean pigment concentration of 0.07 mg/m 3 . Since cloud base height is not a MODIS standard product, a cloud geometrical thickness of 3790 m has been assumed with a homogeneous IWC vertical distribution.
The test is done in 9 MODIS bands located in the SW and NIR, defined in Table 3. The results in the IR range would have been difficult to interpret due to the strong dependence of radiance on the accurate reconstruction of the temperature profile and of cloud top height and thickness.
LBLMS computations are done using same technical choices already discussed, in the spectral ranges that cover the MODIS Terra Relative Spectral Response (RSR) for the selected MODIS bands. These RSR are computed as the average (over all channels) of the individual L1B in-band RSRs, that are available from the MODIS Characterization Support Team (MCST) at the address: http://mcst.gsfc.nasa. gov/l1b/. Figures 23, 24 and 25 summarise the results of the comparisons for the three areas described in Table 2.
In these figures the upper panels are radiance differences and the lower panels are the fractional radiance differences that coincide with fractional reflectance differences. The blue bars denote the (2σ ) variability around the mean measured radiance, that corresponds to same retrieved OT ( and R e ), but slightly different Sun zenith angles and observation angles among the pixels whose radiance is averaged.
Band-2 located at 856.7 nm is currently used for the retrieval of the cloud OT above the ocean (King et al., 1997) and the simulations are very close to measurements at this wavelength. This result would be quite obvious were the Table 2. Relevant parameters and retrieved cloud properties for the three areas selected for the MODIS case studies (the name is defined in column one) over the cirrus cloud of Fig. 22. SUN z.a. is the solar zenith angle; Azi. and Zen. are the azimuth and zenith observation angles; CTP, OT and R e are the MODIS retrieved cloud top pressure, optical thickness and effective radius; OT 1σ is estimated error on OT. IWC is cloud ice water content and CD is assumed cloud layer depth. same procedure used for the forward and inverse computations. In the present case the result is not obvious and in fact it demonstrates that the inverse (MODIS processing) and forward (LBLMS computations) procedures (single-scattering databases and PSDs) are compatible.
As we move to shorter wavelengths, some degradation is observed which is always within the stated 2σ uncertainty level, except for Band-3 in Area Ci1, the case of least opacity, where LBLMS overestimates the measured radiance by 20%. In all the cases under study LBLMS overestimates the radiance measured in Bands-6 and 7, which are used to retrieve the effective radius. Although the relatively large fractional differences (lower panels) are associated to low radiance values, still they indicate that some of the assumptions on which the LBLMS simulations are based are different from the one used in MODIS retrievals; or that, perhaps, R e retrieved from band 7 is representative of the top-to-middle part of the cloud and some of the difference could be related to inhomogeneous vertical distribution of particle sizes. We have no firm explanation of this discrepancy and further work is necessary.

Discussion of results and conclusions
Measurements taken during the 2003 Pacific THORPEX Observing System Test by the MODIS Airborne Simulator (MAS), the Scanning High-resolution Interferometer Sounder (S-HIS) and the Cloud Physics Lidar (CPL) are compared to simulations with a line-by-line and multiple scattering modeling methodology (LBLMS). The extension of LBLMS to the shortwave (SW) and near infrared (NIR) and the treatment of the bi-directional reflectance properties of the marine surface are discussed. LBLMS should provide the best possible simulation (in a plane parallel geometry), albeit applicable only to case studies.
A number of scenes are evaluated: two clear scenes, one with nadir geometry and one cross-track encompassing sun glint, and three cloudy scenes, all with nadir geometry. CPL data is used to estimate the particulate optical depth at 532 nm for the clear and cloudy scenes and the cloud top and bottom heights. Cloud optical depth and D e are retrieved from S-HIS infrared window radiances, and compares well to CPL values. MAS data is simulated convolving high resolution radiances.
The paper discusses, in Sects. 3 and 4, the results of the comparisons for the two clear cases and for the three cloudy cases. The main (problematic) result is that the simulations in cloudy sky conditions, using cloud parameters retrieved from infrared radiances, systematically underestimate the measured radiance in the visible and near infrared by nearly 50%, while cloud optical depths retrieved from infrared data agree, to within natural variability, with those derived from CPL. The change in simulated radiance due to the vertical distribution of ice mass is investigated in Sect. 4.3, and the change in simulated radiance due to the use of diverse PSD is investigated in Sect. 4.4. The conclusion, in both cases, is that the difference in simulated radiances, for the various cases, are much smaller than the difference between measured and simulated radiances described in Sect. 4.2.
In order to understand the cause for the observed discrepancies, MODIS radiances measured from Terra are also compared to LBLMS simulations in cloudy conditions, using retrieved cloud optical depth and effective radius from MODIS. Three case studies are selected corresponding to cloud decks of various opacity and the attention is focused on some SW and NIR MODIS bands. The differences are generally within twice the estimated standard deviation except for band 3 for the least opaque cloud where the difference is +18% and for the near infrared bands (band 6 and 7) where we have the largest relative departures, in presence of relatively low radiance values.
The results presented thus provide evidence that 1) LBLMS simulations in clear conditions are close to (P-TOST) MAS measurements over the oceanic surface, under a number of diverse viewing geometries; 2) LBLMS simulations in cloudy conditions strongly underestimate MAS radiances in the SW and NIR when cloud parameters are derived from infrared retrievals or lidar measurements; 3) LBLMS simulations are in good agreement with MODIS short-wave measurements when cloud parameters derived from same MODIS data are used to define the cloud optical properties; finally 4) the OTs retrieved by RT-RET from S-HIS IR data in 2) agree to within natural variability with OT measured by CPL at 532 nm (see Table 1).
The main difference between the results 2) and 3), presented in Sects. 4.2 and 5, is that in the former the cloud properties used for the radiance simulation are retrieved from hyper spectral S-HIS measurements in the main IR window, while in the latter they are retrieved from MODIS short wave channels. The two retrieval types sense different properties to derive their products: the MODIS retrieval uses scattered radiation and a realistic description of the phase function is thus fundamental. On the other hand RT-RET uses emitted and scattered radiation to infer the extinction OT and D e and scattering is a small fraction of the total signal observed, although it must, and is, fully accounted for by RT-RET. Any retrieval algorithm derives (cloud) parameters that, when used in the same forward model used for the retrieval, reproduce (to within well defined, and usually small, errors) the measured radiances. Point 3) above shows that the adopted cloud modeling (crystal mixture, PSD, databases of optical properties, scattering coefficients and phase function) used with LBLMS is sufficiently close to the cloud modeling used in the MODIS cloud property retrieval. This closeness does not imply that the modeling is correct, since similar properties are used in both the retrieval and the forward model. It is therefore to be expected that best results in the P-TOST cloudy case study are to be obtained in the IR, and that the best results in the MODIS case study are at wavelengths close to the one used for the retrieval, i.e. at 856.7 nm. What is not expected are the large deviations outside the spectral range used for the retrievals.
The retrieval of optical thickness from lidar systems exploits shortwave (scattered) radiation. The CPL OT retrieval methodology is described in McGill et al. (2003) that deals with the use of the airborne CPL in a field campaign. It is stated that the S-ratio (defined as the total scattered and absorbed energy divided by the amount of backscattered energy) can be estimated directly from the lidar data without assumption, under certain favourable circumstances, and among them when the field of view contains an ice cloud above 5 km. Therefore in our experimental conditions, no assumptions on the phase function are required in the CPL retrieval procedure to derive the OT. In our study the good agreement between CPL OT and RT-RET OT (point 4) implies that the angularly integrated properties (coefficients of scattering and extinction) are indeed coherent going from infrared to the short wave. Several studies in fact have shown that infrared retrieved OT well agree with lidar OT (e.g. Maestri and Holz, 2009;Turner and Eloranta, 2008;DeSlover et al., 2003).
The radiance emerging from the cloud layer is proportional, in conditions of single scattering, to the product of direct solar radiance reaching the cloud layer, of the scattering coefficient and of the phase function. Therefore the results outlined above as points 2) and 3) would imply (in condition of single scattering) that the discrepancy is to be ascribed to the phase function adopted, because of the point 4) above. In the actual (multiple scattering) conditions the above conclusion remains a reasonable working hypothesis. Hence the strong radiance (hence reflectance) underestimation in the SW and NIR, obtained in the P-TOST cloud cases, seems likely caused by the structure of the phase function in the short-wave and, to a lesser extent, in the infrared domain. We have taken the utmost care to reconstruct the phase function as originally computed and distributed, therefore the problem is likely present in the original phase function as provided in SSP-SW and SSP-IR. We are aware that the P-TOST cloud case involves only a nadir measurement geometry, and that, in single scattering conditions, our tentative conclusion regarding the phase function would apply only to a limited range of scattering angles.
In Ham et al. (2009), quoted also in the introduction, the RT model DISORT is used to check the consistency of the MODIS cloud products over a broad spectral range, from the IR to the SW. When the MODIS cloud products are used as input parameters in DISORT, there is agreement for the SW channels, a strong overestimation of the simulated reflectance in the 1.38 µm channel and a strong underestimation of the brightness temperature in all IR channels. The discrepancy found in the IR channels, when the cloud optical properties are retrieved in the SW range, resembles the results from our P-TOST case study. In Ham et al. (2009) it is pointed out that a possible explanation (page 1603) could be an incorrect retrieval of the actual cloud top, leading to errors in the positioning in the profile and hence in the emission temperature, but the cloud optical properties could also play a not negligible role. In our P-TOST case study a measurement of cloud top height is available from CPL, while it was decided not to model MODIS IR channels, for reasons explained in the main text.
The causes of our problematic results could be multiple and the authors do not have the practical knowledge required to master the fine details of the computations that were required to generate the SSP-SW and SSP-IR databases. The use, in the two spectral ranges, of different methodologies to compute the volume and projected area, that define the geometry of an ice particle of same maximum dimension, could be one of the reasons for the observed discrepancies. A lot of effort is spent in improving the optical properties of ice crystals: the effect of the inclusion of air bubbles inside ice particles is being studied (Xie et al., 2009) as well as changes to surface texture (Yang et al., 2008b) and roughness (Yang et al., 2008a).
In Zhang et al. (2009) it is shown that ice clouds OT inferred from POLDER (POLarization and Directionality of the Earth's Reflectances) are substantially smaller than the one inferred from collocated MODIS data, and it is stated that this difference is due to the use of different ice particle scattering models and specifically to the different scattering phase functions. Our results are along the same line, but using a different dataset and simulation methodology. In the conclusion section of Zhang et al. (2009) it is suggested that a set of existing or newly developed ice particle models should be used as the common basis to derive climatologies from satellite measurements. Although this unification is an important step, it represents a necessary, but not a sufficient condition. Indeed cirrus clouds pose a fundamental problem because it is the difference in properties among the shortwave and the long-wave that is ultimately important, and not only for climate change research. Therefore the first and most important step is to generate a database of ice particle properties that describes consistently the cloud features that are observed from the SW to the IR.