The dataset used to perform the present overlap analysis is a 4-year (2007–2010) version of the 2D cloud optical depth (COD) field described in O22. COD is resolved vertically (240 m) and along the direction of the satellite track (∼1.1 km), and comes from the 2B-CWC-RVOD CloudSat dataset for liquid phase clouds and from the CALIPSO-enhanced 2C-ICE product for ice phase clouds. Because the first product is constrained by Aqua-Moderate Resolution Imaging Spectroradiometer (MODIS) vertically-integrated (total) optical depth (TAU), which is available only during daytime, the strict daytime availability applies also for the merged COD dataset. Comparisons with the 2B-CLDCLASS-LIDAR product showed missing retrievals for about 20% of liquid phase clouds. Those missing retrievals correspond mainly to thin liquid clouds detected by CALIPSO only in the 2B-CLDCLASS-LIDAR product and missed in the 2B-CWC-RVOD product. O22 describes how these values were filled, namely using coincident MODIS TAUs, and CODs of neighboring cells. We will examine below whether the specific COD filling procedure employed in O22 affects overlap findings.

For cloud occurrence overlap we use the generalized overlap framework first introduced by Hogan and Illingworth (2000), hereafter HI2000. According to this framework, the combined vertically-projected cloud fraction of two layers with cloud fractions C ₁ and C ₂ (expressed in this work either as fractional coverage from 0 to 1 or as percentage between 0 and 100), which are part of a contiguous cloud entity and are separated by a vertical distance Δz, is the linear combination of the combined cloud fractions of the two layers corresponding to the maximum (C ^max ) and random overlap (C ^ran ) assumptions. The relative contribution of C ^max and C ^ran to the combined generalized cloud fraction is regulated by a weighting parameter α: C g e n = α C max + 1 − α C ran (1) where C max = max C 1 , C 2 and C ran = C 1 + C 2 − C 1 C 2 . Because C ^max < C ^ran alpha is positive when C ^gen < C ^ran , with α = 1 for the special case of maximum overlap, and α = 0 for the special case of random overlap. Negative α which occurs when C ^gen > C ^ran indicates that cloud occurrence overlap is less than random, i.e., tending towards minimum overlap. HI2000 found that the overlap of cloud layers that have intervening clear layers is very nearly random, i.e., α ≈ 0 even at small separation distances Δz, a result we also found when analyzing the current dataset (discussed later). For contiguous cloud layers, HI2000 parameterized α as an exponentially decaying function of separation distance with a decorrelation length L _α : α = exp − Δ z / L α (2)

Because this function does not produce negative values, such a parameterization implicitly assumes that cloud overlap cannot be less than random. Maximum overlap corresponds to L _α = ∞, while random overlap to L _α = 0.

For overlap of cloud condensate (cloud optical depth) horizontal distributions, we are interested not as much in how absolute values of members of the distribution are correlated, but rather in how relative values are. Namely, for two cloud layers that may have very different grid-mean optical depth or water contents, we are interested in how the thinner, intermediate, and thicker parts of the layers align vertically, in other words in how their subgrid optical depths overlap in a relative sense. This can be captured by Spearman’s rank correlation coefficient ρ between the two distributions. While α is theoretically unbounded on the negative side, and has maximum value of unity, ρ values remain within ±1. This means that parameterizing its variation as an exponentially decaying function of separation distance with a second decorrelation length L _ρ (distance to the 1/e value), so that ρ = exp − Δ z / L ρ (3) captures only positive rank correlations. When ρ = 1, then the ranks of the COD values that are vertically aligned are correlated maximally which corresponds to L _ρ = ∞, while when ρ = 0 they are correlated randomly which corresponds to L _ρ = 0.

In this work, we calculate the variation of α and ρ as a function of cloud layer separation distance for scenes consisting of 100 consecutive “rays” (subcolumns). The scenes have an approximate length of 110 km (while the along track size of the rays of the merged CloudSat-CALIPSO dataset is about 1.7 km, sampling is performed at 1.1 km scales because of the spatial overlap between adjacent rays). Our overlap results in this work are therefore specific to this scale. Tompkins and Di Giuseppe (2015) conducted an analysis of the dependence of cloud occurrence overlap on scene size (“scale length,”) but such dependence is not investigated here. Our chosen scene size is as in O22, mirroring typical grid sizes of numerical models used for climate simulations.

