Interactive comment on “ Quantifying and correcting the effect of vertical penetration assumptions on droplet concentration retrievals from passive satellite instruments

Abstract. Droplet concentration (N d ) retrievals from passive satellite retrievals of cloud optical depth ( τ ) and effective radius ( r e ) usually assume the model of an idealised cloud in which the liquid water content (LWC) increases linearly between cloud base and cloud top (i.e., at a fixed fraction of the adiabatic LWC) with a constant N d profile. Generally it is assumed that the retrieved r e value is that at the top of the cloud. In reality, barring r e retrieval biases due to cloud heterogeneity, etc., the retrieved r e is representative of that lower down in the cloud due to the vertical penetration of photons at the shortwave infra-red wavelengths used to retrieve r e . This inconsistency will cause an overestimate of N d (referred to here as the penetration depth bias ), which this paper quantifies. Here we estimate penetration depths in terms of optical depth below cloud top ( dτ ) for a range of idealised modelled adiabatic clouds using bispectral retrievals and plane-parallel radiative transfer. We find a tight relationship between dτ and τ and that a 1-D relationship approximates the modelled data well. Using this relationship we find that dτ values and hence N d biases are higher for the 2.1 μm channel r e retrieval ( r e2.1 ) compared to the 3.7 μm one ( r e3.7 ). The theoretical bias in the retrieved N d is likely to be very large for optically thin clouds, nominally approaching infinity for clouds whose τ is close to the penetration depth. The relative N d bias rapidly reduces as cloud thickness increases, although still remains above 20 % for τ 19.8 and τ 7.7 for r e2.1 and r e3.7 , respectively. The magnitude of the N d bias upon climatological N d data sets is estimated globally using one year of daily MODIS (MODerate Imaging Spectroradiometer) data. Screening criteria are applied that are consistent with those required to help ensure accurate N d retrievals. The results show that the SE Atlantic, SE Pacific (where the VOCALS field campaign took place) and Californian stratocumulus regions produce fairly large overestimates due to the penetration depth bias with mean biases of 35–38 % for r e2.1 and 17–20 % for r e3.7 . For the other stratocumulus regions examined the errors are smaller (25–30 % for r e2.1 and 11–14 % for r e3.7 ). Significant time variability in the percentage errors is also found with regional mean standard deviations of 20–40 % of the regional mean percentage error for r e2.1 and 40–60 % for r e3.7 . This shows that it is important to apply a daily correction to N d for the penetration depth error rather than a time-mean correction when examining daily data. We also examine the seasonal variation of the bias and find that the biases in the SE Atlantic, SE Pacific and Californian stratocumulus regions exhibit the most seasonality with the largest errors occurring in the December, January, February (DJF) season. We show that this effect can be corrected for by simply removing dτ from the observed τ and provide a function to allow the calculation of dτ from τ . However, in reality this error will be combined with a number of other errors that affect both the r e and τ , which are potentially larger and may compensate or enhance the bias due to vertical penetration depth.


