Too many exponents

The structure of scale invariant things is determined entirely by the values of various exponents: fractal dimensions, spectral exponents, and more. These exponents are typically defined using some power law function relating a measure of structure and a measure of scale.

But there are myriad ways to define "structure" and "scale". Some definitions are more fundamental within the scope of a theory, others are easier to measure or understand intuitively. Often, particularly before any principled theory has been accepted, these structure-scale relationships are empirical. This is the situation in cloud physics: there is no widely accepted theory, and there are a wide range of empirical power-law relationships found throughout the field. In many cases, it's hard to interpret these relationships physically or relate one exponent to another. During my Ph.D., I've tried to better understand two basic categories of exponents: geometrical exponents for clouds, such as fractal dimensions, and statistical exponents for atmospheric fields, such as spectral exponents for kinetic energy.

Most of these exponents might appear, at first, unrelated. But this apparent diversity is something of an illusion. If a field is scale invariant, there are simply not that many degrees of freedom in how it is constructed, and therefore many of these exponents contain the same information despite their different definitions. This has been proven for at least one class of fields, namely multiplicative cascades, and a subset of exponents, namely moment scaling exponents, by Lovejoy and Schertzer.¹ An apparent infinity of different exponents was shown to contain only three independent degrees of freedom.

This collapsed a large portion of the space of possible exponents and was a big win. But they did not consider all possible exponents, and those they did consider do not encompass all of the possible information about the fields they were interested in. You could also construct fields that are still scale invariant but use an entirely different construction method, and their results may not apply. Most unfortunately, the exponents they did consider, which were statistical, are not the ones that are popular in cloud physics today, which are geometrical.

Their result leaves me with the suspicion that an even larger group of exponents secretly contain the same information even though they have different definitions. In particular, I suspect that many of the geometrical exponents, such as fractal dimensions and size distribution exponents, contain identical information about the field structure and could be related to each other. I also suspect there are simple relationships between many geometrical exponents and many statistical exponents.

If this is true, it would be nice to know what the specific relationships are that relate various exponents. Some exponents are easier to measure than others in observational data, and some are simply more popular in the literature and therefore are more likely to be reported, so it would make it easier to compare reported results across studies. Additionally, since some exponents are easier to understand intuitively, having more relationships between exponents would help us better understand each exponent individually. Exponents that affect the character of turbulent intermittency, for example, are less intuitive than geometrical exponents like fractal dimension.

The problem is that theoretical exponent relationships can be quite difficult to obtain, and small, seemingly subtle definitional changes can invalidate them. In a recent publication titled Toward less subjective metrics for quantifying the shape and organization of clouds, we considered this problem using mostly theoretical arguments for a small number of geometrical exponents that characterize cloud fields. One of the most important points we make is that there is no single "fractal dimension" for cloud fields, even setting aside the question of the threshold used to define cloud. We showed that at least two different dimensions could be usefully defined, and their values are distinct but related through a third exponent defined through a cloud size distribution.

In the time since the bulk of the scientific work was completed for that publication, I've realized that the landscape of possible exponent definitions is far broader than we considered there. So, here, I want to push farther into the space of possible exponents and the relationships between them. I also wanted to record some further findings on the difficulty of estimating these exponents in finite-resolution data. This problem is becoming a persistent thread throughout my work.

Evaluating whether empirical data is scale invariant, or estimating the values of the exponents characterizing the scaling properties, is inescapably and nontrivially biased by finite resolution and domain size. Bias can be partially mitigated by changing the methods by which you compute the exponents, but these methodological changes often result in truly distinct exponents that are not equal in value but are often conflated in the literature. There is no single "fractal dimension", but a whole landscape of distinct dimensions with distinct values even for the same underlying data. Each dimension has different convergence properties in a finite domain, but in nearly all cases, nontrivial discretization error persists to scales that are orders of magnitude larger or smaller than what you might expect.

These issues make empirical analysis of scaling properties very subtle and difficult. However, I am optimistic about the possibility of simulating scale invariant fields with known scaling properties and comparing these to observations. Using simulations, finite domain size and discretization effects can be reproduced in simulation with known scaling properties relatively easily, and in this way one could convincingly isolate discretization bias from the underlying scaling properties of the phenomenon you are trying to measure.

Why exponents?

Firstly, the discussion above of "exponents" implies power-law functions. Is the restriction to power-law functions sufficient when we are considering such a general relation between structure and scale?

The answer is yes for a scale invariant field. The fields I will consider here are random, so consider a single realization $F_i(r)$. This could be one computer generated field with one set of input parameters and one particular random seed $i$, and it would have a unique value at every location $r$. If I changed the value of the seed $i$ but not the input parameters, the generated field would have completely different values everywhere. But for each set of input parameters, the value of any finite statistic of your choice, for example computed across many different seeds, has a single well-defined value. In this sense, the "statistics" of the field are unique for each set of input values but vary if you change the input parameters. Note that I said "finite" statistic. Some statistics are theoretically infinite, which simply means they will not converge to any single value if you tried to compute them over a bunch of different generations/seeds.

Consider taking a single realization and zooming in or out by a factor $\lambda$ to obtain a new field $F_i(\lambda r)$. If the field is scale invariant, then the zoomed field can be rescaled by some factor $c{F}_i(\lambda r)$, to obtain another field that has identical statistics to the original. In other words, there is nothing you could do to distinguish the zoomed and scaled field $c{F}_i(\lambda r)$ from a different generation, with a different seed, of the original: $F_j(r)$. This can be written as equality "in distribution"

$$F(r)\stackrel{d}{=}cF(\lambda r),$$

meaning the two fields are not equal but their statistics are equal. Furthermore, for scale invariance it is important that $c$ is constant regardless of which $\lambda$ you choose. This is what it means to say a field is "scaling". What this implies is that there is nothing you can do to know whether you are looking at the "original" field or a zoomed and rescaled field. And if the zoomed and rescaled fields are statistically identical, then there is no meaningful notion of absolute scale. The field is identical, statistically, at every possible scale.

Now consider that you have some operation $S(l)$ which computes a certain statistic of the field at a certain scale $l$. This could be, for example, the mean area of clouds with width $l$. The same statistic could be computed for a different scale $S(\lambda l)$. How do these statistics differ? Because there is no meaningful notion of absolute scale, the two statistics can only differ by a function of the relative scale $\lambda$:

$$S(\lambda l)=f(\lambda )S(l).$$

The only possible solutions to this functional equation follow a power law form: $f(\lambda )={\lambda }^{\tau }$ where the exponent $\tau$ is a constant whose value depends on your chosen statistic $S$. Thus, it is sufficient to restrict our view to only power-law functions and their exponents.

Exponents in finite domains

The above argument applies in the reverse direction as well. If you measure some statistic of your favorite field, and discover that it is NOT a power law function of scale, then your field is not scale invariant.

This brings up a crucial point that is sometimes missed. All of the exponents I'll discuss are defined assuming a scale invariant field and a power law function. If you observe a non-power law function with a non-constant exponent, the exponent cannot be interpreted in the same way. I've often seen this mistake for the fractal dimension in particular: a study will report that the fractal dimension is "larger for large clouds". If the "fractal dimension" increases with scale, you are observing a non-power-law function, and the concept of a "fractal" does not apply because the thing you are measuring is not scale invariant. It's fine to describe some arbitrary function by saying it has an exponent that increases with scale, but such an exponent should not be interpreted as any kind of fractal dimension.

The question you should then ask is why your observations are not scale invariant. The obvious answer is that the thing you are attempting to measure is not scale invariant. But this is not always the reason. When we measure clouds using satellite imagery, for example, there are two artifacts of the imaging method itself that are decidedly not scale invariant: both your overall image, and each individual pixel, are shaped like rectangles (or sometimes, circles). These Euclidian shapes have a well-defined scale and therefore break scale invariance, even when the underlying cloud you are observing is scale invariant. Correctly identifying and mitigating the problem is difficult and often requires a bespoke methodological correction for each exponent you consider — see, for example, our other paper that attempts to do so for size distribution exponents.

In that paper, we argued that it is far too common to neglect subtle finite domain and resolution effects and to simply assume that the underlying clouds are not scale invariant. This is the inverse error of calling an increasing exponent a "fractal dimension". In the atmosphere, two seemingly incompatible things are often true: the thing you are measuring is scale invariant, and your measurements are not scale invariant.

Does this mean that nothing is ever truly scale invariant, since every image has pixels and Nature itself is quantum-mechanical and discrete? You could take that position, but it's not very useful. Alternatively, we could refine our notion of scale invariance to mean a field that is scaling over only some finite range of scale, or where we find power laws in our observable statistics for some restricted range of sizes. In this case, an image of clouds could very well be scale invariant, but only for scales that are not influenced by the pixel or domain size. And unfortunately, exactly which scales are contaminated depends greatly on which specific power-law function/exponent you are considering.

