This post is image-heavy and may be truncated in the email versions. Consider viewing on Substack.
As a side project, I've been implementing a little-known mathematical method for simulating random fields (below) in Python called the "FIF"1 model. The model was proposed in 1987 by Daniel Schertzer and Shaun Lovejoy as a model for geophysical fields that are shaped by turbulence. Its first applications, for example, were to model rain and cloud fields. It can also make beautiful and almost artistic-looking outputs:
The model never really caught on, which is unfortunate because I’ve become convinced that it is a very good model.2 One of its main advantages is that a simulation can be produced in one go that has the correct statistics for atmospheric fields like rain or clouds. This is very different from traditional hydrodynamic models, which typically require days or even months of model “spin up” where the models slowly transition from a homogeneous initial state to a more realistic and varied state. If the FIF model was further developed for applications, computational efficiency gains would likely be factors of thousands to millions.3 That’s supercomputers-to-laptops.
I wanted to learn how it works and turn Lovejoy’s Mathematica code into a Python package that could be optimized, used more broadly, and maybe even developed for applications. The result is the scaleinvariance package and this visualizer.
I won’t describe the simulation method in detail here. I’ll focus on the results. But as a high-level overview, the FIF model is essentially an opaque mathematical method4 which generates a field with the correct statistics. What are the correct statistics? The basic target is scale invariance: the field has spatial correlations that repeat at every modeled scale, just like the atmosphere. Visually, the field displays structure at all scales like a fractal or a cloud:
On form and pattern in fractal clouds
In fact, the output of the FIF is actually called a “multifractal”, a term I somewhat dislike because it invokes notions of geometric objects. But a field is not a geometric object. A field has a continuous value at every point in space.
So we can think of the FIF as a numerical method to produce any possible scale invariant field. As you might imagine, there is quite a range of possibilities. But incredibly, the full range of possibilities may be described with just three parameters. It is probably easiest to develop intuition about the three parameters by simply looking at examples. For that reason, I put together an interactive website with sliders for the three parameters next to the resulting fields. I recommend playing around with that. But here’s a brief look at some one-dimensional transects of a field:
The Hurst exponent H controls the “smoothness” of the field. More precisely, it answers the question “How much stronger or weaker are larger-scale fluctuations when compared to smaller-scale fluctuations?” Large H implies large-scale fluctuations are dominant, while small H implies the converse. Here it is increasing from H=0 to H=0.4 to H=0.8:
The Levy index α determines how likely extreme jumps are. Lower values (left, α=1.05) create occasional giant spikes, while at higher values (right, α=2) these big spikes are much more rare:
The intermittency exponent C5 determines whether extreme spikes are clustered. Conversely, it determines whether anomalously low regions are clustered. When it is large (right, C=0.2), the field is “intermittent” in the sense that all of the “activity” tends to be concentrated to some small part of the overall field rather than distributed more evenly (left, C=0.001).
It is also possible to specify whether the simulation is “causal” or not. This is useful because the “structure at all scales” requirement implies that each part of the simulation is somewhat correlated with each other part—each point “knows” what the other points are doing. This is great for simulating, for example, the wind field measured along an airplane’s flight path.
But if we want to make temporal simulations, such as the temperature measured at a given station as a function of time, these correlations need to be adjusted. The present can know what happened in the past, but we cannot know what will happen in the future. Compare the below simulations, which have parameters which cause large jumps (low α). The left is acausal while the right is causal:
Moving along the trajectory from left to right, the acausal simulation “anticipates” when the big spike that is about to occur and starts to decrease in advance of the spike (“knowing the future”). But on the right, the spikes are much more sudden. The simulation “doesn’t know” when as spike is about to occur until it happens.
Notes on the simulation accuracy
Again, the target of the simulation is a scale invariant field with structure at all scales. This is a more difficult target than you might think, because any simulation must be performed on a grid with finite grid size. Inevitably, the grid itself imposes structure of its own to the resulting field—a structure that we do not want because it has a single scale: the grid size.
A naive implementation of the FIF model has very bad artifacts from the grid scale. Lovejoy & Schertzer revisited the model in 2010 to address this issue and proposed a method to improve the statistics. I implemented their method, but the simulations are still not perfect. The error is easiest to see by comparing a large simulation, which shows the large-scale characterstics, with a small simulation, which shows the smaller scale statistics. If the simulations were perfect, these would appear identical due to the scale invariant property, which creates the same kind of structure at all scales.
For example, below is a simulation of size 1000 (left) compared to a simulation that was originally of size 100,000, but where I averaged each 100-point chunk to obtain a comparable 1000-point line (right):
As you can see, the line on the left is smoother than the one on the right, even though they should have the same degree of “spikiness” if they were truly scale invariant. Error is also visible in two dimensions (note how much smoother the left image is):
The way to quantify this is to calculate the power spectrum, which measures the amount of structure that exists at each scale. As you’ll see below, Lovejoy & Schertzer’s corrections improve the simulation but it remains far from perfect. Note that, somewhat counterintuitively, it is conventional to plot the power spectrum so that the smaller structures correspond to the right side (not left).
If you use FIF simulations, whenever possible create a larger simulation than necessary (by ~10x or, ideally, ~100x) and then manually coarsen it to your desired size. This helps reduce the numerical error shown above.
Gallery
I’ll highlight some possibilities using my multifractal explorer below.
The “cloud” color scheme is a crude representation of how radiation interacts with water: it jumps quickly from blue to white using a sigmoid function near the median value of the field. Even though this is extremely simplistic, using the correct parameters for clouds (H~0.3, C~0.05, α~1.8) results in a surprisingly realistic image:
Ηere are some different parameters and a different color scale:
Here’s a… moldy mud blob on lava?
Short for “Fractionally Integrated Flux,” but I don’t want to scare anyone with such a technical-sounding term. The point here is not how the model works but what it outputs.
2It is, for example, extremely consistent with all my observational work: this and this and this and this
3If a large eddy simulation with timestep of order 1s requires 1 day of spin up, there are ~100,000 intermediate timesteps which are discarded, none of which are required for the FIF model.
4It is perhaps not opaque to Lovejoy or Shertzer, I don’t know. I say this because the only motivation behind why each step is performed in the algorithm is to obtain the correct statistics. It’s not motivated by a physical argument, and I don’t know what the physical interpretation of, for example, the fractional integrations could be.
5Lovejoy and Schertzer call this the “codimension of the mean” which, while a descriptive name, can be difficult for the uninitiated.