Introduction
Clouds have a major impact on Earth's radiative balance (Hartmann et al., 1992) and small changes in their properties are predicted to have large radiative impacts (e.g., Latham et al., 2008).The amount of shortwave (SW) flux reflected by fully overcast warm (liquid water) clouds for a given sun and scattering angle, or the reflectance of a cloud, is primarily determined by the cloud optical depth (τ ), which in turn can often be characterized by the liquid water path (LWP, the vertical integral of liquid water content) and the cloud droplet number concentration (N d ).For a given cloud updraft, N d is determined by the number concentration and physico-chemical properties of aerosols.Thus, couching cloud reflectance in terms of N d links the cloud albedo to aerosol and microphysical effects via the Twomey (1974) effect making N d a very useful quantity to determine observationally.N d can also influence cloud macrophysical feedbacks via its control on rain formation (Albrecht, 1989;Stevens et al., 1998;Ackerman et al., 2004;Berner et al., 2013;Feingold et al., 2015) and stratocumulus cloud top entrainment (Ackerman et al., 2004;Bretherton et al., 2007;Hill et al., 2009).
Satellite observations of clouds and N d are immensely useful for studying clouds, cloud-aerosol interactions and for model evaluation since they afford large spatial-temporal coverage.A method to obtain N d from passive satellite observations (e.g., from MODIS (MODerate Imaging Spectroradiometer); Salomonson et al. (1998)) of τ and the cloud droplet effective radius (r e ) for stratiform liquid clouds has been previously demonstrated (Han et al., 1998;Brenguier et al., 2000;Nakajima et al., 2001;Szczodrak et al., 2001;Boers et al., 2006;Quaas et al., 2006;Bennartz, 2007;Grosvenor and Wood, 2014;Bennartz and Rausch, 2017) and is described further below.In cloudy environments, aerosol optical depths cannot be retrieved from satellite making cloud property observations such as N d and the cloud droplet effective radius (r e ) the only useful indicator of the influence of aerosol on clouds.An advantage to using N d rather than r e to study cloud-aerosol interactions is that r e is also determined by the cloud water content and thus is a function of cloud macrophysical properties.N d on the other hand is only weakly controlled by cloud macrophysics allowing some separation of microphysical and macrophysical effects.
However, retrievals of N d from space are still somewhat experimental and there is a lack of comprehensive validation of the retrievals and the assumptions required.There is a need to characterize and quantify the associated errors; in this paper we focus on doing this for one source of N d error using a one-year N d data set for stratocumulus clouds from MODIS.
2 The vertical penetration depth N d bias and the adiabatic N d retrieval model N d is retrieved from passive satellite retrievals of r e and τ using an adiabatic cloud model that is described below.However, as shown in Platnick (2000); Bennartz and Rausch (2017), for a retrieval free from other error sources (e.g.those due to cloud heterogeneity), the retrieved r e is representative of the r e value lower down in the cloud due to the vertical penetration of photons at the shortwave infra-red wavelengths used to retrieve r e .In contrast, the retrieved τ is comprised of contributions from the extinction coefficient β ext (h), where h represents height from cloud base, throughout the whole cloud profile :- (1) Here h=0 represents cloud base and h=H is cloud top.
β ext (h) is defined as :- where r is the droplet radius and n(r) is the droplet size number distribution within a cloud unit volume such that N d = ∞ 0 n(r)dr.
Q ext (r) represents the ratio between the extinction and the geometric cross sections of a given droplet and can be approximated by its asymptotic value of 2 (van de Hulst, 1957) since droplet radii are generally much larger than the wavelength of light concerned (typically 0.6 to 0.85 µm) such that the geometric optics limit is almost reached.
r e and liquid water content LW C at a given height are respectively defined as: and where ρ w is the density of liquid water.Combining Eqns. 3 and 4 and inserting into Eqn. 2 gives :- To determine the form of r e (h) in the above equation in terms of L(h) and N d (h) we can utilize the fact that the "k" value, which is a measure of the width of the droplet size distribution, has been shown to be approximately constant in stratocumulus clouds (Martin et al., 1994;Pawlowska et al., 2006;Painemal and Zuidema, 2011).In this study we adopt a value of k = 0.8, which is the most commonly used value in previous studies that perform N d retrievals (e.g., Bennartz, 2007;Painemal and Zuidema, 2011) in stratocumulus.r v is the volume radius, defined as :- where we have used Eqn. 4 to insert LWC and Eqn.6 to write r v as a function of k and r e .Now we utilize the assumptions that N d (h) is constant with height and that LWC(h) is a constant fraction, f ad , of adiabatic.The latter equates to :- where c w is the rate of increase of LWC with height (dLW C/dz, with units kgm −4 ) for a moist adiabatic ascent and is referred to as the "condensation rate" in Brenguier et al. (2000), or the "water content lapse rate" in Painemal and Zuidema (2011).See Ahmad et al. (2013) for a derivation.c w is a constant for a given temperature and pressure.Allowing these assumptions, using Eqn.7 to substitute for r e in Eqn. 5 and combining with Eqns. 1 and 8 we can write :- At this stage, H * is any arbitrary height above cloud base and τ * is thus the optical depth between the cloud base and that height.H * can be expressed as a function of r e (H * ), k, N d and some constants by using Eqns.7 and 8.Then, given r e (H * ) and τ * , N d can be calculated as follows : - Generally, when retrieving N d it is assumed that the r e obtained from satellite is representative of that from cloud top, i.e., r e (H * )=r e (H) (e.g.Bennartz, 2007;Painemal and Zuidema, 2011).This would then mean that τ * is the full cloud optical depth (τ ) as retrieved by the satellite and thus could be used in Eqn. 10 above to obtain N d .However, since the r e obtained by satellite is actually equal to r e (H * ) then τ * <τ and thus τ * should be used in Eqn. 10 instead of the retrieved τ ; the problem lies in the fact that τ * is unknown.However, in this paper we fit a simple function for τ * as a function of τ based on radiative transfer modelling of a variety of idealised clouds.Then we estimate the error introduced in N d retrievals for one year of MODIS data due to the assumptions of r e (H * )=r e (H) and τ * =τ , on the assumption that there are no other biases affecting the r e retrieval.
We label this bias the "vertical penetration bias".
3 Data and Methods

