# Monte Carlo Dropout Ensembles for Robust Illumination Estimation

Firas Laakom1, Jenni Raitoharju2, Alexandros Iosifidis3, Jarno Nikkanen4, and
Moncef Gabbouj1
1Department of Computing Sciences and Communication Sciences, Tampere
University, Finland
2Programme for Environmental Information, Finnish Environment Institute,
Finland
3Department of Engineering, Aarhus University, Denmark
4XIAOMI, Hermiankatu 6-8H, Tampere, Finland
Emails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Computational color constancy is a preprocessing step used in many camera
systems. The main aim is to discount the effect of the illumination on the
colors in the scene and restore the original colors of the objects. Recently,
several deep learning-based approaches have been proposed to solve this
problem and they often led to state-of-the-art performance in terms of average
errors. However, for extreme samples, these methods fail and lead to high
errors. In this paper, we address this limitation by proposing to aggregate
different deep learning methods according to their output uncertainty. We
estimate the relative uncertainty of each approach using Monte Carlo dropout
and the final illumination estimate is obtained as the sum of the different
model estimates weighted by the log-inverse of their corresponding
uncertainties. The proposed framework leads to state-of-the-art performance on
INTEL-TAU dataset.

##  I Introduction

The human eye can effectively adjust to changes in visual conditions and
different illuminations of the scene [1, 2]. The well-acknowledged image
processing problem of computational color constancy tries to mimic this
ability in digital cameras [3]. The aim is to restore the ‘original’ colors of
the objects in a given scene and make it look like it was taken under a
neutral white illumination. The main assumption which is usually made to
reduce the complexity of the problem is that of one global uniform
illumination across the whole scene [3, 4]. Computational color constancy
methods then can be decomposed into two steps: the first step is to estimate
the global illumination of a given scene and the second step is simply to
normalize the pixel values of the scene using the estimated illumination. As
the second step is a straightforward transformation, the problem of color
constancy can be reduced to estimating the illumination and, thus, is usually
referred to as the illumination estimation problem [3]. Over the past few
years, different unsupervised [5, 6, 7, 8, 9, 10, 11, 12, 13] and supervised
[14, 15, 16, 17, 18, 19, 20] techniques have been proposed to tackle this
problem.

Recently, with the advancement of data-driven methods and the availability of
more labeled illumination estimation datasets [21, 22, 23], the direction of
computational color constancy research has shifted towards learning-based
approaches in general and deep learning-based methods in particular [16, 17,
18, 19, 20, 24, 25]. Convolutional Neural Networks have been recently
extensively used to estimate the global illumination of a given scene and have
often outperformed traditional statistical methods [14, 17, 19]. While these
methods usually yield low average error rates, their usage in industrial
applications remains limited [26, 27]. This is due to the inability of a
single model to accurately approximate the illuminations of all types of
scenes and the failure of each of these approaches for some scene content [18,
26, 28].

Figure 1: Combining three different neural networks model outputs
$\textbf{I}_{1}$, $\textbf{I}_{2}$, and $I_{3}$ using their predicted
uncertainties $\sigma_{1}$, $\sigma_{2}$, and $\sigma_{3}$, respectively.

In this paper, we address this limitation by proposing a novel deep learning
approach called Monte Carlo Dropout Ensemble Monte Carlo Dropout Ensembles
(MCDE). In the proposed framework, we aggregate the estimates of different
deep learning methods depending on the scene content and on the confidence of
their predictions, as illustrated in Figure 1. To do so, we need to estimate
the uncertainty of each model alongside its illumination estimation. To this
end, we employ Monte Carlo dropout (MC-dropout) [29, 30] to obtain the output
uncertainty of each CNN model. Given both the illumination estimation and the
model’s uncertainty for each of the models, we propose a confidence score-
based scheme to produce the final global illumination estimation.

The use of MC-dropout for estimating the model’s uncertainty does not require
retraining of the available models. Thus, it can be used to model the
uncertainty of any CNN-based model trained with dropout [30]. With the model
confidence at hand, we can treat uncertain inputs and special cases
explicitly. In this paper, we study as a special case of our approach the
combination of two state-of-the-art CNN-based color constancy models, namely
the Bag of Color Features (BoCF) [19] and the Fully Convolutional Color
Constancy (FC4) [17], by weighting their relative illumination estimations
according to their relative uncertainties. This yields a robust illumination
estimation framework.

The main contributions of this paper can be summarized as follows:

  * •

We propose a novel computational color constancy approach, called Monte Carlo
Dropout Ensembles (MCDE), that aggregates different CNN-based models based on
their confidence in illumination estimation.

  * •

We propose using Monte Carlo dropout to approximate the uncertainty of
different color constancy models.

  * •

To obtain the final estimate, we rely on the log-inverse of the uncertainty as
a confidence score.

  * •

We study as a special case of our approach the combination of Bag of Color
Features (BoCF) and Fully Convolutional Color Constancy (FC4). The proposed
approach yields state-of-the-art results on INTEL-TAU dataset.

The remainder of the paper is structured as follows. Section II provides a
brief review of CNN-based color constancy approaches and Monte Carlo dropout.
Section III presents the proposed approach MCDE. Section IV presents the
experimental setup and reports the evaluation results on INTEL-TAU dataset.
Section V concludes the paper.

##  II Related work

###  II-A CNN-based color constancy

CNNs have been extensively used in the last decade to solve the computational
color constancy problem [16, 17, 18, 19, 20]. The main task is to train a
neural network that approximates the global illumination of the scene. CNN are
usually optimized in an end-to-end manner using back-propagation. We divide
CNN-based approaches used to tackle the illumination estimation problem into
two main categories: methods operating on local patches as input, called
patch-wise methods [16, 18, 31, 32, 33], and methods operating directly on the
full image as input, called single-pass methods [15, 17, 19].

In the first group of methods, the models are trained on small patches of the
scene to overcome the data shortage. Patch-based CNNs were first used for
color constancy by Bianco et al. in [16, 31], where local estimates are
generated from the small patches of the scene. The final estimate is pooled as
the mean or median of these local estimates [16] or passed to a support vector
regressor [31] to estimate the illumination color. In [18], an unsupervised
pretraining phase is proposed to enhance the generalization of the model. A
fully convolutional autoencoder with two branches, one for scene
reconstruction and one for illumination estimation, are initially trained on
the color constancy dataset augmented with ImageNet. Then, the illumination
estimation branch is fine-tuned to approximate the patch-based illumination of
the scene. This resulted in a better generalization capability and a more
robust model. Another patch-based CNN approach is proposed by Bianco and
Schettini [32, 33], which focuses specifically on face regions for color
constancy.

In the second group of approaches, i.e., single-pass methods, networks are
trained on the full image to approximate the illumination of the scene. These
methods are faster as they require a single forward pass in the test phase.
However, they usually operate on high resolution images and thus they have a
higher number of parameters and higher memory requirements. The early works in
this category were proposed in [15], where a CNN operating on the full image
is trained to extract feature hierarchies to achieve robust color constancy.
However, these methods rely on the conventional CNN architecture, i.e.,
convolutional layers followed by max-pooling and a fully connected block in
the end. Thus, in the inference phase, images need to be resized to predefined
dimensions, which may introduce spatial distortions of image content and
degrade the performance of the method. To overcome this limitation, recently
FC4 [17] and BoCF [19] were proposed.

FC4 [17] method relies on a fully convolutional topology using a confidence-
weighted pooling layer. The network is trained end-to-end to incorporate
learning the confidence of each patch of the image in a single pass and
optimizing the weighted sum of the local estimates based on the confidences.
In [19], it is argued that the spatial information is not important in the
color constancy context and, thus, it is discarded using a BoF layer [34]. The
resulting BoCF network is composed of three blocks, namely the feature
extraction, the Bag of Features, and the illumination estimation blocks. In
the first block, a nonlinear transformation of the scene is produced. In the
second block, a histogram representation of this transformation is compiled.
This histogram is used in the third block to approximate the illumination.
Moreover, an attention mechanism is used to adaptively learn to select the
elements of the histogram which best encode the illuminant information.

Both FC4 and BoCF lead to a high performance in terms of average errors.
However, for extreme samples, these methods fail and lead to high errors. In
this paper, we address this limitation by proposing a novel scheme to
aggregate BoCF and FC4 methods to deal with estimation ambiguities. We
estimate the relative uncertainty of each approach using MC-dropout and the
final illumination estimate is obtained as the sum of the different model
estimates weighted by the log-inverse of their corresponding uncertainties.

###  II-B Monte-Carlo dropout

Dropout [35, 36] is a regularization technique employed in many deep learning
methods to avoid over-fitting and improve the robustness of the models.
Dropout is used in the training phase by discounting the output of some
neurons in hidden layers with a fixed probability rate [35]. This technique
efficiently samples from an exponential number of different networks in a
tractable and feasible way [35].

In [30], it was shown that using dropout in CNNs can be interpreted as a
Bayesian approximation of a Gaussian process [37, 38, 39]. Applying dropout in
the test time presents an approach to approximate the model uncertainty in
deep learning without sacrificing neither the computational complexity nor the
test accuracy. For this approach, called MC-dropout, the inference phase is no
longer deterministic as it depends on the randomly dropped-out neurons. For
the same test example, the model yields different predictions each time. The
variance or the standard deviation of the predictions can be interpreted as
the uncertainty of the model. MC-dropout has been applied in many contexts,
such as recurrent neural networks [40], machine translation [41], and medical
diagnostics [42, 43].

In the color constancy context, both BoCF and FC4 contain dropout layers and,
thus, they can be extended with MC-dropout to estimate their confidence along
with illumination output. With the model confidence at hand, we can treat
uncertain inputs and special cases explicitly. However, it should be noted
that using MC-dropout requires multiple forward passes of the model with each
test input. This increases the time complexity of the algorithm.

##  III Proposed framework

In the computational color constancy problem, the main task is to learn a
function $f(\textbf{X})$, which takes an image X as input and outputs an
estimate of its illumination I. The CNN-based approaches approximate the
function $f(\textbf{X})$ with a neural network $f_{\Theta}(\textbf{X})$, where
$\Theta$ is a set of parameters. The weights $\Theta$ are learned in an end-
to-end manner using back-propagation. CNNs have been recently extensively used
in illumination estimation and have often led to state-of-the-art results and
low average error rates. However, for extreme test samples, these methods fail
and lead to high errors and to an undesired distortion of the colors of the
image. This can be explained by the fact that a single model is unable to
efficiently estimate the illumination for all types of scene content and all
environments.

In our framework, i.e., MCDE, we tackle this limitation and propose combining
different CNN-based color constancy models which have dropout layers by
weighting their relative illumination estimations according to their relative
uncertainties. To estimate the uncertainty of each approach, we extend the
original CNN models using MC-dropout. Following this paradigm, we apply
dropout in the test phase. When an unseen sample $\textbf{X}_{1}$ is
presented, $f_{\Theta}(\textbf{X}_{1})$ is non-deterministic and yields
different outputs at each forward pass. For $\nu$ forward passes of $X_{1}$,
we have $\\{\textbf{f}^{1},\textbf{f}^{2},..,\textbf{f}^{\nu}\\}$ different
estimates. The model illumination estimation can be computed as the average of
the estimates:

$\textbf{I}=\frac{1}{\nu}\sum_{i=1}^{\nu}\textbf{f}^{i}.$ (1)