We elected to calculate these two parameters for layers with 0.05 < C < 0.95 and ignore the cloud state of the intervening layers (i.e., whether some or all the layers are clear or cloudy). Layers with very low values of cloud fraction are not used because they contribute little to the combined cloud fraction and make C ^max ≈ C ^ran ; once a layer becomes nearly overcast on the other hand, its complete overlap with another layer is almost assured, rendering quantification of degree of overlap meaningless. The value of the upper C threshold used in cloud overlap calculations was a major focal point in Tompkins and Di Giuseppe (2015) who showed that it is intricately related to scene size. While these authors recommend using a rather low value for the upper threshold to suppress the dependence on scene size, this should not necessarily be in our view the predominant criterion. First, the objective in a climate model application is to parameterize overlap at scale lengths of interest and not at arbitrarily small or large scales. Second, when we examined several randomly chosen individual scenes, we discovered that setting a low C threshold, such as 0.5 suggested by Tompkins and Di Giuseppe (2015), completely distorted their character. Third, since we also wish to calculate COD distribution overlap via rank correlations, preferably from the same data sample, it does not seem prudent to remove the layers with the largest numbers of points of non-zero COD, and therefore the presumably most robust estimates of rank correlation.

Finally, our calculations of α and ρ include all cloudy layers regardless of whether clouds occur in-between. We choose to do so even though differences in overlap between contiguous and non-contiguous clouds are conspicuous (Figure 1). We ignore the distinction between contiguous and non-contiguous clouds in our analysis for several reasons: First, the decorrelation lengths we derive are meant for use in the “Raisanen” generator (Räisänen et al., 2004) which produces subcolumns using a flavor of generalized overlap that ignores the distinction between contiguous and non-contiguous clouds. Second, non-contiguous clouds are quite rare at small separation distances (Figure 1), so including them at these small distances is not impactful. On the flip side, at large separation distances cloud occurrence overlap approaches random behavior anyway regardless of whether the layers are part of contiguous or non-contiguous cloud entities. Third, rank correlation coefficients may be biased in unpredictable ways when making such a contiguous/non-contiguous discrimination, which again, is ignored in the way the Raisanen generator employs rank correlations.

Our analysis process is as follows. We first calculate α and ρ profiles for individual scenes. For two arbitrary cloud layers at levels n and m that are apart by a separation distance Δz _n,m , the scene’s α for these two layers is calculated from: α n , m = C n , m − C n , m ran C n , m m a x − C n , m ran (4) where C _n,m is the true combined cloud fraction of layers n, m, calculated from the binary cloud fraction values (0 or 1) of the portion of the layer contained in subcolumn k. C n , m = ∑ k = 1 K c k , n , m K (5) K is the number of subcolumns in the scene (K = 100, in our case) and c _k,n,m = 1 if either or both of n, m layers of the subcolumn are cloudy and c _k,n,m = 0 when both are clear.

Spearman’s rank correlation coefficient ρ is calculated from: ρ m , n = 1 − 6 ∑ l = 1 L r n , l − r m , l 2 L L 2 − 1 (6) where L is the number of subcolumns where both layers n and m have cloud, and r _n,l and r _m,l are the COD ranks of layers n and m for the l th subcolumn.

The separation distance Δz is resolved at 240 m, thus all separation distances are multiples of 240 m. Values for the same separation distance are averaged regardless of the vertical location of levels n and m, and we therefore drop henceforth layer subscripts from α, ρ, and Δz.

Once α “profiles” (variations vs. Δz) are calculated at the scene level, they are averaged into multi-annual seasonal profiles at 4° latitudinal resolution by appropriate averaging the scene data by latitude φ . Figure 2 compares the zonal profiles of α from the O22 COD fields and from the 2B-CLDCLASS-LIDAR product. We show separately results for two seasons, DJF and JJA. Since our filling procedure by design increases the consistency of the cloud mask (occurrence) implied by the COD field with that in 2B-CLDCLASS-LIDAR (O22), the agreement is not surprising (seen also in Figure 3, discussed later). Still, one should be reminded that the 2B-CLDCLASS-LIDAR dataset does not have a true vertical resolution of 240 m since it resolves clouds as distinct “objects” based on whether horizontally contiguous sequences of rays and vertically contiguous layers belong to the same cloud type. For a consistent calculation of α across datasets we therefore replicate the 2B-CLDCLASS-LIDAR mask to create a resampled 240 m vertical resolution. Despite this oversampling, the α( φ , Δz) distributions of the two datasets look very similar.