Calculation of dτ
In order to calculate we have performed r e retrievals on idealised clouds using a similar algorithm to that used for MODIS retrievals.We produced idealised clouds that span a large range of stratocumlus-like clouds as represented by combinations of N d and LWP.We chose 41 values between N d =10 and 1000 cm −3 that were equally spaced in log space and 91 values between LWP=20 and 200 gm −2 spaced equally in linear space.All of the possible combinations from this sampling were used to sample the 2-D (N d , LWP) phase space.For each combination, discretized adiabatic model profiles following the form of those described in Section 1 (i.e., with a vertically constant N d and LWC that increases linearly with height) were generated using c w =1.81×10 −6 kg m −4 , f ad =0.8 and using a vertical spacing of 1 m.The droplet size distributions at each height were represented by a modified gamma distribution with a k value of 0.72, i.e. representative of an effective variance of 0.1.1-D radiative transfer (RT) calculations, assuming plane-parallel clouds, were performed on these profiles using the DISORT (Discrete Ordinates Radiative Transfer Program; Stamnes et al., 1988) radiation code in order to simulate reflectances at wavelengths of 0.86, 2.1 and 3.7 µm, matching those measured by MODIS to retrieve τ and r e .The RT calculations were performed assuming a black surface, a clear atmosphere (i.e.gaseous absorption is neglected), using a solar zenith angle of 20 o and a nadir viewing angle.These reflectances were then used to retrieve τ and r e values using the Nakajima and King (1990) bi-spectral method, as operationally used by MODIS.To do so, a lookup table was built from reflectances similarly calculated for a range of clouds that were assumed to be plane-parallel in nature, as assumed for the operational MODIS retrievals; i.e. these clouds were uniform in the vertical and horizontal with infinite horizontal extent.Again, a black surface and a k value of 0.72 were assumed along with the same viewing geometry as for the RT calculations on the adiabatic clouds.A fixed depth of 1 km was assumed with cloud base at an altitude of 1 km and cloud top at 2 km, although the cloud depth has no major effect on the reflectances generated.dτ was then calculated by choosing the value of the model profile of τ , as measured from cloud top downwards, that corresponded to the value of the model profile of r e that matched the retrieved r e .
Figure 1a shows a 2-D histogram of dτ values as a function of τ for the 2.1 µm retrieval.It shows that when plotted in this way dτ forms a fairly tight relationship with τ so that for a given τ only a small range of dτ values are possible.This suggests that the relationship can be parameterized based upon a 1-D relationship fitted to this data with little loss of accuracy.The mean value of each τ bin is also plotted (after smoothing over τ windows of 0.2) and this is the relationship used in this paper.
dτ is seen to increase with τ with a gradient that decreases with τ .A 4th order polynomial curve can be fitted (using the least squares method) to this mean value relationship that takes the form:- The coefficients of this fit and for the 3.7 µm retrieval are given in Table 1 along with the maximum absolute τ error for the fit for the range shown.The curve (white line in Fig. 1a) fits the mean data well with a maximum absolute difference of 0.09, although there will be some error when using this relationship (or the mean value relationship) due to the spread in the dτ values seen in the underlying histogram.Figure 1b shows the same results for the 3.7 µm retrieval.Again a tight 1-D relationship is suggested with a dτ value that increases with τ .Here, though, the curve flattens off much more rapidly, so that by τ =7.5 there is little dependence of dτ on τ and dτ saturates at a mean value of ∼2.6.The fit estimate for the mean curve (Eqn.12 and Table 1) again matches the actual curve closely with a maximum absolute error in dτ of 0.14.