As I is an RGB color and has three component
$\textbf{I}=(\textbf{I}_{r},\textbf{I}_{g},\textbf{I}_{b})$, the standard
deviation of each component can be computed as follows:

$\sigma_{m}=\sqrt{\frac{1}{\nu}\sum_{i=1}^{\nu}(\textbf{f}^{i}_{m}-\textbf{I}_{m})^{2}},\hskip
5.69054ptfor\hskip 2.84526ptm=r,g,b.$ (2)

The standard deviation measures the amount of variation or dispersion of the
set of predictions from the mean [44]. A confident model is expected to output
similar values at each forward pass of $\textbf{X}_{1}$ and, thus, its
corresponding standard deviations of
$\\{\textbf{f}^{1},\textbf{f}^{2},..,\textbf{f}^{\nu}\\}$ for each component
are expected to be small, whereas if the input $\textbf{X}_{1}$ is ambiguous
and the model is not confident, the predictions set
$\\{\textbf{f}^{1},\textbf{f}^{2},..,\textbf{f}^{\nu}\\}$ is expected to have
high standard deviations for the components. Thus, the standard deviation can
be used to model the component-wise uncertainty of $f_{\Theta}(\textbf{X})$
for the input $\textbf{X}_{1}$. Moreover, the total uncertainty $\mu$ of the
model $f_{\Theta}(\textbf{X})$ on the prediction I can be approximated as the
product of the standard deviation along the three components, i.e.,
$\mu=\sigma_{r}\sigma_{g}\sigma_{b}$.

Let $\\{\textbf{I}_{1},...,\textbf{I}_{K}\\}$ be the illumination predictions
of K different CNN-based models and $\\{\sigma_{1},...,\sigma_{K}\\}$ their
corresponding uncertainties. In MCDE, we propose using a weighting scheme,
where high weights are given to the more confident models and low weights are
attributed to uncertain models. As a result, we decrease the dependency of the
framework on a single model and provide a straightforward approach to leverage
various models and combine them into a robust color constancy system. Any
function which is inversely proportional to the uncertainty can be used as a
confidence score. The confidence scores of the K models can be computed as
follows:

$c_{i}=g\left(\frac{1}{\sigma_{i}}\right),\hskip 5.69054ptfor\hskip
2.84526pti=1,..,K$ (3)

where $g(.)$ can be a monotonic increasing function, e.g., linear or
logarithm. So, the final estimate can be obtained as as the weighted sum of
$\\{\textbf{I}_{1},...,\textbf{I}_{K}\\}$. The output illuminations of the
models, $\textbf{I}_{i}$ for i=1,..,K , have a norm of 1, i.e.,
$||\textbf{I}_{i}||_{2}=1$. To not distort the norm of the final prediction,
we propose using the angular average. To this end,
$\\{\textbf{I}_{i}\\}^{K}_{i=1}$ are converted to the spherical coordinate
space $\\{\textbf{Is}_{i}=(1,\phi_{i},\varphi_{i})\\}^{K}_{i=1}$ via the
following transformations:

$\phi_{i}=tan^{-1}\left(\frac{{\textbf{I}_{i}}_{g}}{{\textbf{I}_{i}}_{r}}\right),\hskip
5.69054ptfor\hskip 2.84526pti=1,..,K$ (4)
$\varphi_{i}=tan^{-1}\left(\frac{\sqrt{\textbf{I}_{ir}+\textbf{I}_{ig}}}{{\textbf{I}_{i}}_{b}}\right),\hskip
5.69054ptfor\hskip 2.84526pti=1,..,K$ (5)

The spherical coordinates $\textbf{Is}^{est}=(1,\phi,\varphi)$ of the
illumination $\textbf{I}^{est}$ can be computed as follows:

$\phi=\sum_{i=1}^{K}c_{i}\phi_{i}.$ (6)
$\varphi=\sum_{i=1}^{K}c_{i}\varphi_{i}.$ (7)

The estimated illumination $\textbf{I}^{est}$ coordinates can thus be computed
via the following inverse transformations:

$\textbf{I}^{est}_{r}=sin(\varphi)cos(\phi).$ (8)
$\textbf{I}^{est}_{g}=sin(\varphi)sin(\phi).$ (9)
$\textbf{I}^{est}_{b}=cos(\varphi).$ (10)

$\textbf{I}^{est}$ is the output of MCDE and is used to correct the colors of
the input scene. The proposed framework enables us to leverage the benefits of
the different models while dealing with estimation ambiguities. For each
input, higher weights are given to the models with higher confidences and the
contribution of the uncertain models are suppressed. This provides a simple
yet effective way to adapt to the scene content and provide a robust
approximation of the illumination.

Figure 2: Recovery angular errors of FC4 as a function of the Recovery angular
errors of BoCF on the first fold of INTEL-TAU dataset.  Figure 3: Reproduction
angular errors of FC4 as a function of the Reproduction angular errors of BoCF
on the first fold of INTEL-TAU dataset.

##  IV Experimental results and discussion

### Dataset

To evaluate the proposed approach, we use the INTEL-TAU dataset [23]. It is
the largest publicly available high-resolution illumination estimation dataset
with 7022 images. Three different camera models were used, namely Canon 5DSR,
Nikon D810, and Sony IMX135. We evaluate the methods using the 10-folds non-
random cross-validation protocol defined in [23].

### Error metrics

Similar to prior works, we report the mean of the best 25%, the mean, the
median, the trimean, and the mean of the worst 25% of the Recovery angular
error $e_{recovery}$ [45] between the ground truth illuminant and the
estimated illuminant. For better insights on the performance of the methods,
we also provide the same statistics with the recently proposed Reproduction
angular error $e_{reproduction}$ [46]. The error metrics can be computed as
follows:

$\text{$e_{recovery}$}(\textbf{I}^{gt},\textbf{I}^{est})=\cos^{-1}\left({\frac{\textbf{I}^{gt}\textbf{I}^{est}}{\|\textbf{I}^{gt}\|\|\textbf{I}^{est}\|}}\right)$
(11)

$\text{$e_{reproduction}$}(\textbf{I}^{gt},\textbf{I}^{est})=\cos^{-1}\left({\frac{r(\textbf{I}^{gt},\textbf{I}^{est})\hskip
5.69054pt\textbf{n}}{\|r(\textbf{I}^{gt},\textbf{I}^{est})\|}}\right),$ (12)

where $\textbf{I}^{gt}$ is the ground truth illumination, $\textbf{I}^{est}$
is the estimated illumination,
$r(\textbf{I}^{gt},\textbf{I}^{est})=\textbf{I}^{gt}/\textbf{I}^{est}$ is the
element-wise division of $\textbf{I}^{gt}$ by $\textbf{I}^{est}$, and n is the
normalized unit vector, i.e.,
$\textbf{n}=[1/\sqrt{3},1/\sqrt{3},1/\sqrt{3}]^{T}$.

Figure 4: The Recovery angular errors of FC4 and BoCF as a function of the corresponding log-inverse confidence scores on the first fold of INTEL-TAU. TABLE I: Results of benchmark methods on INTEL-TAU dataset using cross-validation protocol. | $e_{recovery}$ | $e_{reproduction}$
---|---|---
Method | Best
25% | Mean | Med. | Tri. | W.
25% | Best
25% | Mean | Med. | Tri. | W.
25%
Grey-World [9] | 1.0 | 4.9 | 3.9 | 4.1 | 10.5 | 1.2 | 6.1 | 4.9 | 5.2 | 13.0
White-Patch [10] | 1.4 | 9.4 | 9.1 | 9.2 | 17.6 | 1.8 | 10.0 | 9.5 | 9.8 | 19.2
Grey Edge [11] | 1.0 | 5.9 | 4.0 | 4.6 | 13.8 | 1.2 | 6.8 | 4.9 | 5.5 | 13.5
Grey Edge 2 [11] | 1.0 | 6.0 | 3.9 | 4.8 | 14.0 | 1.2 | 6.9 | 4.9 | 5.6 | 15.7
Shades-of-Grey [12] | 0.9 | 5.2 | 3.8 | 4.3 | 11.9 | 1.1 | 6.3 | 4.7 | 5.1 | 13.9
Cheng et al. 2014 [21] | 0.7 | 4.5 | 3.2 | 3.5 | 10.6 | 0.9 | 5.5 | 4.0 | 4.4 | 12.7
Weighted GE [13] | 0.8 | 6.1 | 3.7 | 4.6 | 15.1 | 1.1 | 6.9 | 4.5 | 5.4 | 16.5
Yang et al. 2015 [7] | 0.6 | 3.2 | 2.2 | 2.4 | 7.6 | 0.7 | 4.1 | 2.7 | 3.1 | 9.6
Color Tiger [8] | 1.0 | 4.2 | 2.6 | 3.2 | 9.9 | 1.1 | 5.3 | 3.3 | 4.1 | 12.7
Greyness Index [6] | 0.5 | 3.9 | 2.3 | 2.7 | 9.8 | 0.6 | 4.9 | 3.0 | 3.5 | 12.1
PCC_Q2 [47] | 0.6 | 3.9 | 2.4 | 2.8 | 9.6 | 0.7 | 5.1 | 3.5 | 4.0 | 11.9
FFCC[14] | 0.4 | 2.4 | 1.6 | 1.8 | 5.6 | 0.5 | 3.0 | 2.1 | 2.3 | 7.1
Bianco [16] | 0.9 | 3.5 | 2.6 | 2.8 | 7.4 | 1.1 | 4.4 | 3.4 | 3.6 | 9.4
C3AE [18] | 0.9 | 3.4 | 2.7 | 2.8 | 7.0 | 1.1 | 3.9 | 3.3 | 3.5 | 8.8
BoCF [19] | 0.7 | 2.4 | 1.9 | 2.0 | 5.1 | 0.8 | 3.0 | 2.3 | 2.5 | 6.5
FC4 (VGG16) [17] | 0.6 | 2.2 | 1.7 | 1.8 | 4.7 | 0.7 | 2.9 | 2.2 | 2.3 | 6.1
MCDE (linear) | 0.6 | 2.2 | 1.7 | 1.8 | 4.7 | 0.7 | 2.8 | 2.2 | 2.3 | 6.1
MCDE (log) | 0.5 | 2.1 | 1.6 | 1.7 | 4.5 | 0.6 | 2.6 | 2.0 | 2.1 | 5.4

### Performance evaluation

In our experiments, we evaluated the performance of the following state-of-
the-art unsupervised approaches: Grey-World [9], White-Patch [10], Spatial
domain [21], Shades-of-Grey [12], and Weighted Grey-Edge [13], Greyness Index
2019 [6], Color Tiger [8], PCC_Q2 [47], and the method reported in [7]. In
addition, we considered the learning-based approach Fast Fourier Color
Constancy (FFCC) [14] and the four following CNN-based approaches: FC4 [17],
Bianco [16], C3AE [18], and BoCF [19].