The prominent zonal dependence of α and the expected decrease with separation distance are immediately apparent in Figure 2. At the highest separation distance α approaches zero (random overlap) and can even become negative (some degree of minimum overlap). One should keep in mind that sampling becomes progressively poorer the greater the separation distance (cf. Figure 1), so results become noisier, but the tendency towards zero or negative values is clear and consisted with all previous overlap studies. The geographical dependence of α, with the higher values at low latitudes progressively decreasing at high latitudes is also consistent with previous works (Shonk et al., 2010; O12) where it was interpreted as indicative of the greater vertical alignment in cloud types encountered most frequently at low latitudes. Tompkins and Di Giuseppe (2015) casted however doubts on this interpretation arguing instead that the observed behavior rather reflected latitudinal changes in cloud scales (also related to latitudinally-varying frequencies of cloud types) relative to the fixed size of the scene for which α is calculated (∼110 km in our case).

A seasonal dependence of α is also apparent when comparing the DJF (December-January-February) and JJA (June-July-August) panels. The most distinguishing feature is the movement of high α values northward and the expansion of the latitude-height phase space area occupied by positive values of α in the summer hemisphere. Basically, at a given separation distance, there is greater chance to find a higher value of α in the summer hemisphere than the winter hemisphere. This reflects a tendency of greater likelihood of occurrence of vertically developed (and apparently more vertically aligned) clouds in the summer hemisphere than the winter hemisphere. In the tropics, the movement of the Intertropical Convergence Zone (ITCZ) is similarly reflected by higher values of α north of the equator during JJA.

The seasonal variation of α becomes more apparent when expressed in terms of decorrelation length. From the zonal profile of α, the zonal variation of its decorrelation length L _α can be calculated by applying a linear regression fit to lnα vs.–Δz data for which 1/L _α is simply the slope (Figure 3). This type of fit has been previously employed by Mace and Benson-Troth (2002) and Naud et al. (2008), among others, and results in a different L _α than the “effective” decorrelation length of Barker (2008b) and Jing et al. (2016) which is only meaningful at the scene level. Appropriate sampling weights are used in the regression fit to account for the fact that fewer values are available for larger Δz’s. Again, we show separate curves of L _α for DJF and JJA from the two datasets, the filled COD field (left panel) and 2B-CLDCLASS-LIDAR (right panel). To obtain the results shown in Figure 3, we actually calculated an L _α zonal curve separately for each year’s season (i.e., instead of calculating the seasonal average across 4-years) and then averaged the four year-specific L _α curves of each season. We opted to do this in order to suppress the over-smoothing of α profiles caused by extensive averaging (cf. Figure 1). The extreme flipside would have been to calculate L _α from very noisy α profiles of individual scenes (of which there are about four million) and averaging across the very large L _α population. The two methods, namely zonal L _α from fits to zonally-averaged α profiles and zonal L _α from averaging individual scene L _α ’s yield different results. We do not opt for the second method because of the poor regression fits of individual scenes and the distorting effects a few extreme values can have. But even in our method which first conducts extensive averaging to α before calculating L _α , the quality of the regression has some dependence on latitude.

Zonal and seasonal contrasts in cloud occurrence overlap become the most apparent when overlap is expressed in terms of decorrelation length (Figure 3, solid curves). The peak values in L _α reach ∼5 km at the northern (JJA) and southern (DJF) edge of the tropics and drop quite rapidly towards mid and high latitudes. The contrast in L _α between the two seasons is much more pronounced in the NH than the SH. The two seasonal L _α zonal curves, the one coming from O22’s COD field (left panel of Figure 3), and the one from the standard 2B-CLDCLASS-LIDAR product (right panel of Figure 3) are very similar.

The magnitude of peak values of L _α is generally higher than in previous studies, although a direct comparison can only be conducted with O12. The reasons for the discrepancy with O12 are not entirely clear, but the different nature of the underlying dataset likely plays a significant role: α and L _α values in O12 were derived from a dataset based solely on CloudSat measurements, specifically a cloud mask inferred from 2B-GEOPROF reflectivities deemed to come from cells identified as cloudy according to predetermined thresholds. Here we use a daytime COD field, from combined CloudSat and CALIPSO measurements, and with retrievals available only when the algorithm converges to a solution. It is doubtful that the daytime aspect of the COD field explains the rather substantial discrepancy of decorrelation length values. Missing COD values being filled, and the specific method used to achieve this in O22 does not seem to be a factor either (left panel of Figure 4). Actually, had we not filled missing COD values we would have obtained even higher peak L _α values: because the filled values correspond to low (liquid) clouds, the chances of creating occurrences of overlap with high clouds increases; such overlaps tend towards random and correspond to lower L _α .