Here's how this looks in practice for the box dimension of 32 large 2D $8192^2$ simulations of a scale invariant field. I'll define the box dimension and field later more precisely, but for now what's important is that this is a function which is theoretically known to follow a power law for this simulation, which means that departures from power-law behaviour can only be caused by finite pixel and domain geometry. On the left is a plot of the relevant power-law function, and on the right is a plot of the exponent as a function of scale.²

If the function was a power law, the exponent would be constant with scale. You can see that it is almost constant for the middle scales, but near pixel scales and near the domain scale it deviates. With respect to this measure of scale invariance, this field is therefore only moderately scale invariant at scales much larger than a pixel but much smaller than the domain.

If you are actually trying to measure the fractal dimension, you have to make a choice about which scales you use to compute the dimension. You have to choose a threshold for what "much larger" and "much smaller" means, balancing the removal of the worst of the discretization effects while retaining a relatively wide range of scales to fit the exponent to. Although there are some specific heuristics you can use for some specific power-law functions (see again our size distribution paper), there is no general method, and so it usually ends up being a subjective choice, which is unfortunate. Appropriate fitting thresholds can also vary significantly based on which specific power-law function you consider. For example, consider a different power-law exponent, the correlation dimension, which is theoretically equal to the box dimension above:

This exponent is far more constant over the middle range of scales than the box dimension, despite it being computed for the exact same fields. There are still discretization effects at small scales and influence from the square domain shape at large scales, which cause deviations from the theoretical exponent, but both issues are much more contained to the limiting scales. Different power law functions can be better still — the 2D power spectrum for these fields is precisely scale invariant all the way to the exact pixel and domain scales, because the method I used to generate it ensures this by construction.

This example illustrates that the very concept of "scale invariance" can be function-specific when we are considering discrete fields. As discussed above, a theoretically infinite field is scale invariant if and only if all statistics are power law functions, but real images and computer simulations are not theoretically infinite, and so "scale invariance" becomes a much more subtle concept than it already was. The range over which a real field is scaling depends not only on the field itself and the method you used to collect data, but also the method you use to evaluate scale invariance.

The statistical exponents

The exponents I will consider fall into two categories: those that are computed from a continuous field (statistical) and those that are computed from a binary field (geometrical). At least for multiplicative cascades, the continuous exponents are, in some ways, much simpler. There are only three that are independent: the Hurst exponent $H$, the intermittency exponent $C_1$, and the multifractal index $\alpha$. I've written before about how to interpret these exponents, and I think the best thing to do is to simply look at a bunch of simulation output to gain intuition. I have a nice tool for this purpose on my website.

Theoretically, these exponents are nice because they can be simply related to two basic statistical characterizations of a field: first, the moments of the field values, which tell you all there is to know about the probability distribution of the field values, and second, the moments of the two-point statistics, such as autocorrelation functions and power spectra. Two-point statistics contain all of the information about correlations within the field as a function of scale, which is crucial for scale invariant fields which are often correlated at all scales. These two categories include almost everything you could want to know about a given field, but they are fairly abstract and often difficult to measure directly. Thus, one of my main goals here is to look for possible empirical relationships between these field exponents and more commonly measured geometrical exponents.

The three multifractal exponents also serve as the inputs to the main model I'll use here to generate simulations, which is Lovejoy and Schertzer's Fractionally Integrated Flux as implemented in my Python package scaleinvariance. Unfortunately, the model is not perfect and its scaling behavior near grid scales is poor for all measures of scale invariance. To account for this, I generate periodic simulations of shape $16384^2$ and then coarsen them by averaging subsequent $4\times 4$ pixel regions. Although no discrete field is perfectly scaling to grid scales for all metrics, this coarsening method improves scaling for at least some measures of scale invariance. I'll return to this point below.

Each of the resulting $4096^2$ fields represents one statistically independent realization of a scale invariant field that I use to estimate the exponents. For each combination of the input parameters, I generate 128 fields, and unless stated otherwise, all the exponents I report were estimated using the full 128-field ensemble.

A notable special case of Lovejoy and Schertzer's multifractals is when $C_1=0$ (and $\alpha$ no longer influences the result), which reduces to a class of fields that are called "fractional Brownian motion". For these fields, I do not use the Fractionally Integrated Flux, but rather a spectral synthesis method which produces output that is perfectly scale invariant as measured using the radially averaged power spectrum, and so it's not necessary to generate larger simulations and then coarsen them. This allows for larger simulations to be analyzed — specifically, for some experiments I will consider a 32-member ensemble of fBm fields of size $8192^2$.

The main difference between fractional Brownian motion and multifractal fields is intermittency. Highly intermittent fields (with large $C_1$) will have more extreme values, and those extremes will be relatively clustered. Clouds are highly intermittent.

Of course, the parameters $H$, $C_1$, and $\alpha$ are not only input parameters but can also be estimated empirically from the simulation output. This presents another great example of how discretization affects scale invariance, this time for statistical exponents rather than the geometrical exponents I discussed in "Exponents in finite domains" above. The simplest is the Hurst exponent, which can be measured in many ways that theoretically produce the same answer but in practice will not, due to discretization.

Here's an example. The most common method to measure the Hurst exponent is using a structure function, defined through field differences as a function of separation distance:

$$S_{\text{structure}}=⟨|\Delta f(r)|⟩\Delta f(r)=f(x_i+r)-f(x_i)$$

where $f$ is the continuous field, $r$ scale, and the angle brackets denote averaging over many spatial locations $x_i$. Here, I average over all possible spatial locations. We can compare this to Haar fluctuations, which are defined through differences in consecutive field averages as a function of the averaging distance:

$$S_{\text{haar}}=⟨|{\Delta }_{H}f(r)|⟩{\Delta }_{H}f(r)={\overline{f}}_{r/2}(x_i+r/2)-{\overline{f}}_{r/2}(x_i)$$

where the overbars just denote averaging over the first or second halves of the interval $(x_i,x_i+r)$, i.e. ${\overline{f}}_{\ell }(x)=\frac{1}{\ell }{\int }_x^{x+\ell }f(x^{\prime })d{x}^{\prime }$. In a theoretically infinite domain with input Hurst exponent $0\le H\le 1$, both functions scale as the separation or averaging distance raised to the Hurst exponent through

$$S\propto r^H$$

But look at what happens for a field with finite resolution, here a set of $32$ fBm fields of size $8192^2$:

Both are relatively scale invariant for scales between 100 and 1000, but the Haar fluctuation is nearly perfectly scaling to grid scales while the structure function deviates. The discrepancy arises because of the distinction between differences and averages. At any single scale, the Haar fluctuation "picks up" less variability than the structure function, which is why the Haar curve is below the structure function curve in the left plot. Exactly why this happens is interesting to consider, but in any case, I think it's intuitive that there should be at least some difference between the curves given that the definitions are so different.

Consider what happens at the smallest possible scale. The structure function takes the difference between consecutive points. The Haar fluctuation "averages" over the first point, and then "averages" over the second point. But "averaging" doesn't do anything because the average is taken over a single point. So there is, in fact, no difference in these two functions at the smallest possible scale - and as expected, the y-values are equal in the left plot above. A difference only arises at larger scales where there is more than one point to average.

So in a sense, at grid scales structure functions and Haar fluctuation functions are "pulled" to the same value as they approach grid scales. They are no longer parallel and thus one or both are not a power law.

It turns out that simulations produced using the spectral synthesis method are relatively invariant to coarsening via averaging. I used $8192^2$ simulations here, but I could, in principle, have made $65536^2$ simulations and then coarsened them by averaging $8\times 8$ consecutive blocks. This would have made very little difference to the statistics of the resulting field. In particular, the Haar fluctuation would remain almost perfectly scaling to grid scales because it is defined with respect to averages, and the coarsening method also uses averages.

But you can instead coarsen a field via subsampling, not averaging. For example, here I took the set of $8192^2$ fBm fields and coarsened them by simply extracting every eighth pixel in both dimensions to obtain a $1024^2$ grid. The result:

In this case, it is the structure function that is almost perfectly scaling to grid scales. The reason is the coarsening method. Structure functions are defined with respect to values, not averages, and this second coarsening method extracted values, not averages, and so the grid-scale variability is what is expected by the structure function but not Haar fluctuations.

This point can have very practical consequences. For some types of data, "pixels" represent an implicit average, while for others, they represent values at a point. Pixels from an image from a satellite, for example, collect photons from the entire region covered by the pixel and thus represent an average, so a Haar fluctuation would be best suited for these. Financial data, in contrast, typically represent a value at a single moment in time (e.g. close prices) rather than an average over that interval. Structure functions would be better suited here. Picking the correct scaling function for your specific data can mean giving yourself a $10\times$ or $100\times$ wider scaling interval available for fitting a Hurst exponent to. That is a huge deal in many datasets.

The geometrical exponents

To go from the statistical exponents, which are computed for a continuous-valued field, to the geometric exponents, you have to first convert your field to an object. Fields and objects are fundamentally different: for a field, every point in space has a continuous field value, but for an object, every point is simply part of an object or not part of it. An object is a set of points; a field is a set of points plus a field value at each point.