In our experiments, we consider the special case of MCDE composed of K=2
models, namely FC4 and BoCF. For the function $g$ defined in Eq 3, we use the
logarithm function, i.e., $g(x)=log(x)$. This variant of the framework is
called MCDE (log). In this variant, the confidence score are thus obtained as
the log-inverse of the uncertainties. We also experimented with a basic
variant of our framework, called MCDE (linear), where $g$ is defined as the
identity function. In this variant, the confidence scores $c_{1}$ and $c_{2}$
of FC4 and BoCF defined in Eq. 3 are expressed directly as the inverse of the
uncertainty, i.e., $c_{1}=\frac{1}{\sigma_{1}}$ and
$c_{2}=\frac{1}{\sigma_{2}}$. The number of forward passes per model, $\nu$,
is fixed to 30 for both variants of our framework.

In the Figures 2 and 3, we plot the Recovery and the Reproduction angular
errors $e_{recovery}$ and $e_{reproduction}$ of FC4 as a function of the
corresponding angular errors of BoCF on the first fold of INTEL-TAU. As
illustrated in the figure, for a large part of the test images FC4 and BoCF
have different extreme cases, i.e., they succeed and fail with different test
images. This supports our assumption that by combining the methods in an
optimal way, there is potential for improved results. Combining models with
different extreme cases enables MCDE to deal with the different content and to
output a robust estimation for more images.

Figure 4 illustrates the angular errors of FC4 and BoCF as a function of the
corresponding log-inverse confidences, ie., $g(x)=log(x)$. In the optimal
case, the confidence would be high for samples with low errors and vice versa.
We note that this is true in some cases and, especially for FC4, the low
errors tend to get high confidences. This supports our assumption that the
log-inverse of the MC-dropout uncertainty used in MCDE provides a useful tool
to estimate the confidence of the color constancy models. Nevertheless, we see
that the error-confidence correlation is not optimal and, in the future, other
more advanced confidence measures should be developed and compared against the
proposed MC uncertainty.

Figure 5: Visual results on three sample of INTEL-TAU. From left to right: Input image, BoCF output, FC4 output, MCDE output, and ground truth image. Gamma correction was applied for visualization. TABLE II: Results of CNN methods on INTEL-TAU dataset using cross-validation protocol. | $e_{recovery}$ | $e_{reproduction}$
---|---|---
Method | W. 25% | W. 10% | W. 5% | W. 25% | W. 10% | W. 5%
FFCC[14] | 5.6 | 7.9 | 9.8 | 7.1 | 10.2 | 12.4
BoCF [19] | 5.1 | 7.2 | 8.7 | 6.5 | 9.4 | 11.2
FC4 [17] | 4.7 | 6.4 | 7.8 | 6.1 | 8.3 | 10.1
MCDE (linear) | 4.7 | 6.4 | 7.8 | 6.1 | 8.3 | 9.6
MCDE (log) | 4.5 | 6.2 | 7.5 | 5.4 | 8.0 | 9.7
Ideal | 3.8 | 5.3 | 6.6 | 4.7 | 6.7 | 8.2

Table I reports the competitive results of the state-of-the-art approaches
along with the results of our proposed approaches on INTEL-TAU dataset. It can
be seen that the competing CNN-based approaches outperform other methods for
both error functions except FFCC, which is a non-CNN learning-based method and
it outperforms Bianco and C3AE especially in the mean of the worst 25% metric.
Our framework yields the state-of-the-art performance and the lowest error
rates in all error metrics except the best 25%, where FFCC achieves the best
results. For handling the ambiguous input samples, the proposed framework
outperforms the existing methods. This can be seen in terms of the mean of the
worst 25%. In fact, the proposed framework achieves a $0.7^{\circ}$
improvement in this metric for the Reproduction error $e_{reproduction}$
compared to the prior methods.

The proposed framework MCDE aggregates FC4 and BoCF via a weighting schema
using their prediction uncertainty. Figure 5 presents a visualization of the
outputs of the images corrected with the different approaches. Compared to
both of these methods, Table I shows that both variants of the proposed model
achieve at least as low error rates as the best of these two. Thus, we see
that relying on the MC-dropout uncertainty of the model indeed presents an
effective way to model the confidence of the prediction. Moreover, by
comparing both variants of our framework, we see that the log-inverse presents
a more robust method for weighting the contribution of each model. In fact,
this is illustrated especially in terms of the mean of the worst 25%, where
the log-inverse yields a $0.2^{\circ}$ improvement in the Recovery error
$e_{recovery}$ and a $0.7^{\circ}$ improvement in the Reproduction error
$e_{reproduction}$.

To further highlight the robustness of the proposed framework on the extreme
cases, we report in Table II the averages of worst 25%, worst 10%, and worst
5% for $e_{recovery}$ and $e_{reproduction}$. In addition, we report the
results achieved by the hypothetical ideal combination of FC4 and BoCF, i.e.,
the results that would be obtained if always the better of the two methods was
successfully selected. By comparing the errors rates of FC4, BoCF, and our
approach, we note an improvement across all the metrics. MCDE is indeed able
to prune out some of the worst failure cases and combining different
approaches using the proposed confidence scores improves the robustness of the
approach. Moreover, we note that the ideal combination of FC4 and BoCF
achieves robust performance with more than $1^{\circ}$ improvement across all
the metrics compared to each of the combined methods alone. This confirms that
indeed both models succeed and fail in different contents justifies further
the choice of combining these two models in our MCDE framework. However, it
should also be noted that although MCDE (log) yields an improved performance,
there is still room for improvement by further improving the confidence
estimate.

##  V Conclusion

In this paper, we proposed a novel ensemble method for the computational color
constancy problem called MCDE aggregating different CNN models. In the
proposed framework, we estimate the relative uncertainty of each model for a
test sample using Monte Carlo dropout. This enables us to deal with estimation
ambiguities. The final illumination estimation is computed as the sum of the
different models’ estimates weighted by the log-inverse of their corresponding
uncertainties. This yields a robust illumination estimation system and reduces
the dependency of the global framework using a single model. Furthermore, the
proposed approach addresses the limitation of prior works and the inability of
a single model to generalize well for all types of scenes and can handle
extreme cases. We evaluated the special case of MCDE combining FC4 and BoCF on
the INTEL-TAU dataset. The proposed approach led to state-of-the-art
performance.

In future work, extensions of the proposed approach may include incorporating
more CNN-based methods in MCDE, dealing with the overconfidence bias, and
exploring different alternatives to model the uncertainty of deep learning
models in the color constancy context.

## Acknowledgment

This work has been supported by the NSF-Business Finland Center for Visual and
Decision Informatics (CVDI) project AMALIA.

## References

  * [1] M. Ebner, _Color constancy_. John Wiley & Sons, 2007, vol. 7.
  * [2] A. Alsam, “Colour constant image sharpening,” in _2010 20th International Conference on Pattern Recognition_. IEEE, 2010, pp. 4545–4548.
  * [3] K. Barnard, _Practical colour constancy_. Simon Fraser University, 1999.
  * [4] A. Gijsenij, T. Gevers, and J. Van De Weijer, “Computational color constancy: Survey and experiments,” _IEEE Transactions on Image Processing_ , vol. 20, no. 9, pp. 2475–2489, 2011.
  * [5] A. Choudhury and G. Medioni, “Color constancy using standard deviation of color channels,” in _2010 20th International Conference on Pattern Recognition_. IEEE, 2010, pp. 1722–1726.
  * [6] Y. Qian, J.-K. Kamarainen, J. Nikkanen, and J. Matas, “On finding gray pixels,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 8062–8070.
  * [7] K.-F. Yang, S.-B. Gao, and Y.-J. Li, “Efficient illuminant estimation for color constancy using grey pixels,” in _Proceedings of the IEEE conference on computer vision and pattern recognition_ , 2015, pp. 2254–2263.
  * [8] N. Banić, K. Koščević, and S. Lončarić, “Unsupervised learning for color constancy,” _arXiv preprint arXiv:1712.00436_ , 2017.
  * [9] J. Cepeda-Negrete and R. Sanchez-Yanez, “Gray-world assumption on perceptual color spaces,” in _Image and Video Technology_ , 2014, pp. 493–504.
  * [10] A. Rizzi, C. Gatta, and D. Marini, “Color correction between gray world and white patch,” in _The International Society for Optical Engineering_ , 2002\.
  * [11] J. van de Weijer, T. Gevers, and A. Gijsenij, “Edge-based color constancy,” _IEEE Transactions on Image Processing_ , pp. 2207–2214, 2007.
  * [12] G. Finlayson and E. Trezzi, “Shades of gray and colour constancy,” in _Color Imaging Conference_ , 2004, pp. 37–41.
  * [13] A. Gijsenij, T. Gevers, and J. Van De Weijer, “Physics-based edge evaluation for improved color constancy,” in _2009 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2009, pp. 581–588.
  * [14] J. T. Barron and Y.-T. Tsai, “Fast fourier color constancy,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2017.
  * [15] Z. Lou, T. Gevers, N. Hu, M. P. Lucassen _et al._ , “Color constancy by deep learning.” in _BMVC_ , 2015, pp. 76–1.
  * [16] S. Bianco, C. Cusano, and R. Schettini, “Color constancy using CNNs,” in _IEEE Conference on Computer Vision and Pattern Recognition Workshops_ , 2015, pp. 81–89.
  * [17] Y. Hu, B. Wang, and S. Lin, “FC4: Fully convolutional color constancy with confidence-weighted pooling,” in _IEEE Conference on Computer Vision and Pattern Recognition_ , 2017, pp. 4085 – 4094.
  * [18] F. Laakom, J. Raitoharju, A. Iosifidis, J. Nikkanen, and M. Gabbouj, “Color constancy convolutional autoencoder,” in _Symposium Series on Computational Intelligence_ , 2019.
  * [19] F. Laakom, N. Passalis, J. Raitoharju, J. Nikkanen, A. Tefas, A. Iosifidis, and M. Gabbouj, “Bag of color features for color constancy,” _IEEE Transactions on Image Processing, (Early Access) DOI: 10.1109/TIP.2020.3004921_ , 2020.
  * [20] J. T. Barron, “Convolutional color constancy,” _IEEE International Conference on Computer Vision_ , pp. 379–387, 2015.
  * [21] D. Cheng, D. K. Prasad, and M. S. Brown, “Illuminant estimation for color constancy: why spatial-domain methods work and the role of the color distribution,” _Journal of the Optical Society of America. A, Optics, image science, and vision_ , pp. 1049–1058, 2014.
  * [22] G. Hemrit, G. Finlayson, A. Gijsenij, P. Gehler, S. Bianco, B. Funt, M. Drew, and L. Shi, “Rehabilitating the colorchecker dataset for illuminant estimation,” in _Color and Imaging Conference_ , 2018, pp. 350–353.
  * [23] F. Laakom, J. Raitoharju, A. Iosifidis, J. Nikkanen, and M. Gabbouj, “Intel-tau: A color constancy dataset,” _arXiv preprint arXiv:1910.10404_ , 2019.
  * [24] Y. Qian, K. Chen, J.-K. Kämäräinen, J. Nikkanen, and J. Matas, “Deep structured-output regression learning for computational color constancy,” in _2016 23rd International Conference on Pattern Recognition (ICPR)_. IEEE, 2016, pp. 1899–1904.
  * [25] M. Afifi and M. S. Brown, “Deep white-balance editing,” in _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2020, pp. 1397–1406.
  * [26] S.-B. Gao, M. Zhang, C.-Y. Li, and Y.-J. Li, “Improving color constancy by discounting the variation of camera spectral sensitivity,” _JOSA A_ , vol. 34, no. 8, pp. 1448–1462, 2017.
  * [27] M. Afifi, B. Price, S. Cohen, and M. S. Brown, “When color constancy goes wrong: Correcting improperly white-balanced images,” in _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 1535–1544.
  * [28] M. Afifi and M. S. Brown, “Sensor-independent illumination estimation for dnn models,” _arXiv preprint arXiv:1912.06888_ , 2019.
  * [29] Y. Gal and Z. Ghahramani, “Dropout as a bayesian approximation: Insights and applications,” in _Deep Learning Workshop, ICML_ , vol. 1, 2015, p. 2.
  * [30] ——, “Dropout as a bayesian approximation: Representing model uncertainty in deep learning,” in _international conference on machine learning_ , 2016, pp. 1050–1059.
  * [31] S. Bianco, C. Cusano, and R. Schettini, “Single and multiple illuminant estimation using convolutional neural networks,” _IEEE Transactions on Image Processing_ , vol. 26, no. 9, pp. 4347–4362, 2017.
  * [32] S. Bianco and R. Schettini, “Color constancy using faces,” in _2012 IEEE Conference on Computer Vision and Pattern Recognition_. IEEE, 2012, pp. 65–72.
  * [33] ——, “Adaptive color constancy using faces,” _IEEE transactions on pattern analysis and machine intelligence_ , vol. 36, no. 8, pp. 1505–1518, 2014\.
  * [34] N. Passalis and A. Tefas, “Learning bag-of-features pooling for deep convolutional neural networks,” in _IEEE International Conference on Computer Vision_ , 2017.
  * [35] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” _The journal of machine learning research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [36] S. Park and N. Kwak, “Analysis on the dropout effect in convolutional neural networks,” in _Asian conference on computer vision_. Springer, 2016, pp. 189–204.
  * [37] J. M. Bernardo and A. F. Smith, _Bayesian theory_. John Wiley & Sons, 2009, vol. 405.
  * [38] A. Damianou and N. Lawrence, “Deep gaussian processes,” in _Artificial Intelligence and Statistics_ , 2013, pp. 207–215.
  * [39] A. Di Martino, E. Bodin, C. H. Ek, and N. D. Campbell, “Gaussian process deep belief networks: A smooth generative model of shape with uncertainty propagation,” in _Asian Conference on Computer Vision_. Springer, 2018, pp. 3–20.
  * [40] Y. Gal and Z. Ghahramani, “A theoretically grounded application of dropout in recurrent neural networks,” in _Advances in neural information processing systems_ , 2016, pp. 1019–1027.
  * [41] R. Sennrich, B. Haddow, and A. Birch, “Edinburgh neural machine translation systems for wmt 16,” _arXiv preprint arXiv:1606.02891_ , 2016.
  * [42] X. Yang, R. Kwitt, and M. Niethammer, “Fast predictive image registration,” in _Deep Learning and Data Labeling for Medical Applications_. Springer, 2016, pp. 48–57.
  * [43] C. Angermueller and O. Stegle, “Multi-task deep neural network to predict cpg methylation profiles from low-coverage sequencing data,” in _NIPS MLCB workshop_ , 2015.
  * [44] J. M. Bland and D. G. Altman, “Statistics notes: measurement error,” _Bmj_ , vol. 312, no. 7047, p. 1654, 1996.
  * [45] S. Hordley and G. Finlayson, “Re-evaluating colour constancy algorithms,” in _International Conference on Pattern Recognition_ , 2004.
  * [46] G. Finlayson and R. Zakizadeh, “Reproduction angular error: An improved performance metric for illuminant estimation,” _perception_ , 2014.
  * [47] F. Laakom, J. Raitoharju, A. Iosifidis, U. Tuna, J. Nikkanen, and M. Gabbouj, “Probabilistic color constancy,” _arXiv preprint arXiv:2005.02730_ , 2020\.

  *[MCDE]: Monte Carlo Dropout Ensembles
  *[CNN]: Convolutional Neural Network
  *[BoCF]: Bag of Color Features
  *[FC4]: Fully Convolutional Color Constancy
  *[CNNs]: Convolutional Neural Network


