from IPython.display import Image
#---- Plotting
import matplotlib
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
sns.set_context('notebook')
sns.set_style('ticks',
{
'axes.grid': False,
'axes.linewidth': '0.75',
'grid.color': '0.75',
'grid.linestyle': u':',
'legend.frameon': True,
})
xkcd = sns.xkcd_rgb
plt.rc('text', usetex=True)
plt.rc('font', family='Serif')
The data used in this notebook are from a high-resolution, large-eddy simulation model using the System of Atmospheric Modelling (SAM; Khairoutdinov and Randall, 2003). The model was configured to simulate a boundary-layer atmosphere according to BOMEX case data (Holland and Rasmusson, MWR, 1973). The analysis has been performed on two different runs: one at 12.8 km $\times$ 12.8 km $\times$ 3.2 km at 25 m resolution for 12-hours (12-hr run hereafter), and one at 6.4 km $\times$ 6.4 km $\times$ 3.2 km at 12.5 m resolution for 6 hours (12.5-m run hereafter). In both simulations, the model spin-up times were roughly three hours, during which the data are not collected.
The following experiments have been repeated for marine boundary-layer cumulous clouds (CGILS; Blossey et al., 2013) and boundary-layer clouds over a diurnal cycle over land (ARM; Brown et al., 2002) as well. However, seeing that the results do not differ significantly from what is being presented here, they won't be shown for the sake of brevity.
Unfortunately, the analysis had to be repeated multiple times. The suggested project would have been done a model cloud; for example, assuming that we know that the clouds display fractal structure with a known fractal dimension, we could construct artificual shapes with contourted surfaces that will product the right fractal properties according to our assumptions (e.g. Soneira and Peebles, 1978). However, this assumes that (1) individual clouds (as opposed to observed cloud field) behave like fractals and (2) we know exactly their fractal properties beforehand.
Since we performed multiple simulations and found interesting and unexpected results (cf. Section 2.3), I decided to work with the actual measurements taken from the LES simulations. Still, I found that the properties of the simulated cumulus clouds correspond well to earlier studies as well as satellite observations, which assures us that the turbulent properties of the individual simulated clouds are realistic representations of the turbulence during moist convection.
The questions related to this project have been re-arranged slightly, in order to keep this notebook logically organized.
Consider processes as shown in Crum Stull & Eloranta 1987, “Coincident lidar and aircraft observations of entrainment into thermals and mixed layers” (JAS):¶
Compute the fractal dimension of the contorted perimeter of a thermal, and compare it with the fractal dimension of a simple circle.¶
First, we need to ensure that the simulated clouds do indeed have fractal structure. Previous studies using simulated cloud fields (Siebesma and Jonker, 2000) and satellite observations (Lovejoy, 1982; von Savigny et al., 2011; Brinkhoff et al., 2015) have shown that this is indeed the case, but it should be meaningful to show that individual cloud volumes and surfaces have fractal structures.
In order to track the individual cloud histories, we applied the cloud-tracking algorithm initially developed by Dawe and Austin (2012). I have since modified this algorithm extensively, but the heuristics used in determining the tracked volume remains consistent with the original work. Using the cloud-tracking algorithm, we can obtain detailed information about the properties of individual clouds as they develop. For the sake of this work, we focus on the fractal nature of the individually observed clouds.
A generic cloud property $\psi(x, y, z, t)$ obtained from the model output is normally integrated over the horizontal cloud region at every time step, giving a simple two-dimensional variable $\psi(z, t)$. That is, we evaluate the individual cloud properties at every height and every time step. The distribution of $\psi$ are therefore not entirely (and statistically) independent to each other, but the posterier distribution should still be statistically significant given a large enough dataset.
It is surely easier to simply take the projection of the 3D cloud volumes and reproduce satellite measurements (e.g. Lovejoy, 1982). However, this suffers from the discretization effects as well as the overlapping of cloud volumes. When these projections are examined at the individual cloud level, we found that the errors introduced by taking the projection renders the dataset inadequate to represent turbulent properties of moist convection. Still, it works on the scale of the entire cloud field, which you can see in Section 2.3 (Figure 7).
We use the box counting method to calculate the fractal dimension $\mathcal{D}_\mathrm{box} \approx \mathcal{D}_f$ (hereafter $\mathcal{D}_f$ for brevity), also called the Minkowski–Bouligand dimension.
Given a fractal set $\mathcal{S}$, the box-counting dimension can be calculated by counting how many boxes (grid cells in our case) are required to cover the set, and observing how this number changes as the size of the box changes. If $N(l)$ is the number of cells with spatial scale $l$ needed to cover the whole set $\mathcal{S}$, then the box-counting dimension can be defined as:
$$ \mathcal{D}_f = - \lim_{l \to 0} \frac{\log_{10} N(l)}{\log_{10} l}. $$
where $\log$ does not need to be of base $10$, but has been chosen so that it is consistent with the rest of the calculation (and with the literature).
In practice, since we are dealing with the results from large-eddy simulations, the spatial scale $l$ cannot be smaller than the resolution of the LES model $l_0$. In this regard, we can assume that the results are accurate up to the model's resolution ($l_0 = 12.5$ m and $25$ m). I have also tried interpolating the cloud surface to get the scale down to $l_0/2$, but the results were unsatisfactory as the interpolated surface is smoothed out.
For this reason, the measurements of $\mathcal{D}_f$ are performed over the scaling lengths between $l_0$ and $R$, where $R$ is the radius of the parcel. Because the shape of the observed clouds is not cylindrical, however, we need to estimate the radius of the individual clouds. There are two ways to calculate the radius. We can simply calculate the square root of horizontal area, or $r_g = \sqrt{\mathcal{A}}$, called the geometric radius, or we can calculate the average radial distance of all the individual points within the cloud. That is,
$$ r_d = \sqrt{\frac{1}{N} \sum_{i=1}^{N} r_i^2} $$
where $N$ is the number of points in the cloud and $r_i$ is the distance between a point $i$ and the center of the cloud.
The following figure shows the two estimates for cloud radius for a sample cloud:
Image(filename='../png/plot_observed_cloud.png')
Normally, the geometric radius $r_g$ is larger than the average radial distance $r_d$. It is because the former is the radius of a circle with the same area as the observed cloud, whereas the latter is an estimate of the effective radius of a stretched-out region.
Unless otherwise noted, we are using the geometric radius, as we find it works naturally better as a proxy to the cloud size.
Perhaps it might be more accurate to calculate the correlation integral:
$$ \mathcal{C}(r) = \frac{2}{N (N - 1)} \sum_{1 \leq i < j \leq N} H(r - | \pmb{x}_i - \pmb{x}_j |) $$
where $H$ is the Heaviside step function, $N$ is the number of points in the observed cloud, and $\pmb{x}$ is the position of a point. However, while the box counting method is sensitive to the placement of the grid and the scaling of the box size, the particular set of parameters chosen here appear to be accurate enough for our purpose. This seems to be because we are dealing with connected regions of the cloud, not point distributions.
Nonetheless, the following figures show the result of applying the box-counting method to 12-hour and 12.5-m simulation results, respectively.
Image(filename='../png/fdim_histogram_12hr.png')
Image(filename='../png/fdim_histogram_hres.png')
The blue curves denote the number of clouds observed in each bin, and the orange curves represent Gaussian distributions with the corresponding mean $\mu$ and the standard deviation $\sigma$.
The results correspond well to the previous observations. Individual clouds seem to have the fractal dimension $\mathcal{D}_f$ between $1.2$ and $1.4$, with the mean being roughly $1.33$. The exact value of $\mathcal{D}_f$ seems to vary among the cases we have tested, but it remains very close to $1.33$ regardless.
The 12-hour BOMEX simulation has slightly smaller fractal dimension, but interestingly, this seems to be because we have a longer simulation, not because the resolution is relatively lower. When we sample the first six hours of the simulation after a three-hour spin-up, we measure $\mathcal{D}_f=1.325$, which is much closer to the high-resolution case. This is an interesting result, and will be discussed later.
The simulated clouds, as previous observations suggest, appear to be fractals. If the clouds were cylindrical (i.e. if the horizontal shapes were circular), we would observe $\mathcal{D}_f=1$, simply because a circle is a one-dimensional curve.
We should also note that, naturally, this is very close to the fractal dimension of the surface of a Brownian motion. This seems to support an observation that turbulent flows are fractals (Sreenivasan and Meneveau, 1986), which also suggests that we are indeed observing turbulent properties of the individual clouds.
Assume that the aircraft and lidar could see only the larger size lobes and entrained blobs; and assume that the actual contorted shape is fractal-like with similar shaped contortions of ever-smaller spatial scales that exist in the real thermal, but which the lidar and aircraft cannot resolve.¶
- Write or adapt a simple computer program to draw/plot the contorted shape of the thermal that includes the superposition of all the scales of contortion, from the scale of the diameter R of the whole thermal down to some smaller turbulence scale L that you set as a parameter. But try to make your model be consistent with the measurements reported by Crum et al.
As previously mentioned, given the nature of large-eddy simulations, we assume a set of scale lengths $l \in \{l_0, l_0 + 1, \dots, R\}$. The model resolution $l_0$ are $25$ m and $12.5$ m for the two cases we present here. In reality, of course, the scale length $l$ can often become larger than the radius estimate $R$ because the measurements depend largely on the exactly shape of the eddies.
Nonetheless, while the instrumented aircraft and lidar measurements offer better resolutions, with the maximum resolution being 7.5 m for the lidar and smaller for the aircraft, they come with inherent instrumental and observational noise. In that aspect, we consider the $12.5$ m resolution more or less adequte for this study. And as we will see, the fractal properties of clouds observed at $25$ m resolution are consistent to those from the high-resolution run. To study the fractal nature of cumulus clouds, this appears to be a high enough resolution.
Then, what does the observed cloud look like at different scales? The following figure shows two sample clouds from the two model runs, with increasing scale length $l \in \{l_0, l_0 + 1, \dots, R\}$ to make a coarse-resolution observation.
Image(filename='../png/plot_coarse_resolutions_12hr.png')
Image(filename='../png/plot_coarse_resolutions_hres.png')
One unfortunate limitation in making coarse-resolution observations is that the set of scaling lengths $l$ has to be a set of integers $\{l_0, l_0 + 1, \dots, R\}$. There are a few different ways to resolve this. For example, It won't be an issue when the model resolution is higher (so that the set of scaling lengths $l$ is relatively dense) or if the observed clouds are much larger. We do currently have an output from the deep tropical marine convection during GATE on a 90-km domain, and two BOMEX simulations over a 32-km domain, but the output data need to be further processed and that will take a few weeks. It would be interesting to see how the analysis compares to deep convective clouds.
Both clouds shown above are taken at 4 hours into the simulations. It should be noted that the higher-resolution, 12.5-m simulation provides a more detailed representation of the smaller eddies present in the cloud, but that's most likely because in the 25-m simulation, the clouds are more likely to be broken up into smaller regions. This is due to the descrete nature of the LES model run.
Nevertheless, the box counting method is repeated over the entire set of scaling lengths $l \in \{l_0, l_0 + 1, \dots, R\}$ to estimate the fractal dimension (see Figure 4 below). We repeat the coarse-resolution observations above for all the clouds, which allows us to construct a statistically significant distribution of the fractal dimension $\mathcal{D}_f$.
Image(filename='../png/calc_fdim.png')
- Use your computer program to calculate the length P of the perimeter of the contorted shape as a function of L/R.
- Also find the relationship between P and fractal dimension, if any.
We could also measure the perimeter $P$ of the individual cloud slice $\psi(z, t)$ observed over a set of scaling lengths $l \in \{l_0, l_0 + 1, \dots, R\}$. This gives an estimate of a fractal dimension, which we will call the perimeter dimension $\mathcal{D}_p$. This is easy to do, as the perimeter of the fractal surface can be calculated by simply counting the sides of each grid cell between the cloudy region and the dry environment.
Because the grid size increases as the scaling length increases, we examine the relationship between $l$ and the perimeter for each cloud as $l$ increases from $l_0$ to $R$.
This analysis resembles calculating the fractal dimension using the box counting method, as we are interested in the slope between the perimeter $P$ and the length scale $l$ as $l$ changes. One might think of this as a perimeter dimension $\mathcal{D}_p$, where
$$ \mathcal{D}_p= \lim_{l \to 0} \frac{\log_{10} P(l)}{\log_{10} l} $$
which describes how long the perimeter of the cloud remains (i.e. how complex the boundary between the cloud and the environment is) as the resolution of the model becomes coarser. The above relationship transtates to
$$ P(l) \propto l^{- \mathcal{D}_p}. $$
The following figure calculates the slope $\mathcal{D}_p$ for the same cloud as above.
Image(filename='../png/calc_pdim.png')
It should be noted that the calculation of $\mathcal{D}_p$ is much more sensitive to the choice of $l \in \{l_0, l_0 + 1, \dots, R\}$ because once the resolution becomes coarse enough (which is normally around $l \sim R/2$), the perimeter of the observed cloud does not change as the scaling length becomes larger. This is most likely because once the box becomes big enough, we can no longer observe the smaller eddies present in the parcel. This is why in the above figure, the scaling length $\log_{10} l$ only goes up to roughly $1$, as opposed to $1.5$ for the calculation of fractal dimension. We will discuss this further in the next section.
It is not surprising that the perimeter of the observed cloud decreases as the scaling length increases; as $l$ increases and becomes closer to $R$, the smaller contortions are hidden, which results in a smoother boundary between the cloud and the environment.
We can repeat the calculation of $\mathcal{D}_p$ for all the individual clouds. Unfortunately, since we can no longer use the full set of scaling lengths $l \in \{l_0, l_0 + 1, \dots, R\}$, the number of cloud samples used is smaller this time around. We also observed that a portion of our samples are simply too small for this type of analysis.
Image(filename='../png/pdim_histogram_12hr.png')
Image(filename='../png/pdim_histogram_hres.png')
As previously noted, the variability in $\mathcal{D}_p$ is much higher than $\mathcal{D}_f$, especially since the values are much smaller. This seems to be because the calculation is sensitive to the size and the shape of the cloud sample.
It is probably not appropriate to name this a perimeter dimension, but interestingly enough, we have $\mathcal{D}_f \approx \mathcal{D}_p + 1 = 1.33$, which reminds us of the isotropy assumption $\mathcal{D}_s = \mathcal{D}_f + 1$, where the volume-surface dimension $\mathcal{D}_s$ is simply of a higher dimension (Siebesma and Jonker, 2000). Perhaps we are observing the same phenomena here on a lower dimension (that is, this is related to how the smallest patches/points are distributed in the cloud sample), although we could have a better understanding about the perimeter dimension $\mathcal{D}_p$ once we have a larger dataset including larger cloud samples, given the large variability in our $\mathcal{D}_p$ distribution.
We can further look at the perimeter estimates in terms of the fractal nature of cloud samples by examining the area-perimeter relationship. For fractal shapes, it is well-known that the area $A$ and its perimeter $P$ satisfies the following relationship:
$$ A \propto P^{2/\mathcal{D}_2}. $$
Since this is easier to observe, satellite measurement studies have measured this dimension (Lovejoy, 1982; von Savingny et al., 2011; Brinkhoff et al., 2015) to be very consistently close to $4/3 \sim 1.33$, which corresponds well to our estiamte of fractal dimension $\mathcal{D}_f$ shown in Figure 2. We will denote this value as $\mathcal{D}_2$, because it relates specifically and solely to the two-dimensional features of the cloud samples.
We can find $\mathcal{D}_2$ by re-arranging the above equation as:
$$ \begin{align} A &= k P^{2/\mathcal{D}_2} \\ \log_{10} A &= \frac{2}{\mathcal{D}_2} \log_{10} P + \log_{10} k \end{align} $$
where $\mathcal{D}_2$ is simply the slope of the $\log_{10} A$-$\log_{10} p$ distribution.
Interestingly enough, when we perform this calculation for the individually sampled clouds, we obtain $\mathcal{D}_2 \approx 1.2$. This is in fact consistent with Siebesma and Jonker (2000), who found this type of approximation rather inaccurate. Instaead, we repeated the calculation for the volume and the surface area of the cloud, and obtained results consistent with their work (not shown for brevity). This might suggest that estimating turbulent properties of individual cloud slices $\psi(z, t)$ might be inadequate and one might be forced to examining the individual three-dimensional cloud volumes, but it's beyond the scope of this study.
Instead, we replicate the satellite observations by taking a two-dimensional projection of the three-dimensional clouds in our domain. We expect this method to suffer from discretization effects, and the following plot shows the estimated $\log_{10} A$-$\log_{10} p$ over the cloud projection:
Image(filename='../png/plot_area_perimeter.png')
which gives $\mathcal{D}_2 \approx 1.368$.
The result is slightly higher than what we have observed as well as the satellite observations (von Savingny et al., 2011; Brinkhoff et al., 2015), but corresponds well to Lovejoy (1982) who reported $\mathcal{D}_f \approx 1.35$. Given the uncertainties involved in the calculations, and given a standard deviation of $\sigma \approx 0.1$ for the fractal dimensions, we could conclude that the perimeter is a good proxy to the fractal properties of the simulated clouds.
One possible explanation for this is that we are naively taking the projections of the individual cloud volume. Satellte-observed projections, on the other hand, will have multiple clouds overlapping with each other, whereas we solely assume that the individual clouds are independent.
- Define some intermediate scale S such that L < S < R. Further suppose that all scales of contortion smaller than S represent “shell” air (the cylindrical shell around the outside of the thermal that contains a mixture of thermal and environmental air), while all scales of contortion greater than S represent the contorted boundary between thermal air and environmental air. Describe and plot how the relative amount of entrained air varies with S?
- Based on your computer model, is there some radius r where r < R, such that essentially no entrained air is found inside a circle of radius r? If so, do you anticipate that this radius r would decrease with height inside the thermal, while R increases with height?
It's an interesting model of incorporating the idea of moist shells. As we have noted in Section 2.3, we found that there is a certain threshold in the scaling lengths $l \in \{ l_0, l_0+1, \dots, R \}$ where smallest contortions become smoothed out and the perimeter $P(l)$ becomes insensitive to the resolution.
There are a number of ways to define the moist shell surrounding the cloud core in the literature; for example, one could identify a moist downdraft, or simply take the region of moist air with condensed liquid water (cf. Hannah, 2017).
In order to find the relative amount of shell air, we first measure the actual size $\mathcal{A}$ of the cloud (at $l = l_0$) and take the difference between $\mathcal{A}$ and the size of the cloud $\mathcal{A}(l)$ measured at different scaling lengths $l \in \{ l_0, l_0+1, \dots, R \}$. Since the smaller eddies cannot be observed at coarse resolutions, it is not difficult to imagine that the size of the shell $\mathcal{A}_s = \mathcal{A} - \mathcal{A}(l)$ will increase as the scaling length increases (that is, as we take coarser and coarser observations).
I have chosen a few random cloud samples over a range of cloud radii, and plotted the relative size of the observed cloud and the corresponding shell against the relative scale measure $S(l) = l/R$.
Image(filename='../png/find_area_12hr.png')
Image(filename='../png/find_shell_area_12hr.png')
The above figures require a bit of explanation, as it is not exactly easy to visualize how all the observed clouds behave at different resolutions.
It should first be noted that the cloud size $\mathcal{A}(l)$ is observed with a virtual resolution of $l$; that is, we assume that the contortions smaller than this resolution, with a bit of error margin ($5\%$ of grid size), constitute the shell region.
As a result, we see that the relative size of the shell approaches $1$, mainly because the observed cloud size $\mathcal{A}(l)$ becomes smaller (Figure 8.1.) and the size of the shell $\mathcal{A}_s = \mathcal{A} - \mathcal{A}(l)$ gets close to the original size of the cloud $\mathcal{A}$, including the smallest eddies in the thermal (Figure 8.2.).
Perhaps it would be more intuitive to name it total size of eddies smaller than the observational resolution $l$, or simply unresolved eddies (more on this later).
Nonetheless, as expected, the size of the shell increases with increasing $S$, or as the scaling length $l$ approaches the cloud radius, representing the size of the largest eddy in the cloud sample. The relative amount of the entrained air in the model, therefore, increases with the relative scale measure $S(l) = l/R$.
From Figure 8.1. and 8.2., the scaling behaviour can be observed to be oddly consistent over a range of cloud sizes. I have taken two linear regression (using Bayesian regression method) estimates for the shell, one for the relative scale measure less than $0.4$ and one greater. Although this particular parameter is not meant to be an accurate estimate, it is evident that $\mathcal{A}_s/\mathcal{A}$ reaches a plateau as $S$ increases.
We argue that this plateau represents a scale length $l$ where smaller eddies become smoothed out as our observation becomes coarse. Now that there is an inflection point, that might be a great way to pinpoint the scale length $l$ where the smaller eddies (moist shells containing the entrained air) can no longer be identified for the cloud as well as the shell.
How do we identify this inflection point? From Figure 8.1., it appears that as our observation becomes coarser ($l \to R$), the observed size of the cloud goes from $1$ to zero (naturally because we are observing at the model resolution $l = l_0$). The transition occurs somewhat smoothly, so it is not entirely clear if there is a fixed inflection point.
We can think of this model as one where at an observational resolution $l$, we consider any eddies smaller than this resolution as the shell, a mixing region between the thermal and the environment, and eddies larger than $l$ as the cloud.
Let's assume that there is indeed a critical length scale $l_c$ which marks the inflection point where the cloud can no longer be seen (i.e. the observation is too coarse to observe the thermal itself), or where $\mathcal{A}(l) / \mathcal{A} \approx 0$. We can simply conclude that this is roughly at $S(l) \approx 0.4$, but that sounds somewhat unsatisfactory, especially given the amount of data from all the individual cloud samples.
Instead, we investigate if this critical length scale $l_c$ depends on the size of the parcel $R$. To this end, we calculate the changes in the observed cloud size and the relative shell area $\mathcal{A}_s/\mathcal{A}_0$ as a function of $S(l) = l/R$ for the individual clouds, which are shown below:
Image(filename='../png/plot_area_scatter_S.png')
Image(filename='../png/plot_shell_scatter_S.png')
Essentially, Figure 9.1 and 9.2 show all the data points gathered by repeating the same analysis in Figure 8.1 and 8.2 for all the indiviual clouds present in our dataset. Each bin in the two-dimensional histogram in Figure 9.1 and 9.2 contains points that follow the curves observed in Figure 8.1 and 8.2, respectively.
The colors on the above two figures suggest the original size of the cloud $\mathcal{A}$, which is added so that we can track how the observed cloud size (Figure 9.1) and the relative size of the shell (Figure 9.2) change as a function of $S(l)$.
Two things can be noted about the two figures. First, unfortunately, the transition in the observed size of the cloud (Figure 9.1) and of the shell (Figure 9.2) appears to be smooth as $S(l) = l/R$ increases, without an obvious point of inflection. Second, perhaps more importantly, there is a size dependence on the relative measure $S(l) = l/R$ for both the clouds and their shells.
One may also notice that both the size of the cloud $\mathcal{A}$ and of the shell $\mathcal{A}_s$ becomes invisible at $S(l) \approx 0.3$. This is mainly because the two figures only show the mean relative size of the cloud (Figure 9.1) and of the shell (Figure 9.2) within each bin, but does not show the actual frequency distribution. From the figure below, which is the same plot as Figure 9.1 but drawn as a two-dimensional histogram, it is easier to see that as we have observed in Figure 8.1, most clouds become invisible at $S(l) = 0.4 - 0.6$:
Image(filename='../png/plot_area_hist2d.png')
Then, combined with Figure 9.1, we notice that the (relative) critical length scale $S=l_c/R$ decreases with increasing cloud size. This is not surprising, because without any prior knowledge, we are assuming that there will be a fixed length scale $l_c$, inherent to moist convection (and espeically turbulent mixing).
Of course, it is still possible that this critical length scale $l_c$ depends on the size of the cloud. To this end, we plot the relative size of the cloud and of the shell as a function of $l$ instead:
Image(filename='../png/plot_area_scatter_12hr.png')
Image(filename='../png/plot_shell_scatter_12hr.png')
So this is indeed surprising: the larger the cloud, the larger the critical length scale $l_c$ where the primary eddies consisting the thermal become invisible, at least statistically speaking. We can see this in Figure 10.1 by choosing an arbitrary scale at $\mathcal{A}(l)/\mathcal{A} = 0.2$. Given this threshold, it is easy to observe that the larger the cloud size $R$, the larger the size of resolved eddies $l_c$.
We could also observe that, based on the skewness in the distribution of eddies for the largest clouds (purple), the larger the cloud, the larger the primary eddy consisting the parcel. For example, for the largest clouds in Figure 10.1 it is easy to see that a larger portion of the area (measured in $\mathcal{A}(l)/\mathcal{A}$) can be traced back to the largest eddies (with large $l$).
Moreover, this suggests that $S(l)$ depends on the size of the thermal, and it is interesting because it seems to suggest that larger eddies are required to allow the growth and development of larger thermals.
Therefore, in our model where the size of the cloud increases with height, the minimum eddy size of the thermal, which determines the minimum size of the entrained air being observed, will increase accordingly. In our model, this means that the moist shell surrounding the cloud core will increase with height, as the cloud size increases.
Interestinly enough, this is exactly the observation we (as well as Hannah, 2017) made from a series of high-resolution large-eddy simulations, except we do not observe the cloud area increasing with height in these model outputs. However, this is still noteworthy and we will revisit the implications of this model later.
Still, it is rather difficult to conceptualize a model based on the figures we have presented so far, mainly because of the relationship between $l$ and $S(l) = l/R$ in measuring the observed cloud size. Perhaps it is better to step back and make use of the large number of cloud observations (hopefully).
First, we want to make sure that the conclusions from Section 2.4 really holds; that is, the size of the primary eddies increase with increasing cloud size.
In order to make a more direct observation, we measure that relative cloud size at $0.3 R$, where $R$ is the cloud radius. The coefficient $0.3$ isn't particularly significant, as adjusting this parameter does not affect our observations as long as it is not too small (as to be able to resolve the smallest eddies) or too large (where nothing can be observed). Regardless, we take the individual clouds in the cloud field, and at every height where the cloud is defined, measure the relative cloud size $\mathcal{A}(l = 0.3 R)/\mathcal{A}$. The following figure shows the distribution of all the observed clouds:
Image('../png/plot_radius_vs_area_cloud.png')
Fortunately, this is as expected. The observations made at different resolutions (both in $S(l) = l/R$ and $l$) can be reproduced with actual measurements taken from simulated clouds with varying sizes. We could also observe that as the cloud size grows, the relative area covered by an eddy who size increases with the cloud $A(l = 0.3 R)$ decreases.
This is not much different from earlier observations (e.g. from Figure 10.1.), and it further shows that the size of the primary eddy is not exactly a linear function of cloud size. That is, while the primary eddy size increases with cloud radius, the rate at which $l_c$ increases is much slower than the rate of increase in $R$.
Repeating the same measurement as above with a fixed eddy size (e.g. with $l = 125$ [m]) yields the same result, where the relative size occupied by that eddy becomes smaller with increasing cloud size.
If we cannot correctly measure the characteristic eddy size solely using the geometric cloud radius $R_g = \sqrt{\mathcal{A}}$, what else can we use?
This is when I realized that Figure 1 looks more or less like the diagram of the intromission zone (cf. Figure 14 in Crum et al, 1987b), where the average radial distance corresponds to the shape of the cloud core. Because the average radial distance $R_d$ gives an estimate of the largest coherent structure with the shape of a circle (or a cylinder if we consider the entire cloud volume), we could use it to approximate the cloud core size.
There are two ways to use $R_d$ to approximate the cloud core in the individually observed cloud. We could define the cloud core area $\mathcal{A}_c = R_d^2$, since we are assuming a circular cloud core, or we could take the box counting at this resolution $l = R_d$ to measure the observed size of the core eddy. I prefer to do the latter, for reasons that will be described later.
Figure 11.2 below replicates the above figure, but now with the average radial distance $R_d$ as the target eddy size:
Image('../png/plot_radius_vs_area_core_R.png')
Indeed, we could conclude that the size of the largest (primary) eddies of the thermal increases with cloud size.
Now, what about thet shell? Since we defined the shell as the unresolved eddies, the shell is simply the non-core portion of the cloud. Since we could have already shown that the size of the core increases slowly compared to the size of the cloud, it is not difficult to imagine that the shell will become larger as the cloud grows. That is, the larger the cloud, the larger the shell:
Image('../png/plot_radius_vs_area_core_area.png')
where the black line denotes cases where $\mathcal{A}_c = \mathcal{A}$. As the cloud size $\mathcal{A}$ increases, the distribution shifts away from this line, because the cloud core area increases relatively slowly. Naturally, we can define the shell as the difference between the size of the observed cloud size and that of the cloud core, or how far the distribution shifts away from the black curve.
Figure 11.3 describes our model more clearly; as the cloud size increases, the relative size of the core $\mathcal{A}_c/\mathcal{A}$ decreases, while the shell $\mathcal{A}_s$ grows. We imagine that the shell is where the entrainment is occuring, and so our model suggests that (total) entrainment also increases as the cloud grows.
We have a few hypotheses explaining this, and the most relevant one is that the larger clouds consist of multiple primary eddies. It is possible that the clouds can grow larger than the largest eddies in the atmosphere as multiple eddies can merge and form larger clouds, which is why the size of the primary eddy does not depend directly on cloud size (since $R_g = \sqrt{\mathcal{A}}$).
Of course, the fractal nature of clouds is also playing a role as it dictates that when the size of the fractal increases, its perimeter increases according to the power-law relationship ($A \propto P^{2/D_2} \approx P^{1.5}$). We have seen this previously from the area-perimeter relationship in Section 2.3, where the surface of the growing cloud becomes more and more complex, increasing the size of the shell, or the amount of unresolved eddies.
Nonetheless, it must be noted that this is still an approximation; we make the assumption that the cloud cores consist of the largest eddies, whose sizes depend on the average radial distance $R_d$. This is done in order to show the relationship between the primary eddy size and the size of the actual cloud, and fortunately our observations from Section 2.4 still hold.
Therefore we could conclude that for our conceptual model, the size of the largest eddies (which define the cloud core region, without any unresolved eddies) increases with the size of the cloud. As we assume that the size of the thermal increases with height, we could also say that this scale also increases with height.
Perheps we could further simplify our model. As we have used Figure 1 to make an assumption that the size of the cloud core is a function of the average radial distance $R_d$, we could also assume that the outer boundary of the intromission zone can be estimated by the geometric radius $r_g = \sqrt{\mathcal{A}}$. Of course, this won't be an accurate approximation since the observed clouds are largely skewed, elongated, and in many cases, both. But the goal is to have some measure of the size of the shell.
To this end, I have taken measurements of $r_d$ and $r_g$ for all the clouds at every height where they are defined (see Section 2.1 for the calculation of these values) and made a simple plot:
Image('../png/plot_radius_vs_area_radius.png')
where again, the black curve shows where $r_d = r_g$. This also marks an interesting case where the thermal is indeed cylindrical in shape, although it appears that such cases are rare.
It is not surprising that normally, $r_g > r_d$, which would be exactly the case shown in Figure 1 where $r_d$ defines the inner boundary of the intromission zone (cloud core) and $r_g$ the outer boundary (shell). We could see that on this bivariate kernel density estimate (KDE), the density of the Gaussian kernel is higher towards the right-hand side of the black curve.
More importantly, as the size of the cloud increases (which can explicitly be inferred from $r_g = \sqrt{\mathcal{A}}$) the rate of change in $r_d$ seems to be close to the black curve. Therefore, for this simplified model of thermals, we could assume that the size of the shell scales directly with the size of the cloud. Our new model, then, describes how the clouds form and develop. We will discuss this in more detail in the following section.
There are still some cases where $r_g < r_d$ (LHS of the black curve), and it would definitely be interesting to look at what is happening, especially for large clouds. One possibility is that when the horizontal cross-section of the cloud is skewed and elongated, the average circular distance covered by the cloud can be larger than the radius of the equivalent circle of same size. It makes, for example, little sense to calculate $r_d$ for clouds seen in Figure 3.1 and 3.2, and depending on the shape of the observed cloud, simply calculating the core area by $\mathcal{A}_c = r_d^2$ won't yield an acceptable representation of the cloud core.
If we used artificial thermals with a fractal dimension of $1.33$, it is not too difficult to imagine that those thermals will better represent the intromission zones. In such thermals, $r_d$ will exactly correspond to the non-turbulent part of the thermal (i.e. the thermal core) and $r_g$ will yield an estimate of how complex the surface of the thermal is (i.e. the shell representing the turbulent mixing region). Still, the behaviour of $r_d$ and $r_g$ won't change drastically. Thus, we can conclude that the simulated clouds faithfully reproduce the characteristics of fractal clouds.
Regarding Stull & Eloranta 1985 NWD, consider their findings that lidar-observed cloud base is almost exactly equal to the theoretical LCL height. But this theoretical LCL height assumes undiluted air in the thermal. Namely, there is no entrainment into a protected core of the thermal.¶
Discuss how your model gives results that are consistent with the findings of Stull and Eloranta (1985).
Representing the moist shell surrounding the cloud core with unresolved eddies is an interesting concept, that reproduces many characteristics of convective clouds observed in LES modelling studies as we noted above.
The model corresponds well to the observations by Stull and Eloranta (1985) and produces results that are consistent with observational studies.
As we have seen in Section 2.5, both the primary eddies and the cloud core grow as they ascend. The size of the shell (the intromission zone) is proportional to the size of the cloud (Figure 11.3), which implies that the thermal experiences little to no entrainment near LCL. As it ascends above the LCL and grows in size, the turbulent mixing processes occur, forming an intromission zone consisting of unresolved eddies. A simplified diagram of our thermal would look like:
%matplotlib notebook
fig = plt.figure(1, figsize=(6, 6))
plt.plot((-1, 1), (0, 0), '-', c=xkcd['dull orange'], lw=4, label='LCL', zorder=4)
xi = np.linspace(-1, 1, 400)
# Cloud core
yi_c = np.zeros(400)
yi_c[(xi < -0.1)] = -5 * xi[(xi < -0.1)] - 0.5
yi_c[(xi > -0.1) & (xi < 0.1)] = 0
yi_c[(xi > 0.1)] = 5 * xi[(xi > 0.1)] - 0.5
plt.plot(xi, yi_c, '-', lw=2, c=xkcd['bluish grey'], label='Core')
# Shell
yi_s = np.zeros(400)
yi_s[(xi < -0.1)] = -2 * xi[(xi < -0.1)] - 0.2
yi_s[(xi > -0.1) & (xi < 0.1)] = 0
yi_s[(xi > 0.1)] = 2 * xi[(xi > 0.1)] - 0.2
plt.plot(xi, yi_s, '--', c=xkcd['carolina blue'], label='Shell')
plt.fill_betweenx(yi_c, np.linspace(-1, 1, 400), color=xkcd['bluish grey'], alpha=0.5)
plt.fill_between(xi, yi_c, yi_s, color=xkcd['carolina blue'], alpha=0.5)
# Mixing arrows
plt.arrow(0.4, 0.8, 0.2, 0, head_width=0.04, width=0.01, fc=xkcd['dusty rose'], ec=xkcd['dusty rose'])
plt.arrow(0.6, 0.75, -0.2, 0, head_width=0.04, width=0.01, fc=xkcd['azure'], ec=xkcd['azure'])
plt.xlim([-1, 1])
plt.ylim([-0.2, 1])
plt.xticks([])
plt.yticks([])
plt.legend(loc=5)
where we assumed a linear increase in size with height.
The above diagram is an idealized representation of our cloud, but it is still representative of the clouds according to our model where the cloud size is a function of height, and the size of the shell grows as the cloud size increases.
The model describes that the cloud will form at the cloud base (near LCL) and will have little to no shell surrounding the core as it emerges. From Figure 12.1, we made the assumption that the shell scales with the cloud size, although from Figure 11.3 it is obvious that the shell grows faster than the core. Either way, the shell consisting of unresolved eddies increase with height, which implies that the cloud experiences more entrainment as it rises.
This behaviour is consistent with the observation of cloud base height near LCL. When the model thermal ascends, its spatial scale is too small for the lidar to identify or measure any mixing activities. It reaches the LCL, and the cloud forms. From that point, the growing cloud starts to mix with the environment, which can be represented by a growing shell of smaller eddies that contain a mixture of thermal and environmental air. This scale at which no entrained air is found will be slightly different for each cloud, as we do observe a range of (small) clouds with no shell (Figure 11.3). For this simplified model, however, it suffices to say that the cloud emerges without experiencing any significant dilution due to entrainment.
Lastly, we could provide an explanation for the lack of protection for the cloud core at the cloud base, which is that the thermal rising from the surface goes through the same process. As it rises, the thermal sheds the shell and only the protected core of the thermal will be lifted above the LCL and form the cloud. Or, we can think of the shell representing the turbulent mixing region specifically between the cloud and the dry environment. Because mixing near LCL, where the atmosphere is moist enough, does not dilute the thermal. But it is probably an unnecessary addition to our simple model, which already explains the lack of entrainment at the cloud base.
Consider the observations of entrainment into the top of the mixed layer (as given by many of the other readings that I assigned), and how probability distributions (e.g., LCL zone and Entrainment zone EZ) can be found as a function of height.¶
- Use that same concept of probability distributions to describe air entrained laterally into a thermal or cumulus cloud, and thermal air detrained into the environment.
- Discuss how this conceptual model might or might not be different for lateral entrainment into cumulus clouds.
We have defined the turbulent mixing processes between the cloud and the dry environment for our simple cloud model above, but in order to describe entrainment and detrainment processes in more details, it is obligatory to make a few more assumptions.
The precise distribution of total water $q_t$ in cloud cores and in the dry environment is already known for the BOMEX case,
Image('../png/BOMEX_QT_profile.png')
and while it is trivial to assume that the cloud mixes directly with the dry environment following the conventional assumptions, that is not the purpose of our simple cloud model.
Instead, we will describe the properties of the entrained and detrained air in terms of mixtures of mixed layer (ML) air that is being brought up from the cloud base, and the free atmosphere (FA) air that is being entrained from the top of the cloud. This resembles the one-dimensional entraining-detraining plume model used in traditional models, and we believe it is more in line with the goal of this study.
To this end, we need to make a few assumptions regarding the dynamics of turbulent mixing. That is, we will assume that
The intromission into the cloud happens as the cloud ascends, and it depends on the size of the intromission zone.
The last assumption is perhaps the cornerstone of our simple model when it comes to turbulent mixing processes. Then, how does the intromission depend on the size of the the intromission zone? We have previously defined the intromission zone as the shell of moist air surrounding the cloud core, defined by the unresolved eddies. We know that the size of the shell scales with the size of the cloud, which increases with height.
And from what we've observed so far in this study, we can assume that the probability of lateral entrainment depends on the size of the cloud. This goes back to the fractal nature of the clouds, where the small, contorted shapes form around the cloud. It is not too hard to imagine that when the cloud is larger (i.e. the larger the interface between the cloud and the environment), lateral entrainment will be more frequent.
Then, based on the radius of our cloud and the fact that our model clouds have cylindrical shapes, we make the assumption that the likelyhood of entrainment scales with the size of the cloud, or
$$ \begin{align} \mathcal{P_{\epsilon}}(z) &\propto \int^{z_2}_{z_1} A(z) dz\\ &\propto \int^{z_2}_{z_1} R(z)^2 dz \end{align} $$ where the integral is defined betwee the cloud base $z_1$ (or, in this case, LCL) and the cloud top $z_2$. Given the above definition, we plotted our probability density function (PDF) for lateral entrainment below.
%matplotlib notebook
fig = plt.figure(1, figsize=(6, 6))
xi = np.linspace(0, 1, 400)
# Subplot 1: Cloud area as a function of height
plt.subplot(1, 2, 1)
plt.plot((0, 1), (1, 1), 'k-', zorder=4)
plt.plot((0, 1), (0, 0), 'k-', zorder=4)
plt.fill_between(xi, xi, 1, color=xkcd['sky blue'], alpha=0.6)
plt.xlim([0, 1])
plt.xticks([0, 1], ['0', r'R'])
plt.yticks([0, 0.5, 1], [0.6, 1, 1.4])
plt.xlabel('Cloud Radius')
plt.ylabel('Height [km]')
# Subplot 2: PDF of entrainment
plt.subplot(1, 2, 2)
plt.plot((0, 1), (1, 1), 'k-', zorder=4)
plt.plot((0, 1), (0, 0), 'k-', zorder=4)
plt.fill_between(xi**(2), xi, 1, color=xkcd['greyblue'], alpha=0.6)
plt.xlim([0, 1])
plt.xticks([0, 0.5, 1])
plt.yticks([])
plt.xlabel('Entrainment PDF')
plt.show()
We can then use the lateral entrainment PDF to describe the properties of entrained and detrained air. First, it is not difficult to see that at the top of the cloud (the inversion layer) the air being entrained into our cloud will have the properties of the SL air. Likewise, at the cloud base, the air being entrained will have the properties of, or very close to, the ML air.
While the cloud is ascending the cloud layer, then, we can describe the entrained air being a mixture of the cloud air in the intromission zone and the dry environmental air, whose proportion is defined by the entrainment PDF (which describes the likelyhood of finding the dry environmental air in the intromission zone).
For the detrained air, the situation is slightly more complicated; we know that the cloud will detrain ML air near LCL and mostly FA air near the top, but since we are assuming that the entrainment occurs instantly (although it won't make any difference if it does not), the cloud will be detraining a mixture of entrained air from lateral entrainment and the interior cloud air in between.
How do we estimate the properties of the cloud? We could use the entrained air as a proxy to estimate the dilution of the cloud. That is, we assume that the entrainment PDF $\mathcal{P}_\epsilon$ is the cumulative distribution function (CDF) of the detrained air, and that $\mathcal{P}_\epsilon$ now represents the odds of finding the entrained environmental air in our cloud (0 up to $\mathcal{P}_\epsilon$(z)). Conceptually, then, the above PDF represents the total detrainment of our cloud, which dissipates at the free atmosphere. The detrainment PDF, then, is simply a derivative of the entrainment PDF (since CDF is an integral of PDF).
We have plotted the distribution of entrained and detrained air in terms of the properties of ML and FA air below.
%matplotlib notebook
fig = plt.figure(1, figsize=(6, 6))
xi = np.linspace(0, 1, 400)
# Subplot 1: Entrained Air
plt.subplot(1, 1, 1)
plt.plot((0, 1), (1, 1), 'k-', zorder=4)
plt.plot((0, 1), (0, 0), 'k-', zorder=4)
plt.plot(1-xi**(2), xi, color=xkcd['deep rose'], label='Entrained Air PDF')
plt.plot(1-xi**(3), xi, color=xkcd['greyblue'], label='Detrained Air PDF')
plt.xlim([0, 1])
plt.xticks([0, 1], ['FA', 'ML'])
plt.yticks([0, 0.5, 1], [0.6, 1, 1.4])
plt.ylabel('Height [km]')
plt.legend(fontsize=12)
plt.show()
We can make a few interesting observations here. First, there is a bit of inverse relationship between the net entrainment (i.e. net dilution) of the cloud and its size. Our cloud would be the largest near the top, but the net dilution is the smallest there. The dilution is large near $1$ km, where the cloud will still be fairly small. This corresponds well to the observations, where larger clouds dilute relatively slowly. This relationship does not hold near the cloud base, however.
Second, the cloud experiences a net entrainment as it grows. This is not too surprising conceptually, as we expect entrainment to be strictly diluting and detrainment strictly concentrating our cloud.
Lastly, the above plot agrees with Figure 9 in Crum et al. (1987b), as we assume that the air being detrained has the properties of the cloud itself, and we can think of the detrainment PDF as a fraction of ML air in the cloud. The rate of change in our model is a lot slower than the observation made in Crum et al. (1987b), where the cloud loses 80\% of the surface air within the top 20% of the cloud layer. Still, our model, albeit quite simple in nature, reproduces the general tendency observed in instrumented aircraft observations.
Conceptually, this model treats the lateral entrainment as analogous to vertical entrainment. We have intentionally studied the properties of entrained and detrained air as a function of ML and FA air properties. In both cases, our entrainment PDF will look the same regardless of the assumptions.
For vertical entrainment, the only difference would be the eddy size distribution and the direction of entrainment. That is, the clouds are no longer forced to increase in size with height, as most of the entrainment occurs near the top. For this reason, the entrainment PDF will still look the same for vertical entrainment as long as we assume that entrainment occurs less frequently as one moves towards the more turbulent regions (cloud base for vertical entrainment, cloud core for lateral entrainment).
We have looked at the fractal nature of modelled clouds from a high-resolution, large-eddy simulation. The simulated clouds are shown to consistently display fractal nature, both using the box-counting method and observing their area-perimeter relationships.
There are definitely more work to be done for the latter, however. It is already well known that the entrainment rate is a strong function of perimeter (Dawe and Austin, 2013), and taking the horizontal cross-section seems to distort that relationship, as observed by Siebesma and Jonker (2000). Perhaps to study the volume-surface relationship, we cannot rely on cross-sectional cloud properties.
One interesting observation we made, which is not included in the above analysis, is that the fractal dimension $\mathcal{D}_f$ appears to fluctuate over time. There wasn't enough time to study the time series of the fractal dimension, but it was obvious that towards the end of both simulations, the fractal dimension tends to be slightly smaller. It could be because the clouds indeed become less fractal as they grow and merge, or it could simply be because we have more broken-up cloud patches distorting the distribution of fractal dimensions.
As per our conceptual model, there are some interesting properties that relate to the direct entrainment calculations.
For example, considering the shell as the unresolved eddies appears to yield consistent behaviours for the model clouds. This is interesting, because I remember a study where, without directly sampling the cloud core, the researchers would take the spherical region from the cloud to define and track the dynamics of the core, which seems to correspond to the methods used here. Perhaps the dynamics of the cloud core can be reduced to its geometric (as well as thermodynamic) properties.
Also, the shell is normally sampled as the region of condensed liquid water surrounding the (accelerating) core, and the idea that the moist shell plays a significant role in governing the entrainment and detrainment between the cloud and the environment is not new. Although the observed behaviour of the shell is slightly different (e.g. larger clouds would have relatively thinner shells), the idea that the shell contains a mixture of the cloud core and the environment air relates directly to what direct entrainment calculations reveal.
The simple model does not, of course, represent the subtle dynamics of turbulent mixing. For example, we have observed that the properties of the entrained air does not simply correspond to those of the moist shell surrounding the cloud core, and that the cloud preferentially entrain moist air from the shell, which leads to reduced rate of dilution.
Lastly, the model does not differentiate between the entrainment into the cloud and the dilution of the cloud. This is an important topic in studying the turbulent mixing processes using the direct calculation method for entrainment and detrainment rates, but they are used ambiguously in this model. In fact, if we combine the entrainmend PDF above with the fact that both the cloud and environment become drier with height, we will find that the properties of detrained air (which represents the fraction of ML air in the cloud) correspond even better with Figure 9 in Crum et al. (1987b) because entrainment will be concentrated near the top of the cloud layer. Likewise, it would be interesting to see if we can improve our one-dimensional cloud model, by incorporating the idea of dilution due to entrainment, and comapre it to the traditional convective models.