MODIS data
For the MODIS data we use one year (2008) of MODIS Aqua data and follow a similar methodology to that used in Grosvenor and Wood (2014) in order to create a data set akin to the MODIS L3 product (King et al., 1997;Oreopoulos, 2005).We processed MODIS collection 5.1 joint-L2 swaths into 1 o × 1 o grid boxes.Joint-L2 swaths are sub-sampled versions of the full L2 swaths (sampling every 5th 1 km pixel) that also contains fewer parameters.We process the data from L2 to L3 in order to allow the filtering out of data at high solar zenith angles and to provide r e retrievals from both the 2.1µm and 3.7µm MODIS channels, hereafter referred to to as r e2.1 and r e3.7 , respectively.
For this work we relax the screening methodology slightly from that used in Grosvenor and Wood ( 2014) since here we are interested in the effects of the vertical penetration N d bias upon a more general global data set.We applied the following restrictions to each 1 o × 1 o sample that goes into the daily average (since multiple overpasses per day are possible) in order to attempt to remove some artifacts that may cause biases: 1.At least 50 joint-L2 1 km resolution pixels from the MODIS swath that did not suffer from sunglint were required to have been sampled within each gridbox.
2. At least 80 % of the available (non-sunglint) pixels were required to be of liquid phase based upon the "primary cloud retrieval phase flag".Analysis was only performed on these pixels.A high cloud fraction helps to ensure that the clouds are not broken, since broken clouds are known to cause biases in retrieved optical properties due to photon scattering The black line is the mean dτ in each τ bin after smoothing over τ interval windows of 0.2.The white line is the fit to the mean curve using Eqn.12.
through the sides of clouds.Often retrievals of N d are restricted to high cloud fraction fields for this reason (B07; PZ11) and so we focus on such datapoints here.
3. Only the pixels remaining after (2) for which the "cloud mask status" indicated that the cloud mask could be determined, the "cloud mask cloudiness flag" was set to "confident cloudy", successful simultaneous retrievals of both τ and r e  (Cavalieri et al., 1996).
7. Only pixels with τ > 5 were considered for the N d data set due to larger uncertainties from instrument error and other sources of reflectance error for τ and r e retrievals at low τ (Zhang and Plantnick, 2011;Sourdeval et al., 2016).
Figure 2a shows the number of days from the year of data examined in this study (year 2008) that fulfilled the above criteria and thus are likely to produce a good N d retrieval.Regions with high numbers of days where useful N d retrievals can be made have been selected for closer examination in this study; they are listed in Table 2 along with information on the mean and maximum numbers of days of good data and some statistics on the N d biases that will be described later.The permanent marine stratocumulus decks are among those selected, namely those: in the SE Pacific off the western coast of S. America (Region #1); in the SE Atlantic off the western coast of southern Africa (Region #2); off the coast of California and the Baja Peninsula  and r e3.7 .This was done both by using the retrieved τ value instead of τ * in Eqn. 10 (i.e., assuming that re(H * )=r e (H) as is often assumed for N d retrievals) and by estimating τ * using the retrieved τ along with the dτ values that were calculated as described above.This therefore gives N d data sets for the "standard" method and a corrected method, allowing the differences between the two to be examined.

Results
Following Eqn 10, the ratio between the uncorrected and corrected N d values can be shown to be :- The equation shows that, for a constant dτ , the relative N d bias due to an uncorrected τ value would increase with decreasing τ as τ -dτ approaches zero.Figure 3 shows how the relative bias varies as a function of retrieved τ when using the dτ values At τ =5 the relative error is 50% for the r e2.1 retrieval and 30% for the r e3.7 retrieval.At higher τ the errors reduce rapidly, but remain above 10% for the r e2.1 retrieval over the τ range shown.For the r e3.7 retrieval the relative error drops below 10% for τ >∼ 15.Thus, the overall degree of error due to this effect will be determined by the distribution of τ for the regions of interest, which we take into consideration here using MODIS data for a representative N d data set.Figure 2b shows the time-mean τ for the data set as filtered by criteria 1-6 above; i.e. without the τ >5 restriction.Table 2 lists the regional means of these time-mean values along with the regional means of the standard deviations of τ over time.It shows that the mean τ values of the tropical and sub-tropical regions are generally lower than those at higher latitudes.The East China Sea, Barents Sea, NW Atlantic and Southern Ocean regions exhibit the highest mean τ values out of those examined and so should be expected to show the lowest N d biases due to the vertical penetration effect.The SE Atlantic region (and the 5 region to the west of Africa in general) show low τ and can be expected to give high N d biases.Table 2 also lists the fraction of days for which τ ≤ 10 (f τ ≤10 ).τ =10 is the value above which N d biases drop below 33% for the 2.1 µm retrieval and below for that channel.The values in the table indicate that even in the least affected region (Barents Sea) this will occur for 23% of the days.For the SE Atlantic and SE Pacific region the percentages rise to 72% and 55% of the days, respectively.Thus, the vertical penetration depth N d bias is prevalent in all regions for which N d data sets are likely to be used, and particularly so in the sub-tropical stratocumulus regions where N d retrievals have been widely used and studied.
The overall bias is now estimated using one year of actual MODIS data in order to obtain a realistic distribution of τ values.
However, it should be noted that the data set used is deliberately filtered in order to only retain datapoints that are likely to give useful N d data, namely low liquid clouds with extensive 1×1 o cloud fractions; i.e., predominately stratocumulus.This is done in order to assess biases for the types of clouds that N d data sets will typically be used to study.
Table 2. Regional statistics for the various marine stratocumulus regions shown in Fig. 2. Shown are the mean and maximum number of days that fulfill the screening criteria in order to be considered as useful Nd retrievals; the regional means and standard deviations (σ) of the time-averaged optical depths (τ ) for the screened data set; the regional mean of the fraction of days for which τ ≤ 10 (f τ ≤10 ), which is calculated using only data from gridpoints for which the number of days with Nd data was ≥15; regional means of the predicted time-mean percentage biases in Nd due to the vertical penetration depth error; and regional means of the relative (percentage) standard deviations (over time) of the percentage Nd biases (i.e., regional means of the values in Fig. 5).Bias results are shown for both the 2.1 µm and the 3.7 µm re retrievals.Fig. 4 shows a map of the mean percentage biases and Table 2 gives the regional means of the values in the map.Considering firstly the biases for the r e2.1 retrieval, the biases are highest in the tropics and sub-tropics.The regional mean bias is 38%

µm
for the SE Atlantic region (Region #2), which is the stratocumulus region that seems to suffer the most.2. The biases for the r e3.7 retrieval display the same spatial patterns as for r e2.1 , but are significantly lower; the mean value in the region with the maximum bias (SE Atlantic, Region #2) is 20% and that in the region with the lowest bias (Barents Sea, Region #8) is 11%.The left plot shows the results for the re2.1 retrieval and the right for the re3.7 one.
It is also useful to know how variable the biases are from day to day for a given point in space since this will determine how useful the application of a single offset bias correction might be for correcting N d biases for daily data.Figure 5 shows the time variability of the bias in the form of the relative standard deviations (over time) of the percentage N d biases.It reveals that the percentage bias in N d generally has a larger relative standard deviation at latitudes above around 40 o with values typically ranging up to around 30-50% (of the mean percentage N d bias) for the 2.1 µm retrieval.Relative variability is greater for the 3.7 µm retrieval, perhaps due to the much lower mean percentage errors.Some of the selected regions show more variability than others, in particular the Barents Sea and East China Sea regions.
Table 2 gives the regional means of the relative standard deviations revealing values that range from approximately 20 to 40% of the mean percentage biases for the 2.1 µm retrieval and 40-60% for the 3.7 µm one.This shows that the application   DJF for those seasons.The SON season also generally produces higher biases than MAM and JJA for those regions, particularly for SE Pacific.For the East China Sea region the biases are lower in SON and DJF seasons than in the other seasons.We note that there is little data in this region for JJA since there are few low-altitude clouds with large regional liquid cloud fraction there in this season.The other regions either do not show a large amount of seasonal variability, or N d data is only available for part of the year due to a lack of sunlight in the winter months.

Discussion
There are some caveats to the results that we presented here that we now discuss.We have shown that, theoretically, the effect of retrieving a lower r e than the cloud top r e that is assumed in the N d retrieval can be corrected for by simply removing dτ from the observed τ .However, this rests upon dτ having the dependence upon τ shown in Fig. 3 across all of the cloud types relevant for the N d data set.This relationship was based on the retrieved r e for a range of clouds, although only for a nadir viewing angle and a solar zenith angle of 20 o .Platnick (2000) showed that dτ has some dependence on viewing geometry and so the consideration of a wider range of view and solar zenith angles should ideally be made.The modelling of the idealised clouds and the correction rests on the assumption that r e increases monotonically with height within the cloud (following the adiabatic assumption), but there is some suggestion that the development of precipitation-sized droplets might lead to larger droplets being preferentially found below cloud top (Chang and Li, 2002;Nakajima et al., 2010a, b;Suzuki et al., 2010) .However, Zhang et al. (2012) found that MODIS retrievals of r e performed on model generated clouds were not significantly affected by the presence of precipitation.Also, during the VOCALS field campaign in the SE Pacific region, aircraft observations showed that r e generally did increase with height up to cloud top (Painemal and Zuidema, 2011), indicating that this is not a problem at least for the near-coastal clouds tested.Further offshore the likelihood of precipitation increases as clouds become more cumulus-like and so for those clouds the issue may be greater and hence more caution should be exercised when interpreting the results presented here for such regions.
It is also clear that the suggested correction for the vertical penetration effect should only be applied to the retrievals of sub-pixel heterogeneity (Zhang and Plantnick, 2011;Zhang et al., 2012Zhang et al., , 2016)); 3D radiative effects (Marshak et al., 2006); assumptions regarding the degree of cloud adiabaticity (f ad in Eqn.10; Janssen et al., 2011;Merk et al., 2016); the choice of k value (assumed constant; Brenguier et al., 2011;Merk et al., 2016); the assumption of a vertically uniform N d ; the assumed droplet size distribution shape and width (Zhang, 2013); viewing geometry effects (Várnai and Davies, 1999;Horváth, 2004;Varnai and Marshak, 2007;Kato and Marshak, 2009;Liang et al., 2009;Di Girolamo et al., 2010;Maddux et al., 2010;Liang and Girolamo, 2013;Grosvenor and Wood, 2014;Liang et al., 2015;Bennartz and Rausch, 2017); upper level cloud and aerosol layers (Haywood et al., 2004;Bennartz and Harshvardhan, 2007;Davis et al., 2009;Meyer et al., 2013;Adebiyi et al., 2015;Sourdeval et al., 2013Sourdeval et al., , 2016)), etc.These errors have the potential to bias N d in a way that opposes the positive bias expected from the vertical penetration effect such that the overall biases may cancel out.Indeed, the largest source of error in N d is likely that from r e biases given the sensitivity of N d to r e in Eqn.10.MODIS r e has generally been shown to be biased positively compared to aircraft observations (Painemal and Zuidema, 2011;King et al., 2012), which would lead to a negative N d error when taken alone.Thus, the application of the correction described in this paper in isolation has the potential to enhance any negative bias in N d caused by a positive r e bias.Painemal and Zuidema (2011) actually demonstrated that MODIS N d agreed rather well with N d from aircraft for the SE Pacific region despite the fairly large r e bias; this was thought to be due to the fortuitous cancellation of (for N d ) the r e bias with biases in the k parameter and f ad .However, the agreement between aircraft and MODIS N d seen in Painemal and Zuidema (2011) would deteriorate if a correction for the N d bias due to the penetration depth effect discussed here was also applied.
Table 2 indicates that the result would be a MODIS N d underestimate of around 35 % (average for SE Pacific, region#1) for the 2.1 µm retrieval, assuming perfect initial agreement.This indicates that another N d bias may have been operating in order to give the good observed agreement.

Conclusions
We have described and quantified a positive bias in satellite retrievals of cloud droplet concentration (N d ) that make use of the adiabatic cloud assumption to estimate N d from satellite observed cloud optical depth (τ ), effective radius (r e ) and cloud top temperature.We term the N d bias the "vertical penetration bias".The bias is specific to the methodology of the N d retrieval, as opposed to being a bias in the underlying observations.It arises due to the well-documented vertical penetration of photons with wavelengths in the shortwave-infrared range into the upper regions of clouds, so that r e retrievals are representative of values some distance below cloud top (Platnick, 2000;Bennartz and Rausch, 2017) rather than being those at cloud top as assumed by the N d retrieval.Here we quantified the optical depth as measured from cloud top downwards, dτ , at which the retrieved r e equaled the actual r e for adiabatic clouds covering a large range of total cloud optical depths and N d values.We showed that knowledge of dτ allows a corrected N d to be calculated by subtracting dτ from the observed τ and using that in the N d retrieval instead of τ .We characterised dτ as functions of τ for the 2.1 and 3.7 µm r e retrievals (r e2.1 and r e1.6 , respectively) and found that a 1-D relationship approximates the modelled data well.dτ increases with τ and is larger for r e2.1 than for r e3.7 and so the vertical penetration N d bias affects retrievals based on r e2.1 more than those using r e3.7 .
We quantified the vertical penetration N d bias for a one-year N d data set.The new N d formulation presented here suggests that N d errors will increase as the τ value of the cloud scene gets lower.For many regions that are considered trustworthy for N d retrievals (typically stratocumulus regions), there are high frequencies of low τ values and so the N d biases are significant.
For example, for the SE Pacific and SE Atlantic regions clouds with τ ≤ 10 (for which N d errors are expected to be ≥33% for r e2.1 and ≥15% for r e3.7 ) occur, respectively, 56 and 72% of the time on average.The mean r e2.1 vertical penetration N d biases for these regions were 35 and 38%, respectively.Out of the stratocumulus regions examined, these two were the worst affected.
For r e3.7 the N d biases were much smaller; for example, mean biases for the SE Pacific and SE Atlantic regions were 17 and 20%, respectively.N d biases were predicted to be worse for the tropical and sub-tropical regions than for higher latitudes.The time-variability of the biases were also examined and were shown to be significant (regional mean standard deviations of 20-40% and 40-60% for r e2.1 and r e3.7 , respectively).This indicates that long term averages of the vertical penetration N d bias corrections are not useful for correcting N d data over short timescales (e.g.daily N d data).We also examined the seasonality of the N d biases and showed that, for the stratocumulus regions, generally the DJF season was worst affected, followed by SON.
We caution that the correction for the vertical penetration N d bias presented here should only be considered in combination with corrections for other biases that affect N d since otherwise N d biases could be made worse, for example in situations where the fortuitous cancellation of opposing errors leads to initially small N d errors.The latter was suspected to have occurred for the comparison between MODIS N d retrievals and in-situ aircraft observations as presented in Painemal and Zuidema (2011).
We showed that the SE Pacific, which is the region examined in that study, had a mean vertical penetration depth error of 35% suggesting that another unidentified N d bias may have been operating in order to give good agreement.
Previous studies have shown that r e3.7 is less prone to biases due to sub-pixel averaging (Zhang and Plantnick, 2011;Zhang et al., 2012Zhang et al., , 2016)).Thus, combined with the work presented here, this supports the conclusion that r e3.7 likely represents a better choice for use in N d retrievals.
For future work, it is recommended that additional characterization of dτ is performed for a range of viewing geometries in order to ensure that the results presented here are robust for all cloud retrievals.Further investigation into how the presence of precipitation affects our assumptions and results is also warranted.