The simplest way to go from field to object is to apply a threshold such that every point in the field with a value higher than your threshold is part of your object, while any point with a lower value is not part of your set. The resulting object is called the exceedance set.

But this simple definition is almost never used in cloud physics. The closest thing that is used is called the level set, which is the set of locations where the field equals the threshold, rather than being greater than or equal to it. This is used for some fractal dimensions that I'll return to, but it's still uncommon.

What is common is to carve up the field into a large number of distinct objects based on two things: a field threshold and a connectivity definition. This is an attempt to formalize our intuitive notion of "clouds," where clouds are entities with unique identities. I'm not a fan of this idea on physical grounds but it's extremely common and it's going to give you a bunch of new exponents to play with.

If we must have individual clouds, then my favorite rule is to apply a constant threshold, which gives a binary cloudy/not cloudy image, and then define cloudy pixels as "connected" when they are adjacent. A set of connected pixels constitutes a single cloud.

Not all studies use this simple connectivity rule, however. Observational constraints often cause the threshold to effectively vary across the image, for example because the sun angle makes thin clouds more or less visible, even holding their water content fixed. Additionally, many studies apply fancy connectivity thresholds that are designed based on an implicit notion of what "truly" constitutes a cloud, such as dividing clouds if they are connected by only a few pixels. In any case, I will consider the simple method of a constant threshold and adjacent connectivity here, as it's the simplest to understand theoretically.

Because the connectivity definition adds complexity and the simpler level/exceedance sets are better understood theoretically, I'll consider level and exceedance sets first.

Fractal dimensions for level or exceedance sets

The simplest and oldest fractal dimension is probably the box dimension. The box dimension is defined through the procedure used to calculate it: one tiles "boxes" on top of the binarized cloud field of varying size, and for each box size, measured by the box side length $\epsilon$, you count the number $N$ of boxes required to cover the set of interest, such as cloud edges. The box dimension is then defined through

$$N\propto {\epsilon }^{-D_{\text{box}}}.$$

Notice how the box dimension could in principle be obtained for whatever set you want, whether it is the level or exceedance set or even some single, isolated cloud defined with an additional connectivity rule.

I bring up box dimensions here mostly for historical and explanatory reasons. It's the most common way to introduce the fractal dimension and thinking about it can make theoretical arguments easier. However, computing it using the actual method I described (overlaying boxes, counting them, etc.) is not a good idea. As I showed above in "Why exponents?", the box dimension does not converge to its theoretical value even for an $8192^2$ simulation, which is already a far higher resolution than most satellite images.

The convergence issue motivated Grassberger and Procaccia to introduce the "correlation dimension" in 1983.⁹ The correlation dimension is again defined through the method used to calculate it. You first compute all distances between all pairs of points in the set, and then compute a cumulative distribution function $C(l)$, which represents the number of point pairs that are separated by a distance less than $l$. The dimension $D_c$ is defined through

$$C(l)\propto l^{D_c}.$$

This is also sometimes explained by describing the following process: first, draw a circle at each set point and count the number of other set points within the circle as the circle radius $l$ changes. Then, repeat by centering your circle at every other set point and sum the count across all circles. The final summed result is two times the cumulative distribution function $C(l)$, where the factor of two comes from double-counting point pairs. Both of these definitions can be useful to keep in mind.

As I showed above, the correlation dimension converges much more rapidly to this power law form in a finite domain than a box dimension would, and it's therefore much more useful. Its value is not always theoretically identical to the box dimension, however, and it can be computationally expensive given that the number of point pairs is $n(n-1)/2$ for a number of points $n$, which has computational complexity $n^2$. Compare this to the box dimension which has computational complexity $n$, as it simply counts the number of points in a box rather than point pairs. This increased computational complexity is the flipside of the improved convergence of the correlation dimension.

The box and correlation dimensions are, in fact, special cases of a more general algorithm. Consider that the number of point pairs within a given subdomain is $n(n-1)/2$, which for large $n$ scales as $n^2$. Thus, at a given scale, the correlation integral can either be thought of as counting pairwise distances, or alternatively counting set points within the circles and then squaring the result.

This observation suggests a natural extension: rather than raising the count in each circle to a power of $2$, we could in principle raise the count to any power $q$. Conveniently, setting $q=0$ will then give a value of $1$ when a subdomain contains a set point or $0$ otherwise given that, in this case, $0^0=0$,³ which sounds like the box dimension. This resulting family of dimensions are called Rényi dimensions.

For defining subdomains, we could use a "box-dimension-like" subdomain, namely non-overlapping squares, or a "correlation-dimension-like" method, where subdomains are circles centered at any set point. Respectively, the resulting dimensions are called box dimensions $D_{q,r}$ and sandbox dimensions $D_{q,s}$, defined through

$$⟨n^q{⟩}_{\epsilon }\propto {\epsilon }^{{\tau }_{q,r}};{\tau }_{q,r}=D_{q,r}(q-1)$$$$⟨n^{q-1}{⟩}_l\propto l^{{\tau }_{q,s}};{\tau }_{q,s}=D_{q,s}(q-1)$$

where the average $⟨{⟩}_{\epsilon }$ is computed over boxes of the same side length $\epsilon$, i.e. holding $\epsilon$ constant, and $⟨{⟩}_l$ is computed over circles of the same radius $l$. For the sandbox dimension, we require $n^{q-1}$ rather than just $n^q$ to account for the fact that we choose to only center circles at set points, which biases the average relative to the box dimension where we consider a full grid of boxes that cover the domain, regardless of where the cloud edge points are. In objscale, both methods are implemented.

The factors of $(q-1)$ in the exponent $\tau$ are necessary to match the original definitions of the dimensions such that $q=0$ gives the box dimension and $q=2$ gives the correlation dimension. If you take an appropriate limit, even the $q=1$ case gives a named dimension, the "information dimension".

For some objects, like level sets of fBm fields, the Rényi dimensions are not a function of $q$, at least for infinitely large domains. But in general, and particularly for level sets of multifractals, higher values of $q$ will more strongly weigh those portions of the domain that are relatively more dense. This is analogous to how higher-order structure functions weigh rare but intense regions more strongly, and both therefore are related to intermittency.⁴

In the limit of an infinite domain, the Rényi and sandbox dimensions converge, but in finite domains convergence is faster for the sandbox dimension and for higher values of $q$. Both of these effects contributed to our prior observation that the correlation dimension is more scale invariant than the $q=0$ box dimension.

Correlation dimension and domain boundary effects

In a finite domain, unless the effect is corrected for, some of the sandbox circles will inevitably extend beyond the domain boundary and into a no-mans land where there are no set points. In DeWitt et. al 2026, we argued that this should be expected to bias the correlation dimension low, given that there are set points that are "expected" (i.e. would exist in an infinite domain) but not seen or counted. In objscale I implemented a fix for this, toggled on or off via the flag interior_circles_only, which restricts the admissible circle centers to the central portion of the domain such that the circles are never large enough to extend beyond the domain boundary. If this parameter is set to False, any edge point is admissible, regardless of its proximity to the domain boundary, and so some circles will extend beyond that boundary.

This argument still makes a ton of sense to me. However, using the $8192^2$ fBm domains here, my reasoning appears to simply be wrong. Compare the empirical exponent for these domains to the theoretical infinite-domain value represented by the dashed line:

The empirical exponent is clearly closer to its theoretical value over a wider range of scales when we do not restrict circle center locations in almost all cases. It seems that my "fix" actually biases the resulting exponent high in these cases, particularly at larger scales. Honestly, I cannot understand why this is the case. But given this result I'm setting the default in objscale to not apply any circle-center location filtering, and all the other results in this post do not apply such filtering.

Exponents for individual clouds

Let's move on from simple level and exceedance sets to objects defined using both a threshold and a connectivity criterion. The most notable difference is that, with the connectivity criterion, in a given image we no longer consider a single object, as before, but rather a set of distinct objects, namely individual "clouds". This enables us to consider not only their structure as a function of scale but also their number or frequency as a function of scale. These are the size distribution exponents: for any metric of cloud size $s$, scale invariance suggests that the size distribution should follow a power-law:

$$n(s)\propto s^{-(1+\tau )}.$$

Here, I'll consider four options for the object size metric $s$. Two are straightforward: cloud area $a$, defined as the number of connected pixels, and cloud bounding box width $w$, the length along one dimension of the smallest bounding box that covers the cloud. You could also consider the other dimension of the bounding box (height/length), but the fields I consider are isotropic, meaning the statistics along each direction are identical.

The final two variables for $s$ are more subtle, and both represent the perimeter of the cloud. The first, "summed" perimeter $p_s$, is the perimeter along the exterior edge of the cloud plus any perimeter along any edge contained within the cloud. By "contained within" the cloud, I mean a boundary that demarcates a cloud hole. A donut-shaped cloud has a single value for summed perimeter. In contrast, "nested" perimeter $p_n$ treats every individual cloud boundary as a distinct value, whether it constitutes a cloud's exterior edge or the edge of a cloud hole. For this variable, a donut-shaped cloud is responsible for two unique values of $p_n$, one for its exterior and one for its hole.

