Remotely sensed radiances from planetary surfaces or atmospheres reveal fields of immense complexity with structures typically spanning the range of scales from planetary to submillimetric: 10 or more orders of magnitude. The natural framework for analyzing, modeling, and indeed understanding such hierarchies of structures within structures is scale invariance: fractal sets and multifractal fields. A little over 15 years ago, several colleagues and I gave a short review of some of the theories relevant to remote sensing, including several examples (Pecknold et al. 1997; see also Lovejoy et al. 2001a, especially the discussion of correlated scaling processes). This chapter is an update discussing some of the advances that have occurred since then.
Remotely sensed radiances from planetary surfaces or atmospheres reveal fields of immense complexity with structures typically spanning the range of scales from planetary to submillimetric: 10 or more orders of magnitude. The natural framework for analyzing, modeling, and indeed understanding such hierarchies of structures within structures is scale invariance: fractal sets and multifractal fields. A little over 15 years ago, several colleagues and I gave a short review of some of the theories relevant to remote sensing, including several examples (Pecknold et al. 1997; see also Lovejoy et al. 2001a, especially the discussion of correlated scaling processes). This chapter is an update discussing some of the advances that have occurred since then.
Resolutions have improved and channels have multiplied. Today, images are often collected at regular intervals so that today remotely sensed data are very much space–time products. As the quantity of remotely sensed data and the range of scales that they span has grown, the data are increasingly of global extent. In some cases such as topography or cloud imagery, composite data analyses spanning up to eight orders of magnitude are possible. In parallel with that, global numerical models may also show widerange scaling; this includes weather prediction models and reanalysis products. They show that the standard atmospheric state variables are scaling up to planetary scales (typically 5,000 km or larger). This is fortunate because it allows them to be compatible with the scaling of the fields as determined by in situ and remotely sensed measurements. Other significant developments include improvements in some of the older data analysis techniques (notably trace moments, in particular their use to explicitly determine the outer scales) but also in the development of new techniques (Haar fluctuations).
However, many phenomena—especially atmospheric—evolve rapidly enough in time that it is important to understand the space–time behavior. With respect to both atmospheric and oceanic phenomena, there has been important progress, notably in clarifying the notion of weather and climate with the realization that there is a third regime (macroweather) that is intermediate between the two (Lovejoy 2013). A consequence is that it turns out that the climate is not “what you expect”; it is rather macroweather. Whereas in the weather regime (up to 5–10 days, the lifetime of planetary structures), fluctuations on average increase with scale (the exponent H defined below is >0), in macroweather, on the contrary, they decrease with scale (H < 0). In this regime, successive fluctuations tend to cancel each other out: “Macroweather is what you expect.” The climate regime starts at much lower frequencies (≈30 years in the industrial epoch) and again, H > 0. In the section “Space–Time Scaling: Example of MTSAT,” we therefore discuss the space–time development of geostationary thermal IR fields providing a theoretical framework for understanding atmospheric motion vectors (AMVs).
Scale invariance is a symmetry such that when one changes from one scale to another, some aspect is unchanged (conserved). In scaleinvariant systems, fluctuations ΔI over distances Δx have power law dependences: ΔI(Δx) ≈ φ Δx^{H}, where φ is the driving flux whose statistical mean [φ] does not depend on scale. The exponent H characterizing the fluctuations is scale invariant, whereas the fluctuations themselves are power law functions of scale; they are said to be “scaling.” The corresponding (Fourier and power) spectra are of the form E(k) ≈ k^{ −β} where k is a wave number (inverse distance scale) and β is the “spectral exponent.”
Spectra have several useful features: they are traditional, familiar, and applicable to any geophysical signal. They are also very sensitive to breaks in scaling and—when these are due to quasiperiodic nonstationarities (such as daily of annual cycles)—of easily separating these from the otherwise scaling “background.” The main disadvantages are that (1) their interpretation is not as straightforward as for (real space) fluctuations; (2) they are not well adapted to situations with missing data; and (3) they only characterize the secondorder statistical moments so that—unless the signal is quasiGaussian—it only gives a very partial characterization of the statistics. In this section, we give a quick tour of some of examples; in the section “Multifractals, Structure Functions,” we discuss alternative real space analysis techniques.
Let us recall the basics. We will be dealing mostly with twodimensional (2D) remotely sensed radiance fields I(r) (r = (x,y) is a position vector), so that we can define the 2D (Fourier) spectrum P(k) by the following:
where k is the wave vector dual to r, Ĩ(k) is the Fourier transform of I, and δ is the Dirac delta function. For finite data, the delta function is replaced by a finite difference approximation such that at k = −k′, it is equal to the number of degrees of freedom of the system (the number of pixels) N. For a real signal I, we have Ĩ* (k) = Ĩ (−k) (the asterisk [*] indicates a complex conjugate) so that:
where the constant of proportionality is N. If the system is statistically isotropic, then P only depends on the vector norm of k: P(k) = P(k); k = k. Now perform an isotropic Fourier space “zoom” k → λk (i.e., a standard “blowup”) by a factor of λ > 1 so that in physical space there is an inverse blowup: x → λ^{−1}x. If the system is “selfsimilar,” that is, if it is both isotropic and scaling, then the condition that the smaller scales are related to the larger scales without reference to any characteristic size (i.e., that it is scaling) is that the spectra follow power law relations between large wave numbers λk and smaller ones k:
that is, that the form of P is independent of scale. Equation 3.3 is satisfied by the following scaling law for P:
We can now obtain the power spectrum E(k) (with k = k) by integrating over all the directions:
In d = 2, the region of integration is the annuli between the radii k and k + dk. Therefore, if the process is isotropic in 2D, E(k) = 2π kP(k). In terms of data analysis, where one has a finite rather than infinite sample size, this angle integration is advantageous because it reduces the noise. In the following examples, we therefore take the power law dependence of the spectrumml:
as evidence for the scaling of the field f, and the exponent β being the “spectral slope.” Note that in some areas of geophysics, angle averages (P(k)) are used rather than angle integrals (E(k)). This has the disadvantage that the resulting spectral exponents will depend on the dimension of space so that, for example, onedimensional sections will have exponents that differ by 1 from 2D sections. In contrast the angle integrations used here yield the same exponent β in spaces of any dimension (i.e., β is independent of the dimension of space, whereas s is dependent).
When permitted by the data, using the angle integrals over the 2D spectrum P(k_{x}, k_{y}) is advantageous because it reduces the spectrum to a function of a single variable (the modulus of the wave vector) while simultaneously improving the statistics. Often the data are easily amenable to integrating over angles—for example geostationary satellites; others such as orbiting satellites commonly have swaths of limited width while being essentially 20,000 km long (i.e., half a circumference—no two points on the Earth can be further apart) so that they are more one dimensional.
Alternatively, it is possible to estimate the 1D spectra E(k_{x}), E(k_{y}) by integrating out the conjugate coordinates:
and for isotropic statistics, E_{x}(k) = E_{y}(k) = E(k) (to within constant factors). Of course, if the field is not isotropic, for example if it is scaling but with different exponents in different orthogonal directions, $\begin{array}{c}{E}_{x}\left({k}_{x}\right)\propto \end{array}{k}_{x}^{{\beta}_{x}};{E}_{y}\left({k}_{y}\right)\propto {k}_{y}^{{\beta}_{y}}$
with β_{x} ≠ β_{y}, then the isotropic spectrum will in fact display a break (roughly) at the wave number k that satisfies E_{x}(k) = E_{y}(k) with one exponent dominating the behavior at the high frequencies and the other at the low frequencies. This is discussed further in the section on space–time scaling in the context of the space–time analysis of satellite data.Over the last few years, a number of global scale satellite radiances have impressively demonstrated the global scale extent of the scaling. For example, Figure 3.1a and b, taken from over 1,000 orbits of the Tropical Rainfall Measuring Mission (TRMM) satellite, shows that visible, near, and thermal infrared as well as passive microwave radiances have nearly perfect scaling up to the largest scales. This is also demonstrated by Multifunctional Transport Satellite (MTSAT, also known as “Himawari”) geostationary imagery (Figure 3.1c), which shows the results of spectral analysis of a large (roughly 1,000^{3} points) data set in (x,y,t) space showing that the (1D) temporal and (horizontal) spatial statistics are nearly identical. This is an isotropic space–time symmetry that we discuss further in “Space–Time Scaling: Example of MTSAT” (we also analyze and display spectra of the same data over various 2D subspaces). Moving to the opposite extreme of small scales, one can evaluate the spectra at small (“cloud”) scales, this time looking upwards with a handheld largeformat camera (Figure 3.1d), which displays excellent scaling (of the angle integrated) spectrum of individual clouds at visible wavelengths. The theory and numerical modeling of radiative transfer in scaling clouds are reviewed in the section of radiative transfer.
Figure 3.1 (a) Spectra from ≈ 1,000 orbits of the visible and infrared scanner (VIRS) instrument on the Tropical Rainfall Measuring Mission (TRMM) satellite, channels 1–5 (at wavelengths of 0.630, 1.60, 3.75, 10.8, and 12.0 µm, from top to bottom, displaced in the vertical for clarity). The data are for the period January–March 1998 and have nominal resolutions of 2.2 km. The straight regression lines have spectral exponents β = 1.35, 1.29, 1.41, 1.47, and 1.49, respectively, close to the value β = 1.53 corresponding to the spectrum of passive scalars (=5/3 minus intermittency corrections). The units are such that k = 1 is the wave number corresponding to the size of the planet (20,000 km)−1. Channels 1 and 2 are reflected solar radiation so that only the 15,600km sections of orbits with maximum solar radiation were used. The high wave number falloff is due to the finite resolution of the instruments. To understand the figure we note that VIRS Bands 1 and 2 are essentially reflected sunlight (with very little emission and absorption), so that for thin clouds, the signal comes from variations in the surface albedo (influenced by the topography and other factors), whereas for thicker clouds it comes from nearer the cloud top via (multiple) geometric and Mie scattering. As the wavelength increases into the thermal IR, the radiances are increasingly due to black body emission and absorption with very little multiple scattering. Whereas at the visible wavelengths we would expect the signal to be influenced by the statistics of cloud liquid water density, for the thermal IR wavelengths it would rather be dominated by the statistics of temperature variations—themselves also close to those of passive scalars (Adapted from Lovejoy, S., et al., Q. J. Roy. Meteor. Soc., 134, 277–300, 2008). (b) Spectra of radiances from the Thematic Microwave Imager from the TRMM satellite, ≈1,000 orbits from January through March 1998. From bottom to top, the data are from Channels 1, 3, 5, 6, and 8 (vertical polarizations of 2.8, 1.55, 1.41, 0.81, and 0.351 cm) with spectral exponents β = 1.68, 1.65, 1.75, 1.65, and 1.46, respectively, at resolutions of 117, 65, 26, 26, and 13 km (hence the high wave number cutoffs); each are separated by one order of magnitude for clarity. To understand these thermal microwave results, recall that they have contributions from surface reflectance, water vapor, clouds, and rain. Because the particles are smaller than the wavelengths, this is the Rayleigh scattering regime and as the wavelength increases from 3.5 to 2.8 cm the emissivity/absorptivity due to cloud and precipitation decreases so that more and more of the signal originates in the lower reaches of clouds and underlying surface. Moreover, the ratio of scattering to absorption increases with increasing wavelength so that at 2.8 cm multiple scattering can be important in raining regions. The overall result is that the horizontal gradients—which influence the spectrum—increasingly reflect large internal liquid water gradients. (c) Onedimensional spectra of Multifunctional Transport Satellite (MTSAT) thermal IR radiances; the Smith product was developed with similar IR satellite radiance fields. In black: the theoretical spectrum using parameters estimated by regression from the 3D generalization fo equation 3.7 to include the frequency Greek omega. and taking into account the finite space–time sampling volume. The spectra are Ex(kx)≈kx−βx,Ey(ky)≈ky−βy,Et(ω)≈ω−βt with βx ≈ βy ≈ βt ≈ 1.4±0.1. Other parameters are Lw ~ 20,000 km; τw ~ 20±1 days; s ≈ 3.4±0.1. The straight line is a reference line with slope −1.5 (blue). Pink is the zonal spectrum; orange is the meridional spectrum; blue (with the diurnal spike and harmonic prominent) is the temporal spectrum. See the section titled “Space–Time Scaling: Example of MTSAT” and Figure 3.18 for further analysis and discussion. (Reproduced from Pinel, J., et al., Atmos. Res., 140–141, 95–114, 2014. With permission.) (d) The spectra of the 19 (of 38) highest resolution clouds analyzed with a spectral slope β ≈ 2 shown by the reference lines. (Reproduced from Sachs, D., et al., Fractals, 10, 253–265, 2002. With Permission.)
The atmospheric energy budget is essentially determined by the incoming fluxes at visible wavelengths; and the outgoing fluxes at infrared wavelengths, the basic sources and sinks do not introduce characteristic scales. The lower boundary conditions (e.g., the topography and ocean surface) are also scaling (see below), as are the atmospheric dynamical equations. We may therefore expect the state variables (wind, temperature, pressure, humidity) to also be scaling. This is directly demonstrated in Figure 3.2, which shows the spectra of the outputs of atmospheric reanalyses. Reanalyses are datamodel hybrid products in which data of all kinds—increasingly including satellite data—are assimilated into a numerical weather model, which constrains them to satisfy the equations of the atmosphere as embodied in the numerical model.
Figure 3.2 Comparisons of the spectra of different atmospheric fields from the European Centre for MediumRange Weather Forecasts (ECMWF) interim reanalysis. Top (red) is the geopotential (β = 3.35); second from the top (green) is the zonal wind (β = 2.40); third from the top (cyan) is the meridional wind (β = 2.40); fourth from the top (blue) is the temperature (β = 2.40); fifth from the top (orange) is the vertical wind (β = 0.4); at the bottom (purple) is the specific humidity (β = 1.6). All are at 700 mb (roughly 3 km altitude) and between ±45° latitude, every day in 2006 at GMT. The scale at the far left corresponds to 20,000 km in the east–west direction, at the far right to 660 km. Note that for these 2D spectra, Gaussian white noise would yield β = –1 (i.e., a positive slope =+1). (Reproduced from Lovejoy, S., and D. Schertzer, J. Geophys. Res., 116, 2011. With permission.)
The topography is of prime importance for several areas of Earth science, notably as the lower boundary condition for the atmosphere and climate, and also for surface hydrology, oceanography, and hence for hydrosphere–atmosphere interactions. The issue of scaling in topography has an even longer history than it does in atmospheric science, going back over 100 years to when Perrin (1913) eloquently argued that the coast of Brittany was nondifferentiable. Later, Steinhaus (1954) expounded on the nonintegrability of the river Vistula, Richardson (1961) quantified both aspects using scaling exponents, and Mandelbrot (1967) interpreted the exponents in terms of fractal dimensions. Indeed, scaling in the Earth’s surface is so prevalent that there are entire scientific specializations such as river hydrology and geomorphology that abound in scaling laws of all types (for a review see RodriguezIturbe and Rinaldo 1997; see Tchiguirinskaia et al. 2000, for a comparison of multifractal and fractal analysis of basins) and that virtually require the topography to be scaling.
Ever since the pioneering power spectrum of VenigMeinesz (1951) with β ≈ 2, scaling spectra of topography from various regions have routinely been reported with fairly similar exponents (Balmino et al. 1973; Bell 1975 [also with β ≈ 2]; Berkson and Matthews 1983 [β ≈ 1.6−1.8]; Fox and Hayes 1985 [β ≈ 2.5]; Gilbert 1989 [β ≈ 2.1−2.3]; Balmino 1993 [β ≈ 2]; Lavallée et al. 1993; and Gagnon et al. 2006).
Figure 3.3 A log–log plot of the spectral power as a function of wave number for four digital elevation models. From right to left: Lower Saxony, with trees (top), without trees (bottom); United States (in gray), GTOPO30 and ETOPO5. A reference line of slope –2.10 is shown for comparison. The small arrows show the frequency at which the spectra are not well estimated due to their limited dynamical range (for this and scaledependent corrections; Equation 3.8). (Reproduced from Gagnon, J. S., et al., Nonlin. Proc. Geophys., 13, 541–570, 2006. With permission.)
Figure 3.3 shows the global scale spectrum of the ETOPO5 data set (a global gridded topography data set including bathymetry at a 5′ arc, i.e., about 10 km), along with those of other higher resolution but regional digital elevation models (DEMs, i.e., gridded topographic maps), land only. These include GTOPO30 (the continental United States at ≈1 km) as well as two other DEMs: the United States at 90m resolution and part of Saxony in Germany at 50cm resolution. Overall, the spectrum follows a scaling form with β ≈ 2.1 down to at least ≈40 m in scale.
The remarkable thing about the spectra is that the only obvious breaks are near the high wave number (smallscale) end of each data set. In Gagnon et al. (2006), it is theoretically shown that for wave numbers higher than the arrows, the data are corrupted by inadequate dynamical ranges (i.e., the ratio of the largest to the smallest resolvable changes in altitude). The basic idea is straightforward: if the dynamic range is small, then there will be large areas that are nominally at the same altitude, and this leads to spuriously smooth fields. For example, the curve in Figure 3.3 with the largest break in scaling was the DEM at 90m spatial resolution, which only had a 1m altitude resolution. This implies that huge swathes of the country had nominally zero gradients and hence overly smooth spectra. To quantify this, denote the minimum and maximum heights by h_{min} and h_{max}. If they are measured in nondimensional digital counts, then for a spectrum E(k) ≈ k^{−β}, we have the following:
The wave numbers at which this formula predicts that the spectra start to be corrupted are shown with arrows in Figure 3.3. This approximate formula works particularly well for the spectrum of the continental United States at 90m resolution, where it explains the drop in the high frequencies. In fact, the problem of insufficient dynamical range can probably explain many of the scale breaks seen in the literature that are interpreted as characteristic scales of the process. A related cause of spurious high wave number spectral falloffs occurs when data are sampled at a rate higher than the intrinsic sensor resolution (oversampling).
The ocean surface is particularly important due to its exchanges with the atmosphere. Figure 3.4 shows a particularly striking widerange scaling result: a swath over 200 km long at 7m resolution over the St. Lawrence Estuary in eight different narrow visible wavelength channels from the airborne multidetector electrooptical imaging sensor (MEIS). The use of different channels allows one to determine “ocean color,” which itself can be used as a proxy for phytoplankton concentration. For example, the channels fourth and eighth from the top in the figure exhibit nearly perfect scaling over the entire range; these are the channels that are insensitive to the presence of chlorophyll; they give us an indication that over the corresponding range ocean turbulence itself is scaling. In comparison, other channels show a break in the neighborhood of ≈200 m in scale; these are sensitive to phytoplankton. The latter are “active scalars” undergoing both exponential growth phases (“blooms”) as well as being victim to grazing by zooplankton; in Lovejoy et al. (2000), a turbulence theory is developed to explain the break with a zooplankton grazing mechanism.
Figure 3.4 The ocean, Channels 1–8 offset for clarity, eight visible channels characterizing ocean color, 210kmlong swath, 28,500 × 1,024 pixels, 7m resolution. The extreme high wave number is (14 m)−1. (Adapted from Lovejoy, S., et al., Inter. J. Remote Sensing, 22, 1191–1234, 2001a.)
Another important ocean surface field that has been found to be scaling over various ranges is the sea surface temperature (SST). Relevant studies include in situ results such as those of McLeish (1970) (Eulerian, β_{T} ≈ 5/3), Seuront et al. (1996) (Lagrangian, β_{T} ≈ 2), and Lovejoy et al. (2000) (towed instruments, β_{T} ≈ 1.63). Remote sensing using thermal IR images first from aircraft (Saunders 1972) and then from satellites (Deschamps et al. 1981, 1984; Park and Chung 1999) yields, respectively, β_{T} ≈ 5/3, ≈2.2, ≈1.9, ≈2, ≈1.87 ± 0.25 out to distances of ≈100 km. At larger scales (out to at least ≈500 km), Burgert and Hsieh (1989) found β_{T} ≈ 2.1 from “cloud free” satellite data. The satellite data—even if nominally “cloud free”—are somewhat smoothed by atmospheric effects; hence their β_{T} values are probably slightly too high. Monthly averaged in situ SST data were analyzed by Lovejoy and Schertzer (2013), and the literature was reviewed (see also Table 3.1). They concluded that the scaling with β_{T} ≈ 1.8 really does continue up to scales of 5,000 km or more, a conclusion that is bolstered by the corresponding spectral and cascade analyses.
C_{1} 
α 
H 
β 
L_{eff} 
Range of scales^{j} 