Data availability
The data set is built from publically available MODIS Level-2 data (Collection 5).

Figure 1 .
Figure 1.2D histogram of dτ as a function of τ for a range of clouds (see text) for the 2.1 µm re retrieval (a) and the 3.7 µm retrieval (b).
Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License.for the 2.1µm channel were performed and the cloud water path confidence from the MODIS L2 quality flags was designated as "very good confidence" (the highest level possible) were used.This is a little different from the official MODIS L3 product where a set of cloud products are provided that are weighted using the quality assurance flags.Rather than weighting our L3-like product with the QA flags we have simply restricted our analysis to pixels with the highest confidence for water path.4.The mean 1 o × 1 o cloud top height (CTH) is restricted to values lower than 3.2 km.This is done both to avoid deeper clouds for which N d retrievals are likely to be problematic due to the increased likelihood of a breakdown of the assumptions required to estimate N d , such as a constant fraction of the adiabatic value for LWC and vertically constant N d , as well increased retrieval issues due to cloud heterogeneity.CTH is calculated from the MODIS 1 o × 1 o mean cloud top temperature (CTT) and the sea surface temperature (SST) using the method ofZuidema et al. (2009).SST data was obtained from the v2 of the NOAA Optimum Interpolation (OI) Sea Surface Temperature data set (NOAA_OI_SST_V2) that provides weekly SST data at 1 o × 1 o resolution.This was interpolated to daily data on the assumption that SST does not vary significantly over sub-weekly timescales. 5.The mean 1 o × 1 o solar zenith angle (SZA) was restricted to ≤ 65 o following the identification of biases in the retrieved τ , r e and N d at high SZAs (Grosvenor and Wood, 2014).6. 1 o × 1 o grid-boxes were rejected if the maximum sea-ice areal coverage over a moving two week window exceeded 0.001 %.The sea-ice data used was the daily 1 o × 1 o version of the "Sea Ice Concentrations from Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Data, Version 1" data set