In the last section, we mainly considered fractal dimensions of cloud fields overall, defined by a level or exceedance set, but we also mentioned they could be defined for individual clouds too. We could even define a box dimension for a single, isolated edge contour, such that the boundary of any cloud hole would not contribute to the fractal dimension. The simplest way to do this operationally is to "fill" any cloud hole before computing the dimension, which simply removes boundaries associated with holes so that they do not contribute to the dimension.

In DeWitt et. al 2026, we introduced the term "ensemble" dimension to indicate a level set calculation done over all clouds in the cloud field at once — an ensemble of clouds — and "individual" to indicate a dimension computed for a single isolated cloud defined through a level set plus a connectivity rule. We now have to further distinguish between "filled individual", indicating a dimension computed for an isolated, hole-filled cloud, and "unfilled individual", indicating a dimension computed for an isolated cloud that does not have its holes filled.

Some people might argue that filling holes "misses an important physical property defining the cloud field". But if you are already isolating out a single cloud and ignoring its neighbors for the purposes of computing a fractal dimension, is it really that much less physical to ignore the cloud holes? In my view it's already physically unrealistic to turn clouds into objects in the first place. But in the interest of simply characterizing the geometry of a thing, it makes sense to me to either consider the fractal dimension of all edge contours in an image, i.e. the ensemble dimension of the level set, or a single edge contour, i.e. the filled dimension, not some arbitrary mix of selection of some but not all edge contours, which is what the unfilled dimension measures. Regardless, as I'll show, any of these choices results in different numerical values for each of the dimensions, so they should be considered as distinct things. Too many exponents!

For individual box dimensions, there is an alternative methodology which is theoretically equivalent to the method described above but is less strongly affected by discretization effects. I won't derive it here, but there is a derivation in DeWitt et. al 2026. For the individual unfilled dimension, the method starts by computing summed perimeters and bounding box widths for each cloud. For the individual filled dimension, you would fill cloud holes first and then compute filled widths $w_f$ and perimeters $p_f$. Note that filling holes does not change the bounding box width so $w_f=w$. Once you have a list of widths and perimeters, the exponent can be fit using the form

$$⟨p_f⟩\propto w^{D_{i,f}}$$$$⟨p_s⟩\propto w^{D_{i,uf}}$$

where $p_f$ are perimeters of clouds after their holes are filled, $p_s$ is the sum of hole and exterior perimeter for individual unfilled clouds. Averages are taken over narrow bins of bounding box width. The "individual filled dimension" $\left(D_{i,f}\right)$ and "individual unfilled dimension" $\left(D_{i,uf}\right)$ can then be estimated from a linear regression to a scatterplot of mean filled/summed perimeter as a function of width.

If you've made it this far you are probably more familiar with a different relationship, namely one relating perimeters and areas. But as we argue in DeWitt et. al 2026, for the perimeter-area or perimeter-width exponent to be interpretable as a box dimension of the cloud's boundary, the cloud length metric (either width or the square root of area) must be resolution invariant. We argued that this requirement is satisfied by cloud bounding box width and filled area but not unfilled area. Thus, in the paper we recommended sticking with the individual filled dimension, which can equivalently be computed using either of these equations:

$$⟨p_f⟩\propto w^{D_{i,f}}$$$$⟨p_f⟩\propto a_f^{D_{i,f}/2}$$

where, here, each variable is computed for hole-filled clouds. These are equivalent if we assume $w\propto a_f^{1/2}$, which is implicitely an assumption about the dimension of the set of cloudy points for an individual cloud. If that set has dimension 2, then the area is resolution invariant as required. If that set is fractal with dimension less than 2, the exponent relating unfilled perimeters and areas would be some mix of the dimension of cloudy points and cloud edge points. Filling holes is a simple way to enforce the cloudy-point dimension to be equal to 2, fulfilling the requirement of the resolution independence of cloud area.

The argument was inspired by Lovejoy & Schertzer 1990¹⁰ who argued $w\propto a^{1/2}$ holds only for fBm but not a multifractal. They argued that an unfilled perimeter-area dimension would in fact be some mix of the dimension of the cloudy points of an individual cloud and the dimension of the perimeter of the cloud.

Empirically, their argument seems correct but the effect is subtle. To evaluate this I compared cloud width $w$ to both $⟨\sqrt{a}⟩$ and $⟨\sqrt{a_f}⟩$, where an exponent of 1 would indicate that the set of individual cloudy points are non-fractal:

The cleanest case is comparing the $H=0$ fBm to the multifractal with $\alpha =1.3,C_1=0.1$. Those cases show that the filled vs. unfilled exponents differ by 0.029 for fBm but 0.036 for the multifractal. Overall this is a tiny effect that may be indistinguishable from noise, but it is in the direction you would expect, with the fBm cloudy points showing a dimension closer to 2 relative to the multifractal. This empirical distinction disappears completely when you consider the $H=0.3$ case, where the two exponents are nearly identical. The implications are that Lovejoy & Schertzer's argument may depend on the value of $H$, the distinction may not matter in practice for finite datasets, and overall the problem is more subtle than we suggested it was in DeWitt et. al 2026.

Regardless of the mechanism, the key empirical result of DeWitt et. al 2026 was that the unfilled perimeter-area relationship implies a larger dimension than the filled relation, and moreover the unfilled relation is less scale invariant. This is also true in both fBm and multifractal fields. Note how much flatter the local dimension curves are for both the filled area and width relationships, when compared to the unfilled relationship:

Based on a similar observation in cloud data, where the filled dimension is more scale invariant than the unfilled dimension, we argued that the filled dimension is more interpretable as a true "dimension" and is therefore preferrable to the more common unfilled dimension. I still agree with this result but for a slightly different reason. Based on the analysis here, it appears that the primary reason that the unfilled and filled dimensions differ is not that the dimension of cloudy points themselves is affected, as we suggested in DeWitt et. al 2026. Rather, filling holes affects the dimension of the cloud edge points by omitting hole perimeters from the analysis. The set of cloud edge points is a different set, with a different dimension, based on whether you include edge points at hole boundaries.

If one is measuring fractal dimensions, it makes sense to either consider the dimension of the level set of the entire field, i.e. the ensemble dimension, or the dimension of a single closed contour defined by the level set, i.e. the individual dimension. The fractal dimension of an isolated cloud containing holes is some strange mix, where you consider more than one closed contour (one for the exterior cloud boundary, and one for each hole), but not all closed contours together. Calculating the individual fractal dimension using hole-filled clouds still appears to be the best method to compute the fractal dimension of isolated contours.

General methods for exponent estimation

By far the most common method for estimating all of these exponents is to perform a least-squares linear regression to the logarithm of the empirical scaling function vs. the logarithm of the measure of scale. But this requires deciding which scales to include in the regression. It's best to fit over a relatively wide range of scales, but you must also make sure to exclude the regions most affected by discretization at the small end or finite domain effects at the large end.

Unfortunately this usually comes down to a subjective judgement call. There are sometimes principled ways to choose the regression limits — Tim Garrett and I wrote a whole paper about one method for size distributions — but each specific exponent is likely to need a specific method for choosing limits, and the method may need to be tailored to multifractals in finite domains. It does seem to me that it would be worthwhile to derive better estimators for some of the more important exponents like the Hurst exponent.

Here, mainly because of the variety of exponents I'm considering, we are stuck with general methods like linear regression. But there still remains a question of how to pick the relevant range of scale. Luckily, some of the exponents have known theoretical values, and I can use this fact to test the accuracy of various estimators.

The simplest general way to define a restricted range of scale is to use some fraction of the total available range, centered logarithmically. I'll compare the middle 90%, 75%, 50%, 25%, and the central point. For example, a domain with minimum scale 1 and maximum scale 100 has a span of two orders of magnitude. For the 50% threshold, I would perform a regression to the central one order of magnitude, corresponding to scales between $\sim 3$ and $\sim 30$. For the 75% threshold, it would be the central 1.5 orders of magnitude, and so on. Finally, by "central point" I mean the smallest possible "linear regression" to the two data points closest to the center, again logarithmically spaced. This is equivalent to $d\log \varphi /d\log x$, evaluated using a finite difference at the central point, for some scaling function $\varphi$ and some measure of scale $x$.

Here I'll consider a set of metrics which have known theoretical values. For FIF fields, these are the structure function exponent ${\xi }_{sf}(q)$ and Haar fluctuation function exponent ${\xi }_h(q)$ evaluated for various orders $q$, and the exponent of the power spectrum $\beta$. For the fBm fields, I consider the Rényi dimensions, evalued for several orders $q$ and using both box and sandbox methods. All Rényi level set dimensions equal $2-H$ for fBm fields, regardless of threshold.

As a very coarse overall measure of the accuracy of the fitting methods for all these metrics, below is a table of the mean absolute difference between the theoretical and empirical exponents, aggregated across all of the exponents above.