# The development of near-vent volcanic ash cloud layers due to inhomogeneous
atmospheric turbulence and relationship to wind shear

Marcus Bursik 1,*ORCID: 0000-0002-9312-5202, Qingyuan Yang 2
Adele Bear-Crozier 3, Michael Pavolonis 4 and Andrew Tupper 3
1 Center for Geohazards Studies,
University at Buffalo, Buffalo NY USA
2 Earth Observatory of Singapore,
Nanyang Technological University, Singapore
3 Bureau of Meteorology, Melbourne, Australia
4 NOAA Cooperative Institute for
Meteorological Satellite Studies University of Wisconsin, Madison WI USA
Correspondence<EMAIL_ADDRESS>Tel.: +1-716-645-4265

###### Abstract

Volcanic ash clouds often become multilayered and thin with distance from the
vent. We explore one mechanism for development of this layered structure. We
review data on the characteristics of turbulence layering in the free
atmosphere, as well as examples of observations of layered clouds both near-
vent and distally. We then explore and contrast the output of volcanic ash
transport and dispersal models with models that explicitly use the observed
layered structure of atmospheric turbulence. The results suggest that the
alternation of turbulent and quiescent atmospheric layers provides one
mechanism for development of multilayered ash clouds by modulating the manner
in which settling occurs.

##### Keywords:

ash cloud; volcanic cloud; Pinatubo

##  1 Introduction

Volcanic ash is a multi-billion dollar economic hazard to aviation, as shown
during the 2010 eruptions of Eyjafjallajökull, Iceland (Casadevall, 1994;
Mazzocchi et al., 2010). It is also a risk to flight safety, with hundreds of
encounters of varying severity recorded, and several instances of multiple
engine flame-out in flight. The International Airways Volcano Watch (IAVW),
which seeks to safely separate aircraft from volcanic ash in flight, relies on
detecting areas of ash and forecasting its future movement (Tupper et al.,
2007). However, the forecasting of ash presence and concentration is generally
poorly resolved vertically, although there is some progress in this direction,
e.g., (Heinold et al., 2012; Kristiansen et al., 2015). Aircraft flying in a
supposedly ash-contaminated region at a particular altitude may encounter no
ash or significant and potentially damaging amounts, due to the high degree of
ash stratification with altitude. The improved understanding and forecasting
of stratification would assist enormously in managing the hazard and support
the continuing development of the IAVW.

Photography and satellite imagery of numerous volcanic eruptions show that
stratification or layering of volcanic clouds is a fundamental aspect of
volcanic cloud development (Fig. 1). Lidar backscatter data have been key in
defining this layered structure in distal regions (Fig. 2). Volcanic layers
can be stratospheric as well as tropospheric.

Figure 1: Photograph and tracing of eruption of Redoubt volcano, AK, 21 Apr
1990, showing plume overshoot, main umbrella intrusion and particle rich
secondary intrusion. Tracings show growth of cloud layers through different
time steps (numbered solid lines), long dashed where partly obscured. White
dashed line, tropopause. Modified from Woods and Kienle (1994) Figure 2: (a)
CALIOP nadir LIDAR total attenuated backscatter (along track shown in (b))
showing complex layering of Eyjafjallajokull ash cloud on 8 May, 2010,
isotherms (kelvins, black dotted lines) and tropopause (black solid line). (b)
Aqua-MODIS RGB false color image (Pavolonis et al., 2013) of North Atlantic
capturing this ash cloud (pink hues).

The layers form and separate by numerous processes, some unique to volcanic
clouds. The primary volcanic cloud layer near the vent, such as the volcanic
umbrella cloud or anvil cloud, arises from the driving of hot eruptive gas and
ash parcels outward around their equilibrium level, or neutral buoyancy height
(Sparks et al., 1997). Ash accretion, ash re-entrainment, source variability –
injection of ash at different altitudes with changing eruption rate and wind
fields, and gas-ash separation cause development of multiple layers (Holasek
et al., 1996a; Tupper et al., 2004; Thorsteinsson et al., 2012). Double
diffusion and convective sediment flux to a single (Woods and Kienle, 1994;
Bursik, 1998; Hoyal et al., 1999a, b; Carazzo and Jellinek, 2012) and multiple
(Carazzo and Jellinek, 2013) levels by descending fingers that intrude below
the level of a major volcanic cloud layer have been observed, and recreated in
the laboratory. As suggested by brightness temperatures over the surface of
near-vent clouds and ground based photography, umbrella clouds can be
solitary, accompanied by a single lower intrusion resulting from re-
entrainment and column-edge downflow (Woods and Kienle, 1994; Bursik, 1998) or
accompanied by lower level skirt clouds (Barr, 1982), which may or may not
contain ash. Mechanical unmixing of particulate-laden and gaseous volcanic
cloud components has been noted as a further cause of volcanic layer formation
(Holasek et al., 1996b; Fero et al., 2009), perhaps enhanced by gravity
current slumping of the particle laden component (Prata et al., 2017).

