JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 109, D04104, doi:10.1029/2003JD003742, 2004 Spatial statistics of marine boundary layer clouds Gregory M. Lewis1 The Fields Institute, Toronto, Ontario, Canada Philip H. Austin Atmospheric Sciences Programme, University of British Columbia, Vancouver, British Columbia, Canada Malgorzata Szczodrak Rosentiel School of Marine and Atmospheric Science, University of Miami, Miami, Florida, USA Received 1 May 2003; revised 17 December 2003; accepted 9 January 2004; published 25 February 2004. [1] An analysis is presented of the structure functions and scalar spectra for 25 satellite- derived marine stratocumulus cloud optical depth fields. The scenes, which cover a horizontal domain of 58 Â 58 km at a resolution of 28.5 m, are partitioned into two ensembles on the basis of cloud fraction. For the fully cloudy scenes, although there is wide scene-to-scene variability, both the average isotropic scalar spectrum and the average isotropic second-order structure function exhibit power law behavior over approximately two decades, with scale-invariant exponents equal to those expected for inertial-subrange passive tracer fluctuations. Higher-order structure functions show anomalous scaling that closely matches that observed for wind tunnel temperature fluctuations and for other fully cloudy observations. The partly cloudy scenes, while scaling, show different behavior. The average isotropic second-order structure function and average isotropic scalar spectrum have scale-invariant exponents that are significantly smaller than those of the fully cloudy scenes, and the analysis of the higher-order structure functions indicates that the field has much more intermittent fluctuations than the fully cloudy scenes. Fits to random cascade models for the fully cloudy scenes show that the increment statistics are consistent with an underlying log normal distribution. For the partly cloudy scenes a divergence of higherorder moments is predicted, indicating that the field fluctuations are necessarily derived from fat-tailed distributions and that there will be significant realization dependence of the measured statistics. In addition, the presence of long-range correlations in all the data predicts that single-point histograms of the field values will have significant scene-to-scene variability, or equivalently, the use of spatial averages in the approximation of the parameters of the single-point probability density function of the field will result in random INDEX TERMS: 0320 Atmospheric Composition and fluctuations of the estimated parameters. Structure: Cloud physics and chemistry; 3250 Mathematical Geophysics: Fractals and multifractals; 3307 Meteorology and Atmospheric Dynamics: Boundary layer processes; 3359 Meteorology and Atmospheric Dynamics: Radiative processes; 3360 Meteorology and Atmospheric Dynamics: Remote sensing; KEYWORDS: cloud radiative structure, structure functions and scalar spectra, cloud parameterization Citation: Lewis, G. M., P. H. Austin, and M. Szczodrak (2004), Spatial statistics of marine boundary layer clouds, J. Geophys. Res., 109, D04104, doi:10.1029/2003JD003742. 1. Introduction [2] Satellite measurements of cloud reflectivity, and inferred cloud properties such as optical depth (t) and liquid water path (lwp), contain detailed information about cloud structure. In particular, the one-point statistics (e.g., the mean and the variance) and the two-point statistics (e.g., the scalar spectrum and the pth-order structure functions) of 1 Now at Faculty of Science, University of Ontario Institute of Technology, Oshawa, Ontario, Canada. Copyright 2004 by the American Geophysical Union. 0148-0227/04/2003JD003742$09.00 these fields give a statistical description of the cloud layer that can provide insight into underlying cloud processes and introduce strong constraints on the results of cloud resolving models. [3] Estimates of these cloud statistics, and their dependence on the state of the atmosphere, are also needed for many subgrid-scale cloud parameterizations being developed for use in large-scale models. Such ‘‘statistical’’ cloud schemes typically specify the subgrid-scale fluctuations of the excess saturation e = qv À qs + s, where qv is the mean total water content in the cell prior to condensation, qs is the grid cell mean saturation vapor density, and s is a random variable with an assumed distribution. Liquid water or ice D04104 1 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS occurs in a cell when e > 0; the resulting vertical profile of optical depth is then combined with a cloud-overlap assumption to calculate the radiative flux profiles in the column [e.g., Sommeria and Deardorff, 1977; Barker, 1996; Pincus and Klein, 2000; Larson et al., 2001; Smith and Del Genio, 2002]. [4] Subgrid-scale models specify the form of the s probability density function (PDF) and either diagnose or predict the PDF parameters from the resolved dynamic and thermodynamic fields. Example distributions include two parameter PDFs such as the symmetric triangle function of Smith [1990] or three parameter distributions that incorporate the skewness and/or the possibility of bimodality [e.g., Lewellen and Yoh, 1993; Bony and Emanuel, 2001; Tompkins, 2002]. The diagnostic or prognostic equations for the distribution variance and skewness may use information from subgrid-scale turbulence and convection schemes or the resolved variance in a region surrounding the grid cell [Cusack et al., 1999; Tompkins, 2002]. [5] Single-point histograms of satellite, radar, and lidar measurements of optical depth or liquid water/ice path can, in principal, be used to constrain the e PDF at the 10 to 200 km spatial scales typical of mesoscale and climate model grids [Smith and Del Genio, 2002; Carlin et al., 2002; Jeffery and Austin, 2003]. However, across this range of scales, the two-point statistics (in particular, the scalar spectra and the structure functions) of cloud radiance or optical depth measurements typically satisfy a power law relation with a scale-invariant exponent which indicates the presence of long-range correlations [e.g., Cahalan and Snider, 1989; Lovejoy et al., 1993; Barker and Davies, 1992; Gollmer et al., 1995; Davis et al., 1999]. Such a power law relation implies a statistical symmetry that is usually referred to simply as ‘‘scaling.’’ Determining the nature and extent of this scaling is important both because the prognostic equations for the PDF moments require estimates of the unresolved statistics and because establishing the extent of long-range correlations is a prerequisite to obtaining reliable estimates of the one-point statistics. [6] An examination of departures from scaling can also provide information about cloud structure and cloud-radiation interactions. For example, Oreopoulos et al. [2000] used Landsat Thematic Mapper images of two cloud scenes to show that variations in the shape of the scalar spectra of radiance fields at scales from $1 to 5 km were consistent with fluctuations caused by variations in cloud-top height, while Wood and Hartmann [2001] used local maxima in the scalar spectrum of satellite-retrieved liquid water path to study the climatology of organized boundary layer convection. Similarly, Davis et al. [1997] used the scalar spectrum and second-order structure function from Landsat radiances to show that a break in scaling at small scales (%300 m) is consistent with radiative smoothing, while von Savigny et al. [2002] argued that the break in scaling at larger (%20 km) scales, that is observed in the same Landsat scene, is consistent with saturation in the radiance field. [7] Additional information about the optical depth PDF can be obtained using higher-order two-point statistics. Cloud radiance or optical depth fields typically have skewed histograms that indicate non-Gaussian tails in the optical depth PDFs. Such tails are a signature of intermittency, which can be described by multifractal analyses using D04104 higher-order structure functions or other techniques such as the wavelet maximum modulus transform (WMMT) and singular measures. Pierrehumbert [1996] examined higherorder structure functions from 30 days of satellite measurements of emitted longwave radiance from high clouds in the tropics and found that the multifractal properties were very similar to those of temperature fluctuations observed in wind tunnel experiments. Arneodo et al. [1999] found nearly identical scaling behavior in the visible reflectivity field of a single Landsat image. Lovejoy et al. [2001] determined the ‘‘universal multifractal’’ parameters for an ensemble of 909 satellite visible and infrared radiance measurements with resolutions varying from 1 to 5 km and found that the statistics were consistent with generation by an anisotropic cascade spanning scales from 1 to 5000 km. [8] In this paper, we analyze the two-point statistics of the optical depth of 12 fully cloudy and 13 partly cloudy marine stratocumulus boundary layers using the Landsat-retrieved optical depth data of Barker et al. [1996]. We treat the fully cloudy scenes and the partly cloudy scenes as realizations from two separate ensembles. We compare the character of the scaling and the scene-to-scene variability using both the scalar spectrum and the second-order structure function. These two-point second-order statistics, which are related by an integral transform, give complimentary estimates of the spatial variability and are affected differently by cloudfree or saturated pixels and by the finite extent of the satellite scenes. These differences lead to discrepancies in the analysis which we investigate in order to gain further insight into the character of the variability of the cloud layers. We also use higher-order structure functions to compare the intermittency of the two data sets with laboratory experiments, in situ aircraft measurements, and other satellite data sets. [9] In section 2 we discuss two-point statistics, paying particular attention to the second-order statistical measures: the second-order structure function and scalar spectrum. We discuss the satellite data set in section 3 and present results of the second-order statistics in section 4. In section 5 the observed discrepancies in the second-order statistics are discussed, and results of statistics of order other than second are presented in section 6. A discussion follows in section 7. 2. Two-Point Statistics: The Scalar Spectrum and Structure Functions 2.1. Scalar Spectrum [10] The spectral density, É(k), of a statistically homogeneous (but possibly anisotropic) two-dimensional random field, t(x), is defined by ÉðkÞ ¼ ð2pÞÀ2 Z htð0ÞtðrÞieikÁr d 2 r Z 2 ¼ tðxÞeikÁx d 2 x ; ð1Þ ð2Þ where x = [x, y] is the position, ht(0)t(r)i is the (auto) covariance function at lag r = x0 À x = [rx, ry], the 2 of 16 LEWIS ET AL.: SPATIAL STATISTICS D04104 angled braces, hÁi, indicate an ensemble average (i.e., an average over all possible realizations of the field), k = [kx, ky] is the wave number, and the integrals are over all of the respective space. Here we use the term statistically ‘‘homogeneous’’ as the spatial equivalent of statistically ‘‘stationary,’’ which is generally used for time-dependent processes. Equation (1) is known as the Wiener-Khinchin formula, and equation (2) indicates that the spectral density É(k) is the modulus squared of the Fourier transform of the field t(x) itself. [11] If the field is statistically isotropic, then É(k) = É(k), where k = jkj, and an annular integration can be performed to obtain the isotropic scalar spectrum E(k) = 2p kÉ(k), which is also referred to as the energy spectrum or the power spectrum. For a given (anisotropic) realization of an isotropic field or for a marginally anisotropic field, E(k) is still meaningful and is defined as Z 2p EðkÞ ¼ ÉðkÞ k dq; ð3Þ 0 where q = tanÀ1 (ky/kx) is the polar angle. [12] If the increments t(x + r) À t(x) of an isotropic field are homogeneous, then the isotropic scalar spectrum E(k) of the field, as defined by equations (2) and (3), satisfies Z k0 k 2 EðkÞdk < 1 ð4aÞ 0 Z 1 EðkÞdk < 1 ð4bÞ k0 for any k0 [e.g., Monin and Yaglom, 1975, p. 88]. Scalar spectra of in situ cloud liquid water content, cloud reflectivity fields, and passive scalar fluctuations in fully developed isotropic turbulence are all observed to satisfy (4) [Barker and Davies, 1992; Davis et al., 1996, 1999; Oreopoulos et al., 2000]. [13] If, in addition to equation (4), the scalar spectrum satisfies D04104 second-order structure function instead of the covariance function. 2.2. Structure Function [15] The structure function Sp(x, r) of order p of a field t(x) is given by: Sp ðx; rÞ ¼ hjtðx þ rÞ À tðxÞjp i; ð6Þ where x is the position, r is a lag with respect to x, and the angled braces, hÁi, indicate an ensemble average. The second-order structure function S2 is obtained by setting p = 2 in equation (6). Structure functions of high order (large values of p) provide information regarding the spatial distribution of the large fluctuations of the field, while the smaller fluctuations generally dominate the lower-order structure functions. A comparison of the lower-order with the higher-order statistics can be used to characterize the intermittency of the field as described in section 6 below. [16] If increments of the field t(x + r) À t(x) are assumed to be statistically homogeneous, Sp(x, r) is only a function of the lag r (is therefore still multivariate) and may be written as Sp(r). In this case the ergodic theorem indicates that a spatial average (over x) may be used to approximate the ensemble average in equation (6). [17] If the field is statistically isotropic as well as having homogeneous increments, Sp(r) is only a function of r = jrj, the distance between the locations of the field values. We can then write Sp = Sp(r), where Sp(r) will be called the isotropic structure function of order p. For such a field, Sp(r) will have contours that are concentric circles, and thus an average can be taken about annuli of constant r to obtain a more reliable approximation of Sp(r). [18] The scalar spectrum E(k) and the second-order structure function S2(r) of an isotropic two-dimensional field are related by the extended version of the Wiener-Khinchin formula Z S2 ðrÞ ¼ 2 ½1 À cosðk Á rÞ É ðkÞ dk Z 1 ½1 À J0 ðkrÞ EðkÞ dk; ¼2 ð7Þ 0 Z 1 EðkÞdk < 1; ð5Þ 0 then the field is homogeneous. In this case the integral on the left hand side of equation (5) gives the variance of the field. [14] Real fields will have finite variance and thus will automatically satisfy equation (5). However, for many geophysical fields the scale of the largest fluctuation exceeds the scale of the observation, which prevents an accurate estimate of the variance or autocovariance function of the field [e.g., Monin and Yaglom, 1975, p. 85]. When this occurs and, in addition, equation (4) is satisfied (as is the case for the spectra presented in this paper), it is useful to assume that the increments of the field are homogeneous without assuming the homogeneity of the field itself. Under these assumptions the Wiener-Khinchin formula (1) is not valid. However, as described in the next section, an extended version of equation (1) relates the scalar spectrum to the provided that E(k) satisfies equation (4), where J0 is a Bessel function of the first kind of order zero, r = jrj, and k = jkj [e.g., Monin and Yaglom, 1975, p. 87]. 2.3. Scaling [19] A field t(x) is said to be scaling (in an isotropic sense) if the structure function S2(r) satisfies S2 ðrÞ ¼ lz2 S2 ðr1 Þ ð8Þ over a range of scales, where z2 is a scale-invariant exponent, r = l r1, l is a scaling factor, and r1 is an arbitrary unit scale. If S2(r) satisfies equation (8), then it satisfies the power law relation S2 ðrÞ / rz2 ; ð9Þ and it can be described, over the given range of scales, by only two parameters: the exponent z2 and a prefactor (a multiplicative constant). 3 of 16 LEWIS ET AL.: SPATIAL STATISTICS D04104 Table 1. Solar Zenith Angle q0, Cloud Fraction Ac, Average Optical Depth t ", Maximum Optical Depth tmax, and Fraction of Saturated Pixels for the 25 Scenes Image q0 Ac t " tmax Pixels at tmax, % A6 A8 A9 A10 A11 A12 A13 A14 A15 A16 A17 A18 B1 B2 B3 B4 B5 B6 B7 B8 B9 B10 B11 B12 B13 32 30 30 30 31 28 28 32 32 32 32 32 30 30 47 28 29 29 29 29 32 37 69 58 58 1.000 1.000 0.999 0.982 1.000 0.845 1.000 1.000 1.000 1.000 0.998 1.000 0.931 0.644 0.523 0.890 0.737 0.951 0.914 0.814 0.916 0.409 0.601 0.809 0.974 17.0 17.5 10.3 9.5 16.7 14.8 13.6 13.4 17.0 18.6 14.6 20.4 9.2 3.4 2.4 6.8 4.4 14.9 11.6 4.6 5.7 3.7 15.6 3.7 6.9 25.8 25.1 20.2 25.1 25.8 39.8 39.8 46.8 48.1 48.1 48.1 48.1 25.1 20.8 39.8 39.8 42.0 40.9 42.0 42.0 22.5 37.7 100.0 100.0 100.0 5 5 0 7 11 1 0.2 2 0.5 1 3 3 5 0 0 1 0.01 11 5 0 0 0.01 4 0 0 [20] If the relation (7) holds and S2 (r) satisfies equation (9) for all r, then the scalar spectrum will also satisfy a power law relation E(k) / kÀb for all k, where b is the scale-invariant spectral exponent. Equation (4) will be satisfied for a power law E(k) as long as 1 < b < 3; however, equation (5) will not be satisfied for any b. Thus the condition 1 < b < 3 implies the validity of equation (7) and implies that the corresponding field is nonhomogeneous with homogeneous increments. In addition, given a power law E(k) with 1 < b < 3, equation (7) implies the following relationship between the spectral exponent and the exponent of the second-order structure function: z2 ¼ b À 1 ð10Þ [e.g., Frisch, 1995, p. 54]. This relation is derived under the assumption that the second-order structure function and scalar spectrum satisfy power law relations for all scales. When this assumption does not hold, equation (10) only holds approximately, and equation (7) may still be valid even when b is not between 1 and 3 (see section 5 below). For passive scalars in fully developed isotropic turbulence, the similarity argument of Kolmogorov-Obukhov-Corrsin (KOC) predicts an inertial subrange (limited range of scaling) with z2 = 2/3 corresponding to b = 5/3; for Brownian motion we expect b = 2 corresponding to z2 = 1. 3. Data Analysis [21] The cloud optical depth t data set is derived from 25 scenes of marine boundary layer clouds using radiance measurements from the Landsat 5 Thematic Mapper. Table 1 lists the solar zenith angle, scene cloud fraction, mean and maximum optical depth, and the percentage of pixels at maximum for each scene, using the scene numbers D04104 of Barker et al. [1996]. The clouds have mean optical depth 2 < t " < 20.5 and are located in the Atlantic and Pacific Oceans between the latitudes 20°N and 35°N. All scenes span an area of 58 Â 58 km and contain 2048 Â 2048 pixels at 28.5 m resolution. The optical depth values are those of Harshvardhan et al. [1994], who inferred them using channel 4 (0.83 m m) reflectance measurements and assuming a droplet size distribution with effective radius 10 mm and variance 0.1 mm2. This constant radius assumption means that the statistics for the retrieved optical depth are very similar to those of the 0.83 mm raw radiance. Variability in the radiance field may be decoupled from variability in cloud properties such as liquid water path by cloud-top shadowing at intermediate scales and photon diffusion at small scales as discussed in the next section. [22] The scenes are sorted into two classes based on contiguous cloud fraction and the histogram of the optical depth: class A, overcast stratocumulus or ‘‘fully cloudy,’’ with a cloud fraction between 0.98 and 1.0, (the one exception A12 contains a small cloud-free region at one edge of the scene); and class B, broken stratocumulus or ‘‘partly cloudy,’’ with cloud fraction between $0.4 and 0.975. The optical depth histograms, which are presented by Barker et al. [1996], are positively skewed for all scenes. While the histograms of the fully cloudy A scenes all possess a single mode, the histograms of the partly cloudy B scenes decrease monotonically from low to high t values, with the thinnest resolvable cloud layers being the most probable. Considine et al. [1997] point out that the absence of a t mode in the partly cloudy scenes is consistent with an underlying excess supersaturation distribution which has a mode em < 0. In this case, only the positive portion of the e distribution will be mirrored in the cloud liquid water distribution, and hence detectable by satellite. [23] Because the data are discrete, the two-dimensional structure function Sp(r) of order p is approximated by À Á X X À Á À Á t xk;l þ ri;j À t xk;l p ; Sp ri;j ¼ k ð11Þ l where xk,l is the pixel with position x = (x, y) = (k, l) (likewise for ri,j), xk,l + ri,j = xk+i,l+j, the image is of size N Â N, and the limits of the sum are chosen such that 1 k + i N, 1 l + j N, 1 k N, 1 l N, ÀN i N, and ÀN j N. The computation scales as N4 but parallelizes easily. The approximation to the isotropic structure function, S2(rt), is obtained by averaging S2(ri,j) over concentric annuli with radii rt, the (integer) distance from the origin. [24] Cloud-free pixels, saturated pixels or other data not representative of the underlying field, can simply be ignored in the calculation of the structure function. This possibility does not exist for the computation of the scalar spectrum. We examined the bias introduced by including saturated or cloud-free pixels in the calculations by comparing the second-order structure functions computed with these pixels included to the second-order structure functions computed when these pixels are ignored. We find that when the percentage of saturated or cloud-free pixels is less than 30%, the resulting bias is relatively small. Therefore all pixels are used in the structure function calculations that are presented below. [25] The two-dimensional spectral density is approximated from the discrete Fourier transform, which is computed 4 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS D04104 Figure 1. Second-order two-point statistics for the optical depth of the fully cloudy scene A17. (a) Optical depth; (b) two-dimensional second-order structure function; (c) two-dimensional discrete spectral density; (d) isotropic structure function; and (e) isotropic scalar spectrum. Figures 1d and 1e also show straight line fits (in the log-log coordinates) to data and two reference lines with (z2 = 2/3 and b = 5/3) corresponding to values expected for passive tracer fluctuations in the inertial subrange. using the two-dimensional fast Fourier transform (FFT). A Parzen window is applied to the image prior to applying the transform. 4. Second-Order Statistics for Landsat Imagery [26] The two-dimensional second-order structure function S2(r), the isotropic second-order structure function S2(r), the two-dimensional spectral density É(k), and the isotropic scalar spectrum E(k) are calculated for all the fully cloudy and partly cloudy scenes of Table 1. Figure 1 shows all four functions for the fully cloudy scene A17, which has a moderately thick t = 14.6 and 3% of pixels at saturation. A comparison of Figures 1a and 1b shows that S2(r) provides a much clearer view of the A17 scene anisotropy than É(k), primarily because the Fourier transform has poor resolution at large wavelengths, but also because neighboring Fourier coefficients are less correlated than adjacent values of the structure function. [27] For S2(r) and E(k), power law behavior appears as linear behavior in log-log coordinates. Thus Figures 1d and 1e show a scaling regime that begins at r % 150 m and spans 2 decades for S2(r) and about 3 decades (to the largest resolved scale) for E(k). Figure 1d includes r increments larger than 30 km which are not shown in the full twodimensional plot of Figure 1b. Scene anisotropy is removed by the annular averaging; this would not be the case for the one-dimensional scalar spectrum E1(k1) which is calculated using one-dimensional slices along an arbitrary direction through the scene, where k1 indicates the wave number along the chosen direction of the slices. [28] At smaller length scales, both S2(r) and E(k) exhibit a fall-off that is consistent with the photon diffusion smoothing discussed by Davis et al. [1997]. The least squares best fit (in log-log coordinates) to the scaling regimes for S2(r) and E(k) is given by the light lines in each figure, where the slope of the line gives the value of the scale-invariant exponent. For this scene the exponent values of z2 = 0.63 and b = 1.66 are very close to the (z2, b = 2/3, 5/3) values from KOC theory. [ 29 ] Figure 2 shows the optical depth field and corresponding second-order statistics for scene B2, which 5 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS D04104 Figure 2. As in Figure 1, but for the partly cloudy scene B2. has a cloud fraction of 64% and a mean optical depth t = 3.4. As Figures 2c and 2d show, the scaling regime in scene B2 is limited to a single decade for the structure function and to 2 decades for the scalar spectrum, with smaller values of (z2, b) = (0.43, 1.39). Once again, E(k) exhibits power law behavior to the largest resolved scale, while S2(r) indicates a break in the scaling at r % 10 km. For both Figures 2c and 2d, there is a break in the scaling at a scale of $1 km. [30] In Figure 3 the isotropic structure functions S2(r) and scalar spectra E(k) for all the fully cloudy scenes are presented together to facilitate a comparison of the qualitative features of the plots. In order to determine the variability of the scale-invariant exponents and the extent of the scaling regimes, lines are fit to all S2(r) and E(k) in log-log coordinates. As seen in Figure 3, it is generally not possible to fit a line over the whole range of scales, and so ranges for the fit are subjectively chosen over which the S2(r) and E(k) showed reasonable power law behavior. [31] For the annularly averaged structure functions, three different regimes are generally observed: (1) An inner regime consistent with photon smoothing reported by Davis et al. [1997] (most prominent in scenes A10 and A13 – A18). In this regime z2 is between 1.0 and 1.2 across scales ranging from the smallest resolved scale of 28.5 m to between 200 and 400 m. (2) A middle regime beginning at a lower scale that ranges from %50 to %400 m and ends at an outer scale that ranged from as small as 3 km (A9), to as large as 40 km (A10). In this regime, z2 has a mean of 0.68 ± 0.08, and 7 of the 12 scenes have z2 between 0.6 and 0.7. (3) An outer regime which is seen in most cases at scales larger than $10 km. In this regime, z2 is much smaller than 2/3 and is often close to zero. Scene A12 is a moderate exception to this general form with a z2 of 0.89 between 50 and 6 km and a z2 of 0.41 between 8 and 30 km. [32] As seen in Figures 1 and 2, a scene-by-scene comparison of the isotropic structure functions and the isotropic scalar spectra, shown in Figure 3b, reveals that, although an outer regime is clearly observed in the structure function, none is observed in the scalar spectra. For example, in Figure 3b, the scalar spectrum of scene A11 clearly increases with scale (corresponding to a decrease in wave number) up to the largest resolved scaled, while a roll-off in the structure function is observed at scales larger than $10 km. [33] There is also a noticeable change in slope in many of the scalar spectra starting at 1/k % 0.1 km. Examination of the images (not shown) indicates that the scenes with the strongest departures from scaling (A9, A15, and A16; see also B3 below) had clearly developed cellular structures. 6 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS Figure 3. (a) Isotropic second-order structure functions for all fully cloudy scenes. (b) Isotropic scalar spectra for all fully cloudy scenes. Curves are offset for clarity. Also shown are reference lines with slopes z2 = 2/3 and b = 5/3. The presence of a local break in scaling owing to organized mesoscale structure has been seen in both AVHRR and MODIS images [Barker and Davies, 1992; Wood and Hartmann, 2001] and in several subscenes of a Landsat image by Oreopoulos et al. [2000], who demonstrated that this ‘‘radiative roughening’’ of the radiance field is particularly apparent at low Sun angles when shadowing occurs owing to cloud-top variability. [34] Figure 4 presents all isotropic structure functions S2(r) and scalar spectra E(k) for the partly cloudy scenes. The isotropic structure functions, shown in Figure 4a, show a variety of qualitative differences, but generally two regimes are observed. The location of the transition between the two regimes varies from separations as small as 300 m (scene B13) to as large as 30 km (scene B6). Scenes B5 and B10 – B13 show a transition at $1 km, while most other scenes show a transition above 5 km. Values of z2 at small scales lie between 1 and 1.2. At large scales, z2 exponents are observed to range from near zero (scene B11) to 0.67 (scene B6) but are generally less than the values of z2 for the fully cloudy scenes. [35] The relation (10) between z2 and b is exactly satisfied by few individual scenes in either the fully cloudy or the partly cloudy scenes. Of the 12 A scenes for which scaling ranges exist, the ratio Á = jb À 1 À z2 j/b is less that 20% for only four, with a mean value Á = 29%. For the partly cloudy B scenes, Á <20% for only 2 of 12 scenes, with Á = 33%. However, the range of scales over which the D04104 scale-invariant exponents are measured is not chosen identically for the structure function and scalar spectrum of individual scenes; this is expected to influence these results. [36] Summarizing the individual features of Figures 3 and 4: [37] 1. Individual scenes exhibit wide variability in the location and extent of their scaling regimes and in the values of their scale-invariant exponents z2 and b. In addition, the exponents satisfy equation (10) only approximately. [38] 2. At the largest scales the isotropic structure function S2(r) and scalar spectra E(k) exhibit qualitatively different behavior. The structure functions plateau at large scales, while the isotropic scalar spectra increase up to the largest resolved scales. [39] 3. The z2 values for the slopes in the partly cloudy scenes are smaller than for the fully cloudy scenes. [40] The significant scene-to-scene variability in the structure functions and scalar spectra is not an indication that the optical depth field is not scaling. The structure function, and so scaling, is defined over an ensemble average (likewise for the scalar spectrum). In fact, each realization (cloud scene) is expected to break the scaling in some way; the scaling is excepted to be recovered only if sufficient data are averaged. We assume that the fully cloudy and partly cloudy scenes are realizations from different ensembles, and thus to approximate the two ensemble averages, we compute, separately, an average over all of the fully cloudy scenes and an average over all of the partly cloudy scenes. Specifically, we compute a point-wise average of the two-dimensional statistics and then compute annular averages to obtain two sets of Figure 4. As in Figure 3, but for all partly cloudy scenes. 7 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS D04104 Figure 5. Average of data over all fully cloud scenes. (a) Average for the two-dimensional second-order structure function; (b) average for the two-dimensional discrete spectral density; (c) average for the isotropic second-order structure function; and (d) average for the isotropic scalar spectrum. In Figures 5c and 5d, vertical lines mark reference scales, between which the structure function/scalar spectrum show clear scaling behavior. Also shown are reference lines with slopes z2 = 2/3 and b = 5/3. isotropic statistics. For such an average the individual structure functions and scalar spectra that are of higher average value will contribute more to the combined statistics. The results are presented in Figures 5 and 6 for the fully cloudy and the partly cloudy scenes, respectively. [41] As Figure 5 shows, the average structure function for the fully cloudy scenes is isotropic and displays clear power law behavior between 200 and 10 km. Below 200 m, there is a transition to the photon smoothing regime. The scalar spectrum is approximately power law Figure 6. As in Figure 5, but for averages over all partly cloudy scenes. 8 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS between %300 m and 60 km (the largest resolved scale), and the scale-invariant exponents are measured to be (z2, b) = (0.66, 1.66). Furthermore, in contrast to many of the individual realizations, the average second-order statistics of the A scenes satisfy equation (10) almost exactly. The discrepancy between the isotropic structure functions and the scalar spectra that is seen at the large scales of the individual realizations carries over to the average statistics: Figure 5 shows a clear break in scaling of S2(r) at 10 km which is not mirrored in E(k). The average scalar spectrum also shows slight deviations from scaling, while the scaling regime for S2(r) is nearly exact. [42] The statistics averaged over all B scenes are shown in Figure 6. The average two-dimensional statistics are approximately isotropic, although some anisotropy is still observed at large scales. As with the fully cloudy averaged statistics, unambiguous power law behavior is observed in the isotropic structure function between 1 and 20 km, while the isotropic scalar spectrum shows growth to the largest resolved scales. Values of (z2, b) = (0.22, 1.02) are both smaller than for their fully cloudy counterparts, and in contrast to the average of the fully cloudy scenes, the two scaling exponents do not satisfy equation (10); in particular, b predicts a rougher field than z. 5. Investigation of the Discrepancies Between S2(r) and E(k) [43] The second-order structure function S2(r) and the scalar spectrum E(k) are related by the integral transform equation (7) and should, in principle, supply redundant information. However, for the fully cloudy scenes, the average of the second-order structure functions exhibits a break in scaling at $10 km, while the average of the scalar spectra indicates no such break. Similarly, the partly cloudy scenes exhibit a scale break in the average structure function at $20 km, with no corresponding break in the average scalar spectrum. A second discrepancy is also seen in the statistics of the partly cloudy scenes; the z2 = b À1 relation that is expected for scaling fields does not hold. In this section we discuss the possible origins of these discrepancies. [44] One explanation for the discrepancy between z2 and b in the partly cloudy scenes is that the Fourier basis functions that are used in the computation of the scalar spectrum have difficulty in resolving sharp cloud-clear transitions, leading to estimates of b smaller than the true value. To test this, we have simulated partly cloudy layers with specified scaling behavior by applying a threshold to scaling log normal random fields [Lewis and Austin, 2002]. Comparing the structure function and scalar spectrum for these fields (not shown) indicates that the structure function does in fact produce a more accurate estimate of the scale-dependent variability than does the scalar spectrum. [45] We next consider the discrepancy that is observed at large scale. A break in the power law behavior of the structure function may arise from several possible causes. As von Savigny et al. [2002] argue, a scale break may be induced when the field has a limited dynamic range, i.e., the field is bounded. In this case the structure function will have a maximum value, and therefore power law growth cannot D04104 be maintained at all scales. A limited dynamic range in the optical depth field could be due to a physical mechanism such as precipitation thinning of thicker clouds or from retrieval issues such as the insensitivity of cloud reflectivity to cloud thickness at large optical depths. The Landsat Thematic Mapper also has a limited dynamic range and cannot measure the reflectivity of the brightest pixels in thick stratocumulus clouds. [46] A scale break may also be introduced if a field does not have a large-scale trend or if a trend is removed from a field. That is, if r is of the order of the image size L, and t(x + r) % t(x), then the increment t(x + r) À t(x) will be close to zero. Thus, as in the case of a limited dynamic range, the structure function will be unable to sustain power law growth at large r, and a break will be observed. In contrast, the calculation of the scalar spectrum assumes periodic boundary conditions, i.e., t(x + L) = t(x), and errors (aliasing) will be introduced if the field violates this assumption. This divergence between the large-scale behavior of the structure function and scalar spectrum is similar to that observed in both the fully and partly cloudy scenes. [47] Results on simulations (not shown) suggest that both a limited dynamic range and a lack of trend can induce a scale break in the structure function at a scale of approximately half the image size (L/2), without a corresponding break in the scalar spectrum. However, a break in the structure function at scales smaller than L/2 could not be induced by limiting the dynamic range or by trend removal. Thus the scale break observed in the fully cloudy scenes appears to occur at scales too small to result from a limited dynamic range or a lack of trend alone. [48] In order to explore further possible circumstances under which a scale break may be observed in S2(r) and not in E(k), we investigate the direct theoretical relation between the two statistics, i.e., the extended version of the Wiener-Khinchin formula given by equation (7). Before we discuss the consequences of a scale break, it is instructive to give a short derivation of the z2 = b À 1 relation. Therefore assume that for all wave numbers k the scalar spectrum satisfies the power law relation EðkÞ ¼ E0 jkjÀb ; ð12Þ where b is the scale-invariant spectral exponent, and E0 is a constant (prefactor). We substitute this form of E(k) into the integral relation (7), and, with r fixed, make a change of variables s = kr to obtain Z 1 S2 ðrÞ ¼ 2E0 ! ½1 À J0 ðsÞsÀb ds rbÀ1 ¼ Crz2 ; ð13Þ 0 where z2 = b À 1, C = 2 E0A is a constant when the integral Z 1 A¼ ½1 À J0 ðsÞsÀb ds ð14Þ 0 converges. Using standard techniques, it can be shown that equation (14) converges whenever 1 < b < 3. This ends the derivation. [49] Now consider the case when the power law relation does not hold over all k, but instead the scalar spectrum has 9 of 16 LEWIS ET AL.: SPATIAL STATISTICS D04104 D04104 two scaling regimes with an abrupt transition at the wave number k1. That is, the scalar spectrum is given by EðkÞ ¼ Ea jkjÀba ; 0 < jkj < k 1 ; Eb jkjÀbb ; k 1 < jkj < 1; ð15Þ where ba, bb, Ea and Eb are constants, and to ensure continuity of E(k) at k = k1, we must have Eb = Eak1Àba+bb. [50] As above, we substitute E(k) given by equation (15) into the integral relation (7) and obtain Z k1 S2 ðrÞ ¼ 2Ea ½1 À J0 ðkrÞk Àba dk 0 Z 1 þ 2Eb ½1 À J0 ðkrÞk Àbb dk: ð16Þ k1 We make a change of variables s = kr in all the integrals and obtain S2 ðrÞ ¼ Ia ðrÞrza þ Ib ðrÞrzb ; ð17aÞ where za = ba À1, zb = bb À1, Z k1r Ia ðrÞ ¼ 2Ea ½1 À J0 ðsÞsÀba ds; ð17bÞ ½1 À J0 ðsÞsÀbb ds; ð17cÞ 0 and Z 1 Ib ðrÞ ¼ 2Eb k1r whenever the integrals converge. The limits of integration now depend on r. [51] If Ia and Ib were Heaviside functions which turn on and off, respectively, at some number, then S2(r) would exhibit a scale break (at this number) that separates two different scaling regimes (as the scalar spectrum does). This is not the case, and, technically, the scaling and the relation z2 = b À 1 is only recovered asymptotically in the limits r ! 1 and r ! 0. However, for practical purposes, two scaling regimes, which are separated by a transition region, can be considered to exist. [52] In Figure 7a an example is plotted in logarithmic coordinates (with k1 = 5, ba = 4/3, and bb = 2), where, for each r, the integrals Ia and Ib are approximated numerically using the trapezoidal rule. The graph covers four orders of magnitude with the outer radius L = 10; the behavior for different L is essentially identical. By a property of Ia and Ib, a break in S2(r) seems to occur near r1 = 1/k1, as expected. However, the figure also shows a transition region between a small-scale regime and a large-scale regime, both of which display approximate scaling. The main part of this transition region covers $1/2 order of magnitude from scales larger than r = 1/(2 k1) to scales smaller than r = 1/k1. Thus the scale break in the small scale regime is observed at approximately r = 1/(2 k1). [53] Furthermore, although the transition appears to cover about 1/2 order of magnitude, the graphs on either side are not actually straight lines, i.e., a power law relation is not satisfied. It is necessary to go at least an order of magnitude Figure 7. (a) Structure function calculated using equation (17) with k1 = 5, ba = 4/3 and bb = 2. Lines are fit to regimes before and after the break. The corresponding slopes are zb = 1.10 and za = 0.47. Reference lines with these slopes are drawn with an offset above the curve. Also shown are vertical lines at r = r1 = 1/k1 and r = r1/2 = 1/(2k1) bracketing the transition between scaling regimes. (b) The corresponding integrals Ia and Ib from equation (17). (on each side) before the graphs become approximately linear, and scaling is approximately recovered. This can be seen in Figure 7, where lines with slopes corresponding to estimated exponents are drawn offset from the plot. The estimates of the exponents are za = 0.47 and zb = 1.10 (where za = 1/3 and zb = 1 would be expected from the asymptotic relation z2 = b À 1). Other examples (not shown) indicate that the transition is wider and errors larger as the difference between ba and bb grows. [54] In Figure 7b the integrals Ia and Ib are plotted together. Here the difference between the integrals and step functions, which results in the deviation from power law behavior, can be seen. The r axis is in logarithmic units so that it can more easily be related to Figure 7a. Although the basic features of Ia and Ib are similar to Heaviside functions, it does not seem that Ia has reached the ‘‘asymptotic range’’ 10 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS even for scales that are almost two orders of magnitude larger than the scale of the break. For small r, it seems that this range has been reached. [55] Thus, in the presence of a break in the scaling, the relation between the scale-invariant exponents of the scalar spectrum and the structure function z2 = b À 1 no longer holds for all scales, although the relation is approximately recovered for scales far from the scale of the break. In addition, if the scale break in the scalar spectrum occurs at wave number k1, then the corresponding break in the structure function occurs at r1 = 1/(2k1). This indicates that if no break is seen in the scalar spectrum of a field of length L, then a break in the structure function should not be seen at scales smaller than L/2. This is consistent with the statistics of the partly cloudy scenes, where the break observed in the structure function was close to half the image size. However, the results alone do not explain the discrepancy between the range of scaling of the scalar spectrum and the range of scaling of the structure function that is observed in the fully cloudy scenes. [56] Thus we consider two more examples. In the first we assume an integral scale at 20 km, and we consider similar scales as resolved by the data (L % 60 km). That is, we take ba = 0, bb = 5/3, and k1 = 1/(20 km). In this case, even though ba = 0, the integral Ia still converges for all r, and za tends asymptotically to zero; i.e., the relation z2 = b À 1 does not hold even in the limit r ! 1. The corresponding S2(r) and E(k) are presented in Figure 8, where we can see that the break in the scaling of the structure function occurs at a scale similar to the observed break in the structure function of the fully cloudy scenes of Figure 3. Although an integral scale is not observed in the average scalar spectrum, the break could be masked in the observed spectra by the aliasing that can occur owing to the assumption of periodicity in the calculation of the scalar spectra; i.e., power from wave numbers of scales larger than the image size could obscure the break by introducing power in the largest resolved wave numbers. Of course, the presence of power in wave numbers larger than the image size implies that the scale break does not correspond to an integral scale. [57] This leads to the proposition that another scaling regime exists for scales larger than that of the data. That is, there is a ba = 5/3 for scales ranging from below 1 km to 10 km, a bb = 0 for scales close to 10 km, and at scales larger than $20 km there is another (unknown) exponent bc. Making the assumption that the scalar spectrum has three scaling regimes with breaks at k1 and k2, the same analysis as above can be performed to investigate the effects this would have on the structure function. [58] Figure 9 shows the results of a third example, where k1 = 1/30, k2 = 1/20, ba = 5/3, bb = 0, and bc = 5/3. The figure shows qualitatively that the corresponding scalar spectrum, with a flattening over a relatively small range of scales, can produce a second-order structure function with characteristics similar to those of the structure function calculated for the fully cloudy scenes. Such a restricted deviation from scaling may not be observable at this scale in the scalar spectrum owing to the aliasing that would occur in this case. A more detailed comparison to the observed statistics reveals that this example produces an estimate of the scale-invariant exponent for the structure function that is too large at scales just below 10 km and produces a D04104 Figure 8. (a) The structure function that would be observed if (b) a field had the scalar spectrum plotted, i.e., k1 = 1/20 kmÀ1, ba = 0 and bb = 5/3. A scale break at $20 km in the scalar spectrum produces a scale break in the structure function similar to that seen in the data. Such a break in the scalar spectrum is not observed the Landsat data. structure function which decreases for scales just above 10 km. Smoothing the idealized scale break would bring the structure function of this example into closer agreement with that calculated from the fully cloudy scenes. These results indicate that an integral scale does not occur at scales that are resolved by the data, and thus long-range correlations must be present. 6. Statistics of Order Other Than Second: Multifractal Analysis [59] Two random fields that have the same second-order two-point statistics may still display very different variability. In particular, the intermittency of a field cannot be found from second-order statistics alone and must be determined by evaluating the structure function defined by equation (6) for a range of moments p. For scaling fields the intermit- 11 of 16 LEWIS ET AL.: SPATIAL STATISTICS D04104 D04104 of the increments t(x + r) À t(x) for separation r, provided the tails of the PDF decrease at least exponentially [Vainshtein et al., 1994]. Although observational constraints typically limit the determination of z(p) to a restricted range of p, a fit to measured z(p) values can be compared to laboratory measurements, other atmospheric observations, and to random cascade models which can predict z(p) for all p. [61] Figure 10 shows z(p) data for the fully cloudy (FC) and the partly cloudy (PC) scenes analyzed in section 4, together with a combination of other laboratory, aircraft, and satellite data sets summarized in Table 2. The symbols mark the individual data points which are the exponents z(pj) approximated from the isotropic average structure function of order pj, where j = (1, 2, 3, 4, 5). As in section 4, the isotropic average structure functions are calculated separately for the FC and PC scenes by first averaging the twodimensional statistics over all scenes and then computing annular averages. The z(pj) values are determined by linear regression on log-log plots over a range of scales where the average structure functions exhibit power law behavior; for both scene types and for all orders pj, clear power law behavior is observed over approximately the same range of scales as for the second-order average structure functions (see Figures 5 and 6). [62] The curves on Figure 10 are fits to z(p) data using either an empirical two-parameter function introduced by Pierrehumbert [1996] or the three-parameter ‘‘universal multifractal’’ model of Schertzer and Lovejoy [1991]. The Pierrehumbert [1996] function has the form zðpÞ ¼ z00 p : 1 þ z00 p=z1 ð19Þ 0 The coefficient z 0 = dz/dp|p=0 is a measure of the average size of the fluctuations; the most probable fluctuations on Figure 9. (a) The structure function that would be observed if (b) a field had the scalar spectrum plotted, with three scaling regimes with k1 = 1/30 kmÀ1 (r1 = 30 km), k2 = 1/20 kmÀ1 (r2 = 20 km), ba = 5/3, bb = 0, and bc = 5/3. As in Figure 8, this scalar spectrum produces a scale break in the structure function similar to that seen in Figure 3a. The scale breaks in Figure 9b are qualitatively consistent with the observed cloudy spectra of Figure 3b. tency can be characterized from the variation of the scaleinvariant exponents z(p) with p, where z(p) is defined by a generalized version of (8): Sp ðrÞ / rzðpÞ : ð18Þ A smooth field viewed at scale r has finite gradients and exponents z(p) / p. This is known as simple scaling; fields with structure functions that deviate from this relationship are intermittent and are said to be multifractal or to exhibit anomalous scaling [e.g., Frisch, 1995, p. 143]. The degree of intermittency can be characterized in terms of the deviation from simple scaling, where larger deviation from simple scaling indicates increased intermittency of the field. [60] Knowledge of the exponents z(p) for all p, with corresponding prefactors, is equivalent to knowing the PDF Figure 10. Multifractal analysis. The scale-invariant exponent z(p) is plotted versus the moment p. See Table 2 plus text for legend descriptions. 12 of 16 LEWIS ET AL.: SPATIAL STATISTICS D04104 D04104 Table 2. Data Sets Used for the Multifractal Analysis of Figure 10 Abbreviation Description SIM FC FCf P96 A84 R96 LWP VAPOR ASTEX SOCEX FIRE LO1 PC PCf simple scaling fully cloudy optical depth fit to FC using (20) cloud longwave radiances temperature fluctuations temperature fluctuations cloud liquid water path water vapor fluctuations, tropical boundary layer cloud liquid water concentration (LWC) fluctuations, Azores LWC fluctuations, Southern Ocean LWC fluctuations, eastern Pacific cloud long and short wave radiances partly cloudy optical depth fit to PC using equation (20) 0 Data Source the image scale as rz0. The large p asymptote z1 is related to the intermittency of the field. Specifically, z1 = 1 represents simple scaling, with decreasing values of z1 corresponding to larger probabilities that the field contains jump discontinuities [Pierrehumbert, 1996; Davis et al., 1999]. For a two-dimensional field, D = 2 À z1 represents the dimension of the subset of the field that contains the 0 discontinuities. For the FC data the best fit values for (z 0, z1) are (0.43, 2.78); the fact that D < 0 for this fit implies that while the fully cloudy scenes exhibit anomalous scaling they contain very few true discontinuities. [63] The line labeled P96 lies closest to the FC data on Figure 10; it is a fit using equation (19) to a representative z( p) taken from 30 days of infrared radiance data from the central Pacific [Pierrehumbert, 1996]. Most of the cloud cover in this region consisted of high0 anvils formed from deep convection. Although we use z 0 = 0.5, z1 = 2.25 for P96, Pierrehumbert [1996] found that z1 fluctuated 0 between 2 and 2.5 throughout the month, while z 0 was nearly constant. Also shown for comparison are two sets of wind tunnel measurements of temperature fluctuations: A84 [Antonia et al., 1984] and R96 [Ruiz-Chavarria et al., 1996]. These were taken at Taylor microscale Reynolds numbers of Rl % 850 (A84) and 130 < Rl < 470 (R96). [64] The curve labeled LWP is a fit to ground-based measurements of integrated liquid water path taken at the Atmospheric Radiation Monitoring Site in Oklahoma for a range of cloud types [Davis et al., 1999]. The fitting 0 coefficients for LWP are z0 = 0.49, z1 = 1.52. This field has similar intermittency characteristics as measurements of water vapor fluctuations in a cloud-free tropical boundary layer (represented by the curve labeled VAPOR) and is less intermittent than measurements of liquid water content fluctuations measured in horizontal flight legs in stratocumulus layers off of California (FIRE), the Azores (ASTEX), and the Tasman Sea (SOCEX) [Cho et al., 2001; Davis et al., 1999]. [65] The empirical relation (19) is limited to z( p) data that can be fit by a branch of a hyperbola. An alternative approach is to use forms of the z( p) function that are derived from random cascade models with specified generators. A principal advantage of such a fit is that the estimated parameters of the z( p) function can be used to simulate random fields with ensemble-averaged statistics that match those of the observations. Such simulations can be used, for instance, to study the realization-to-realization variability in the single-point Reference satellite Frisch [1995] this paper satellite wind tunnel wind tunnel ground-based radiometer aircraft aircraft aircraft aircraft satellite satellite Pierrehumbert [1996] Antonia et al. [1984] Ruiz-Chavarria et al. [1996] Davis et al. [1999] Cho et al. [2001] Davis et al. [1999] Davis et al. [1999] Davis et al. [1999] Lovejoy et al. [2001] this paper histograms for comparison with observed data. Furthermore, the form of the z( p) of certain cascade models can be more flexible than equation (19) permitting, for example, the description of z( p) which have dz/dp < 0, for some p, e.g., PCf and L01 in Figure 10. Curve L01 is derived from 909 infrared and visible satellite measurements spanning a wide variety of both cloud types and cloud fractions [Lovejoy et al., 2001]. Lovejoy et al. [2001] fit L01 using the ‘‘universal multifractal’’ model of Schertzer and Lovejoy [1991], which is a continuous cascade model with log Le´vy statistics. For this model, zð pÞ ¼ Hp À C1 ð pa À pÞ; aÀ1 a 6¼ 1 ð20Þ with parameters a, C1, and H. The parameter 0 < a < 2 is the Le´vy index, which characterizes the degree of multifractality (as a ! 0 a simple scaling model is approached, while a = 2 leads to a log normal model), C1 is the codimension of the mean of the field and characterizes the sparseness of the mean (the larger the value of C1, the more clustered the field), and the parameter H = z(1) is the Holder exponent, or Hurst’s roughness exponent, which characterizes the smoothness of the field (larger H corresponds to smoother fields, with larger jumps less likely between closely spaced points) [e.g., Marshak et al., 1997]. Following Lovejoy et al. [2001], we choose the refined similarity hypothesis parameter a = 1 [see also Schertzer et al., 1997; Frisch, 1995, p. 164]. [66] The fit to the FC data using this model yields parameter values a = 2, C1 = 0.029, and H = 0.354; the corresponding curve is labeled as FCf in Figure 10. The value a = 2 indicates that the statistics of the fully cloudy optical depth field are consistent with a log normal cascade model. Although the best fit was obtained for a = 2.0, reasonably good fits can also be obtained by setting a as low as 1.5, indicating a degree of degeneracy in the model for the p moments that are calculated from this data set. The parameter values of a = 1.91 and C1 = 0.077 measured by Lovejoy et al. [2001] also indicate near-log normal statistics, although the value of C1 is closer to the value measured for the partly cloudy data (see below). The curve L01 is drawn with these values of a and C1, and a representative H = 0.3; Lovejoy et al. [2001] observed that although the values of a and C1 did not have large scene-to-scene variability, the value of H fluctuated between 0.2 and 0.6. 13 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS [67] Both the partly cloudy PCf curve and the L01 curve have scaling exponents z( p) which decrease with p for approximately p > 3. Fields that have a decreasing z( p), for sufficiently small p, are extremely intermittent and contain singularities when the field has infinite resolution and is scaling down to the smallest scales. In practice, such fields will always contain some values which are very large relative to the mean. [68] Decreasing z( p) also indicates the possibility of a divergence of higher-order moments [Schertzer et al., 1997]. Such a divergence would imply that the increments of the field have algebraic tails (e.g., those of log Le´vy distributions), which are much fatter than, for instance, the exponential tails which are often used to describe turbulence statistics [e.g., Frisch, 1995, p. 128]. In such a case the random field is so intermittent that no single realization will contain the extreme values of the field that are necessary for an accurate approximation of the high order moments. Thus the approximation of the moments above some critical value pc will depend on the number of realizations used in the approximation. Furthermore, as more realizations are averaged, the estimate of the moments higher than a second critical moment pd will grow without bound (hence the term divergence of moments) [Schertzer et al., 1997]. [69] In a single realization the critical pc value is marked by a transition of z( p) to a linear asymptote of the form z( p) = C À g p, for p > pc. If many realizations are averaged, then the transition to linear behavior occurs at the moment, ps, where pc < ps < pd. This behavior, which indicates a divergence of moments, is seen in measurements of the structure functions of transverse vertical velocity fluctuations in the atmospheric boundary layer [Cho et al., 2001], as well as in the results of Lovejoy et al. [2001]. For curve PCf a transition appears to lie between p = 2 and p = 3, and thus moments p > 3 are not used in the approximation of the universal multifractal parameters. Estimates of a = 1.37, C1 = 0.089, and H = 0.18 are obtained for the partly cloudy data. The value of a = 1.37 indicates that a log normal model cannot adequately describe the statistics. However, the lack of data points used in the fit should be considered when making conclusions regarding these quantitative results. 7. Discussion [70] We have studied the spatial variability of marine boundary layer clouds using structure functions and scalar spectra. The results indicate that, although significant sceneto-scene variability is observed, the averaged two-point statistics for both the fully cloudy and partly cloudy data clearly satisfy a power law relation over a range of scales. [71] For each scene the second-order structure function and the scalar spectrum, which present complimentary views of the optical depth fluctuations, are computed and compared. In the fully cloudy scenes, both of these secondorder statistics display a break in scaling at a length scale of %200 m, consistent with the diffusion regime produced by photon smoothing of the radiance field. At separations between %200 m and 10 km the behavior varies widely between realizations, but the averaged statistics for the fully cloudy scenes show well-defined scaling, with the scale- D04104 invariant exponents, z2 and b, closely following the b = z2 + 1 relationship expected for scaling fields and having values similar to those predicted by KOC theory. Although we might expect the anomalous scaling of the radiance field to deviate from that produced by passive scalar transport in isotropic turbulence, we find, as did Arneodo et al. [1999], that there is close agreement between these cloud fluctuations and wind tunnel temperature measurements. [72] At separations larger than 10 km the variability for the fully cloudy scenes as measured by the isotropic secondorder structure function S2(r) and the isotropic scalar spectrum E(k) diverges, with S2(r) showing a break in scaling for nearly every realization, while E(k) does not. By investigating the integral transform between S2(r) and E(k), we find that this discrepancy is consistent with the difference with which the two statistics resolve a deviation from scaling over a restricted range of scales. An important consequence of this result is that the observed break in the structure function does not indicate the existence of an integral scale within the range of scales resolved by the data; i.e., there are long-range correlations over all these scales. We find no evidence of measurement effects that can cause this deviation, which supports the conclusion that it is of physical origin. These results do not contradict Lovejoy et al. [2001], who observed no break in scaling, because, given the resolution of their data, a restricted deviation from scaling at $10 km would be easily obscured. In addition, Lovejoy et al. [2001] averaged the statistics of a wide range of cloud types, while in this paper we have restricted our samples to marine boundary layer clouds. [73] From the integral transform we also find that if a break in the scaling exists in the scalar spectrum, the scaling of the structure function is only recovered asymptotically. Loosely, approximate scaling is recovered with a transition between scaling regimes (that spans one or more orders of magnitude). Furthermore, if a scale break is observed in the scalar spectrum at wave number k = k1, then the structure function will show a corresponding break in the small scale regime at approximately r = 1/(2 k1). [74] The average statistics of the B scenes also show welldefined scaling between approximately $1 and 20 km. For the B scenes the exponents z2 and b are significantly smaller than those of the A scenes, indicating the importance of considering the A scenes and B scenes as realizations from separate ensembles. Also, the exponents for the B scenes do not exactly satisfy the expected relationship, i.e., b 6¼ z2 + 1. This discrepancy may be explained by the differing ways in which the structure function and scalar spectrum respond to the cloud-clear transitions that are present in the partly cloudy scenes. [75] The comparison of the scaling of the higher-order structure functions (see Figure 10) shows that the intermittency is much larger for the B scenes than either the A scene measurements or in situ measurements of liquid water content or water vapor fluctuations. This is consistent with the structure of the broken clouds, in which portions of the field with zero value are randomly interspersed with relatively large activity (the cloud). Furthermore, for the broken liquid water clouds, the multifractal analysis of section 6 predicts a divergence of moments, which indicates that the PDF of the field increments is fat-tailed, as the log Le´vy distribution, and that no individual scene can capture all of 14 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS the fluctuations, i.e., there is significant realization dependence of the measured statistics. [76] As mentioned in section 1, the A and B ensembles also have differing single-point histograms: all B scene histograms lack a mode value, so that the smallest values of optical depth are most probable in these layers. Both models and observations indicate that underlying this modeless PDF in liquid water clouds is a distribution of excess supersaturation fluctuations, e, which when thresholded produces the observed distribution of integrated liquid water path [Sommeria and Deardorff, 1977; Bechtold and Siebesma, 1998]. Much current work on subgrid-scale parameterizations for both ice and water clouds is focused on finding relationships between quantities such as the optical depth and cloud fraction at resolved scales and the PDF of the subgrid-scale fluctuations [Barker and Fu, 2000; Carlin et al., 2002; Smith and Del Genio, 2002]. The results presented here suggest that it is particularly difficult to determine the underlying single-point statistics for these partly cloudy scenes, both because the thresholding process eliminates most of the fluctuation distribution from the remotely sensed sample of the cloud and because the essentially inhomogeneous field exhibits large realizationdependent variability. The high degree of intermittency observed in these partly cloudy fields is also an issue for efforts to characterize the impact of unresolved cloud on satellite-derived estimates of albedo, cloud fraction, and cloud liquid water path [e.g., Oreopoulos and Davies, 1998; Wood and Hartmann, 2001]. [77] Given the assumption that a field contains no fluctuations at scales greater than some ‘‘integral’’ scale and the assumption that all realizations are derived from a single ensemble, all realizations will have the same underlying single-point PDFs. In theory, the parameters of this PDF can be found by an ensemble average; however, in practice, owing to the presence of long-range correlations, the PDF parameters approximated by spatial averages over a given scale, smaller than the integral scale, will be random variables rather than a single constant, i.e., the single-point histograms of the fields, and thus the estimated parameters of the model PDF will exhibit significant scene-to-scene variations. Indeed, this realization dependence of the histogram is observed for the data presented here [see Barker et al., 1996, Table 2]. Furthermore, accurate predictions of the statistics of a single subsample of a cloud field cannot be found deterministically from, e.g., the mean of the sampled data, if the sample size is smaller than the scale of the largest fluctuations. [78] In section 6 we fit the observed structure function exponents to parameters of multifractal cascade models. Given the model parameters H, C1, and a, it is possible to generate scaling random fields that can be used to establish, for example, confidence intervals about the PDF parameter estimates needed for subgrid-scale parameterizations [e.g., Wilson et al., 1991; Gupta and Waymire, 1993]. In addition, for a given set of model parameters, we can obtain the ensemble-averaged probability that the optical depth will exceed a particular threshold, i.e., tail probability estimates [e.g., Lovejoy and Schertzer, 1992; Gupta and Waymire, 1993]. [79] More data are required for a better understanding of the variability of the full and broken cloud layers at large D04104 scales. We plan to extend this analysis to satellite observations that span larger spatial domains and that are stratified over narrower ranges of cloud fraction. Colocated comparisons of the two-point statistics of broken clouds with subcloud aircraft measurements would be a particularly useful addition to the results of section 6. [80] Acknowledgments. This work has benefited from conversations with Chris Jeffery, the helpful comments of reviewers, and the programming assistance of Charles Leung and Chris Seymour. It was supported through funding of the Modeling of Clouds and Climate Proposal by the Canadian Foundation for Climate and Atmospheric Sciences, the Meteorological Service of Canada, and the Natural Sciences and Engineering Research Council. References Antonia, R. A., E. J. Hopfinger, Y. Gagne, and F. Anselmet (1984), Temperature fluctuations in turbulent shear flows, Phys. Rev. A, 30, 2704 – 2707. Arneodo, A., N. Decoster, and S. G. Roux (1999), Intermittency, log-normal statistics and multifractal cascade process in high-resolution satellite images of cloud structure, Phys. Rev. Lett., 83(6), 1255 – 1258. Barker, H. W. (1996), A parameterization for computing grid-averaged solar fluxes for inhomogeneous marine boundary layer clouds. Part I: Methodology and homogeneous biases, J. Atmos. Sci., 53, 2289 – 2302. Barker, H. W., and J. A. Davies (1992), Cumulus cloud radiative properties and the characteristics of satellite radiance wavenumber spectra, Remote Sens. Environ., 42, 51 – 64. Barker, H. W., and Q. Fu (2000), Assessment and optimization of the gamma-weighted two-stream approximation, J. Atmos. Sci., 57, 1181 – 1188. Barker, H. W., B. Wielicki, and L. Parker (1996), A parameterization for computing grid-averaged solar fluxes for inhomogeneous marine boundary layer clouds. part II: Validation using satellite data, J. Atmos. Sci., 53, 2304 – 2316. Bechtold, P., and P. Siebesma (1998), Organization and representation of boundary layer clouds, J. Atmos. Sci., 55, 888 – 895. Bony, S., and K. A. Emanuel (2001), A parameterization of the cloudiness associated with cumulus convection: Evaluation using TOGA COARE data, J. Atmos. Sci., 58, 3158 – 3183. Cahalan, R. F., and J. B. Snider (1989), Marine stratocumulus structure, Remote Sens. Environ., 28, 95 – 107. Carlin, B., Q. Fu, U. Lohmann, G. G. Mace, K. Sassen, and J. M. Comstock (2002), High cloud horizontal inhomogeneity and solar albedo bias, J. Clim., 15, 2321 – 2339. Cho, J. Y. N., B. E. Anderson, J. D. W. Barrick, and K. L. Thornhill (2001), Aircraft observations of boundary layer turbulence: Intermittency and the cascade of energy and passive scalar variance, J. Geophys. Res., 106(23), 32,469 – 32,479. Considine, G., J. A. Curry, and B. Wielicki (1997), Modeling cloud fraction and horizontal variability in marine boundary layer clouds, J. Geophys. Res., 102(D12), 13,517 – 13,525. Cusack, S., J. M. Edwards, and R. Kershaw (1999), Estimating the subgrid variance of saturation, and its parameterization for use in a GCM cloud scheme, Q. J. R. Meteorol. Soc., 125, 3057 – 3076. Davis, A., A. Marshak, W. Wiscombe, and R. Cahalan (1996), Scale invariance of liquid water distribution in marine stratocumulus. Part I: Spectral properties and stationarity issues, J. Atmos. Sci., 53, 1538 – 1558. Davis, A., A. Marshak, R. Cahalan, and W. Wiscombe (1997), The Landsat scale break in stratocumulus as a three-dimensional radiative transfer effect: Implications for cloud remote sensing, J. Atmos. Sci., 54, 241 – 260. Davis, A. B., A. Marshak, H. Gerber, and W. J. Wiscombe (1999), Horizontal structure of marine boundary layer clouds from centimeter to kilometer scales, J. Geophys. Res., 104(D6), 6123 – 6144. Frisch, U. (1995), Turbulence, Cambridge Univ. Press, New York. Gollmer, S., M. Harshvardhan, R. F. Cahalan, and J. B. Snider (1995), Windowed and wavelet analysis of marine stratocumulus cloud inhomogeneity, J. Atmos. Sci., 52(16), 3013 – 3030. Gupta, V. K., and E. C. Waymire (1993), A statistical analysis of mesoscale rainfall as a random cascade, J. Appl. Meteorol., 32, 251 – 267. Harshvardhan, B., A. Wielicki, and K. M. Ginger (1994), The interpretation of remotely sensed cloud properties from a model parameterization perspective, J. Clim., 7(12), 1987 – 1998. Jeffery, C. A., and P. H. Austin (2003), Unified treatment of thermodynamic and optical variability in a simple model of unresolved low clouds, J. Atmos. Sci., 60, 1621 – 1631. 15 of 16 D04104 LEWIS ET AL.: SPATIAL STATISTICS Larson, V. E., R. Wood, P. R. Field, J.-C. Golaz, T. H. V. Haar, and W. R. Cotton (2001), Systematic biases in the microphysics and thermodynamics of numerical models that ignore subgrid-scale variability, J. Atmos. Sci., 58, 1117 – 1128. Lewellen, W. W., and S. Yoh (1993), Binormal model of ensemble partial cloudiness, J. Atmos. Sci., 50, 1228 – 1237. Lewis, G. M., and P. H. Austin (2002), An iterative method for generating log-normal scaling simulations, in Proceedings of the 11th Conference on Atmospheric Radiation, Am. Meteorol. Soc., Boston, Mass. (Available as http://ams.confex.com/ams/11AR11CP/11cldphy/abstracts/42772.htm) Lovejoy, S., and D. Schertzer (1992), Multifactals and rain, in New Uncertainty Concepts in Hydrology and Water Resources, edited by Z. W. Kunzewicz, Cambridge Univ. Press, New York. Lovejoy, S., D. Schertzer, P. Silas, Y. Tessier, and D. Lavalle´e (1993), The unified scaling model of atmospheric dynamics and systematic analysis of scale invariance in cloud radiances, Ann. Geophys., 11, 119 – 127. Lovejoy, S., D. Schertzer, and J. D. Stanway (2001), Direct evidence of multifractal atmospheric cascades from planetary scales down to 1 km, Phys. Rev. Lett., 86(22), 5200 – 5203. Marshak, A., A. Davis, W. Wiscombe, and R. Cahalan (1997), Scale invariance of liquid water distribution in marine stratocumulus. Part II: Multifractal properties and intermittency issues, J. Atmos. Sci., 54, 1423 – 1444. Monin, A. S., and A. M. Yaglom (1975), Statistical Fluid Mechanics: Mechanics of Turbulence, vol. II, MIT Press, Cambridge, Mass. Oreopoulos, L., and R. Davies (1998), Plane parallel albedo biases from satellite observations. Part I: Dependence on resolution and other factors, J. Clim., 11, 919 – 932. Oreopoulos, L., A. Marshak, R. F. Cahalan, and G. Wen (2000), Cloud three-dimensional effects evidenced in Landsat spatial power spectra and autocorrelation functions, J. Geophys. Res., 105, 14,777 – 14,788. Pierrehumbert, R. T. (1996), Anomalous scaling of high cloud variability in the tropical Pacific, Geophys. Res. Lett., 23(10), 1095 – 1098. Pincus, R., and S. A. Klein (2000), Unresolved spatial variability and microphysical process rates in large-scale models, J. Geophys. Res., 105(22), 27,059 – 27,065. Ruiz-Chavarria, G., C. Baudet, and S. Ciliberto (1996), Scaling laws and dissipation scale of a passive scalar in fully developed turbulence, Phys. D, 99, 369 – 380. Schertzer, D., and S. Lovejoy (1991), Non-linear geodynamical variability: Multiple singularities, universality and observables, in Non-linear Variability in Geophysics: Scaling and Fractals, edited by D. Schertzer and S. Lovejoy, pp. 41 – 82, Kluwer Acad., Norwell, Mass. D04104 Schertzer, D., S. Lovejoy, F. Schmitt, Y. Chigirinskaya, and D. Marsan (1997), Multifractal cascade dynamics and turbulent intermittency, Fractals, 5, 427 – 471. Smith, R. N. B. (1990), A scheme for prediciting layer clouds and their water content in a general circulation model, Q. J. R. Meteorol. Soc., 116, 435 – 460. Smith, S. A., and A. D. Del Genio (2002), A simple conceptual model of cirrus horizontal inhomogeneity and cloud fraction, Q. J. R. Meteorol. Soc., 128, 149 – 171. Sommeria, G., and J. W. Deardorff (1977), Subgrid-scale condensation in models of nonprecipitating clouds, J. Atmos. Sci., 34, 344 – 355. Tompkins, A. M. (2002), A prognostic parameterization for the subgridscale variability of water vapor and clouds in large-scale models and its use to diagnose cloud cover, J. Atmos. Sci., 58, 1917 – 1942. Vainshtein, S. I., K. R. Sreenivasan, R. T. Pierrehumbert, V. Kashyap, and A. Juneja (1994), Scaling exponents for turbulence and other random processes and their relationships with multifractal structure, Phys. Rev. E, 50, 1823 – 1835. von Savigny, C., A. B. Davis, O. Funk, and K. Pfeilsticker (2002), Timeseries of zenith radiance and surface flux under cloudy skies: Radiative smoothing, optical thickness retrievals and large-scale stationarity, Geophys. Res. Lett., 29(17), 1825, doi:10.1029/2001GL014153. Wilson, J., D. Schertzer, and S. Lovejoy (1991), Continuous multiplicative cascade models of rain and clouds, in Non-Linear Variability in Geophysics: Scaling and Fractals, edited by D. Schertzer and S. Lovejoy, pp. 185 – 207, Kluwer Acad., Norwell, Mass. Wood, R., and D. L. Hartmann (2001), Investigations of liquid water path spatial variability using MODIS, in Proceedings of the 11th Conference on Atmospheric Radiation, Am. Meteorol. Soc., Boston, Mass. (Available as http://ams.confex.com/ams/11satellite/11satellite/abstracts/23925. htm). ÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀÀ P. H. Austin, Department of Earth and Ocean Sciences, University of British Columbia, 6339 Stores Road, Vancouver, BC V6T 1Z4, Canada. (paustin@eos.ubc.ca) G. M. Lewis, Faculty of Science, University of Ontario Institute of Technology, 2000 Simcoe Street North, Oshawa, Ontario L1H 7K4, Canada. (greg.lewis@uoit.ca) M. Szczodrak, Rosentiel School of Marine and Atmospheric Science, University of Miami, Miami, FL 33149-1098, USA. (goshka@rrsl.rsmas. miami.edu) 16 of 16
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Faculty Research and Publications /
- Spatial statistics of marine boundary layer clouds
Open Collections
UBC Faculty Research and Publications
Spatial statistics of marine boundary layer clouds Lewis, Gregory M.; Austin, Philip H.; Szczodrak, Malgorzata Feb 25, 2004
pdf
Page Metadata
Item Metadata
Title | Spatial statistics of marine boundary layer clouds |
Creator |
Lewis, Gregory M. Austin, Philip H. Szczodrak, Malgorzata |
Publisher | American Geophysical Union |
Date Issued | 2004-02-25 |
Description | An analysis is presented of the structure functions and scalar spectra for 25 satellite-derived marine stratocumulus cloud optical depth fields. The scenes, which cover a horizontal domain of 58 × 58 km at a resolution of 28.5 m, are partitioned into two ensembles on the basis of cloud fraction. For the fully cloudy scenes, although there is wide scene-to-scene variability, both the average isotropic scalar spectrum and the average isotropic second-order structure function exhibit power law behavior over approximately two decades, with scale-invariant exponents equal to those expected for inertial-subrange passive tracer fluctuations. Higher-order structure functions show anomalous scaling that closely matches that observed for wind tunnel temperature fluctuations and for other fully cloudy observations. The partly cloudy scenes, while scaling, show different behavior. The average isotropic second-order structure function and average isotropic scalar spectrum have scale-invariant exponents that are significantly smaller than those of the fully cloudy scenes, and the analysis of the higher-order structure functions indicates that the field has much more intermittent fluctuations than the fully cloudy scenes. Fits to random cascade models for the fully cloudy scenes show that the increment statistics are consistent with an underlying log normal distribution. For the partly cloudy scenes a divergence of higher-order moments is predicted, indicating that the field fluctuations are necessarily derived from fat-tailed distributions and that there will be significant realization dependence of the measured statistics. In addition, the presence of long-range correlations in all the data predicts that single-point histograms of the field values will have significant scene-to-scene variability, or equivalently, the use of spatial averages in the approximation of the parameters of the single-point probability density function of the field will result in random fluctuations of the estimated parameters. An edited version of this paper was published by AGU. Copyright 2004 American Geophysical Union. |
Genre |
Article |
Type |
Text |
Language | eng |
Date Available | 2016-11-07 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
DOI | 10.14288/1.0041775 |
URI | http://hdl.handle.net/2429/32765 |
Affiliation |
Science, Faculty of Earth and Ocean Sciences, Department of |
Citation | Lewis, Gregory M. , Philip. H. Austin, Szczodrak, Malgorzata (2004) Spatial statistics of marine boundary layer clouds. Journal of Geophysical Research Atmospheres 109, D04104, |
Publisher DOI | 10.1029/2003JD003742. |
Peer Review Status | Reviewed |
Scholarly Level | Faculty |
Copyright Holder | Ausin, Philip H. |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 52383-Austin_AGU_2004_JD003742.pdf [ 1.59MB ]
- Metadata
- JSON: 52383-1.0041775.json
- JSON-LD: 52383-1.0041775-ld.json
- RDF/XML (Pretty): 52383-1.0041775-rdf.xml
- RDF/JSON: 52383-1.0041775-rdf.json
- Turtle: 52383-1.0041775-turtle.txt
- N-Triples: 52383-1.0041775-rdf-ntriples.txt
- Original Record: 52383-1.0041775-source.json
- Full Text
- 52383-1.0041775-fulltext.txt
- Citation
- 52383-1.0041775.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.52383.1-0041775/manifest