Cloudy Example

Overview

This example is based on Cloudy, a microphysics model that simulates how cloud droplets collide and coalesce into larger drops. Collision-coalescence is a crucial process for the formation of rain.

Cloudy is initialized with a mass distribution of the cloud droplets; this distribution is then evolved in time, with more and more droplets colliding and combining into bigger drops according to the droplet-droplet interactions specified by a collision-coalescence kernel. The evolution is completely determined by the shape of the initial distribution and the form of the kernel.

This example shows how ensemble Kalman methods can be used to learn the parameters of the initial cloud droplet mass distribution from observations of the moments of that mass distribution at a later time. The collision-coalescence kernel is assumed to be known, but one could also learn the parameters of the kernel instead of the parameters of the droplet distribution (or both).

Cloudy is used here in a "perfect model" (aka "known truth") setting, which means that the "observations" are generated by Cloudy itself, by running it with the true parameter values. In more realistic applications, this parameter estimation procedure will use actual measurements of cloud properties to obtain an estimated droplet mass distribution at a previous time.

Prerequisites

In order to run this example, you need to install Cloudy.jl (the "#master" lets you install the current master branch):

pkg > add Cloudy#master

Structure

The file Cloudy_example_eki.jl sets up the inverse problem and solves it using ensemble Kalman inversion, and the file Cloudy_example_uki.jl does the same using unscented Kalman inversion. The file DynamicalModel.jl provides the functionality to run the dynamical model $\Psi$, which in this example is Cloudy.

Running the Example

Once Cloudy is installed, the examples can be run from the julia REPL:

# Solve inverse problem using ensemble Kalman inversion
include("Cloudy_example_eki.jl")

# Solve inverse problem using unscented Kalman inversion
include("Cloudy_example_uki.jl")

What Does Cloudy Do?

The mathematical starting point of Cloudy is the stochastic collection equation (SCE; sometimes also called Smoluchowski equation after Marian Smoluchowski), which describes the time rate of change of $f = f(m, t)$, the mass distribution function of liquid water droplets, due to the process of collision and coalescence. The distribution function $f$ depends on droplet mass $m$ and time $t$ and is defined such that $f(m) \text{ d}m$ denotes the number of droplets with masses in the interval $[m, m + dm]$ per unit volume.

The stochastic collection equation is an integro-differential equation that can be written as