In contrast to the well-defined and complex layers discussed above, volcanic
ash transport and dispersal models (VATDs), used to forecast ash cloud motion,
display layering of only a single ash cloud, if derived from wind shear; or
after the fact from a numerical inversion for variable source height, which is
otherwise assumed a fixed parameter ((Kristiansen et al., 2015). Although
horizontal, planview resolution can be good, VATDs thus have difficulty in
reproducing the vertical thickness and multilayering of distal volcanic clouds
(Devenish et al., 2012; Folch et al., 2012; Heinold et al., 2012). In VATDs,
dispersal in the vertical direction is commonly described by a single vertical
diffusivity, $\kappa_{z}$, which in many models is taken to have the same
value as the horizontal diffusivity, $\kappa_{h}$; the result is then uniform
or isotropic ash dispersion.

###  1.1 Problem Statement

The formation mechanism and morphology of distal ash cloud layers are poorly
understood, and one potential mechanism for their formation and
characteristics is the subject of the present contribution. Our working
hypothesis is that, with particle settling, the structure of the atmosphere
itself can cause layer formation, through the process of enhanced suspension
in vertically restricted regions of high turbulence. Given that the mechanism
explored herein results from the properties of the atmosphere itself, it may
be the primary, distal layer-forming mechanism. Understanding the reason for
the occurrence and morphology of specific ash-rich cloud layers is critical
for correctly characterizing extended persistence of ash in the atmosphere,
especially when multiple separate layers occur together. Satellite sensors
penetrate only partially into the highest layer of an optically thick cloud,
and satellite remote sensing algorithms are most sensitive to the column
integrated ash properties, when multiple optically thin layers are present. A
correct understanding of layer formation and morphology is also critical for
any attempts to construct VATD models capable of producing output consistent
with observations of vertical dispersion.

##  2 Materials and Methods

###  2.1 Background Data

####  2.1.1 Atmosphere

Volcanic clouds both create turbulence and are subject to ambient atmospheric
turbulence. In volcanic clouds near the vent, turbulence is created by both
the Rayleigh-Taylor and Kelvin-Helmholtz mechanisms, as the clouds intrude
into the atmosphere as gravity currents. Kelvin-Helmholtz instability is
driven by the shear between the intruding cloud and the atmosphere (Britter
and Simpson, 1981). Rayleigh- Taylor instability is driven by convective
sedimentation, fingering and local, eddy scale density reversal (Woods and
Kienle, 1994; Holasek et al., 1996b; Chakraborty et al., 2006).

Information of sufficient resolution in the vertical direction to discover and
characterize the layered structure of the atmosphere is obtained from airborne
measurement campaigns or rawinsonde balloon releases (Dehghan et al., 2014;
Cho et al., 2003; Pavelin et al., 2002). Several methods have been developed
to derive turbulence from rawinsonde and other high-resolution data. Vasseur
and Vanhoenacker (1998) measured changes in the refractive index structure
parameter for radio waves, as turbulence causes changes in the refractive
index, based on rawinsonde pressure, temperature, humidity, wind speed and
wind direction data. Clayson and Kantha (2008) used variations in the
potential temperature profile from an idealized profile to calculate the
Thorpe scale, and derive turbulent dissipation and diffusivity. These
estimates can be made for single rawinsonde profiles with simple calculations.
There are drawbacks of course to extrapolating such high-resolution or point
data to a regional scale because of spatial and temporal inhomogeneity; the
troposphere is highly transient and spatially variable (Clayson and Kantha,
2008; Thouret et al., 2000). Tropospheric isobaric surfaces are not
necessarily parallel to the earth’s surface, especially at fronts and in
mountain waves Sharman et al. (2012). Fronts are associated with tropopause
folds, non-horizontal isobaric surfaces separating cold from warm air, and
turbulence in folds is generated by local dynamic and convective
instabilities. Mountain waves form as the density stratified atmosphere flows
past the lee side of a mountain or mountain range. These waves can break,
resulting in local turbulence concentrated in non-horizontal layers.

In the free atmosphere, parameters such as moisture content and temperature do
not change monotonically with height; there are regions of relatively
homogeneous, convecting or turbulent air, separated by regions in which
parameters vary rapidly (Vasseur and Vanhoenacker, 1998; Sharman et al.,
2012). Both the stratosphere and the troposphere are layered on scales of
$\mathcal{O}[0.1-1km]$ (Maekawa et al., 1993; Wilson et al., 2014), but the
layers in the troposphere tend to be more transient and discontinuous (Gage et
al., 1980). In the troposphere, layers of high turbulence, including high
vertical turbulence, can be indicated by constant relative humidity (RH) or
mixing ratio, $q$ (Cho et al., 2003). RH is thus an important proxy for
turbulence intensity. (Atmospheric moisture is also important in aiding plume
lift, especially in plumes from weak sources, or at low latitude, where the
moisture content is high (Sparks et al., 1997; Tupper et al., 2009)). As a
result of the layering, large-volume or bulk turbulence is highly anisotropic
($\kappa_{h}\gg\kappa_{z}$), and only within thin, well-defined layers is it
approximately isotropic ($\kappa_{h,local}\approx\kappa_{z,local}$) (Gage et
al., 1980).

####  2.1.2 Ash Clouds

Data are available from a number of sources on the shape and structure of both
near-vent and distal ash clouds. In the near-vent region, data tend to be more
limited, due to the greater optical depth. Nevertheless, cloud brightness
temperature (BT) as measured from nadir-looking geostationary and low-earth
orbiting satellites provides much useful information, as the topography of the
top of the near-vent clouds can be quite variable, with a distinct high point
or swell above the central vent that might be many kilometers above the top of
the main umbrella cloud (Fero et al., 2009). In addition, airborne and ground-
based photography and videography have provided extensive data on the features
at the base of the main umbrella or anvil, and within the underlying cloud
layers. Visible satellite imaging of the near vent cloud top consistently
reveals strong, well-defined three-dimensional vortex structures above the
vent, which evolve to smooth, somewhat more diffuse structures in the umbrella
cloud (Pouget et al., 2016).

Although the air can be choked with opaque, diffuse ash bodies that extend to
ground level near vent (Fig. 3a, b), and although gravitational intrusions,
such as umbrella clouds, are wedge-shaped by nature and therefore of variable
depth (Fig. 3c), measurements have been made of the $\Delta$BT between cloud
top and edge of the main cloud. The results suggest that umbrella clouds at
the vent typically encompass depths of $\mathcal{O}[5]$ km (Fig. 1; Table 1),
which makes a large mass of ash available for transport at sometimes high
levels. Some of the lower near-vent clouds also become distal clouds (Fig.
3d). In some cases, therefore, it might be possible to find the entire
troposphere and even lower stratosphere charged with ash, or only a distinct
layer or two of $\mathcal{O}[5]$ km depth.

Figure 3: (a) Eruption of Pinatubo, June 15, 1991. Emission from vent as well
as pyroclastic flows results in ash injecting into the atmosphere at all
heights up to 40 km. In the public domain. (b) Eruption of Grimsvotn, May 22,
2011. Ash injected in umbrella cloud (above), and in lower intrusion from
secondary, ash-rich cloud (below). Umbrella at 18-20 km. Modified from
https://www.dailymail.co.uk/news/article-1389846/Iceland-volcano-
eruption-2011-Grimsvotn-hurles-ash-plume-12-miles-sky.html. (c) Eruption of
Calbuco, April 22, 2015. Ash injected in single umbrella cloud at c. 15 km.
Modified from https://volcano.si.edu/volcano.cfm?vn=358020. (d) Eruption of
Eyjafjallajökull, April 24 and May 1, 2010. Ash injected at multiple heights
in relatively quiescent atmosphere on April 24. On May 1, under windy
conditions, ash mostly injected from single downwind plume. Modified from
photographs by M. Gudmundsson.  Table 1: Measured main, near-vent, upper
volcanic cloud depth from top to cloud edge. Measured from geostationary
imagery in first scene after eruption start (photography for Redoubt). Data
from Bear-Crozier et al. (2020); and Pouget et al. (2013) and Holasek et al.
(1996a) (below divider)

Volcano |  Eruption start date |  Height, km ASL |  Depth, km |  No. layers
---|---|---|---|---
Tinakula |  Oct 20, 2017, 2350 UT |  16.6 |  4.9 |  1
Tinakula |  Oct 20, 2017, 1930 UT |  15.1 |  3.4 |  1
Rinjani |  Aug 1, 2016, 0345 UT |  5.5 |  4.0 |  1
Manam |  Jul 31, 2015, 0132 UT |  13.7 |  9.2$\dagger$ |  2
Sangeang Api |  May 11, 2014, 0832 UT |  15.4 |  12.0$\dagger$ |  2
Kelud |  Feb 13, 2014, 1632 UT |  15.3 |  2.8 |  1
Manam |  Jan 27, 2005, 1400 UT |  24.0 |  3.0 |  1
Manam |  Oct 24, 2004, 2325 UT |  18.5 |  1.5 |  1
Pinatubo |  Jun 15, 1991, 2241 UT |  23.6 |  4.7 |  1
Redoubt |  Apr 21, 1990, 1412 UT |  12.0 |  4.9 |  2
Mount St. Helens |  May 18, 1980, 2020 UT |  13.0 |  1.0 |  1

$\dagger$ Depth possibly overestimated.

In the distal region, airborne lidar, EARLINET-AERONET and CALIOP data
typically show much thinner, more discontinuous cloud structures (Fig. 2).
Three separate tabulations of distal ash cloud layer data for the
Eyjafjallajökull plume have been published (Jonsson et al., 1996; Ansmann et
al., 2011; Winker et al., 2012). These data suggest that distal clouds from
this tropospheric eruption were typically 0.3-3 km thick, made up of 2-3
layers, with individual layers of 0.3-1.4 km depth, and maximum age of 129 h
(¡ 1 week) (Table 2). Vernier et al. (2013) using CALIOP data, discerned two
or more well-defined layers in the cloud from Puyehue-Cordon Caulle three
weeks after the eruption. Some of the layers showed fold or wrap-around
structures (Figure 1a in Vernier et al. (2013)), perhaps related to vertical-
plane chaotic mixing (Pierce and Fairlie, 1993). Clouds were up to 3 km thick,
with individual layers of 0.1-2 km depth, in the upper troposphere-lower
stratosphere (UTLS), centered on the tropopause at 8-14 km altitude. On July
12-13, 1991, 26 days after the last major eruption of Pinatubo, a lidar flight
noted numerous stratospheric layers (Winker and Osborn, 1992). The data showed
a number of well-defined layers of 0.5-1 km depth between about 14 and 25 km
altitude (Figure 1 in Winker and Osborn (1992)). Along much of the line of
flight there were two layers, but in places there were up to five. In all
studies cited above, distal ash layers were horizontal or tilted relative to
the horizon, and had extents in the cross-transport direction of hundreds of
km in the troposphere (Eyjafjallajökull), to thousands of km in the
stratosphere (Puyehue, Pinatubo). Ash is retained longer in the stratosphere
than in the troposphere as suggested by the data cited herein.

Once the particles have propagated far from vent (¿11 hr in the case of a
large eruption, e.g., Pinatubo), they no longer retain significant memory of
source conditions (Fero et al., 2009). Back trajectories of distal ash clouds
for Eyjafjallajökull and Puyehue-Cordon Caullé are generally consistent with
theoretically possible cloud heights at the source (Winker et al., 2012;
Vernier et al., 2013). However, it is not clear that ash was injected at these
altitudes at the source, given uncertainties in vertical parcel motion and
settling speed, or lack of incorporation thereof in the models (Madankan et
al., 2014; Vernier et al., 2013). It is at least possible that some cloud
layers were generated at distance from the volcanic source.

Table 2: Measured distal Eyjafjallajökull cloud depths from CALIOP lidar. Data from Winker et al. (2012) (top), Marenco et al. (2011) (middle) and Schumann et al. (2011) (lower section). Note that number of layers varies with spatial position. Date |  Cloud |  Height range |  Depth |  Age, hr |  No. layers
---|---|---|---|---|---
|  |  km ASL |  km |  |
Apr 15 |  20100415 |  $1.41-3.23$ |  0.51 |  $<6$ |  –
Apr 16 |  20100416-a |  $3.77-5.50$ |  0.58 |  30 |  $>1$
Apr 16 |  20100416-b |  $1.97-7.27$ |  0.67 |  24 |  $>1$
Apr 17 |  20100417-a |  $0.20-6.28$ |  0.76 |  42 |  1
Apr 17 |  20100417-b |  $0.05-4.00$ |  0.61 |  42 |  1
Apr 18 |  20100418-a |  $3.14-5.59$ |  0.81 |  66 |  –
Apr 18 |  20100418-b |  $3.75-6.49$ |  0.86 |  66 |  –
Apr 19 |  20100419-a |  $3.20-5.26$ |  1.06 |  71 |  –
Apr 19 |  20100419-c |  $2.48-3.94$ |  0.45 |  30 |  –
Apr 19 |  20100419-d |  $4.63-5.20$ |  0.41 |  $114–126$ |  –
Apr 20 |  20100420 |  $0.05-1.88$ |  1.08 |  $20–24$ |  –
May 4 |  20100504 |  $2.3-5.5$ |  0.5 |  – |  $1-2$
May 5 |  20100505 |  $2.4-4.5$ |  0.9 |  – |  $1-2$
May 14 |  20100514 |  $5.1-8.1$ |  1.1 |  – |  $1-3$
May 16 |  20100516 |  $3.4-5.5$ |  1.2 |  – |  $1-3$
May 17 |  20100517 |  $3.5-5.6$ |  1.3 |  – |  $1-3$
May 18 |  20100518 |  $2.5-4.9$ |  0.9 |  – |  $1-3$
Apr 19 |  20100419-1 |  $3.9-5.6$ |  1.7 |  $105-111$ |  $>1$
Apr 19 |  20100419-2 |  $3.5-3.8$ |  0.3 |  $104-108$ |  1
Apr 19 |  20100419-3 |  $3.9-4.2$ |  0.3 |  $105-108$ |  1
Apr 22 |  20100422-4 |  $0.7-5.5$ |  – |  $49-50$ |  diffuse
Apr 23 |  20100423-5 |  $2.1-3.4$ |  1.3 |  $40-58$ |  $>1$
May 2 |  20100502-6 |  $1.6-3.7$ |  2.1 |  $7.1-12$ |  $>1$
May 9 |  20100509-7 |  $3.5-4.9$ |  1.4 |  $97-129$ |  1
May 13 |  20100513-8 |  $2.8-5.4$ |  $0.4-0.7$ |  $71-78$ |  1 tilted
May 16 |  20100516-9 |  $3.6-7.0$ |  3.4 |  $58-66$ |  $>1$
May 17 |  20100517-10 |  $3.2-6.3$ |  3.1 |  $66-88$ |  $>1$
May 18 |  20100518-11 |  $2.8-3.4$ |  0.6 |  $81-100$ |  1
May 18 |  20100518-12 |  $4.0-5.7$ |  1.7 |  $66-78$ |  $>1$

The data presented herein, in text, Figures 2, 3, and Tables 1, 2, suggest
that the number of layers increases with time, and the depth decreases (see
also Dacre et al. (2015)). Volcanic source conditions are eventually lost
after an eruption, meaning that the atmosphere completely controls cloud
shape. Layers are more transient in the troposphere than in the stratosphere,
as mentioned for certain eruptions (also, cf. data in Table 1 with 2; (Thouret
et al., 2000)). In summary, the features of volcanic clouds include the
following:

  1. 1.

Near-vent clouds often evolve from sharply defined with clear eddy or vortex
structure to diffuse or diaphanous, into distal, thin well-defined layers with
sharp, smooth boundaries

  2. 2.

There are often multiple distal layers

  3. 3.

Collocation in position in planview is common (stacking), but not pervasive

  4. 4.

Sometimes distal cloud forms are horizontal; sometimes sloping or tilted
relative to horizontal

  5. 5.

Separate particle clouds exist, not only separate gas and ash clouds

  6. 6.

Horizontal extent $\gg$ vertical extent

  7. 7.

Vertical extent of near-vent layers is $\mathcal{O}[5]$ km, which with time
results in

  8. 8.

Vertical extent of single, distal layers being $\mathcal{O}[0.1-1]$ km

The present contribution seeks to provide explanations for some of these
features.

###  2.2 Model

The advection-diffusion equation forms the basis for all VATD models. To
illustrate differences and potential problems in implementation of VATD and
the underlying physics, in this section, we introduce two zeroth-order
simplifications. Our goal is to contrast the basic behavior of an ash cloud
under conditions of isotropic (or constant $\kappa_{h}$ and $\kappa_{z}$)
turbulence, as assumed in the models, and layered turbulence, as we find in
the atmosphere. We base our modeling on generation of synthetic atmospheres
with and without multiple turbulent layers separated by relatively quiescent
air.

In the case of isotropic turbulence, we begin by assuming Cartesian
coordinates, $(x,y,z)$, with velocity components, $(u,v,w)$. The three
components of the turbulent diffusivity, $(\kappa_{x},\kappa_{y},\kappa_{z})$
are the same, $\kappa$. The concentration of particles in the $i$-size
fraction, $C_{i}$, varies in time, $t$ and space as:

$\frac{\partial C_{i}}{\partial t}+\frac{\partial}{\partial
x}\left(uC_{i}\right)+\frac{\partial}{\partial
y}\left(vC_{i}\right)+\frac{\partial}{\partial
z}\left(wC_{i}\right)=\frac{\partial^{2}}{\partial x^{2}}\left(\kappa
C_{i}\right)+\frac{\partial^{2}}{\partial y^{2}}\left(\kappa
C_{i}\right)+\frac{\partial^{2}}{\partial z^{2}}\left(\kappa
C_{i}\right)+\Phi$ (1)

where $\Phi$ represents the source/sink function, which in the case of ash
clouds is mostly represented by aggregation and disaggregation of small
particles. In the present case, such processes are set to zero. We assume a
two-dimensional system with a point-source in time and space, $w=w_{s}$, the
settling speed, and that, following a streamtube, the motion of the volcanic
cloud can be characterized by a single downwind coordinate direction $s$ – for
which the axis is everywhere tangent to the plume centerline, e.g., (Wright,
1977; Hopkins and Bridgman, 1985) – and speed $U$ in that direction. Under
these assumptions, the advection-diffusion equation becomes:

$\frac{\partial C_{i}}{\partial t}+\frac{\partial}{\partial
s}\left(UC_{i}\right)+\frac{\partial}{\partial
z}\left(w_{s}C_{i}\right)=\frac{\partial^{2}}{\partial s^{2}}\left(\kappa
C_{i}\right)+\frac{\partial^{2}}{\partial z^{2}}\left(\kappa C_{i}\right)$ (2)

with the well-known solution for the impulse initial condition (Csanady, 1980;
Roberts and Webster, 2002):

$C_{i}(s,z,t)=\frac{C_{i0}}{4\pi\kappa
t}\exp\left[-\frac{(s-s_{0}-Ut)^{2}+(z-z_{0}-w_{s}t)^{2}}{4\kappa t}\right].$
(3)

It is reasonably clear that the solution is a Gaussian in $(s,z)$, in which
ash spreads, settles and is blown downwind with time.

In the second case, of layered turbulence, more realistic for the free
atmosphere, we assume particles filling a layer of finite vertical extent. Due
to its internal turbulence, concentration varies in $t$ but not in $s$ or $z$
within the layer, $C_{i}(s,z,t)=C_{i}(t)$. There is no flux at the upper
boundary of such a layer, as any particles injected upwards by eddies will
settle back down into the layer. There is a flux boundary condition at the
lower boundary where $\kappa_{z}\rightarrow 0$, and in this case, the
advection-diffusion equation becomes, at the lower boundary:

$\frac{\partial}{\partial t}\left(C_{i}\right)+\frac{\partial}{\partial
z}\left(w_{s}C_{i}\right)=0$ (4)

If furthermore it can be assumed that some particles near the cloud base fall
from the turbulent layer when their weight overcomes the internal turbulence
at the lower cloud edge, thus developing a step-like concentration gradient at
the base of the streamtube, then:

$C_{i}(t,z)=H(z)C_{i}(t)$ (5)

where $H(z)$ is the Heaviside step function, then:

$\frac{\partial}{\partial t}\left(C_{i}\right)+\frac{\partial}{\partial
z}\left(w_{s}C_{i}\right)=\frac{\partial C_{i}}{\partial
t}+w_{s}C_{i}\frac{\partial{H}(z)}{\partial z}=0$ (6)

Integrating through the layer depth, $h$:

$\frac{\partial C_{i}}{\partial
t}\int_{0}^{h}dz=-w_{s}C_{i}\int_{0}^{h}\delta(z)dz$ (7)

we obtain:

$\frac{dC_{i}}{dt}=-\frac{w_{s}}{h}C_{i}$ (8)

which has solution:

$C_{i}=C_{i0}\exp\left(-\frac{w_{s}(t-t_{0})}{h}\right)$ (9)

In the quiescent layer below the boundary, particles only settle and are
advected downwind, there is no turbulence mechanism to enhance persistence
within the layer. Thus, turbulent layers can retain particles longer than do
quiescent layers because of continuing re-entrainment in eddies. Particles
fall relatively rapidly through the quiescent layers because of uninhibited
settling, sometimes even enhanced by the effects of convective sedimentation
(Hoyal et al., 1999a), which is not included in the present model.

Based on similarity theory, the timescales for the processes under the
different particle transport conditions arising in different layers can be
used to examine the conditions under which diffusion or settling dominates.
From Eq 3, the timescale of vertical diffusion, $\tau_{1}$, through a layer of
depth, $h$, is given as $\tau_{1}=\frac{h^{2}}{\kappa}$. From Eq 9, the
timescale of settling through the same layer, $\tau_{2}$, is
$\tau_{2}=\frac{h}{w_{s}}$. The ratio of the two timescales indicates
domination of particle transport by settling or dispersion in the vertical
direction. The ratio is given by the dimensionless group,$\Pi_{1}$:

$\Pi_{1}=\frac{\tau_{1}}{\tau_{2}}=\frac{hw_{s}}{\kappa}$ (10)

In the following section, we explore results from these simplifications and
the similarity analysis, as well as numerical solutions to more complicated
cases. Numerical solutions are provided for the Ash3D VATD model (Schwaiger et
al., 2012), as well as both Eulerian and Lagrangian model codes (Table 3).

##  3 Results

Following from Eq 3, spread from a point source in a VATD model, with
isotropic turbulence and a wind of constant speed with height, is shown in
Figure 4a, b. Ash diffuses and progressively spreads from the source as the
center of mass descends at the settling speed. Using a higher settling speed,
the rate at which the center of mass descends increases, but the rate at which
the particles disperse from the center of mass remains constant. Thus, at any
one height below the source, particles with a higher settling speed should be
spread less distance from the source.

Figure 4: Cross sections through Ash3D output, with: (a) instantaneous source, showing advection, settling and isotropic dispersal of ash. 10 min after release. (b) instantaneous source, showing advection, settling and isotropic dispersal of ash. 120 min after release. (c) maintained source, showing advection with wind shear, settling under isotropic turbulence. 10 min after release. (d) maintained source, showing advection with wind shear, settling in isotropic turbulence. 120 min after release. Note that maintained source and wind shear result together in elongated dispersal pattern, and layer development. Red cross, source location. Table 3: Simulation parameters. Duration refers to emission from vent.  Parameter |  $\kappa=$const |  $\kappa=$const |  Layered Atmosphere |  |
---|---|---|---|---|---
|  wind=const |  windshear |  |  |
Simulation type |  Ash3D |  Ash3D |  Lagrangian-1 |  Lagrangian-2 |  Eulerian
Source type |  Point |  Point |  Point |  Point |  Layer
Source height, km |  24.75 |  24.75 |  24.75 |  4.2 |  $3.8-4.2$
Particle size, $\mu$m |  1000 |  1000 |  1000 |  40 |  40
Settling speed, m/s |  3 |  3 |  3 |  0.02 |  0.02
Amount |  0.01 km3 |  0.01 km3 |  1000 parcels |  1000 parcels |  10000 particles
Duration |  0.2 hr |  0.2 hr |  Instantaneous |  Instantaneous |  Instantaneous
Turbulent layer heights |  – |  – |  $24-25$, $14.5-19$ km |  $3.8-4.2$, $2-2.7$ km |  $3.8-4.2$, $2-2.7$ km
Wind speed |  5 m/s all elev |  Below 19.75 km: 0 m/s |  5 m/s |  10 m/s |  0 m/s
|  |  19.75 to 21.75 km: 10 increasing to 50 m/s |  |  |
|  |  21.75 to 23.75 km: 50 decreasing to 10 m/s |  |  |
|  |  Above 23.75 km: 10 m/s |  |  |

Dispersion in the presence of wind shear is shown in Figure 4c, d. The shear
distorts the dispersal pattern from the idealized, spherically symmetric
pattern seen in the constant wind field, causing an elongation in the
dispersal pattern centered at the wind speed maximum. Because the deformation
is by simple shear, cloud thinning does not occur. Wind-shear produced
elongation thus creates a volcanic cloud layer that continues to deepen by
diffusion.

In a layered atmosphere, we refer to Eq 10 to explore asymptotic behavior. In
layers for which $\Pi_{1}>1$, the diffusivity is low relative to the settling
speed, the timescale of diffusion is therefore long, hence motion is
controlled by settling. In layers for which $\Pi_{1}<1$, the diffusivity is
high relative to the settling speed, the timescale of settling is long, hence
motion is dominated by diffusion. Note also that as $h$ increases, the
timescale of diffusion increases faster than does that for settling, meaning
it becomes more likely the particles will exit a layer by settling than by
diffusion. For a typical diffusivity of $\kappa=1$ m2/s in the UTLS (Wilson,
2004) at 10 km altitude, and layer of depth $h=1$ km, the critical settling
speed, $w_{s,crit}$, dividing settling from diffusion dominated motion is c. 1
m/s, which would correspond to a pumice particle of diameter c. 100 $\mu$m at
about 2400 kg/cu m (e.g., (Scollo et al., 2005)).

Results for a simple layered system are shown in Figures 5 and 6. These
figures are simplified from the observations of Cho et al. (2003), who point
out two layers of especially striking turbulence from 2-2.7 km and 3.8-4.2 km
(Fig. 5a), but do not give specific values for turbulence intensity in any
layers. We therefore apply turbulence in these two layers through a Lagrangian
random walk, and in other layers, no turbulence. Particles in the turbulent
layers, then, initiate the random walk, being “stuck” within the eddies of the
turbulent layer (Fig. 5b). Consider the lower boundary of a layer with strong
turbulence. All the particles there are subject to a random walk. They have a
50% percent probability of going up, and a 50% probability of going down.
Those particles sent above the lower boundary due to turbulence are sent to a
position higher above the boundary than their original position at the
boundary. This will give them a greater chance to spend a longer time in the
turbulent layer, whether or not one considers settling. Thus, for those layers
dominated by the random walk and dispersion, $\Pi_{1}<1$, Eq 9 holds, and the
behavior seen in Figures 5 (purple layers) and 6 (red and blue lines) occurs.
Particles accumulate in the lowermost turbulent layers (e.g., blue curve in
Figure 6), and once reaching a peak, settle out only slowly. Thus, the lower,
turbulent layers are the most likely ones in which to observe particles on
longer timescales.

Figure 5: (a) Turbulent tropospheric layers (shaded) indicated by constant
mixing ratio ($q$), high turbulent energy dissipation rate
($\log\varepsilon$), and bounded by high shear ($(dU/dz)^{2}$). Modified from
Cho et al. (2003). Shaded layers are used in simplified layer models
Lagrangian-2 and Eulerian (Table 3) (b) Lagrangian paths of two volcanic
particles settling through the turbulent layers, in which suspension is
enhanced. This is output from model Lagrangian-2.  Figure 6: Eulerian model
(Table 3) of number of particles in layered atmosphere. Middle layer is not
turbulent (green dashed curve); upper and lower layers are turbulent (red and
blue solid curves)

The two modes of behavior, isotropic diffusion and layered diffusion, contrast
starkly, as seen in explicitly comparing VATD output with that from the
Lagrangian-1 model (Table 3; Fig. 7). Isotropic diffusion alone cannot lead to
generation of ash layers away from source, but layered diffusion can, even in
a constant windfield.

##  4 Discussion and Conclusions

In the present work, we have presented data and models on near-vent and distal
volcanic cloud morphology and loading. We have performed numerical experiments
comparing dispersal in an atmosphere with constant $\kappa$ and wind, constant
$\kappa$ and wind shear, and variable $\kappa_{z}$ with height. The
observational data suggest that more distal clouds of depth $\mathcal{O}$[0.1
to 1] km develop from near-vent clouds of depth generally $1-5$ km. The
downwind clouds occur at heights consistent with the original eruption column
heights for both tropospheric and stratospheric eruptions. The depth range of
the distal layers, being markedly less than the near-vent depth range, and the
common stacking of more distal ash cloud layers, suggest that their
development is controlled by atmospheric processes. The observations are
consistent with the working hypothesis that the layering of the atmosphere in
turbulence intensity, causing alternating suspension and settling dominated
behavior of particles, is a cause of distal layer morphology.

The present model outputs (Figs. 4, 7) are consistent with those of VATD and
backtrajectory models that layers can be produced under conditions of certain
source characteristics or wind shear (Devenish et al., 2012; Folch et al.,
2012; Heinold et al., 2012; Winker et al., 2012; Vernier et al., 2013). The
model outputs assuming the atmosphere has multiple turbulent layers support
the working hypothesis regarding the effects of an atmosphere layered with
respect to turbulence intensity. Thus, in addition to the near-vent, volcanic
and sedimentation processes causing formation of volcanic layers, processes
associated with atmospheric layered turbulence also produce layering in ash
clouds. These are often multilayered due to multiple, alternating layers of
turbulent and quiescent air. The distal ash layers furthermore scale to the
depth of the alternating turbulent and quiescent layers in the atmosphere,
which are $\mathcal{O}[0.1-5]$ km deep.

Figure 7: (a) Cross section through Ash3D output (Table 3), with instantaneous
source and no wind shear, showing advection, settling and isotropic dispersal
of ash. (b) Cross section through Lagrangian-1 dispersion model (Table 3 with
turbulence layered atmosphere, showing advection, settling and non-isotropic
ash dispersal. Red $\times$ is point of origin for particles, and color
gradient is scaled to concentration in both.

Because of the alternating turbulent and quiescent structure in the
troposphere and stratosphere, volcanic ash clouds tend to separate vertically
over time, lending to them the distinct layering or banded appearance in
imagery. The turbulent layers retain particles longer than the quiescent
layers because the turbulence retains particles in suspension. Particles fall
more rapidly through the relatively quiescent layers (lower $\kappa_{z}$) by
single particle settling, or because of convective sedimentation.

The results suggest that to better model the position and morphology of ash
clouds for aviation safety and other purposes in VATDs, the vertical
characteristics of the atmosphere need to be better resolved than is typical
at present. Because of the importance of turbulence and moisture to layer
formation, it is critical that these two parameters especially be estimated
well, and at as high a vertical resolution as possible.

##### Abbreviations.

The following abbreviations are used in this manuscript:

AERONET | AErosol RObotic NETwork
---|---
BT | Brightness Temperature
CALIOP | Cloud-Aerosol Lidar with Orthogonal Polarization
EARLINET | European Aerosol Research Lidar Network
IAVW | International Airways Volcano Watch
RH | Relative Humidity
UTLS | Upper Troposphere - Lower Stratosphere
VATD | Volcanic Ash Transport and Dispersal

## References

  * Ansmann et al. [2011] Ansmann, A., et al., Ash and fine-mode particle mass profiles from earlinet-aeronet observations over central europe after the eruptions of the eyjafjallajökull volcano in 2010, Journal of Geophysical Research: Atmospheres, 116(D20), doi:10.1029/2010JD015567, 2011.
  * Barr [1982] Barr, S., Skirt clouds associated with the soufriere eruption of 17 april 1979, Science, 216(4550), 1111–1112, 1982.
  * Bear-Crozier et al. [2020] Bear-Crozier, A., S. Pouget, M. Bursik, E. Jansons, J. Denman, A. Tupper, and R. Rustowicz, Automated detection and measurement of volcanic cloud growth: towards a robust estimate of mass flux, mass loading and eruption duration, Natural Hazards, pp. 1–38, 2020.
  * Britter and Simpson [1981] Britter, R. E., and J. E. Simpson, A note on the structure of the head of an intrusive gravity current, Journal of Fluid Mechanics, 112, 459–466, doi:10.1017/S0022112081000517, 1981.
  * Bursik [1998] Bursik, M., Tephra dispersal, in The Physics of Explosive Volcanic Eruptions, Special Publications, vol. 145, edited by J. S. Gilbert and R. S. J. Sparks, pp. 117–146, Geological Society of London, illus. incl. sketch maps, 1998.
  * Carazzo and Jellinek [2012] Carazzo, G., and A. M. Jellinek, A new view of the dynamics, stability and longevity of volcanic clouds, Earth and Planetary Science Letters, 325, 39–51, 2012.
  * Carazzo and Jellinek [2013] Carazzo, G., and A. M. Jellinek, Particle sedimentation and diffusive convection in volcanic ash-clouds, Journal of Geophysical Research, 118, 1420–1437, 2013.
  * Casadevall [1994] Casadevall, T. J., Volcanic ash and aviation safety: Proceedings of the first international symposium on volcanic ash and aviation safety, U.S. Geological Survey Bulletin, 2047, 450, 1994.
  * Chakraborty et al. [2006] Chakraborty, P., G. Gioia, and S. Kieffer, Volcán reventador’s unusual umbrella, Geophysical research letters, 33(5), 2006.
  * Cho et al. [2003] Cho, J. Y., R. E. Newell, B. E. Anderson, J. D. Barrick, and K. L. Thornhill, Characterizations of tropospheric turbulence and stability layers from aircraft observations, Journal of Geophysical Research: Atmospheres, 108(D20), 2003.
  * Clayson and Kantha [2008] Clayson, C. A., and L. Kantha, On turbulence and mixing in the free atmosphere inferred from high-resolution soundings, Journal of Atmospheric and Oceanic Technology, 25(6), 833–852, 2008.
  * Csanady [1980] Csanady, D. T., Turbulent diffusion in the environment, D.Reidel Publishing Co, p. 248p, dordrecht, 1980.
  * Dacre et al. [2015] Dacre, H., A. Grant, N. Harvey, D. Thomson, H. Webster, and F. Marenco, Volcanic ash layer depth: Processes and mechanisms, Geophysical Research Letters, 42(2), 637–645, 2015.
  * Dehghan et al. [2014] Dehghan, A., W. K. Hocking, and R. Srinivasan, Comparisons between multiple in-situ aircraft turbulence measurements and radar in the troposphere, Journal of Atmospheric and Solar-Terrestrial Physics, 118, 64–77, 2014.
  * Devenish et al. [2012] Devenish, B., D. Thomson, F. Marenco, S. Leadbetter, H. Ricketts, and H. Dacre, A study of the arrival over the united kingdom in april 2010 of the eyjafjallajökull ash cloud using ground-based lidar and numerical simulations, Atmospheric Environment, 48, 152 – 164, doi:https://doi.org/10.1016/j.atmosenv.2011.06.033, volcanic ash over Europe during the eruption of Eyjafjallajöekull on Iceland, April-May 2010, 2012\.
  * Fero et al. [2009] Fero, J., S. N. Carey, and J. T. Merrill, Simulating the dispersal of tephra from the 1991 pinatubo eruption: implications for the formation of widespread ash layers, Journal of Volcanology and Geothermal Research, 186(1-2), 120–131, 2009.
  * Folch et al. [2012] Folch, A., A. Costa, and S. Basart, Validation of the fall3d ash dispersion model using observations of the 2010 eyjafjallajökull volcanic ash clouds, Atmospheric Environment, 48, 165 – 183, doi:https://doi.org/10.1016/j.atmosenv.2011.06.072, volcanic ash over Europe during the eruption of Eyjafjallajöekull on Iceland, April-May 2010, 2012\.
  * Gage et al. [1980] Gage, K., J. Green, and T. VanZandt, Use of doppler radar for the measurement of atmospheric turbulence parameters from the intensity of clear-air echoes, Radio Science, 15(2), 407–416, 1980.
  * Heinold et al. [2012] Heinold, B., I. Tegen, R. Wolke, A. Ansmann, I. Mattis, A. Minikin, U. Schumann, and B. Weinzierl, Simulations of the 2010 eyjafjallajökull volcanic ash dispersal over europe using cosmo–muscat, Atmospheric Environment, 48, 195 – 204, doi:https://doi.org/10.1016/j.atmosenv.2011.05.021, volcanic ash over Europe during the eruption of Eyjafjallajöekull on Iceland, April-May 2010, 2012\.
  * Holasek et al. [1996a] Holasek, R. E., S. Self, and A. W. Woods, Satellite observations and interpretation of the 1991 mount pinatubo eruption plumes, Journal of Geophysical Research, 12(27), 635–27, 1996a.
  * Holasek et al. [1996b] Holasek, R. E., A. W. Woods, and S. Self, Experiments on gas-ash separation processes in volcanic umbrella clouds, Journal of Volcanology Geothermal Research, 70, 169–181, 1996b.
  * Hopkins and Bridgman [1985] Hopkins, A. T., and C. J. Bridgman, A volcanic ash transport model and analysis of mount st. helens ashfall, Journal of Geophysical Research, 90, 10,620–10,630, 1985.
  * Hoyal et al. [1999a] Hoyal, D. C. J. D., M. I. Bursik, and J. F. Atkinson, Setting-driven convection; a mechanism of sedimentation from stratified fluids, Journal of Geophysical Research, C, Oceans, 104(4), 7953–7966, illus. incl. 1 table, 1999a.
  * Hoyal et al. [1999b] Hoyal, D. C. J. D., M. I. Bursik, and J. F. Atkinson, The influence of diffusive convection on sedimentation from buoyant plumes, Marine Geology, 159(1-4), 205–220, 1999b.
  * Jonsson et al. [1996] Jonsson, H. H., J. C. Wilson, C. A. Brock, J. E. Dye, G. V. Ferry, and K. R. Chan, Evolution of the stratospheric aerosol in the northern hemisphere following the june 1991 volcanic eruption of mount pinatubo: Role of tropospheric-stratospheric exchange and transport, Journal of Geophysical Research, 101, 1553–1570, 1996.
  * Kristiansen et al. [2015] Kristiansen, N. I., A. Prata, A. Stohl, and S. A. Carn, Stratospheric volcanic ash emissions from the 13 february 2014 kelut eruption, Geophysical Research Letters, 42(2), 588–596, 2015.
  * Madankan et al. [2014] Madankan, R., et al., Computation of probabilistic hazard maps and source parameter estimation for volcanic ash transport and dispersion, Journal of Computational Physics, 271(0), 39 – 59, doi:http://dx.doi.org/10.1016/j.jcp.2013.11.032, frontiers in Computational Physics Modeling the Earth System, 2014.
  * Maekawa et al. [1993] Maekawa, Y., S. Fukao, M. Yamamoto, M. D. Yamanaka, T. Tsuda, S. Kato, and R. F. Woodman, First observation of the upper stratospheric vertical wind velocities using the jicamarca vhf radar, Geophysical research letters, 20(20), 2235–2238, 1993.
  * Marenco et al. [2011] Marenco, F., B. Johnson, K. Turnbull, S. Newman, J. Haywood, H. Webster, and H. Ricketts, Airborne lidar observations of the 2010 eyjafjallajökull volcanic ash plume, Journal of Geophysical Research: Atmospheres, 116(D20), 2011.
  * Mazzocchi et al. [2010] Mazzocchi, M., F. Hansstein, and M. Ragona, The 2010 volcanic ash cloud and its financial impact on the european airline industry, in CESifo Forum, vol. 11, pp. 92–100, München: ifo Institut für Wirtschaftsforschung an der Universität München, 2010.
  * Pavelin et al. [2002] Pavelin, E., J. A. Whiteway, R. Busen, and J. Hacker, Airborne observations of turbulence, mixing, and gravity waves in the tropopause region, Journal of Geophysical Research: Atmospheres, 107(D10), ACL–8, 2002.
  * Pavolonis et al. [2013] Pavolonis, M. J., A. K. Heidinger, and J. Sieglaff, Automated retrievals of volcanic ash and dust cloud properties from upwelling infrared measurements, Journal of Geophysical Research: Atmospheres, 118(3), 1436–1458, doi:10.1002/jgrd.50173, 2013.
  * Pierce and Fairlie [1993] Pierce, R. B., and T. D. A. Fairlie, Chaotic advection in the stratosphere: Implications for the dispersal of chemically perturbed air from the polar vortex, Journal of Geophysical Research: Atmospheres, 98(D10), 18,589–18,595, 1993.
  * Pouget et al. [2013] Pouget, S., M. Bursik, P. Webley, J. Dehn, and M. Pavolonis, Estimation of eruption source parameters from umbrella cloud or downwind plume growth rate, Journal of Volcanology and Geothermal Research, 258, 100–112, 2013.
  * Pouget et al. [2016] Pouget, S., M. Bursik, C. G. Johnson, A. J. Hogg, J. C. Phillips, and R. S. J. Sparks, Interpretation of umbrella cloud growth and morphology: implications for flow regimes of short-lived and long-lived eruptions, Bulletin of Volcanology, 78(1), 1–19, doi:10.1007/s00445-015-0993-0, 2016.
  * Prata et al. [2017] Prata, F., M. Woodhouse, H. E. Huppert, A. Prata, T. Thordarson, and S. Carn, Atmospheric processes affecting the separation of volcanic ash and so2 in volcanic eruptions: inferences from the may 2011 grímsvötn eruption, Atmospheric Chemistry and Physics, 17(17), 10,709–10,732, 2017\.
  * Roberts and Webster [2002] Roberts, P. J., and D. R. Webster, Turbulent diffusion, in Environmental Fluid Mechanics: Theories and Applications, edited by H. H. Shen, A. H. Cheng, K.-H. Wang, M. H. Teng, and C. C. Liu, pp. 7–45, ASCE Press, Reston, Virginia, 2002.
  * Schumann et al. [2011] Schumann, U., et al., Airborne observations of the eyjafjalla volcano ash cloud over europe during air space closure in april and may 2010, Atmos. Chem. Phys., 11, 2245–2279, doi:10.5194/acp-11-2245-2011, 2011.
  * Schwaiger et al. [2012] Schwaiger, H. F., R. P. Denlinger, and L. G. Mastin, Ash3d: A finite-volume, conservative numerical model for ash transport and tephra deposition, Journal of Geophysical Research: Solid Earth, 117(B4), 2012\.
  * Scollo et al. [2005] Scollo, S., M. Coltelli, F. Prodi, M. Folegani, and S. Natali, Terminal settling velocity measurements of volcanic ash during the 2002–2003 etna eruption by an x-band microwave rain gauge disdrometer, Geophysical Research Letters, 32(10), 2005.
  * Sharman et al. [2012] Sharman, R. D., S. B. Trier, T. P. Lane, and J. D. Doyle, Sources and dynamics of turbulence in the upper troposphere and lower stratosphere: A review,, Geophys Res Lett, 39, doi:10.1029/2012GL051996, 2012.
  * Sparks et al. [1997] Sparks, R. S. J., M. I. Bursik, S. N. Carey, J. S. Gilbert, L. S. Glaze, H. Sigurdsson, and A. W. Woods, Volcanic Plumes, John Wiley & Sons, London, 574p., 1997.
  * Thorsteinsson et al. [2012] Thorsteinsson, T., T. Jóhannsson, A. Stohl, and N. I. Kristiansen, High levels of particulate matter in iceland due to direct ash emissions by the eyjafjallajökull eruption and resuspension of deposited ash, Journal of Geophysical Research: Solid Earth, 117(B9), 2012\.
  * Thouret et al. [2000] Thouret, V., J. Y. Cho, R. E. Newell, A. Marenco, and H. G. Smit, General characteristics of tropospheric trace constituent layers observed in the mozaic program, Journal of Geophysical Research: Atmospheres, 105(D13), 17,379–17,392, 2000.
  * Tupper et al. [2004] Tupper, A., S. Carn, J. Davey, Y. Kamada, R. Potts, F. Prata, and M. Tokuno, An evaluation of volcanic cloud detection techniques during recent significant eruptions in the western ’ring of fire’, Remote Sensing of Environment, 91(1), 27–46, 2004.
  * Tupper et al. [2007] Tupper, A., I. Itikarai, M. Richards, F. Prata, S. Carn, and D. Rosenfeld, Facing the challenges of the international airways volcano watch: the 2004/05 eruptions of manam, papua new guinea, Weather and Forecasting, 22(1), 175–191, 2007.
  * Tupper et al. [2009] Tupper, A., C. Textor, M. Herzog, and H.-F. Graf, Tall clouds from small eruptions: modelling the sensitivity of eruption height and fine ash fallout to tropospheric instability, Natural Hazards, doi:10.1007/s11069-009-9433-9, 2009.
  * Vasseur and Vanhoenacker [1998] Vasseur, H., and D. Vanhoenacker, Characteristics of tropospheric turbulent layers from radiosonde data, Electronic Letters, 34, 318–319, 1998.
  * Vernier et al. [2013] Vernier, J.-P., et al., An advanced system to monitor the 3d structure of diffuse volcanic ash clouds, Journal of applied meteorology and climatology, 52(9), 2125–2138, 2013.
  * Wilson [2004] Wilson, R., Turbulent diffusivity in the free atmosphere inferred from mst radar measurements: a review, Annales Geophysicae, 22, 3869–3887, 2004.
  * Wilson et al. [2014] Wilson, R., H. Luce, H. Hashiguchi, N. Nishi, and Y. Yabuki, Energetics of persistent turbulent layers underneath mid-level clouds estimated from concurrent radar and radiosonde data, Journal of Atmospheric and Solar-Terrestrial Physics, 118, 78–89, 2014.
  * Winker and Osborn [1992] Winker, D., and M. Osborn, Preliminary analysis of observations of the pinatubo volcanic plume with a polarization-sensitive lidar, Geophysical research letters, 19(2), 171–174, 1992.
  * Winker et al. [2012] Winker, D., Z. Liu, A. Omar, J. Tackett, and D. Fairlie, Caliop observations of the transport of ash from the eyjafjallajökull volcano in april 2010, Journal of Geophysical Research: Atmospheres, 117(D20), 2012\.
  * Woods and Kienle [1994] Woods, A. W., and J. Kienle, The dynamics and thermodynamics of volcanic clouds: Theory and observations from the april 15 and april 21, 1990 eruptions of redoubt volcano, alaska, Journal of Volcanology and Geothermal Research, 62(1-4), 273–299, 1994.
  * Wright [1977] Wright, S. J., Effects of ambient crossflows and density stratification on the characteristic behavior of round turbulent buoyant jets, W.M. Keck Laboratory, Caltech, Pasadena, Report KH-R-36, 1977.