Remote Sensing of the Surface 

Topography^{a} (Earth) 
Altitude 
0.12 
1.8 
0.7 
2.1 
20,000 
40 m–20,000 km 
Topography^{b} (Mars) 
Altitude 
0.10 
1.7 
0.5 
1.8 
8,000 
10 km–10,000 km 
Sea surface temperature^{c} 
SST 
0.12 
1.9 
0.50 
1.8 
16,000 
500 m–50,000 km 
Soil moisture index^{d} 
MODIS 
0.05 
2.0 
0.14 
1.2 

512 m–25 km 
Vegetation index^{d} 
MODIS 
0.06 
2.0 
0.16 
1.2 


Surface magnetic field^{e} 
(Aircraft) <10 km 
0.14 
1.9 
0.6 
2.0 

250 m– 20 km 
(Aircraft) >10 km 
0.08 
2 
0.1 
1 

20 km–1,000 km 

Mostly Atmospheric^{k} 

State variables^{f} 
u, v 
0.09 
1.9 
1/3 (0.77) 
1.6 (2.4) 
(14,000) 
280 m–≈5,000 km 

w 
(0.12) 
(1.9) 
(–0.14) 
(0.4) 
(15,000) 


T 
0.11 (0.08) 
1.8 
0.50 (0.77) 
1.9 (2.4) 
5,000 (19,000) 