Family	90%	75%	50%	25%	center	N
Haar	0.1032	0.0992	0.0961	0.0957	0.0947	162
Rényi (box)	0.1300	0.1231	0.1173	0.1097	0.1097	420
Rényi (sandbox)	0.1203	0.1151	0.1095	0.1057	0.1019	420
Spectrum	0.0520	0.0537	0.0650	0.0647	0.0741	27
Structure Function	0.1544	0.1494	0.1437	0.1411	0.1402	162

As might be expected, estimators generally get better as we consider scales more and more removed from grid and domain scales. Given that the observed fitting bias for the central 25% and the center two points are very similar, and that fitting to only two points is more brittle, going forward I will fit all exponents using the central 25% of the logarithmically-spaced fitting window.

I should note that, here, I am in a somewhat unusual situation where all of the fields I consider are known to be theoretically scale invariant, so I'm free to pick the exponent estimator by simply minimizing error. In a more realistic case, you may not know for certain whether the system you consider is scale invariant in the first place, even without discretization effects. You would then want to check that a power law function is appropriate, which typically requires checking linearity over a relatively wide range of scales, in which case you may want to fit over a wider range of the available scales than I do.

Additionally, it's worth noting that the overall error shown in the table above is quite high in an absolute sense. I'm using theoretically perfect simulations of size $4096^2$, much higher resolution than most geophysical data, and yet the departure from theory is large enough that it will be difficult to distinguish many of the exponents such as individual fractal dimensions of various types.

A thought I had was that, perhaps, the spread of the slope plot, defined as the standard deviation of $d\log F/d\log x$ for scaling function $F$ and size metric $x$ over the fitting interval, could be used as an estimate of the bias in the measured exponent. I tried this, and it basically does not work. Too bad.

Statistical error and parameter uncertainty

The overall magnitude of the departure from theory I showed above makes one wonder how large a domain would have to be to reliably estimate the exponents. Two closely related questions in this vein are, first, how many realizations of a given multifractal parameter set are required for the empirical parameter to converge to its infinte-ensemble value? And second, how one might estimate statistical uncertainty arising from analyzing a finite sample size?

All of the analysis methods in both packages I'm using, scaleinvariance and objscale, return estimates of uncertainty, which are simply derived as $2\times$ the standard error of the linear regression. The truth is, I feel somewhat conflicted about this, because these standard statistical packages assume that each data point in the scaling function is statistically independent of the rest. This is also true of fancier methods such as Maximum Likelihood Estimation (MLE) which are commonly recommended.⁵

The data points that make up a scaling function are not statistically independent if they are derived from a scale invariant field. The simplest example is the presence of one very large cloud within a small domain. If the cloud is large enough, it could make it impossible for another cloud to fit in the same domain, which violates statistical independence.

Violation of statistical independence is fairly straightforward to show empirically. Take a size distribution exponent for unfilled cloud areas as an example. If I compute this exponent for a single realization of a multifractal, there will be some statistical noise in the probability distribution because, say, clouds of a given size happened to be more common in that realization than they would be in expectation. If such noise was independent, the uncertainty calculation in linear regression packages would be well calibrated. Specifically, it would be calibrated in the sense that the spread represented by the reported uncertainty would be close to the spread in actual estimates of the exponents if I were to run the experiment a large number of times.

To show this, for three subsets of the FIF parameter space, I generated 1024 realizations of the $4096^2$ fields. Each of these realizations is truly statistically independent. For a varying sample size $n$, I randomly selected 32 samples of $n$ fields out of the 1024 population with replacement, representing a standard bootstrapping procedure. I then computed three different exponents for each of the 32 samples and computed the standard deviation of the 32 exponents for varying $n$ between 1 and 128. This represents a truly well calibrated estimate of the uncertainty, which I can then compare to the uncertainty returned by the objscale or scaleinvariance function for one random selection of $n$ fields, representing the (inaccurate) linear-regression-derived uncertainty. Here is the result:

As you can see, in almost all cases spanning the three exponents and three input fields, the bootstrap uncertainty, which represents the "true" uncertainty, is much higher for small $n$, i.e. small dataset sizes. In some cases, the linear regression-derived uncertainty is actually an overestimate relative to the true uncertainty for large dataset sizes. Interestingly, the linear regression uncertainty for the area distribution exponent appears relatively close to the bootstrapped uncertainty, although the overal magnitude of the uncertainty is much higher than the other two exponents.

Bias due to discretization

A crucial point is that the statistical error discussed above is entirely distinct from bias arising from discretization. On the one hand, I could have, say, a billion fields of size $100^2$ and compute exponents using all of the fields, which would result in very low statistical error but very high discretization bias. On the other hand, I could in principle generate a single $1,000,000,000^2$ field, where a computed exponent could have very little discretization bias but relatively higher statistical error since it's a single array. For most of the experiments I'm reporting here, unless stated otherwise I chose 128 fields of size $4096^2$ as an attempt to balance statistical error with discretization bias given the computational constraint of me being impatient.

Let's consider the second effect, discretization bias, more closely for fBm fields, which are less memory intensive to simulate and therefore have a higher upper bound for what I can consider with my 64 GB of available RAM. These also have theoretically known values for the correlation dimension and Hurst exponent. Below, I computed three exponents for an ensemble of fields of varying sizes, ranging from $256^2$ to $16384^2$. As an attempt to keep statistical error constant, I decreased the number of arrays used as the grid size increased while holding the total number of grid points constant.⁶