For a practical use of the observed L _α values, the latitudinal and seasonal variability should ideally be captured by a parameterization. We follow here on the footsteps of O12 who used time-varying Gaussian functions whose parameters vary with Julian date (J), namely: L = μ 1 + μ 2 exp − φ − μ 3 2 2 μ 4 2 (7) (we drop the subscript from L because this parameterization is also applied for COD overlap) where μ ₁, μ ₂, μ ₄ are constants, while μ ₃ varies with day of the year according to the equations below controlling the latitude at which the decorrelation length peaks: μ 3 = − 4 μ 3,0 J − 272 / 365   w h e n   J > 181 μ 3 = 4 μ 3,0 J − 91 / 365   w h e n   J ≤ 181 (8)

Such a parameterization can be easily incorporated in a climate model that has a cloud subcolumn generator (Räisänen et al., 2004), as shown in O12. Figure 3 (dotted curves) shows the Gaussian curves for January 1st and July 1st and includes the parameters of the Gaussian fit. The present parameterization has the exact same functional form as the one introduced by O12 and encapsulates the migration of the most vertically aligned (and presumably developed) clouds (those with the highest value of L _α ) northward (southward) during boreal summer (winter). The new parameters produce higher values of L _α (more maximum overlap) at low latitudes and smaller values (more random overlap) at high latitudes, compared to O12. In Section 4 we present an implementation of this parameterization to assess scene and grand-average simulated cloud properties and average cloud radiative effects and to evaluate the impact of transitioning from O12’s to the current Gaussian fits.

The “profiles” (Δz dependences) of rank correlation coefficient ρ of COD distributions are calculated and averaged similarly to profiles of α. Figure 5 provides a visualization of the zonal dependence of ρ profiles for DJF and JJA as in Figure 2. We see a similar pattern for ρ( φ , Δz) as for α( φ , Δz), namely decrease with height and latitude, and shift toward higher values in the summer hemisphere. Negative values below −0.1 are virtually non-existent for ρ values that have been averaged extensively over multiple scenes; this suggests that the inability of the inverse exponential to yield negative values is likely inconsequential in estimates of L _ρ . Values greater than 0.7 are quite rare and limited to Δz ≤ 480 m.

Estimates of L _ρ are obtained by regressing ln(ρ) against -Δz, similarly to how we calculated L _α ; namely by performing again a regression on each year’s seasonal zonal profiles, and averaging seasonal values across the 4 years. The zonal variation of L _ρ for DJF and JJA is shown in Figure 6 (solid curves). The zonal pattern is very similar to that of L _α even though the cloud fraction overlap parameter and the rank correlation of COD are distinct physical parameters describing different aspects of cloud vertical structure. A model-oriented parameterization vs. latitude and day of the year is again accomplished with a Gaussian function of the same type as for L _α (Eq. 7), with parameter values included in Figure 6 which shows parameterized curves for January 1st and July 1st (dotted curves).

Because our COD filling procedure assigns equal values of liquid COD to cells when they belong to rays (cloudy subcolumns) where TAU from MODIS is available to be used as constraint, a tendency to overestimate ρ might have been expected. Such a potential overestimate implies higher values of L _ρ which is however not seen; in reality the unfilled COD fields imply somewhat higher rank correlations (right panel of Figure 4). First of all, one should keep in mind that equal COD values in two different layers do not necessarily have the same rank. Second, the lower rank correlation of the filled COD field may be explained by the fact that many cells are filled by a secondary procedure that utilizes the available COD of neighboring liquid cells, and this process may actually increase randomness and therefore decrease rank correlations. Third, lower rank correlations can also result from an increase in the frequency of distant cloud layer pairs after filling, and such pairs are expected to have COD distributions that are more uncorrelated.

The peak magnitudes of L _ρ at low latitudes are somewhat smaller than those of L _α , but at some higher latitudes the values of L _ρ surpass those of L _α . This behavior differs from that in O12 where L _ρ < L _α was universal, but with L _ρ calculated directly from radar reflectivities of cloudy cells, and not actual cloud retrievals, and in Oreopoulos and Norris (2011) who used cloud condensate retrievals from ground-based radar. In an analysis of cloud resolving model (CRM) fields Pincus et al. (2005) found that L _ρ > L _α is possible depending on the subset of clouds used and the method of calculation. Räisänen et al. (2004) in their own analysis of different CRM fields find L _ρ < L _α consistently. O22 found that halving the magnitude L _α while at the same doubling the magnitude of L _ρ compared to the original parameterization of O12 did not have serious undesired consequences on the radiative effects of cloud fields constructed by the Raisanen cloud subcolumn generator, which implies that L _ρ > L _α is plausible and not necessarily unphysical.