(
Region #3); in the Bering Sea off the SW coast of Alaska (Region #6); and in the Barents Sea to the north of Scandinavia (Region #8).These regions are where the highest numbers of selected days occur with values ranging up to a maximum of 141 days (for the Bering Sea region).The Barents Sea region has the lowest maximum number of days out of this group, reflecting the fact that N d retrievals cannot be made during a lot of the winter season in this region due to a lack of sunlight.The Southern Ocean (Region #5) and the NW Atlantic (Region #7) regions frequently produce stratocumulus, although it is often associated with the cold sectors of cyclones and so its location from day to day is more transient.These regions are also affected by high solar zenith angles in the winter seasons, which also restricts the number of retrievals possible there.The East China Sea region Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License.(Region #4) produces the lowest mean and maximum numbers of days since the stratocumulus areas are mostly restricted to near the coast and occur mostly in the winter season.

Figure 2 .
Figure 2. Left: Number of days in 2008 that fulfilled the criteria required to be counted as a valid Nd retrieval.See the text for details on the criteria.Various regions of interest are also denoted by the boxes and numbers.Right: Mean optical depth for data set with only filtering criteria 1-6 applied (see text); i.e., the τ > 5 restriction was not applied.Following this screening, the 1 o × 1 o gridboxes associated with each MODIS Aqua overpass were averaged into daily mean values.N d was then calculated using Eqn. 10 from the 1 o × 1 o daily mean τ , r e and CTT.N d was calculated for both r e2.1 10 based on those from Fig. 1 (mean curve, black line).9 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 3 .
Figure 3.The ratios of Nd values from the standard MODIS calculation (using τ * =τ , see Eqn. 10) to those from the corrected calculation (using τ * =τ -dτ ; the dτ values used are those shown by the black line in Fig. 1) vs retrieved τ .
The biases are a little lower for the other major stratocumulus regions; e.g. for the SE Pacific region (Region #1) and the Californian region (Region #3) the mean biases are 35%, although the biases increase further west where the dominant cloud regime tends to shift towards trade cumulus clouds.The remaining regions all have mean biases of 25-30%.The Barents Sea region (Region #8) has a value Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License. of only 25%, representing the stratocumulus region with lowest mean bias.These results indicate higher τ values for the clouds in the East China Sea, Southern Ocean, Bering Sea, NW Atlantic and Barents Sea regions relative to the Californian and S.E.Pacific stratocumulus regions, with the SE Atlantic region exhibiting the lowest τ values.This is confirmed by the mean τ values shown in Table

Figure 4 .
Figure 4. Maps of the annual mean percentage error for uncorrected Nd retrievals using a year (2008) of daily MODIS data that has been filtered to select data points in which Nd retrievals are favourable and therefore most likely to be used for Nd data sets (see text for details).
of a single annual mean offset bias correction is likely to lead to fairly large biases for the N d estimates for individual days for regions where the mean N d errors are significant.If daily data is used to determine relationships between cloud properties and 12 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License.N d without correcting for the biases examined here then significant variability in N d might be introduced that may affect those relationships via non-linear effects.

Figure 5 .
Figure 5.As for Fig. 4 except showing the relative (as a percentage) standard deviation of the percentage Nd bias over time.

Fig. 6
Fig. 6 shows how the percentage N d biases change with season for the r e2.1 retrieval only.Interestingly, the highest biases tend to occur in the DJF season for the SE Pacific and SE Atlantic stratocumulus regions, indicating that τ values are lower in

Figure 6 .
Figure 6.Seasonal mean percentage Nd biases for the re2.1 retrieval only.
N d with consideration of other bias sources.These other potential error sources are numerous and include r e biases due to 14 Atmos.Meas.Tech.Discuss., https://doi.org/10.5194/amt-2017-477Manuscript under review for journal Atmos.Meas.Tech.Discussion started: 4 January 2018 c Author(s) 2018.CC BY 4.0 License.

Table 1 .
Coefficients for the fitting curve (Eqn.12) to estimate the mean dτ value as a function of τ .The maximum absolute error from the fit is also shown.