The Hurst exponent shows the clearest convergence behavior, smoothly approaching its theoretical value over the full range of array sizes. But even at $16384^2$, the estimated value remains below the theoretical value by about 0.005. The area distribution exponent is next, which shows possible convergence around size $4096^2$ (although in this case I don't know the theoretical value), and the correlation dimension is surprisingly consistent across array sizes, suggesting it converges even for very small arrays. There seems to be a running theme of the correlation dimension being one of the most accurate and well-converged metrics I'm considering here.

Bias due to an imperfect numerical method

A third, independent source of error arises from the accuracy of the method by which I am generating fields in the first place. I mentioned above that fields generated using the FIF algorithm benefit from being initially constructed at finer resolution, before being upscaled by averaging consecutive blocks. In this article, I generated $16384^2$ FIF fields and coarsened by averaging consecutive $4\times 4$ blocks. You can see that this is independent of discretization error as I mentioned above, because a native resolution $4096^2$ field has the same domain geometry as a coarsened $16384^2$ field, yet the coarsened field has superior scaling properties as I'll show below. In other words, any scaling function computed from a $4096^2$ FIF field is not scale invariant for two reasons: discretization error, which is a fixed property of the domain geometry, and numerical error, which can be reduced by creating larger simulations and then coarsening them.

This may remind you of the discussion in The statistical exponents where I showed that certain coarsening methods pair better with certain scaling functions, with coarsening by subsampling and coarsening by averaging being two limiting cases. This ambiguity cannot account for numerical error in FIF simulations, however, because at native resolution there is no scaling function that is perfectly scaling to grid scales for the FIF. This is not true for the fBm simulations, which have a perfectly scale invariant 2D power spectrum to grid scales by construction, and in this sense have no numerical error at all.

So how much can this numerical error be reduced? Is coarsening by $4\times$ sufficient? The answer, again, depends on the particular scaling function used. This is not only because some scaling functions are implicitely averaging while others are implicitely subsampling, but is also because some scaling functions weigh rare, extreme field values differently. Specifically, let's consider Haar fluctuations and explore how coarsening via averaging reduces numerical error. Recall that the Haar fluctuation function is defined with respect to differences of respective field averages as defined above:

$$S_{\text{haar}}=⟨|{\Delta }_{H}f(r)|⟩{\Delta }_{H}f(r)={\overline{f}}_{r/2}(x_i+r/2)-{\overline{f}}_{r/2}(x_i)$$

This is implicitely a "first order" fluctuation function. A more general definition involves some arbitrary order $q$:

$$S_q=⟨|{\Delta }_{H}f(r){|}^q⟩$$

where the simpler version above is just the $q=1$ case. Considering other values of $q$ applies a different weighting for unusually large or small differences, and this is one way to measure intermittency: a highly intermittent field contains relatively more small-scale but large-amplitude jumps. In any case, there is a theoretical relationship between the scaling exponent ${\xi }_q$ and the moment $q$ for multifractals:

$$S_q\propto r^{{\xi }_q}$$$${\xi }_q=qH-K(q),K(q)=\frac{C_1}{\alpha -1}\left(q^{\alpha }-q\right)$$

These provide a means of evaluating numerical error for FIF simulations as a function of coarsening amount. Here, I'll consider one-dimensional rather than two-dimensional fields given that $4\times$ is the ceiling for what I can consider numerically in 2D. Below, I consider Haar fluctuation functions, and the computed exponents, for an ensemble of 4096 one-dimensional simulations, each of size 4096 and with $\alpha =2$. I will compare native-resolution simulations to simulations that were simulated initially $u$ times larger and then coarsened via averaging consecutive regions of length $u$.

As you can see, $4\times$ averaging is not sufficient, while averaging by $128\times$ is essentially identical to $1024\times$. Emprical exponents for almost all cases are very close to theory, but somewhat worryingly the highly intermittent $C_1=0.3$ case does not seem to converge to the theoretical value for $q=3$. I don't know why.⁷ In any case, numerical error is worse in higher dimensional simulations, and so the 2D results below should be expected to be somewhat problematic, particularly for higher-order statistics.

Empirical relationships between scaling exponents

With all of the methodological details out of the way, we are finally ready to consider empirical relationships between the exponents computed from the ensemble of synthetic fields. I generated 128 fBm and FIF synthetic fields for each combination of $(H,C_1,\alpha )$ values including $H\in \{0,0.3,0.5,0.7,1\}$, $C_1\in \{0,0.03,0.1,0.2\}$, and $\alpha \in \{1.3,1.8,2\}$. To limit numerical overflow errors I avoided some low-$\alpha$ and high-$C_1$ combinations. The full list of parameter combinations is below.

The 27 parameter combinations

$H$	$\alpha$	$C_1$	Type
0	—	0	fBm
0	1.3	0.03	FIF
0	1.8	0.03	FIF
0	2	0.03	FIF
0	1.3	0.1	FIF
0	1.8	0.1	FIF
0	2	0.1	FIF
0	2	0.2	FIF
0.3	—	0	fBm
0.3	1.3	0.03	FIF
0.3	1.8	0.03	FIF
0.3	2	0.03	FIF
0.3	1.8	0.05	FIF
0.3	2	0.05	FIF
0.3	1.3	0.1	FIF
0.3	1.8	0.1	FIF
0.3	2	0.1	FIF
0.3	2	0.2	FIF
0.5	2	0.1	FIF
0.7	—	0	fBm
0.7	2	0.03	FIF
0.7	2	0.1	FIF
0.7	2	0.2	FIF
1	—	0	fBm
1	2	0.03	FIF
1	2	0.1	FIF
1	2	0.2	FIF

Once the fields were generated I computed all the geometrical and statistical exponents mentioned so far for the full set of fields. For the geometrical exponents, I computed them for a threshold equal to the field mean, and also for five percentile-based thresholds: 5th, 20th, 50th, 80th, and 95th. Percentiles were computed across the full 128-realization ensemble.

Level set dimension as a function of threshold and intermittency

The first thing I'd like to consider is how the level set ensemble dimensions vary as a function of the threshold used to binarize the field, as well as the strength of the intermittency of the field. Theoretically, for non-intermittent fBm fields, the value of the level set Rényi dimensions should not depend on either the threshold or the order of the dimension, but for intermittent FIF multifractals the value depends both on the threshold and dimension order.

I am not aware of any explicit theoretical relationship between threshold, Rényi order, and the multifractal parameters $C_1$ and $\alpha$, although I believe one is likely to exist. One point that is important to keep in mind is that the scaling between Rényi dimension and order $q$ is distinct from the moment scaling of field values. Both represent moment scaling functions, but for moments computed for different things: either the count of thresholded set points or field values.

This first plot is only for fields thresholded at their mean:

As might be expected, higher intermittency translates to smaller dimensions for a fixed threshold and Rényi order. Interestingly, this is cleaner when intermittency is measured by $C_1$ rather than $K(q)$. The behavior at a fixed level of intermittency, and across Rényi orders, is more surprising: a higher order does not uniformly increase or decrease the measured dimension. Rather, the value of the Hurst exponent matters: fields with $H=0.3$ show $D_q$ decrease with $q$ while fields with $H=0.7$ show the reverse.

The next plot holds Rényi order $q$ constant while varying the threshold used to binarize the field. I plot the difference between the dimension computed at a 50th percentile threshold to the threshold in question, rather than the dimension itself, in order to isolate the trend between dimension and threshold. As a reminder, theoretically fBm should show no such trend, while multifractals should show decreasing $D_e$ for extreme exceedance fractions — we should expect the curves to curve downward at the left and right extremes, and for this dependence to be stronger for higher values of $C_1$.

For the 20th, 50th, and 80th percentile thresholds, Rényi dimensions are relatively stable, typically differing by less than 0.05. The extreme 5th and 95th percentile thresholds show more variability, with typical difference less than 0.2, but this variation remains minor compared to variability due to changing the Hurst exponent, which changes the dimensions by up to 1. Strangely, there seems to be no consistent trend between increasing $C_1$ and increasing the strength of the threshold dependence, which contradicts the expectation for multifractals. I can only assume that the reason is, again, discretization and numerical errors. We should expect relatively more bias for the relatively extreme thresholds that result in a smaller sample size of set points. Still, I'm surprised we do not see a stronger trend between threshold dependence and the level set dimensions.

Other relationships

I performed a massive brute-force fit to all combinations of the computed exponents to three functional forms:

$$a=x_{1}b+x_{2}c+x_3$$$$a=x_{1}bc+x_2$$$$a=x_{1}b+x_2$$

where $a$, $b$, and $c$ are arbitrary geometrical or statistical exponents. I chose these forms because they are some of the simplest combinations of three exponents that seem nontrivial and would recover several theoretical relationships that I'm aware of. For each unique combination of exponents $a,b$, and $c$, I performed an ordinary least squares fit to determine the constants $x_1$ and $x_2$, and ranked candidates by $R^2$. There are 66 exponents in total, implying nearly 300,000 individual regressions.

The 66 computed and input exponents

Input parameters (5)

Math	Description
$H$	Hurst exponent
$C_1$	Codimension of the mean (intermittency parameter)
$\alpha$	Multifractal index (Lévy index)
$K(q)$, $q\in \{2,3\}$	Theoretical moment scaling exponent (computed analytically from $C_1,\alpha$)

Statistical exponents (7)

Math	Description
$\beta$	Spectral slope of the 1D power spectrum
${\xi }_h(q)$, $q\in \{0.5,1,1.5,2,2.5,3\}$	Order-$q$ Haar fluctuation scaling exponent

Geometrical exponents (54)

Math	Description
$D_{q,s}^{\mathrm{lev}}$, $q\in \{0,0.25,\dots ,5\}$	Sandbox Rényi dimension of the level set (21 values)
$D_{q,s}^{\mathrm{exc}}$, $q\in \{0,0.25,\dots ,5\}$	Sandbox Rényi dimension of the exceedance set (21 values)
$D_{i,f}^{p,a}$, $D_{i,uf}^{p,a}$	Individual fractal dimension from perimeter–area, holes filled / unfilled
$D_{i,f}^{p,w}$, $D_{i,uf}^{p,w}$	Individual fractal dimension from perimeter–width, holes filled / unfilled
$D_{i,f}^{p,h}$, $D_{i,uf}^{p,h}$	Individual fractal dimension from perimeter–height, holes filled / unfilled
$D_c^{\mathrm{f}}$, $D_c^{\mathrm{uf}}$	Correlation dimension of the largest cloud, holes filled / unfilled
${\tau }_{\mathrm{area}}$	Cloud area size-distribution exponent
${\tau }_{\mathrm{perim}}$	Cloud summed-perimeter size-distribution exponent
${\tau }_{\mathrm{width}}$	Cloud bounding-box-width size-distribution exponent
${\tau }_{\mathrm{nested}}$	Cloud nested-perimeter size-distribution exponent

Total: 66 exponents (5 input + 7 statistical + 54 geometrical).

I should note that I am intentionally using a very simple statistical technique here. Since the expressions above have different numbers of free parameters, $R^2$ values are not directly comparable across the different relationship forms. Additionally, even if $a$, $b,$ and $c$ were completely independent and random, some relatively high $R^2$ values would still be seen because we are sampling 300,000 individual combinations. All of this is on top of statistical, numerical, and discretization errors I discussed above. The reality is that I am overall quite skeptical about these results so any interpretation needs to come with a substantial level of caution.

I've included a data explorer here that makes it much easier to browse the exponent relation landscape. The data explorer is divided into three datasets from which exponents are computed. The first is the full set of $4096^2$ fBm and FIF fields shown in the table above where the mean and percentile thresholds used to binarize the fields contribute unique exponent values for the regression. For the second, I still consider all the $4096^2$ fBm and FIF fields, but only the mean threshold, resulting in a smaller number of values for each geometrical exponent. Finally, I consider the $8192^2$ fBm fields in isolation, thresholded at five different places: 5th, 20th, 50th, 80th, and 95th percentiles. This last set is a useful point of comparison because we know that the fBm fields have zero numerical error by construction.

Interactive regression explorer

Analytical relationships between scaling exponents

A recurring theme in this article has been the difficulty of eliminating finite-size, discretization, and statistical uncertainties and biases, even considering the very large simulations I have considered. I have primarily pursued an empirical approach here to emphasize these problems, which are often ignored in more analytical and theoretical analyses but are centrally important when it comes to connecting theory to the real world.

Having said that, it is still useful to have theoretical results that would apply to an infinitely well resolved domain, and so for completeness I will detail those relationships I am aware of here. We can get quite far for nonintermittent fields, and in particular, I will show how any of the above-discussed geometrical exponents can be expressed purely in terms of the Hurst exponent. Unfortunately, few of these results can be expected to hold in the more realistic intermittent case, but we can partly address this gap using the above empirical results.

Size distribution exponents and individual fractal dimensions

To begin, there are some straightforward relationships worth mentioning. The size distribution can be written as

$$\frac{dn}{ds}\propto s^{-(1+\tau )}$$

for any metric of cloud size $s$. If you integrate this you obtain $n\propto s^{-\tau }$, implying

$$\frac{\log n}{\log s}=-\tau +c$$

where $c$ is some constant that may be removed by differentiating:

$$\frac{d\log n}{d\log s}=-\tau .$$

This allows two different size distribution exponents to be related by cancelling $d\log n$. For example:

$$\frac{d\log p}{d\log w}=\frac{{\tau }_w}{{\tau }_p}$$

which, conveniently, relates to our definition for the individual fractal dimension for either filled or unfilled clouds. Those expressions were of the form $p\propto w^{D_i}$, and we can take the logarithm of this equation and differentiate it to obtain

$$\frac{d\log p}{d\log w}=D_i.$$

So for filled or unfilled perimeters, respectively, we have

$$D_{i,uf}={\tau }_w/{\tau }_p$$$$D_{i,f}={\tau }_w/{\tau }_{p_f}.$$

There are also a number of analogous relationships you can derive following a very similar procedure.

Ensemble dimension and Hurst exponent

This next relationship is a well known result for fBm and I'll include it here for completeness.

Consider a 2D fractal field as a surface in three dimensions. At resolution $\epsilon$, the field on average fluctuates a vertical distance ${\epsilon }^H$. This follows directly from the structure function definition, for example. To cover this 3D surface, we need at least one box at every column, and the number of columns is proportional to $(L/\epsilon {)}^2$. We must also stack boxes within a column if the surface height fluctuates more than the box height. Since its average vertical fluctuation of the surface is ${\epsilon }^H$, we need on average $\sim {\epsilon }^H/\epsilon$ boxes stacked vertically, since the surface fluctuates by ${\epsilon }^H$ but each box covers only $\epsilon$. The total number of boxes is proportional to ${\epsilon }^{-2+H-1}$. So the box dimension of this three-dimensional surface is $3-H$. Now, since this is the fractal dimension of a 3D surface, but we want the fractal dimension of a level set (a binary threshold), we cut a plane through the surface to get a fractal object in 2D space. Going from 3D space to 2D space like this reduces the fractal dimension by one, so the box dimension of the perimeters is $2-H$. So:

$$D_e=2-H.$$

Ensemble and individual dimensions

In DeWitt et. al 2026, at great effort we derived a relationship between ensemble and individual fractal dimensions, namely

$$D_e=D_i{\tau }_{p,n}$$

We argued that in this case ${\tau }_{p,n}$ should be for nested perimeters if $D_i$ was computed for filled areas/perimeters. To derive this, we used an advection argument, where changing the resolution could be thought of as "advecting" clouds in perimeter space, i.e. moving each cloud to a different location on the perimeter axis.

Here's a simpler way to think about it. The two dimensions can be defined through a box dimension covering two different sets. The first set is an individual closed contour, such as a filled cloud's boundary, with perimeter $p$:

$$p\propto {\epsilon }^{1-D_i}$$

and second, the level set including contributions from all clouds, giving the total perimeter $P$:

$$P\propto {\epsilon }^{1-D_e}$$

The question is, how do $D_i$ and $D_e$ relate? The first thought is to just integrate the size distribution of the individual contours, each of variable length $p$:

$$P={\int }_{p_{\text{min}}}^{p_{\text{max}}}p\frac{dn}{dp}dp$$

The problem with directly calculating this integral is that perimeter is a function of resolution $\epsilon$, making the integral conceptually subtle, because we ultimately want to allow $\epsilon$ vary to give the box dimension meaning through the box definitions above. This is why we introduced the idea of "advection in perimeter space" in DeWitt et. al 2026. But it is simpler to change variables to a resolution-independent metric of cloud size, for which a good choice is cloud width $w$, as discussed in Exponents for individual clouds. Width is related to perimeter through $p\propto {\epsilon }^{1-D_i}{w}^{D_i}$ (Eqn. 4 in DeWitt et. al 2026). Using this equation, we can change variables to

$$P\propto {\epsilon }^{1-D_i}{\int }_{w_{\text{min}}}^{w_{\text{max}}}{w}^{D_i}\frac{dn}{dw}dw$$

If we assume that this integral is dominated by small clouds, then

$$P\propto {\epsilon }^{1-D_i}{w}_{\text{min}}^{D_i-{\tau }_w}$$

where ${\tau }_w$ is the distribution exponent for cloud widths. Since the smallest resolvable cloud width $w_{\text{min}}$ is proportional to the box size $\epsilon$, this simplifies to

$$P\propto {\epsilon }^{1-{\tau }_w}$$

Comparing this to the defintion of the ensemble fractal dimension

$$P\propto {\epsilon }^{1-D_e}$$

we can see $D_e={\tau }_w$. This relationship is interesting in itself, as it relates a fractal dimension to a size distribution exponent, which might initially seem like two very different things given their different definitions. An analogous result is Korcak's law,¹¹ which relates the area distribution exponent to the fractal dimension, which could be obtained by assuming ${\tau }_w={\tau }_a/2$, using an analogous argument discussed in Size distribution exponents and individual fractal dimensions, implying $D_e={\tau }_a/2$. As I noted in Exponents for individual clouds, this corresponds to assuming the dimension of the set of cloudy points is equal to 2 for an individual cloud. This holds for fBm but its status for multifractals is unclear.

In Size distribution exponents and individual fractal dimensions, we derived ${\tau }_w={\tau }_p{D}_i$, which implies $D_e={\tau }_p{D}_i$. This is the expression proposed in DeWitt et. al 2026.

One advantage of this argument is that the distinction between nested and filled perimeter is more clear. Regardless of which is used, the definition must be consistent in the two key equations $p\propto w^{D_i}$ and $n(p)\propto p^{-(\tau +1)}$. This suggests that there could actually be two distinct relationships:

$$D_e={\tau }_{p,n}{D}_{i,n},p_n\propto w_n^{D_{i,n}},n(p)\propto p_n^{-({\tau }_{p,n}+1)},$$$$D_e={\tau }_{p,s}{D}_{i,s},p_s\propto w_s^{D_{i,s}},n(p)\propto p_s^{-({\tau }_{p,s}+1)},$$

where, in the first equation, all perimeters are defined as nested perimeters, and in the second, all perimeters are defined as summed perimeters. Note that we even have to distinguish between "summed cloud width" and "nested cloud width". This is because, if we are considering cloud hole perimeters as a distinct value in the distribution as in the nested case, then we must also include a value for the width of each corresponding cloud hole. Values corresponding to the width of a cloud hole would not appear in the summed case but are required for the nested case.

Notably, "filled" perimeters do not appear here because filled perimeters omit hole perimeters from contributing to the total $P$. In DeWitt et. al 2026 we implicitely assumed that a dimension computed for nested widths and nested perimeters would equal a dimension computed for filled widths and filled perimeters. It is worth being more explicit about this assumption. In both cases, you are simply using a relationship between the bounding box width and perimeter of individual closed contours, but the "nested" case aggregates statistics using additional closed contours corresponding to cloud holes compared to the summed case. So the equivalence assumption comes down to assuming that closed contours defining cloud holes have the same dimension as closed contours defining clouds themselves. Is this true?

Consider the fBm case. Since positive and negative excursions have identical statistics for fBm, the field is statistically invariant to flipping the sign of the field values. Thus, if we defined "cloud" by regions of positive field values, any dimension or size distribution exponent must have the same value if we defined "cloud" by regions of negative field values, where objects that were previously cloud holes are now clouds and vice versa. The implication is that closed contours defining hole boundaries have the same dimension as closed contours defining clouds. Therefore:

$$D_{i,n}=D_{i,f}$$

This argument breaks down for multfractals, because the statistics are not the same for above-threshold and below-threshold excursions. Some of the other above arguments may also break down for multifractals, particularly because multifractals have fractal dimensions that vary spatially, and this is not true for fBm. A multifractal cloud field will therefore imply different clouds, in different locations, each have different fractal dimensions. But I assumed above that the individual fractal dimension is the same for all clouds. I do not know how to make appropriate corrections to account for intermittency, but it should be noted that some such corrections are likely to be necessary.

Additional results by Claude

I also asked Claude to investigate whether it was able to derive any additional interesting analytical relationships, given a draft of this blog post and a substantial amount of additional context on my work. I also gave it the option to iterate after several rounds of feedback from me and let it decide what to cover and how much effort to invest. The final report took Claude several hours to produce over a few different Claude Code sessions.

Expand Claude's full derivation report

Claude derived a few additional useful relationships. It filled in the size distribution exponent relationships I mentioned above as being "analogous" but did not derive myself. It also provided a relationship between three individual cloud dimensions: namely, the dimension of the set of cloudy points, the dimension computed for unfilled clouds using a perimeter-area relation $D_{i,uf}^{p,a}$, and one computed using a perimeter-width relation $D_{i,uf}^{p,w}$. As discussed in Exponents for individual clouds, when the dimension of an individual cloud's cloudy points is equal to 2, such as in the fBm case, Claude's relation reduces to $D_{i,uf}^{p,w}=D_{i,uf}^{p,a}$. Note that all of these dimensions remain distinct from a dimension for filled clouds.

Claude also derived $D_e=D_i^f{\tau }_{p,n}$, the relation I considered in the last section, in two other ways. It also extended this result in a way I quite like to both cases where small clouds dominate the width distribution or when large clouds dominate the distribution. I considered only the former. Its result is nicely packaged through

$$D_e=\text{max}\left(D_i^f,{\tau }_w\right)$$

where ${\tau }_w=D_i^f{\tau }_{p,n}$ as I had above. In DeWitt et. al 2026 we did not consider the case $D_e=D_{i,f}$ because it does not apply to cloud size distributions as seen from space, but overall this is a nice extension. Intuitively, when a single large cloud dominates the distribution, that cloud's fractal dimension dominates the fractal dimension of the entire field. Claude notes that the crossover point between $D_e=D_{i,f}$ and $D_e={\tau }_w$ is not known analytically, although with the conjecture below it would occur at $H=1$, with $D_e={\tau }_w$ for $H<1$. Claude also derived some statistical exponent relations, such as relating the spectral exponent to the Hurst exponent.

Most interestingly to me, Claude remembered a paper that gives an explicit relationship (that I verified) between individual filled fractal dimension and the Hurst exponent as a conjecture, which I will call the Kondev-Henley-Salinas (KHS) conjecture. This is the final missing link that allows us to represent all the geometrical exponents as a function of the Hurst exponent alone, at least for fBm. As I noted in the introduction, this must be possible in principle given that the statistics of fBm are fully set by the value of the Hurst exponent, but the exact relationships are not trivial to derive, as you can see by now.

Claude also flagged a number of open questions that it attempted but could not close, most of which are intermittency corrections to geometrical/statistical exponent links. Yep, it's a difficult problem.

Geometrical exponents expressed as a function of the Hurst exponent, with empirical multifractal corrections

With the analytical relationships in hand, I'd like to see how well they are respected empirically. Let me begin by summarizing the known theoretical relationships for fBm, starting with the relation derived in Ensemble dimension and Hurst exponent:

$$D_e=2-H$$

Although I only derived this for the box dimension, it also holds for all Rényi dimensions, including the correlation dimension, for fBm. This equation is reproduced relatively well empirically, with some departure near the limits $H\to 1$ and $H\to 0$:

Based on messing with the empirical intermittency correction slider in the above plot, I propose as a guess the following intermittency-corrected relationship:

$$D_e=2-H-\frac{1}{2}{C}_1.$$

In Claude's derivation work, it pointed out that intermittency corrections should be a function of both $C_1$ and $\alpha$, not just $C_1$ as I proposed here. I agree with this, although I do suspect that numerical and discretization error are far too large for this to matter in practice. I did look into this a bit but am not showing it here.

Since we also showed that $D_e={\tau }_w$ in Ensemble and individual dimensions, we have

$$\begin{array}{rl}{\tau }_w& =2-H\\ {\tau }_a& =1-H/2.\end{array}$$

The KHS conjecture⁸ (their 4.3) relates the individual fractal dimension and nested perimeter size distribution exponent through:

$$\begin{array}{rl}{D}_{i,n}& =\frac{3-H}{2}\\ {\tau }_{p,n}& =\frac{4-2H}{3-H}.\end{array}$$

Combining this with $D_e={\tau }_{p,n}{D}_{i,n}$ and $D_{i,n}=D_{i,f}$ from Ensemble and individual dimensions we get $D_e=2-H$ as expected.

In that paper, they also tested these equations empirically, but the paper was published in 2000 and so they could only compute empirical values for arrays up to $512^2$. They also did not consider multifractals.

Overall, I think these empirical results broadly support the KHS conjecture, although a simplistic intermittency correction involving only $C_1$ seems less compelling than the above relation I suggested for $D_e$. There are some substantial departures from their conjecture at the extreme 5th and 95th percentiles, but I would not read too much into these.

Conclusion

Condensed water in the atmosphere is best thought of as a continuous, scale invariant field, but most humans, scientists or not, tend to binarize this field into "cloudy" regions and "clear" regions. Such a binarization forces us to trade statistical exponents, which characterize the properties of a field, for geometrical exponents, which characterize the properties of an object. Further separating out the "cloudy" region into distinct, unconnected objects -- individual clouds -- opens up a wide landscape of possible geometric exponents to choose from.

For the simplest non-intermittent fields, some progress can be made in relating the geometrical exponents to each other and to the statistical exponents characterizing the field, as I show here. But these relationships break down for the more realistic multifractal case, and worse, multifractals allow for an even wider range of geometrical exponents with distinct values, such as dimensions that depend on threshold or Rényi order. Theoretical relationships for multifractals are also much more challenging analytically.

These are purely theoretical challenges. But additional complications arise when it comes to estimating any of these exponents in a domain of finite size, measured at finite resolution, and spanning a finite sample size. These issues bias empirical measurements in a non-negligible way, even for massive datasets like the $8192^2$ fBm I consider here, a dataset that is much larger than typical observational data such as satellite images which tend to have pixel counts in the low thousands. I worry about these issues in particular because they are usually ignored by theoretical analyses.

My aim was to show the importance of considering finite-size effects and to give a sense of the vastness and arbitrariness of the zoo of geometrical exponents that can be defined. There is far more than just "the fractal dimension". At the end of "Toward less subjective metrics for quantifying the shape and organization of clouds" we put it like this:

Although understanding clouds as objects might seem simple and intuitive at first, the hidden subtleties identified here prove surprisingly problematic ... One might reasonably wonder whether the subtleties described here are worth the care required, or whether some continuous field-based approaches might ultimately prove superior (Lovejoy and Schertzer, 2010).

This is why.

---

All code and data used to generate the simulations, figures, and exponent relationships in this post are available on GitHub.

---

1. S. Lovejoy & D. Schertzer, The Weather and Climate: Emergent Laws and Multifractal Cascades, Cambridge University Press (2013). The reduction of multifractal moment scaling to three parameters ($H$, $C_1$, $\alpha$) for universal multifractals is developed there.

   ↩︎

2. Specifically, the exponent is computed from a linear regression to the logarithm of the scaling function to the logarithm of the scale measure. Each point represents a regression computed to scales spanning a half order of magnitude, logarithmically centered at the plotted scale value. This represents the value of the power-law exponent at a local range of scales.

   ↩︎

3. It's not always true that $0^0=0$ because it depends how the limit is reached. In this case, consider taking $n\to 0$ first with finite positive $q$, then let $q\to 0$. Since $0^q=0$ for any positive $q$, the limit $n=0,q=0$, $n^q=0$. If the order of the limits were reversed, then $0^0=1$, in which case the box dimension would not be recovered.

   ↩︎

4. I suspect there is an analytical relationship here, but I haven't been able to figure out what it is. Boxes 5.4/5.5 in Lovejoy and Schertzer 2013 seem to come close, but do not ultimately provide any clear relationship for the level set of a multifractal field. A particular issue for multifractals is that the dimensions depend on the threshold used to binarize the field.

   ↩︎

5. and are truly superior for the independent-data case.

   ↩︎

6. For example, I used 40960 arrays of size $256^2$ and 10 arrays of size $16384^2$, since $40960\times 256^2=10\times 16384^2$.

   ↩︎

7. you may wonder about the divergence of statistical moments, but in this case the critical moment should be $10/3$ for $C_1=0.3$ and $\alpha =2$, which is (slightly) larger than 3. This might be close enough to account for the departure, but I'm not sure. It could also be a result of an error in my implementation, though this would be somewhat surprising given the agreement with theory for all the other cases. This would be worth looking into further.

   ↩︎

8. J. Kondev, C. L. Henley & D. G. Salinas, Nonlinear measures for characterizing rough surface morphologies, Physical Review E 61, 104 (2000).

   ↩︎

9. P. Grassberger & I. Procaccia, Measuring the strangeness of strange attractors, Physica D (1983). This is where the correlation integral $C(l)$ corresponding dimension were introduced.

   ↩︎

10. S. Lovejoy & D. Schertzer, Multifractals, universality classes and satellite and radar measurements of cloud and rain fields, Journal of Geophysical Research (1990).

    ↩︎

11. Korčák's law is an empirical size–frequency relation for geographical objects (islands, lakes) dating to the 1930s. See A. R. Imre & J. Novotný, Fractals and the Korcak-law: a history and a correction, European Physical Journal H (2016).

    ↩︎