\[ \frac{\partial f(m, t)}{\partial t} = \frac{1}{2} \int_{m'=0}^{\infty} f(m', t) f(m-m', t) \mathcal{C}(m', m-m')\text{d}m' - f(m, t) \int_{m'=0}^\infty f(m', t)\mathcal{C}(m, m') \text{d}m', \]

where $\mathcal{C}(m', m'')$ is the collision-coalescence kernel, which encapsulates the physics of droplet-droplet interactions – it describes the rate at which two drops of masses $m'$ and $m''$ come into contact and coalesce into a drop of mass $m' + m''$. The first term on the right-hand side of the SCE describes the rate of increase of the number of drops having a mass $m$ due to collision and coalescence of drops of masses $m'$ and $m-m'$ (where the factor $\frac{1}{2}$ avoids double counting), while the second term describes the rate of reduction of drops of mass $m$ due to collision and coalescence of drops having a mass $m$ with other drops.

We can rewrite the SCE in terms of the moments $M_k$ of $f$, which are the prognostic variables in Cloudy. They are defined by

\[ M_k = \int_0^\infty m^k f(m, t) \text{d}m\]

The time rate of change of the k-th moment of $f$ is obtained by multiplying the SCE by $m^k$ and integrating over the entire range of droplet masses (from $m=0$ to $\infty$), which yields

\[ \frac{\partial M_k(t)}{\partial t} = \frac{1}{2}\int_0^\infty \left((m+m')^k - m^k - {m'}^k\right) \mathcal{C}(m, m')f(m, t)f(m', t) \, \text{d}m\, \text{d}m' ~~~~~~~~ (1)\]

In this example, the kernel is set to be constant – $\mathcal{C}(m', m'') = B = \text{const}$ – and the cloud droplet mass distribution is assumed to be a $\text{Gamma}(k_t, \theta_t)$ distribution, scaled by a factor $N_{0,t}$ which denotes the droplet number concentration:

\[f(m, t) = \frac{N_{0,t}}{\Gamma(k_t)\theta_t^k} m^{k_t-1} \exp{(-m/\theta_t)}\]

The parameter vector $\phi_t= [N_{0,t}, k_t, \theta_t]$ changes over time (as indicated by the subscript $t$), as the shape of the distribution evolves. In fact, there is a priori no reason to assume that the distribution would retain its Gamma shape over time, but this is a common assumption that is made in order to solve the closure problem (without this assumption, one would have to keep track of infinitely many moments of the mass distribution in order to uniquely identify the distribution $f$ at each time step, which is obviously not practicable).

For Gamma mass distribution functions, specifying the first three moments ($M_0$, $M_1$, and $M_2$) is sufficient to uniquely determine the parameter vector $\phi_t$, hence Cloudy solves equation (1) for $k = 0, 1, 2$. This mapping of the parameters of the initial cloud droplet mass distribution to the (zeroth-, first-, and second-order) moments of the distribution at a specified end time is done by DynamicalModel.jl.

Setting up the Inverse Problem

The goal is to learn the distribution parameters at time $t = 0$, $\phi_0 = [N_{0,0}, k_0, \theta_0]$, from observations $y = [M_0(t_{end}), M_1(t_{end}), M_2(t_{end})]$ of the zeroth-, first-, and second-order moments of the distribution at time $t_{end} > 0$ (where t_end = 1.0 in this example). This is a known truth experiment, in which the true parameters $\phi_{0, \text{true}}$ are defined to be:

N0_true = 300.0  # number of particles (scaling factor for Gamma distribution)
θ_true = 1.5597  # scale parameter of Gamma distribution
k_true = 0.0817  # shape parameter of Gamma distribution

Priors

In the code, the priors are constructed as follows:

par_names = ["N0", "θ", "k"]
# constrained_gaussian("name", desired_mean, desired_std, lower_bd, upper_bd)
prior_N0 = constrained_gaussian(par_names[1], 400, 300, 0.4 * N0_true, Inf)
prior_θ = constrained_gaussian(par_names[2], 1.0, 5.0, 1e-1, Inf)
prior_k = constrained_gaussian(par_names[3], 0.2, 1.0, 1e-4, Inf)
priors = combine_distributions([prior_N0, prior_θ, prior_k])

We use the recommended constrained_gaussian to add the desired scale and bounds to the prior distribution, in particular we place lower bounds to preserve positivity (and numerical stability).

Observational Noise

Cloudy produces output $y = [M_0(t_{end}), M_1(t_{end}), M_2(t_{end})]$, which is assumed to be related to the parameter vector $\theta$ according to:

\[ y = \mathcal{G}(\theta) + \eta,\]

where $\mathcal{G} = \Psi \circ \mathcal{T}^{-1}$ is the forward map, and the observational noise $\eta$ is assumed to be drawn from a 3-dimensional Gaussian with distribution $\mathcal{N}(0, \Gamma_y)$. In a perfect model setting, the observational noise represents the internal model variability. Since Cloudy is a purely deterministic model, there is no straightforward way of coming up with a covariance $\Gamma_y$ for this internal noise. We decide to use a diagonal covariance with the following entries (variances):

Γy = convert(Array, Diagonal([100.0, 5.0, 30.0]))

Artificial observations ("truth samples") are then generated by adding random samples from $\eta$ to G_t, the forward map evaluated for the true parameters:

for i in 1:n_samples
    y_t[:, i] = G_t .+ rand(MvNormal(μ, Γy))
end

truth = Observations.Observation(y_t, Γy, data_names)

Solution and Output

Cloudy_example_eki.jl: The optimal parameter vector determined by the ensemble Kalman inversion is the ensemble mean of the particles after the last iteration, which is printed to standard output. An output directory is created, where two files are stored: parameter_storage_eki.jld2 and data_storage_eki.jld2, which contain all parameters and model output from the ensemble Kalman iterations, respectively (both as DataContainers.DataContainer objects). In addition, an animation is produced that shows the evolution of the ensemble of particles over subsequent iterations of the optimization, both in the computational (unconstrained) and physical (constrained) spaces.
Cloudy_example_uki.jl: In addition to a point estimate of the optimal parameter (which is again given by the ensemble mean of the last iteration and printed to standard output), unscented Kalman inversion also provides a covariance approximation of the posterior distribution. Together, the mean and covariance allow for the reconstruction of a Gaussian approximation of the posterior distribution. The evolution of this Gaussian approximation over subsequent iterations is shown as an animation over the computational (unconstrained) space. All parameters as well as the model output from the unscented Kalman inversion are stored in an output directory, as parameter_storage_uki.jld2 and data_storage_uki.jld2.

Playing Around

If you want to play around with the Cloudy examples, you can e.g. change the type or the parameters of the initial cloud droplet mass distribution (see Cloudy.ParticleDistributions for the available distributions), by modifying these lines:

ϕ_true = [N0_true, θ_true, k_true]
dist_true = ParticleDistributions.GammaPrimitiveParticleDistribution(ϕ_true...)

(Don't forget to also change dist_type accordingly).

You can also experiment with different noise covariances (Γy), priors, vary the number of iterations (N_iter) or ensemble particles (N_ens), etc.