h 
0.09 
1.8 
0.51 
1.9 
10,000 


z 
(0.09) 
(1.9) 
(1.26) 
(3.3) 
(60,000) 

Precipitation^{g} 
R 
0.4 
1.5 
0.00 
0.2 
32,000 
4–12,000 km 
Passive scalars^{h} 
Aerosol concentration 
0.08 
1.8 
0.33 
1.6 
25,000 
100 m–125 km 
Radiances^{i} 
Infrared 
0.08 
1.5 
0.3 
1.5 
15,000 
2–15,000 km 

Visible 
0.08 
1.5 
0.2 
1.5 
10,000 
2–15,000 km 

Passive microwave 
0.1–0.26 
1.5 
0.25–0.5 
1.3–1.6 
5,000–15,000 
20–15,000 km 
Notes: When available (and when reliable), the aircraft data were used in precedence over the reanalysis values. In those cases where there was no comparable in situ value or when reanalysis was significantly different from the in situ value, the latter is given in parentheses. For the estimate of the effective outer scale, L_{eff}, where the anisotropy is significant, the geometric mean of the north–south and east–west estimates is given (the average ratio is 1.6:1 EW/NS, although for the precipitation rate the alongtrack TRMM estimate was used). Finally, the topography estimate of L_{eff} is based on a single realization (one Earth, one Mars); in both cases, the values pertain to the largescale scaling regimes: on Mars, >10 km; on Earth, >40 m. See Lovejoy and Schertzer (2013) for an extensive review. Note that the half circumference on Mars (the largest Martian distance) is 10,600 km, compared with 20,000 km on Earth.
From Gagnon et al. (2006); data from GTOPO30 (≈1 km resolution), ETOPO5 (≈10 km resolution), and the continental United States at 90 m and Saxony (Germany) at 50 cm.
From Landais et al. (2015); at scales smaller than 10 km there is also a (nearly) monofractal regime with H ≈ 0.75.
These values are from both SST data at ≈500 km resolution and also satellite data at ≈500 m resolution; see Table 8.2 of Lovejoy and Schertzer (2013).
From Lovejoy et al. (2007a).
From airborne magnetic field anomaly measurements at 800 m altitude over Canada (two regions). The two different scale regions correspond to a break at the horizontal scale corresponding to the Curie depth (see Lovejoy et al. 2001b; Pecknold et al. 2001). The inner scale (approximately indicated here as 20 km) is thus only a rough average; the Curie depth varies from ≈10 to 50 km.
From instrumented (Gulfstream) aircraft 280 m to 1,100 km (Lovejoy et al. 2010), the values in parentheses are from reanalyses from ≈100 to ≈5,000 km by Lovejoy and Schertzer (2011).
From gauge networks from 200 to 5,000 km, reanalyses (≈100–8,000 km) and from TRMM satellite radar data, 4–12,000 km (Lovejoy et al. 2012).
From airborne lidar backscatter ratios from 100 m to 125 km (horizontal) and 3 m to 5 km (vertical) (Lilley et al. 2004, 2008; Lovejoy et al. 2008).
From >1,000 orbits of TRMM satellite data with visible, nearIR, IR, and passive microwave at varying resolutions (see Figure 3.1a and b). Also from MTSAT (geostationary) satellite data (Figure 3.3c) in the IR, 30–5,000 km, and GOES (Geostationary Operational Environmental Satellite; IR, visible geostationary) satellite data (Lovejoy et al. 1993b; Lovejoy et al. 2001c); there are many other scaling analyses of cloud radiances that used other analysis techniques that give exponents that cannot be used to determine the parameters in this table.
This column gives the range of the scaling observed with the cited references. If the range covered by the cited data covers a welldefined scale break, then the range covered by the data may be larger. Similarly, the scaling range may cover scales beyond those that were empirically studied.
With the exception of the passive scalars from lidar and precipitation (from satellite radar), the radiances are not purely dependent on atmospheric conditions; they depend also on the surface emissivities and temperatures (IR, microwave) and albedos (visible).
Many surface fields are scaling over wide ranges, particularly as revealed by remote sensing. Figure 3.5a shows six moderateresolution imaging spectroradiometer (MODIS) channels at 250m resolution over Spain (each a 512 × 512 pixel “scene”). The scaling is again excellent except for the single lowest wave number, which is probably an artefact of the contrast enhancement algorithm that was applied to each image before analysis. These channels are used to yield vegetation and surface moisture indices by dividing channel pair differences by their means, so that the scaling is evidence that both vegetation and soil moisture are also scaling.
The vegetation and soil surface moisture indices are standard products described by Rouse et al. (1973) and Lampkin and Yool (2004):
Figure 3.5 (a) Spectra of six bands of MODIS radiances over a 512 × 512 pixel region of Spain (at 250m resolution; k = 1 corresponds to 128 km): E(k) as a function of the modulus of the wave vector. In order from top to bottom at the point log10k = 0.7, the curves are as follows: purple = Band 6, black = Band 1, magenta = Band 7, light green = Band 2, cyan = Band 4, dark green = Band 3. Reference lines have slopes –1.3. The band wavelengths are (in nm): Channel 1: 620–670, Channel 2: 841–876, Channel 3: 459–479, Channel 4: 545–565, Channel 5: 1,230–1,250, Channel 6: 1,628–1,652, Channel 7: 2,105–2,155. These data are used for determining both vegetation and soil moisture indices. Adapted from Figure 3.3a in Lovejoy et al. (2007a). (b) The spectra of Bands 1 and 2 and the vegetation index (Equation 3.9) compensated by k1.17 so that a spectrum E(k) ≈ k1.17 is flat. Note that the variations are less than ±0.2, whereas the actual spectra vary by almost five orders of magnitude (Figure 3.5a). (c) The spectra of Bands 6 and 7, and the soil moisture index (Equation 3.9) compensated by k1.17 so that a spectrum E(k) ≈ k−1.17 is flat.
where VI and SM are the vegetation and soil surface moisture indices and I is the radiance (the band number is indicated in the subscript); Equation 3.9 is usually stated without reference to any particular scale. We may immediately note that, although the Terra MODIS data have a resolution of 500 m (Bands 1 and 2 were degraded from 250 m) and the variability of the radiances and surface features continues to much smaller scales, the surrogates are defined at a single (subjective) resolution equal to that of the sensor. One of the applications of our analyses was to investigate how the relations between the surrogates and the bands used to define them change with scale: we anticipate that since the scaling properties are different, making the surrogates with data at different resolutions would produce fields with different properties. This is investigated in detail in the section “Soil Moisture and Vegetation Indices.”
In spite of the fact that gravity acts strongly at all scales, the classical theories of atmospheric turbulence have all been quasiisotropic in either two or three dimensions. Whereas a few models tentatively predict possible transitions in the horizontal (for example between k^{−5/3} and k^{−3} spectra for a transition from 3D to 2D isotropic turbulence for the wind), in contrast, in the vertical any 2D/3D “dimensional transition” would be even more drastic (Schertzer and Lovejoy 1985a; this is also true for passive scalars such as chaff; see e.g., Lesieur 1987).
The fact that these predicted drastic transitions in horizontal spectra have not been observed was the starting point for the anisotropic scaling model with the “in between” dimension D_{el} = 23/9 = 2.5555 D proposed by Schertzer and Lovejoy (1985a), in which the horizontal and vertical have power law spectra but with different exponents β_{h} and β_{v}. When—as observed—β_{v} > β_{h} (and it is not obvious), this implies that structures with horizontal extent L have vertical extents L^{Hz} with H_{z} = (β_{v} − 1)/(β_{h} − 1) so that they become flatter and flatter at larger and larger scales (see Figure 3.7b for an illustration). In addition, the volume of typical structures is LLL^{Hz} = L^{Del}, where D_{el} = 2 + H_{z} is the “elliptical” dimension characterizing the effective dimension of the (nonintermittent, i.e., spacefilling) structures.
Although the existence of different scaling exponents in the horizontal and vertical was confirmed by in situ measurements of many different atmospheric variables (see notably Lilley et al. 2004, 2008; Lovejoy et al. 2007a, 2007b, 2008; Radkevitch et al. 2008; Lovejoy and Schertzer, 2010a; and especially Chapter 6 in Lovejoy and Schertzer 2013), the overall situation took a long time to clarify. It turned out that there was a problem with the interpretation of aircraft data in the horizontal: spurious breaks in the scaling were caused by aircraft following isobars rather than isoheights (Lovejoy et al. 2009a). This had the effect of leading to a break in the horizontal scaling where the vertical displacement of the aircraft became large enough that the fluctuations in wind were dominated by the much stronger vertical shears (and larger vertical exponents). This was recently directly demonstrated by the first joint vertical–horizontal structure function analyses of the horizontal wind from 14,500 aircraft flights (Pinel et al. 2012).
Although it was important to clarify the status of the in situ wind measurements, the basic result—anisotropic scaling—had in fact already been very clearly established through the use of remotely sensed aerosol backscatter from airborne lidar (Lilley et al. 2004, 2008; Lovejoy et al. 2007a, 2007b, 2008; Radkevitch et al. 2008; see box 6.1 in Lovejoy and Schertzer 2013, for analogous results using CloudSat data and cloud liquid water radar reflectivities). The key differential stratification can be observed directly by eye in Figure 3.6a and b, which shows vertical cross sections of lidar aerosol backscatter fields with resolutions down to 3 m in the vertical. Starting at the lowresolution image (Figure 3.6a), we can see that structures are generally highly stratified. However, zooming in (Figure 3.6b) we can already make out waves and other vertically (rather than horizontally) oriented structures. In Figure 3.6c, we confirm—by direct spectral analysis—that the fields are scaling in both the horizontal and vertical directions and that the exponents are indeed different in both directions: the critical exponent ratio (β_{h} − 1)/(β_{v} − 1) = H_{z} is quite near the theoretical value 5/9, which is the ratio of the (horizontal) Kolmogorov value 1/3 and the (vertical) Bolgiano–Obukhov value 3/5. The wave number at which the spectra in Figure 3.6c typically cross is in the range from (10 cm)^{−1} to (100 cm)^{−1}; this is the scale at which structures are roughly “roundish,” called the “spheroscale.”
Figure 3.6 (a) Typical vertical–horizontal lidar backscatter crosssection acquired on August 14, 2001. The scale (bottom) is logarithmic: darker is for smaller backscatter (aerosol density surrogate), lighter is for larger backscatter. The black shapes along the bottom are mountains in the British Columbia region. The line at 4.6 km altitude shows the aircraft trajectory. (b) Enlarged content of the (700–1,600 m) box in (a). Note that small structures become more vertically aligned, whereas large structures are fairly flat. The aspect ratio is 1:96. Zoom of the previous image, showing that at the small scales, structures are beginning to show vertical (rather than horizontal) “stratification” (even though the visual impression is magnified by the 1:40 aspect ratio, the change in stratification at smaller and smaller scales is visually obvious). (Reproduced from Lilley, M., et al., Phys. Rev. E, 70, 036307, 2004. With permission.) (c) The lower curve is the power spectrum for the fluctuations in the lidar backscatter ratio, a surrogate for the aerosol density (B) as a function of horizontal wave number k (in m−1) with a line of best fit with slope βh = 1.61. The upper trace is the power spectrum for the fluctuations in B as a function of vertical number k with a line of best fit with slope βv = 2.15; hence Hz = 0.61/1.15 ≈ 0.53. (Adapted from Lilley, M., et al., Phys. Rev. E, 70, 036307, 2004.)
We have presented a series of striking widerange scaling spectra covering many significant atmospheric and surface fields. In this “tour” of the scaling, we have exclusively used a common statistical analysis technique (the power spectrum). The conclusion that scaling is a fundamental symmetry principle of wide applicability is hard to escape, although it is not universally embraced. The main difficulty seems to be that we are inculcated with the “scale bound” (Mandelbrot 1981) idea that every time we zoom into structures (by magnifying them), we expect to see qualitatively different structures and different processes. When trying to understand and model such phenomena, we are taught to search for mechanistic explanations involving processes acting over only narrow ranges of scale, and when we find them, we typically assume that they are incompatible with statistical theories and explanations. In atmospheric science, a particularly extreme form of this was recently highlighted by Lovejoy (2014), who compared a (still frequently cited) 1970s “mental picture” of atmospheric processes that proposed a spectrum dominated by narrow range spikes (such as the diurnal and annual cycles), with the rest being essentially uninteresting whitenoise “background” processes (Mitchell 1976). By comparing this speculation with modern data, it was found that the mental picture was in error by a factor of ≈ 10^{15} and that almost all of the variance was in the nontrivial scaling background. Similarly, publication of the nearly perfect space–time scaling of thermal infrared data over the Pacific (figure 3c, Pinel et al. 2014), was delayed for over a year due to strenuous referee objections based on a perceived incompatibility between the scaling statistics and the usual mechanistic phenomenology of tropical meteorology (e.g., the dominance of waves). There was a similar reaction to the scaling statistics of satellitebased Martian reanalyses in Lovejoy et al. (2014) and Chen et al. (2016).
Figure 3.7 (a) A selfsimilar (isotropic) multifractal cloud simulation. Each image is enlarged by a factor of 1.7 (the areas enlarged are shown in yellow and red rectangles for the first few enlargements, top rows). (b) A sequence “zooming” into a vertical crosssection of an anisotropic multifractal cloud with Hz = 5/9. Starting at the upper left corner, moving from left to right and from top to bottom, we progressively zoom in by factors of 1.21 (total factor ≈ 1,000). Note that while at large scales, the clouds are strongly horizontally stratified, when viewed close up they show structures in the opposite direction. The spheroscale is equal to the vertical scale in the leftmost simulation on the bottom row. The film version of this (and other anisotropic space–time multifractal simulations) can be found at http://www.physics.mcgill.ca/~gang/multifrac/index.htm. (Adapted from Lovejoy, S., and D. Schertzer, The weather and climate: Emergent laws and multifractal cascades, Cambridge University Press, Cambridge, 2013.)
Such reactions illustrate the “phenomenological fallacy” (Lovejoy and Schertzer 2007a), which arises when phenomenological approaches are only based on morphologies rather than underlying dynamics. The fallacy is to identify phenomenologically defined forms, structures, or morphologies with distinct dynamical mechanisms when in actual fact a unique dynamical mechanism acting over a wide range of scales can also lead to structures that change with scale. The fallacy thus has two aspects. In the first, form and mechanism are confounded so that different morphologies are taken as prima facie evidence for the existence of different dynamical mechanisms. In the second, scaling is reduced to its special isotropic “selfsimilar” special case in which small and large scales are statistically related by an isotropic “zoom” or blowup (Figure 3.7a). In fact, scaling is a much more general symmetry: it suffices for small and large scales to be related in a way that does not introduce a characteristic scale. Moreover, the relation between scales can involve differential squashing, rotation, and so on, so that small and large scales can share the same dynamical mechanism yet nevertheless have quite different appearances. Compare Figure 3.7a with b, which is the same process but with different horizontal and vertical exponents corresponding to differential (scaledependent) “squashing.”
Figure 3.8 This selfaffine simulation illustrates the “phenomenological fallacy,” because both the top and bottom look quite different while having the same anisotropic mechanism at scales differing by a factor of 64 (top and bottom blowup). The figure shows the proverbial geologist’s lens cap at two resolutions differing by a factor of 64. Seen from afar (top), the structures seem to be composed of lefttoright ridges; however, closer inspection (bottom) shows that in fact this is not the case at the smaller scales. (Reproduced from Lovejoy, S., and D. Schertzer, Nonlin. Processes Geophys., 14, 1–38, 2007b. With permission.)
In order to illustrate how morphologies can change with scale when the scaling is anisotropic, consider Figure 3.8. This is a multifractal simulation of a rough surface with the parameters estimated for the topography; its anisotropy is in fact rather simple in the framework of the generalized scale invariance (GSI) (Schertzer and Lovejoy 1985b; for a review, see Chapter 7 of Lovejoy and Schertzer 2013). More precisely, it is an example of linear GSI (with a diagonal generator) or “selfaffine” scaling. The technical complexity with respect to selfsimilarity is that the exponents are different in orthogonal directions, which are the eigenspaces of the generator, so that structures are systematically “squashed” (stratified) at larger and larger scales. The underlying epistemological difficulty, which was not simple to overcome and which still puzzles phenomenologists, corresponds to a deep change in the underlying symmetries. The top figure illustrates the morphology at a “geologist’s scale” as indicated by the traditional lens cap reference. If this were the only data available, one might invoke a mechanism capable of producing strong left–right striations. However, if one only had the bottom image available (at a scale 64 times larger), then the explanation (even “model”) of this would probably be rather different. In actual fact, we know by construction that there is a unique mechanism responsible for the morphology over the entire range.
Figure 3.9 Examples of continuous in scale anisotropic multifractals in 3D (256 × 256 × 64), showing the effect of changing the spheroscale (ls) on multifractal models of clouds with Hz = 5/9. The cloud statistical parameters are as follows: α = 1.8, C1 = 0.1, and H = 1/3 (similar to CloudSat and aerosols; see Table 3.1). From left to right, we decrease ls (corresponding to zooming out by factors of 4) so that we see the initially vertically aligned structures (bottom left) becoming quite flat at scales 64 times larger (right). At the same time, the horizontal structures have scaling anisotropies so that they too change orientation and elongation (the horizontal spheroscale starts at 1 pixel, far left; for a review, see Lovejoy and Schertzer [2013, chapters 6 and 7] for this generalized scale invariance). The middle row is a falsecolor rendition of the liquid water density field; the bottom row is the corresponding vertical sections (side view); the top row is the corresponding single scatter visible radiation; the mean optical thickness is 2; isotropic scattering phase function; sun incident at 45° to the right. (Reproduced from Lovejoy, S., et al., Atmos. Chem. Phys., 9, 1–19, 2009a. With permission.)
Figure 3.9 gives another example of the phenomenological fallacy, this time with the help of multifractal simulations of clouds. Again (roughly) the observed cascade parameters were used, yet each with a vertical “spheroscale” (this is the scale where structures have roundish vertical cross sections) decreasing by factors of 4, corresponding to zooming out at random locations. It is evident from the vertical crosssection (bottom row) that the degree of vertical stratification increases from left to right. These passive scalar cloud simulations (liquid water density, bottom two rows; single scattering radiative transfer, top row) show that by zooming out (left to right) diverse morphologies appear. Although a phenomenologist might be tempted to introduce more than one mechanism to explain the morphologies at different scales, in the figure we are simply seeing the consequence of single underlying mechanism repeating scale after scale. The phenomenological fallacy can undermine many classical ideas. For example, Lovejoy and Schertzer (2013, Box 6.1) argued with the help of CloudSat analyses that the classical twoscale theories of convection are incompatible with data that are scaling—that division into qualitatively distinct small and large regimes is unwarranted.
Most wavelengths used for remote sensing are sensitive to the radiative transfer properties of clouds, and we have seen (Figure 3.1a through d) that the associated radiances are scaling over a wide range of scales. It is thus reasonable to simulate clouds using scaling models for the distribution of scatterers and absorbers and then to model the radiative transfer through them at different wavelengths. Such models will be needed either to understand the clouds themselves—of fundamental importance in atmospheric science—or to understand the role of clouds in modulating the transmission/absorption characteristics of the underlying surface fields. An understanding of cloud and radiation variability and their interrelations over wide is also a challenging problem in the physics of disordered media. Figure 3.10 shows some examples for a purely scattering atmosphere (albeit with single scattering only). Figure 3.10a and b shows the external “faces” of the 3D liquid water concentration field for four realization of an anisotropic (stratified) multifractal cloud with exponents near to those observed; the difference in the figures is due to the different degrees of stratification, as can be seen best from the side views (Figure 3.10b). Figure 3.10c and d shows the corresponding visible radiation fields above and below the cloud using single scattering only and an isotropic phase function (no absorption; numerical details on the multifractal simulations can be found in Lovejoy and Schertzer [2010c] and Lovejoy and Schertzer [2010b]; for some theory and numerics for multiple scattered radiative transfer in multifractal clouds, see Lovejoy et al. [2009b] and Watson et al. [2009]).
The classical theory of radiative transfer is elegant (Chandrasekhar 1950) but is only relevant in 1D (“plane parallel,” horizontally homogeneous) media, yet the use of 1D models has long dominated the field. This is because when we turn to horizontally inhomogeneous media, there is no consensus on the appropriate model of heterogeneity, nor is the transport problem analytically tractable. As a consequence, the effect of horizontal variability was underestimated and usually reduced to the problem of inhomogeneity of the external cloud/medium boundaries (e.g., cubes, spheres, and cylinders; Busygin et al. 1973; McKee and Cox 1976; Preisendorfer and Stephens 1984) with the internal cloud and radiance fields still being considered smoothly varying if not completely homogeneous. When stronger internal horizontal inhomogeneity was considered, it was typically confined to narrow ranges of scale so that various transfer approximations could be justified (Weinman and Swartzrauber 1968; Welch et al. 1980).
When the problem of transfer in inhomogeneous media finally came to the fore, the mainstream approaches were heavily technical (see Gabriel [1993] for a review), with emphasis on intercomparisons of general purpose numerical radiative transfer codes (see the C^{3} initiative; Cahalan et al. 2005) and the application to large eddy simulation cloud models (Mechem 2002). At a more theoretical level, the general problem of the consequences of smallscale cloud variability on the largescale radiation field has been considered using wavelets (Ferlay and Isaka 2006) but has only been applied to numerical modelling. As a consequence, these “3D radiative transfer approaches” have generally shed little light on the scalebyscale statistical relations between cloud and radiation fields in realistic scaling clouds (e.g., see the collection in Marshak 2005). Overall, there has been far too much emphasis on techniques and applications with little regard for understanding the basic scientific issues.
Figure 3.10 (a) The top layers of threedimensional cloud liquid water density simulations (false colors); all have anisotropic scaling with horizontal exponents d = 1, c = 0.05, e = 0.02, f = 0 (see Lovejoy and Schertzer, 2013, for a review and definition of these parameters); stratification exponent Hz = 0.555; and multifractal statistical exponents α = 1.8, C1 = 0.1, and H = 0.333. They are simulated on a 256 × 256 × 128 point grid. The simulations in the top row have horizontal spheroscales of ls = 8 pixels (left column), 64 pixels (right column), with additional (trivial) anisotropy exponents ξ = 0 (top row), ξ = 3/4 (bottom row); they all have ξz = 0.25 (the ξ parameters—not to be confused with the structure function exponent ξ(q)—are explained in Lovejoy and Schertzer [2013], from which all these figures were reproduced). Note that in these simulations, ls = 8 and 64 applies to both the vertical and horizontal crosssections (i.e., ls = lsz). (b) A side view. (c) The top view with single scattering radiative transfer; incident solar radiation at 45° from the right; mean vertical optical thickness = 50. (d) The same as Figure 3.10c except viewed from the bottom.
The simplest interesting transport model is diffusion (on fractals see the reviews by Bouchaud and Georges [1990] and Havlin and BenAvraham [1987] and on multifractals see Meakin [1987], Weissman [1988], Lovejoy et al. [1993a], Lovejoy et al. [1998], Marguerite et al. [1997]). However—except in 1D (Lovejoy et al. 1993a, 1995)—diffusion is not in the same universality class as radiative transport (Lovejoy et al. 1990).
The first studies of radiative transport on fractal clouds (with a constant density on the support) were by Barker and Davies (1992), Cahalan (1994), Cahalan (1989), Cahalan et al. (1994), Davis et al. (1990), Gabriel et al. (1990), Gabriel (1986), Lovejoy et al. (1990), and Lovejoy (1989). These works used various essentially academic fractal models and focused on the (spatial) mean (i.e., bulk) transmission and reflectance. They clearly showed that (1) fractality generally leads to nonclassical (“anomalous”) thick cloud scaling exponents; (2) the latter were strongly dependent on the type of scaling of the medium; and (3) the exponents are generally independent of the phase function (Lovejoy et al. 1990).
Some general theoretical results exist for conservative cascades (H = 0; see the section “Multifractals, Structure Functions” for a definition of H and α and the section “Multifractals and Trace Moment Analysis” for cascades), single scattering for lognormal clouds (Lovejoy et al. 1995) (α = 2), and for more general (i.e., α < 2), “universal” multifractal clouds dominated by low density “Lévy holes” (frequent lowdensity regions where most of the transport occurs; Watson et al. 2009). The latter shows how to “renormalize” cloud density, that is, to relate the mean transmission statistics to those of an equivalent homogeneous cloud. Lovejoy et al. (2009b) extended these (numerically) to H > 0 and with multiple scattering including the case of very thick clouds. By considering the (fractal) path of the multiply scattered photons, it was found that due to longrange correlations in the cloud, the photon paths are “subdiffusive,” and that the corresponding fractal dimensions of the paths tend to increase slowly with mean optical thickness. Reasonably accurate statistical relations between N scatter statistics in thick clouds and single scatter statistics in thin clouds were developed, showing that the renormalized single scatter result is remarkably effective. This is because of two complicating effects acting in contrary directions: the “holes,” which lead to increased single scatter transmission, and the tendency for multiply scattered photons to become “trapped” in optically dense regions, thus decreasing the overall transmission.
All results to date are for statistically isotropic media; for more realism, future work must consider scaling stratification as well as the statistical properties of the radiation fields and their (scaling) interrelations with cloud density fluctuations.
Spectral analysis is often convenient, but it is not always easy to interpret. In addition, it is only a secondorder statistic so that—unless the process is quasiGaussian—it characterizes neither the intermittency nor the extremes. The real space alternatives are based on fluctuations of various sorts. Wavelets provide a general formalism for defining and handling fluctuations; indeed, it is so general that the choice of wavelets is often made on the basis of mathematical convenience or elegance. Although such wavelets may be useful in localizing singularities (if they are indeed localized!) or in estimating statistical scaling exponents, the interpretation of the fluctuations is often opaque. An exception is the Haar fluctuation, which was in fact the first wavelet (Haar 1910; before the full theory). Haar fluctuations are simple to understand; the wavelet formalism is not needed in order to be able to easily apply them in remote sensing and geophysical applications. We could also mention the detrended fluctuation analysis (Peng et al. 1994; Kantelhardt et al. 2002) method, which is nearly a wavelet method but, unfortunately with fluctuations that are not simple to interpret.
The Haar fluctuation ΔI(Δx) at timescale Δx of the intensity field I is simply the difference of the mean of I over the first and second halves of the interval Δx:
where the subscript Haar was added to distinguish it from other common definitions of fluctuation, and the x dependence was suppressed because we assume that the fluctuations are statistically homogeneous. With an appropriate “calibration” constant (a factor of 2 is used below and is canonical), in scale regions where H > 0, the Haar fluctuations are nearly equal to the differences; in scale regions where H < 0, they are nearly equal to the anomalies:
where Ī is the mean over the entire series. The Haar fluctuation is essentially the difference fluctuation of the anomaly; it is also equal to the anomaly fluctuation of the difference.
Now that we have defined the fluctuations, we need to characterize them; the simplest way is through (generalized) structure functions (generalized to fluctuations other than the usual differences, and generalized to moments of order other than the usual value of 2): ⟨ΔI (Δx)^{q}⟩ where “[·]” indicates statistical (ensemble) averaging; this is the qth order structure function.
Physically, if the system is scaling, then the fluctuations are related to the driving flux φ by the following:
where the subscript Δx on φ indicates that it is the flux at resolution Δx. The structure function is
Turbulent fluxes φ are conserved from scale to scale so that ⟨φ_{Δx}⟩ = constant (independent of scale), implying that ⟨ΔI(Δx)⟩ ∝ Δx^{H}, so that H is the mean fluctuation exponent. Beyond the simplicity of interpretation, the Haar fluctuations give a good characterization of the variability for stochastic processes with –1 < H < 1. In contrast, fluctuations defined as differences or as anomalies are only valid over the narrower ranges 0 < H < 1 and –1 < H < 0, respectively (see Lovejoy and Schertzer 2012; Lovejoy et al. 2013). Outside these ranges in H, the fluctuation at scale Δx is no longer dominated by wave numbers ≈Δx^{−1}, so that the fluctuations depend spuriously on details of the finite data sample, typically either the highest or the lowest frequencies that happen to be present in the sample. For H outside these ranges, other definitions of fluctuations (wavelets) must be used. However, empirically, almost all geophysical fields fall into this range.
In the spatial domain, remotely sensed fields typically have 0 < H < 1 (see Table 3.1), so that there is not much to be gained by using Haar fluctuations instead of differences. However, in the time domain, –1 < H < 0 is quite general for atmospheric fields with resolutions of 5–10 days or more; for oceanic fields, of resolutions ≈1 year or more (see chapter 10 of Lovejoy and Schertzer 2013), so that Haar fluctuations have the advantage of being able to handle both temporal and spatial analyses, whereas differences are often not appropriate in the corresponding temporal domains.
The generic scaling process is multifractal so that, in general, φ has the following statistics:
where K(q) is a convex function (the constant of proportionality depends on the external scale of the the process; see Equation 3.20). Substituting this into Equation 3.13, we obtain
where ξ(q) is the “structure function exponent.” Although we return to this in more detail below, for the moment note that the mean (q = 1) flux <φ_{Δx}> is independent of Δx, so that K(1) = 0 and hence ξ(1) = H. Note also that for quasiGaussian processes, none of the moments of φ_{Δx} have any scale dependence, so that K(q) = 0 and ξ(q) = qH (all the scale dependencies are characterized by H). Finally, the RMS fluctuation [⟨ΔI(Δx)^{2}⟩^{1/2}] has the exponent ξ(2)/2 so that the error in using the quasiGaussian approximation for the variance (i.e., ξ(2) = 2H) is ξ(2)/2 − H = K(2)/2.
In the section titled “The Horizontal,” we displayed the excellent spectral scaling of the MODIS bands used to determine soil moisture and vegetation indices, noting that these surrogates are defined (Equation 3.9) in an ad hoc way using the (subjective and finest available) resolution. Figure 3.11a shows the corresponding (difference fluctuation based) structure functions for the MODIS bands contributing to the vegetation and soil moisture indices. As expected from the spectra, except for the smallest scales where there are sampling and sensor smoothing artefacts, the scaling is excellent. The slopes (estimates of ξ(q)) are shown in Figure 3.11b.
Figure 3.11 (a) Structure functions based on difference fluctuations for the MODIS bands discussed in the text, plotted as a function of lags log10Δx (with distance Δx measured in pixels = 512 m). The regression lines (slope ξ(q), fit for the range Δx = 8–256 pixels) are also shown; the linearity is an indication of the quality with which the scaling is respected. Adapted from Lovejoy et al. (2007a). (b) The structure function exponent ξ(q) for MODIS bands 7, 1, 6, and 2 (top to bottom, the slopes from Figure 3.11a). The figure is derived from Lovejoy et al. (2007a).
The fact that the bands have nontrivial scaling implies that if the surrogates are defined at different resolutions, their statistical properties will be different; in other words, the singlescale surrogate (sss) can be (at most) correct at a single resolution. To see this more clearly, these sss (at scale ratio λ = L/Δx, where L is the large [image scale] and Δx is the resolution) can be contrasted with the corresponding multiscale surrogates (mss, at scale ratio λ). Mathematically, the difference can be expressed as follows:
where Λ = L/l is the maximum scale ratio (satellite image scale L)/(single pixel scale l) and the notation [I_{i, Λ}]_{λ} and ${\left[{\sigma}_{\mathrm{\Lambda}}^{\left(s\right)}\right]}_{\lambda}$
denotes averaging from this finest resolution up to the intermediate ratio λ < Λ (i.e., Δx > l). The mss is the surrogate that would be obtained by applying an identical algorithm (Equation 3.9) to satellite data at the lower resolution, whereas the sss is the surrogate at the same resolution but based instead on the finest scale available.From the single and multiplescale definitions of the surrogates, we can define the difference fluctuations:
where for the multiplescale definition (second line), we have Δx = L/λ. In Figure 3.12a (bottom), we show that the resulting mss obtained across the same range of scales as the sss (Figure 3.12a, top) is quite different (it is statistically significantly larger). The comparison of the corresponding exponents (ξ(q)) is given in Figure 3.12b.
Figure 3.12 (a) A comparison of the single and multiplescale indices (top and bottom rows, respectively) defined by Equation 3.17, that is, by degrading indices defined at the finest resolution to an intermediate resolution (single scale index) and by first degrading the resolution of the bands and then determining the indices (multiplescale index). The moments from top to bottom are in order q = 1.9, 1.7, 1.5, 1.3, 1.1, 0.9, 0.7, 0.5, 0.3, and 0.2. Adapted from Lovejoy et al. (2007a). (b) Comparison of the structure function exponents for the vegetation and soil moisture indices (red and black, respectively), estimated using the singlescale and multiplescale index definitions (solid and dashed lines, respectively; the slopes in Figure 3.12a). Derived from Lovejoy et al. (2007a).
An impressive application of the structure function method with Haar fluctuations was recently published by Landais (2014) and Landais et al. (2015), using data from the Mars Orbiter Laser Altimeter (MOLA). Figure 3.13a shows an analysis over more than four orders of magnitude of scale, displaying two scaling regimes with a break at ≈10 km. In this case, the altimeter has numerous holes (some large) caused notably by obstruction of the signal by dust clouds. The usual method of filling such holes is by interpolation, but if there is a significant amount of interpolation, then the results can be badly corrupted. The reason is that even the lowest order interpolation (linear interpolation) is once differentiable it has H = 1, whereas the topography is only differentiable of order H ≈ 0.52 (see Figure 3.13b). Therefore, linear interpolation will mix segments of H = 1 with segments of H = 0.52, generally yielding a spurious result. Thus, in order to estimate the Haar fluctuations it is best to use an algorithm that does not require interpolation. One such simple algorithm is described in the appendix of Lovejoy (2014) and is freely available at: http://www.physics.mcgill.ca/~gang/software/index.html; it was used to generate the figure.
Figure 3.13 (a) The Haar fluctuations from 5 × 105 Mars Orbiter Laser Altimeter altitude estimates for moments of order 0.1, 0.2, …, 1.9, and 2 (bottom to top). Δx is in meters; S(Δx) is the corresponding structure function for the q = 1 moment (the second red line from the bottom); the units are also in meters; and S(Δx) gives the mean vertical change over a distance Δx. (b) The structure function exponent ξ(q) estimated from the regression slopes in Figure 3.13a for the Martian topography (the blue curve). The bottom red curve is for the region >10 km and the top blue curve is for scales <10 km. While the latter is nearly linear and hence monofractal (with ξ(1) = H ≈ 0.75), the former is curved and multifractal. For this curve, we also show the graphical estimates of the mean fluctuation exponent H using H = ξ(1) ≈ 0.52 and (from the tangent at q = 1) the intermittency correction near the mean, C1 = ξ(1) − ξ’(1) ≈ 0.1. (Adapted from Landais, F., et al., Nonlin. Processes Geophys. Discuss., 2, 1007–1031, 2015.)
The previous examples estimating ξ(q) from MOLA and MODIS data underline a basic problemml: that the scaling generally involves an entire nonlinear convex function ξ(q) (equivalently the concave function K(q) = qH − ξ(q)) for its characterization, and empirically this is the equivalent of determining an infinite number of parameters. In order to make the problem manageable (to reduce the parameters to a finite number), we can make use of a multiplicative analogue of the usual (additive) central limit theorem for random variables. This shows that under fairly general circumstances, K(q) is determined by only two parameters that define multifractal “universality classes”:
where α is the Levy index and C_{1} is the codimension of the mean (Schertzer and Lovejoy 1987). From Equation 3.18, we see that C_{1} = K’(1); this provides a convenient way of estimating the parameters. Figure 3.13b shows how it can be graphically estimated from ξ(q); for α, it is also possible to use α = K″(1)/K′(1). Table 3.1 shows the resulting parameter estimates for many geophysical fields, including several that were remotely sensed. In addition, note that the difference between the exponents of the RMS and mean (important for interpreting the slopes in the RMS Haar graphs) ξ(2)/2 − H = K(2)/2 = A_{α}C_{1}, where A_{α} = (2^{α−1} − 1)/(α − 1). Because empirically 1.5 > α > 2 (Table 3.1), we find 0.83 > A_{α} > 1 (near 1) so that often C_{1} provides a good estimate of the error in using the RMS exponent ξ(2)/2 in place of H. From the table it is evident that, in space, the difference ξ(2)/2 − H can readily be ≈ 0.1, which is significant. The intermittency quantified by C_{1}, α leads to much larger variability than would be expected from classical (quasiGaussian) processes. Additional difficulties in the interpretation of data analyses arise when the scaling is anisotropic (see Lovejoy and Schertzer [2007b] or section 7.1.6 of Lovejoy and Schertzer [2013]).
Note that with the help of remote sensing, many solid earth fields have also been shown to be scaling over various ranges. These include the rock density, magnetic susceptibilities, surface magnetic fields, and surface gravity fields; see the review by Lovejoy and Schertzer (2007b) for examples, references, and theory for these geopotential fields. Visible and infrared emissions from geothermal regions have been studied by Laferrière and Gaonacߡh (1999), Harvey et al. (2002), Gaonacߡh et al. (2003), and Beaulieu et al. (2007).
We investigated the nonlinearity of ξ(q), K(q) (intermittency) directly from the structure function, estimating ξ(q) from the exponents of the qth order moments and then estimating K(q) as ξ(q) + qξ(1) (see Equation 3.15, recalling that K(1) = 0). However, frequently, the H value is quite a bit larger than the C_{1} value, so that the linear contribution to ξ(q) is much larger than the nonlinear contribution, making C_{1}, α difficult to accurately estimate. However, by estimating φ and hence K(q) directly, we can obtain improved estimates of C_{1}, α. We are also able to estimate the outer scale of the underlying cascade process.
To see how this works on gridded data, we can conveniently estimate ΔI as the absolute finite difference along the transect. The normalized highresolution flux is then obtained by dividing Equation 3.12 by its ensemble average to obtain the following:
Finally, this highresolution flux can be systematically degraded to lower resolution by averaging. The generic multifractal process is a multiplicative cascade; if such a process starts at outer scale L_{eff}, then the statistics follow:
where λ′ is the scale ratio of the outer scale to the resolution scale of the degraded flux Δx and φ_{λ′} is the normalized flux at scale ratio λߡ. When $\u27e8{\phi}_{{\lambda}^{\prime}}^{q}\u27e9$
is estimated in this way (by taking q powers of the successively degraded flux), it is called a trace moment. In empirical analyses, L_{eff} is not known a priori; it has to be estimated from the data. Here, we replace it by a convenient reference scale, L_{ref} = 20,000 km, which is the largest distance on Earth (half the circumference), and use λ = L_{ref}/Δx. If Equation 3.20 holds, then for all q the lines of $\text{log}\u27e8{\phi}_{\lambda}^{q}\u27e9$ against log λ will cross at a scale corresponding to λ = λ_{eff} = L_{ref}/L_{eff}.Figure 3.14a shows the results of estimating the various moments of order 2 ≥ q ≥ 0 for the TRMM orbiting weather radar data analyzed along the orbit direction. It is evident from the log–log linearity that the scaling is excellent up to nearplanetary scales of 10,000 km, and the convergence of the lines to a common outer scale, λ_{eff}, is a direct confirmation of the multiplicative (“cascade”) nature of the statistics (Equation 3.20). The lines plausibly cross at a scale of the order of the size of the planet (see Table 3.1). The fact that L_{eff} = L_{ref}/λ_{eff} can be a little larger than the size of Earth is because, even at planetary scales (20,000 km), there is some residual variability due to the interaction of the precipitation field with other atmospheric fields—L_{eff} is simply the “effective” scale at which the cascade would have had to start in order to explain the statistics over the observed range.
Precipitation is particularly interesting because, according to Table 3.1, it is the most intermittent of the common geophysical fields (the value of the intermittency parameter C_{1} is the highest, ≈0.4) and therefore is presumably the most difficult to measure. Figure 3.14b shows the corresponding results for other precipitation products (both gauge and reanalysisbased) as well as for a lowerresolution TRMM product, all analyzed in the east–west direction. It is evident that the qualitative behaviors are very similar. More detailed analysis (Lovejoy et al. 2012) shows that there are nevertheless significant differences in parameters, notably C_{1}. Similarly, in the macroweather regime (for precipitation, ≈2–4 days up to 40 years in the industrial epoch), this is also true (de Lima and Lovejoy 2015; Lovejoy and de Lima 2015), leading to the conclusion that the problem of accurately estimating areal precipitation is still unsolved.
Figure 3.14 (a) The TRMM reflectivities (4.3km resolution). The moments M (see Equation 3.20) are for q = 0, 0.1, 0.2, …, 2, taken along the satellite track. The poor scaling (curvature) for the low q values (bottom with negative slopes) can be explained as artefacts of the fairly high minimum detectable signal. Lref = 20,000 km so that λ = 1 corresponds to 20,000 km; the lines cross at the effective outer scale ≈32,000 km, C1 ≈ 0.63. (Reproduced from Lovejoy, S., et al., Geophys. Resear. Lett., 36, L01801, 2009c. With permission.) (b) East–west analyses of the gridded precipitation products discussed in the text. Upper left: The TRMM 100 × 100 km, 4day averaged product. Upper right: The ECMWF interim stratiform rain product (all latitudes were used). Note that the data were degraded in constant angle bins so that the outer scale is 180°. To compare with the other analyses, a mean map factor of 0.69 has been applied (the mean east–west outer scale was ≈14,000 km). Lower left: The CPC hourly gridded rainfall product (US only). (Reproduced from Lovejoy, S., et al., Adv. Water Res., 45, 37–50, 2012. With permission.) (c) The temporal analyses of the various precipitation products analyzed spatially in Figure 3.14b (for the TRMM satellite radar at 4 days to 1 year; the ECMWF stratiform rain product, 3 hours to 3 months; and the CPC network, 1 hour to 29 years; the 20CR reanalysis was for 45°N, every 6 hours for 138 years). (Reproduced from Lovejoy, S., et al., Adv. Water Res., 45, 37–50, 2012. With permission.)
Interestingly, the same data can be analyzed in the temporal domain. The reason for expecting scaling in time as well as space is that the wind field that connects the two is scaling (see, e.g., the spectrum in Figure 3.1c). Figure 3.14c shows that the straight lines on the log–log plot converge—the behavior expected for multiplicative cascades (Equation 3.20) and the signature of multifractality—but at timescales a bit longer than the outer scale of the weather regime (the point timescale at which the scaling becomes poor, closer to ≈10 days). The main exception is the TRMM data, which seem to have excellent scaling out to the limit of the 5,900 orbits (about 1 year). It is possible that the reason for this significantly larger outer timescale is that the TRMM data are mostly captured over the tropical oceans and the oceans have a significantly longer outer scale (about 1 year) due to the much lower energy rate densities in the ocean (see chapter 8 of Lovejoy and Schertzer, 2013, for a review of the theory and data).
Other significant examples of the use of trace moments are for characterizing the Earth’s energy budget as a function of scale. Figure 3.15a and b shows the results for the TRMM satellite (from two of the channels whose spectra were analyzed in Figure 3.1a and b) at visible and thermal infrared wavelengths, respectively. Upon inspection, it is evident that the multiplicative cascade Equation 3.20 holds very accurately (to about ±0.5%) up to about 5,000 km. Figure 3.15c shows the analogous result for the (geostationary) MTSAT thermal IR radiances over the Pacific (same as the data in Figure 3.1c). Finally, Figure 3.16 shows the results of trace moment analysis on meteorological state variables estimated through a reanalysis (see Figure 3.2). It is therefore not surprising that numerical weather prediction models also show excellent scaling and similar multiplicative trace moments (see Stolle et al. [2009] and Stolle et al. [2012] for analyses of several forecast models). In situ aircraft data of meteorological variables of state also show very similar cascade characteristics (Lovejoy et al. 2010).
Figure 3.15 (a) TRMM visible data (0.63 mm) from the VIRS instrument, Channel 1, with fluxes estimated at 8.8 km. Only the welllit 15,000km orbit sections were used. Lref = 20,000 km, so that λ = 1 corresponds to 20,000 km, the lines cross at Leff ≈ 9,800 km. (Reproduced from Lovejoy, S., et al., Geophys. Resear. Lett., 36, L01801, 2009c. With permission.) (b) Same as Figure 3.15a except for VIRS thermal IR (Channel 5, 12.0 µm), Leff ≈ 15,800 km (Reproduced from Lovejoy, S., et al., Geophys. Resear. Lett., 36, L01801, 2009c. With permission.) (c) Logs of normalized moments M (Equation 3.20) vs. log10λ for 2 months (1,440 images) of MTSAT, thermal IR, 30km resolution over the region 40oN to 30oS, 130o east–west over the western Pacific, the average of east–west and north–south analyses. Lref = 20,000 km so that λ = 1 corresponds to 20,000 km; the lines cross at the effective outer scale ≈ 32,000 km (from Pinel 2012) and C1 ≈ 0.074 (close to the TRMM thermal IR results, Table 4.7a, VIRS 4, 5).
Finally, we can return to the topography data whose spectra were analyzed in Figure 3.3. Figure 3.17 shows the cascade structure of the topographic gradients obtained by combining the four different data sets used in Figure 3.3 spanning the range 20,000 km down to submetric scales. As for the spectrum, the scaling holds quite well until around 40 m. Gagnon et al. (2006) argued that this break is due to the presence of trees (for the high resolution data set used over Germany, 40 m is roughly the horizontal scale at which typical vertical fluctuations in the topography are of the order of the height of a tree). Over the range of planetary scales down to ≈40 m, it was estimated that the mean residue of the universal scaling form with parameters C_{1} = 0.12, α = 1.79 (for all moments q ≤ 2) was ±45% over this range of nearly 10^{5} in scale (this error estimate was for the “reduced” moments ⟨φ^{q}⟩^{1/q}, e.g., for q = 2, root mean square moments).
Figure 3.16 The trace moment analysis of the 700mb ECMWF interim reanalysis products analyzed in the east–west direction, data at 700 mb (about 3 km altitude) between ±45o latitude. The fields are specific humidity (hs), temperature (T), east–west wind (u), north–south wind (v), vertical wind (w), and geopotential height (z); fields at 24hour intervals for the year 2006 were used. The moments are for orders q = 2, 1.9, and 1.8 (descending from top to bottom); λ = 1 corresponds to 20,000 km. See Table 3.1 for some parameter estimates. (Reproduced from Lovejoy, S., and D. Schertzer, The weather and climate: Emergent laws and multifractal cascades, Cambridge University Press, Cambridge, 2013. With permission.)
Remotely sensed products remapped into convenient coordinate systems are increasingly available at regular time intervals. This allows for the systematic study and assessments of trends. It also requires us to understand the full (joint) space–time variability. As an example, in this section we study global scale thermal infrared data from MTSAT (see Figure 3.1c for the 1D spectra and Figure 3.18a through c for three 2D subspaces) and show how knowledge of its spectral density in (k_{x}, k_{y}, ω) space (horizontal wave number and frequency) can be used to determine AMVs, which are used as surrogates for winds.
Figure 3.17 The trace moments of the data sets whose spectra were analyzed in Figure 3.4: ETOPO5 (circles), continental United States (Xs), and lower Saxony (squares). For the latter, a subsample was also analyzed (light circles) that was (mostly) free of trees, the difference indicating the effect of trees. The regression lines distinguish between the q values of the associated points: q = 2.18, 1.77, 1.44, 1.17, 0.04, 0.12, and 0.51 (top to bottom). (Reproduced from Gagnon, J. S., et al., Nonlin. Proc. Geophys., 13, 541–570, 2006. With permission.)
The MTSAT data set used for this study is comprised of nearly 1,300 hours of hourly geostationary MTSAT data (at 30km resolution) from an 8,000 × 13,000 km^{2} region centered on the equatorial Pacific (see Pinel et al. 2014, for details). It is convenient to use Fourier techniques. Introduce the nondimensional space–time lag ΔR and nondimensional wave vector K:
The nondimensionalization can be achieved by using the pixel and sampling periods or the planet scale (space) and the corresponding lifetime of planetary structures (time, roughly 10 days). We can now estimate the space–time (“st”) spectral density P_{st} (K) ∝ ⟨Ĩ(K)^{2}⟩ (see Equation 3.2). We can determine the three 1D spectra by successively integrating out the conjugate variables:
Figure 3.18 (a) Contours of log P(kx, ky), the spatial spectral density. Black represents empirical contours; colored lines are from the fit using Equation 3.32 with a = 1.2 ± 0.1. (b) Contours of log P(kx, ω), the zonal wave number/frequency subspace. Black represents empirical contours; colored lines are from the fit using Equation 3.32. The orientation is a consequence of the mean zonal wind, –3.4 m/s. (c) Contours of log P(ky, ω), the meridional wave number/frequency subspace. There is very little if any “tilting” of structures because the mean meridional wind was small: 1.1 m/s. Black represents empirical contours; colored lines are from the fit using Equation 3.32. (Reproduced from Pinel, J., et al. Atmos. Resear., 140–141, 95–114, 2014. With permission.)
where the righthand side equalities are consequences of assuming space–time scaling; the exponents can in principle all be different. Figure 3.1c shows the result: it is evident that the east–west, north–south, and temporal spectra are nearly identical up to ≈(5–10 days)^{−1} and are nearly perfect power laws (most of the deviations from linearity in the figure can be accounted for by the finite resolution and finite data “window”; see the black line that theoretically takes these limits into account). The main exceptions are the two small spectral “bumps” at (1 day)^{−1}, (12 hours)^{−1}. The fact that the 1D spectral densities are nearly identical shows that β_{x} = β_{y} = β_{t} = β and that the spectrum satisfies the isotropic scaling symmetry:
where empirically, s = β + 2 ≈ 3.4. Note that at any given scale P_{s} will generally display anisotropy; Equation 3.22 simply implies that the anisotropy doesnߡt change with scale (see Pinel et al. 2014, for the full analysis).
By necessity, the scale symmetry (Equation 3.22) can only hold up to planetary scales (L_{e} = L_{ref} = 20,000 km); this implies a breakdown in the time domain at scales τ, which is interpreted as the lifetime of planetary scale structures. In Figure 3.1c, we find τ ≈ 5–10 days; Lovejoy and Schertzer (2010a) describe how this is determined by the turbulent energy flux ε(power/mass): $\tau ={\epsilon}^{1/3}{L}_{e}^{2/3}$
, where ε itself is determined by the solar flux. This theory well describes the spectrum of the many atmospheric variables—ocean temperatures as well as the Martian weather and macroweather (Lovejoy et al. 2014), see also Chen et al. 2016.In order to solve the functional Equation 3.22 expressing the scaling symmetry on P, it is convenient to introduce real and Fourier space–time scale functions 〚ΔR〛, 〚K〛_{F}. The scale functions are generalizations of the usual notion of vector norm (distance). They generally satisfy a scale function equation (for this GSI, see Chapter 7 in Lovejoy and Schertzer 2013). In the general case (needed, for example, to deal with vertical stratification), we have the following:
where G is the generator of the anisotropy (in the basic linear case, G is a matrix so that the anisotropy depends on scale, but not position). In the special isotropic case in Equation 3.23, we simply have G = 1 = the identity we have:
With this, the empirical spectra (Figure 3.18, Equation 3.22) can be written as follows:
where P_{0} is a constant. A basic result (a corollary of the Wiener–Khinchin theorem) relates the secondorder structure function and the spectral density: it states that ⟨ΔI^{2}⟩ and 2(1−P) are Fourier transform pairs. With this, we obtain the following:
with s = 2 + β as before.
Theoretical considerations based on turbulence theory (Lovejoy and Schertzer 2010a; Pinel et al. 2014) lead to the convenient result for the nondimensional scale functions:
with the matrices B, B^{−1} given by
µ_{x} and µ_{y} are the components of the nondimensional average advection velocity:
where V_{w} is the largescale turbulent velocity and a is the north–south/east–west aspect ratio (taking into account the fact that the largescale gradients are typically a times larger in the north–south direction than in the east–west direction). If we now introduce the modified nondimensional speeds and frequency:
we find:
The corresponding spectral density is as follows:
Equation 3.32 has the interpretation as a mean advection by µ′ with a spatial squashing by factor a in the east–west (x) direction.
In order to test this in detail, Pinel et al. (2014) considered the spectral densities in the three 2D subspaces (k_{x}, ω), (k_{y}, ω), and (k_{x}, k_{y}) obtained by integrating out the (single) conjugate variables from P_{st}(k_{x}, k_{y}, ω) (Figure 3.18a through c). This enabled the theory to be tested and the parameters to be estimated. A further analysis (Pinel and Lovejoy 2014) showed that the (small) residuals could be interpreted in terms of multiplicative wavelike perturbations to P_{st}. This analysis is theoretically motived by improvement of the wellknown “turbulence background” of space–time IR radiances analyzed in the traditional wave context by Wheeler and Kiladis (1999) and Hendon and Wheeler (2008).
Finally, Pinel et al. (2014) showed how to use the real space structure function (based on the quadratic form determined by the matrix B) as a theoretical justification for a series of essentially ad hoc techniques starting with Hubert and Whitney (1971) for measuring “satellite winds,” now more accurately called AMVs (Atmospheric Motion Vectors). Using the real space parts of Equations 3.27, 3.28 and 3.29, we obtain an equation for the (secondorder) structure function (Equation 3.26) in terms of the winds. The structure function is linearly related to the autocorrelation function and hence gives useful formulae for relating the maximum correlation (which turns out to be the minimum structure function) and the winds. Although there are variants, this maximum crosscorrelation method is currently used operationally with MTSAT and GOES geostationary IR, visible satellite imagery (see, e.g., Szantai and Sèze [2008] for an overview and comparison).
As resolutions improve, remote sensing allows us to take a veritable voyage through scales, effectively zooming further and further into complex geofields. The journey potentially takes us over 10 orders of magnitude in scale—from planetary to submillimetric. In order to understand, characterize, and model such geocomplexity, we need an appropriate theoretical framework. The traditional framework is scale bound; it assumes a priori that as we change scales we uncover qualitatively different processes and phenomena. It postulates simple, usually deterministic models, at best explaining the variability over a narrow range of space–time scales. It is usually justified by an appeal to phenomenology, to the fact that structures at different scales often look different.
In this paper, we argued that on the contrary most of the geovariability (for example, most of the spectral variance) is in the widerange scaling background processes. With the help of numerical multifractal simulations, we pointed out the phenomenological fallacy: that scaling is generally anisotropic, in which case structures do change appearance with scale, even though the underlying mechanism simply repeats scale after scale—it is scale invariant. If in these cases, the morphologies (rather than statistics) were used to infer the dynamics, then several distinct mechanisms might be invoked, whereas there was in fact only one—hence the fallacy. Although it is true that we have discussed natural geosystems, which often respect a scaleinvariance symmetry, it is possible that the same considerations might also apply to at least some humancreated morphologies such as the spatial distribution of urban areas.
This paper is an update of Pecknold et al. (1997), which concentrated on developing the basic framework of multifractal processes and GSI needed for handling extremely intermittent scaling variability and anisotropic scaling, respectively. We discussed some striking new global scale examples: visible, infrared, passive, and active microwave sensing of the atmosphere; meteorological reanalyses (meteorological state variables) of ocean color, soil moisture, and vegetation indices; topography (Earth and Mars). We focused on improvements in (real space) fluctuation analysis methods (especially the Haar fluctuations), as well as on trace moment analysis, which can now be routinely used to estimate the external scale of the process and to directly test the multiplicative nature of the intermittency. Finally—reflecting the increasing availability of regular temporal series of images—we reviewed some recent work on space–time analyses from geostationary infrared satellite data. It showed, remarkably, that the (horizontal) space–time scaling has an isotropic exponent, while at the same time having “trivial” anisotropy that reflects the advection and turbulent velocities and can be used as theoretical bases for AMV determinations.