Model/Package Overview¶

This is an implementation of a simple, 0D adiabatic cloud parcel model tool (following Nenes et al, 2001 and Pruppacher and Klett, 1997). It allows flexible descriptions of an initial aerosol population, and simulates the evolution of a proto-cloud droplet population as the parcel ascends adiabatically at either a constant or time/height-dependent updraft speed. Droplet growth within the parcel is tracked on a Lagrangian grid.

You are invited to use the model (in accordance with the licensing), but please get in touch with the author via e-mail or on twitter. Up-to-date versions can be obtained through the model’s github repository or directly from the author. If you use the model for research, please cite this journal article which details the original model formulation:

Daniel Rothenberg and Chien Wang, 2016: Metamodeling of Droplet Activation for Global Climate Models. J. Atmos. Sci., 73, 1255–1272. doi: http://dx.doi.org/10.1175/JAS-D-15-0223.1

Documentation Outline¶

Scientific Description¶

The simplest tools available for describing the growth and evolution of a cloud droplet spectrum from a given population of aerosols are based on zero-dimensional, adiabatic cloud parcel models. By employing a detailed description of the condensation of ambient water vapor onto the growing droplets, these models can accurately describe the activation of a subset of the aerosol population by predicting how the presence of the aerosols in the updraft modify the maximum supersaturation achieved as the parcel rises. Furthermore, these models serve as the theoretical basis or reference for parameterizations of droplet activation which are included in modern general circulation models ([Ghan2011]) .

The complexity of these models varies with the range of physical processes one wishes to study. At the most complex end of the spectrum, one might wish to accurately resolve chemical transfer between the gas and aqueous phase in addition to physical transformations such as collision/coalescence. One could also add ice-phase processes to such a model.

Model Formulation¶

The adiabatic cloud parcel model implemented here is based on models described in the literature ([Nenes2001], [SP2006],) with some modifications and improvements. For a full description of the parcel model, please see ([Rothenberg2016]) The conservation of heat in a parcel of air rising at constant velocity \(V\) without entrainment can be written as

(1)\[\frac{dT}{dt} = -\frac{gV}{c_p} - \frac{L}{c_p}\frac{d w_v}{dt}\]

where \(T\) is the parcel air temperature. Assuming adiabaticity and neglecting entrainment is suitable for studying cloud droplet formation near the cloud base, where the majority of droplet activation occurs. Because the mass of water must be conserved as it changes from the vapor to liquid phase, the relationship

(2)\[\frac{d w_v}{dt} = - \frac{dw_c}{dt}\]

must hold, where \(w_v\) and \(w_c\) are the mass mixing ratios of water vapor and liquid water (condensed in droplets) in the parcel. The rate of change of water in the liquid phase in the parcel is governed solely by condensation onto the existing droplet population. For a population of \(N_i\) droplets of radius \(r_i\), where \(i=1,\dots,n\), the total condensation rate is given by

(3)\[\frac{dw_c}{dt} = \frac{4\pi \rho_w}{\rho_a}\sum\limits_{i=1}^nN_ir_i^2\frac{dr_i}{dt}\]

Here, the particle growth rate, \(\frac{dr_i}{dt}\) is calculated as

(4)\[\frac{dr_i}{dt} = \frac{G}{r_i}(S-S_{eq})\]

where \(G\) is a growth coefficient which is a function of the physical and chemical properties of the particle receiving condensate, given by

(5)\[G = \left(\frac{\rho_w R T}{e_s D'_v M_w} + \frac{L\rho_w[(LM_w/RT) - 1]}{k'_a T}\right)^{-1}\]

Droplet growth via condensation is modulated by the difference between the environmental supersaturation, \(S\), and the droplet equilibrium supersaturation, \(S_{eq}\), predicted from Kohler theory. To account for differences in aerosol chemical properties which could affect the ability for particles to uptake water, the \(\kappa\)-Köhler theory parameterization ([PK2007]) is employed in the model. \(\kappa\)-Kohler theory utilizes a single parameter to describe aerosol hygroscopicity, and is widely employed in modeling of aerosol processes. The hygroscopicity parameter \(\kappa\) is related to the water activity of an aqueous aerosol solution by

\[\frac{1}{a_w} = 1 + \kappa\frac{V_s}{V_w}\]

where \(V_s\) and \(V_w\) are the volumes of dy particulate matter and water in the aerosol solution. With this parameter, the full \(\kappa\)-Kohler theory may be expressed as

(6)\[S_{eq} = \frac{r_i^3 - r_{d,i}^3}{r_i^3 - r_{d,i}^3(1-\kappa_i)}\exp\left( \frac{2M_w\sigma_w}{RT\rho_w r_i} \right) - 1\]

where \(r_d\) and \(r\) are the dry aerosol particle size and the total radius of the wetted aerosol. The surface tension of water, \(\sigma_w\), is dependent on the temperature of the parcel such that \(\sigma_w = 0.0761 - 1.55\times 10^{-4}(T-273.15)\) J/m\(^2\) . Both the diffusivity and thermal conductivity of air have been modified in the growth coefficient equation to account for non-continuum effects as droplets grow, and are given by the expressions

\[D'_v = D_v\bigg/\left(1 + \frac{D_v}{a_c r}\sqrt{\frac{2\pi M_w}{RT}}\right)\]

and

\[k'_a = k_a\bigg/\left(1 + \frac{k_a}{a_T r \rho_a c_p}\sqrt{\frac{2\pi M_a}{RT}} \right)\]

In these expressions, the thermal accommodation coefficient, \(a_T\), is assumed to be \(0.96\) and the condensation coefficient, \(a_c\) is taken as unity (see Constants). Under the adiabatic assumption, the evolution of the parcel’s supersaturation is governed by the balance between condensational heating as water vapor condenses onto droplets and cooling induced by the parcel’s vertical motion,

(7)\[\frac{dS}{dt} = \alpha V - \gamma\frac{w_c}{dt}\]

where \(\alpha\) and \(\gamma\) are functions which are weakly dependent on temperature and pressure :

\[\alpha = \frac{gM_wL}{c_pRT^2} - \frac{gM_a}{RT}\]

\[\gamma = \frac{PM_a}{e_sM_w} + \frac{M_wL^2}{c_pRT^2}\]

The parcel’s pressure is predicted using the hydrostatic relationship, accounting for moisture by using virtual temperature (which can always be diagnosed as the model tracks the specific humidity through the mass mixing ratio of water vapor),

(8)\[\frac{dP}{dt} = \frac{-g P V}{R_d T_v}\]

The equations (8), (7), (3), (2), and (1) provide a simple, closed system of ordinary differential equations which can be numerically integrated forward in time. Furthermore, this model formulation allows the simulation of an arbitrary configuration of initial aerosols, in terms of size, number concentration, and hygroscopicity. Adding additional aerosol size bins is simply accomplished by tracking one additional size bin in the system of ODE’s. The important application of this feature is that the model can be configured to simulate both internal or external mixtures of aerosols, or some combination thereof.

Model Implementation and Procedure¶

The parcel model described in the previous section was implemented using a modern modular and object-oriented software engineering framework. This framework allows the model to be simply configured with myriad initial conditions and aerosol populations. It also enables model components - such as the numerical solver or condensation parameterization - to be swapped and replaced. Most importantly, the use of object-oriented techniques allows the model to be incorporated into frameworks which grossly accelerate the speed at which the model can be evaluated. For instance, although models like the one developed here are relatively cheap to execute, large ensembles of model runs have been limited in scope to several hundred or a thousand runs. However, the framework of this particular parcel model implementation was designed such that it could be run as a black box as part of a massively-parallel ensemble driver.

To run the model, a set of initial conditions needs to be specified, which includes the updraft speed, the parcel’s initial temperature, pressure, and supersaturation, and the aerosol population. Given these parameters, the model calculates an initial equilibrium droplet spectrum by computing the equilibrium wet radii of each aerosol. This is calculated numerically from the Kohler equation for each aerosol/proto-droplet, or numerically by employing the typical Kohler theory approximation

\[S \approx \frac{A}{r} - \kappa\frac{r_d^3}{r^3}\]

These wet radii are used as the initial droplet radii in the simulation.

Once the initial conditions have been configured, the model is integrated forward in time with a numerical solver (see parcel_model().run for more details). The available solvers wrapped here are:

LSODA(R)
LSODE
(C)VODE

During the model integration, the size representing each aerosol bin is allowed to grow via condensation, producing something akin to a moving grid. In the future, a fixed Eulerian grid will likely be implemented in the model for comparison.

Aerosol Population Specification¶

The model may be supplied with any arbitrary population of aerosols, providing the population can be approximated with a sectional representation. Most commonly, aerosol size distributions are represented with a continuous lognormal distribution,

(9)\[n_N(r) = \frac{dN}{d \ln r} = \frac{N_t}{\sqrt{2\pi}\ln \sigma_g}\exp\left(-\frac{ \ln^2(r/\mu_g)}{2\ln^2\sigma_g}\right)\]

which can be summarized with the set of three parameters, \((N_t, \mu_g, \sigma_g)\) and correspond, respectively, to the total aerosol number concentration, the geometric mean or number mode radius, and the geometric standard deviation. Complicated multi-modal aerosol distributions can often be represented as the sum of several lognormal distributions. Since the parcel model describes the evolution of a discrete aerosol size spectrum, can be broken into \(n\) bins, and the continuous aerosol size distribution approximated by taking the number concentration and size at the geometric mean value in each bin, such that the discrete approximation to the aerosol size distribution becomes

\[n_{N,i}(r_i) = \sum\limits_{i=1}^n\frac{N_i}{\sqrt{2\pi}\ln\sigma_g}\exp\left(-\frac{\ln^2(r_i/\mu_g)}{2\ln^2\sigma_g}\right)\]

If no bounds on the size range of \(r_i\) is specified, then the model pre-computes \(n\) equally-spaced bins over the logarithm of \(r\), and covers the size range \(\mu_g/10\sigma_g\) to \(10\sigma_g\mu_g\). It is typical to run the model with \(200\) size bins per aerosol mode. Neither this model nor similar ones exhibit much sensitivity towards the density of the sectional discretization .

Typically, a single value for hygroscopicity, \(\kappa\) is prescribed for each aerosol mode. However, the model tracks a hygroscopicity parameter for each individual size bin, so size-dependent aerosol composition can be incorporated into the aerosol population. This representation of the aerosol population is similar to the external mixing state assumption. An advantage to using this representation is that complex mixing states can be represented by adding various size bins, each with their own number concentration and hygroscopicity.

References¶

[Nenes2001]

Nenes, A., Ghan, S., Abdul-Razzak, H., Chuang, P. Y. & Seinfeld, J. H. Kinetic limitations on cloud droplet formation and impact on cloud albedo. Tellus 53, 133–149 (2001).

[SP2006]

Seinfeld, J. H. & Pandis, S. N. Atmospheric Chemistry and Physics: From Air Pollution to Climate Change. Atmos. Chem. Phys. 2nd, 1203 (Wiley, 2006).

[Rothenberg2016]

Daniel Rothenberg and Chien Wang, 2016: Metamodeling of Droplet Activation for Global Climate Models. J. Atmos. Sci., 73, 1255–1272. doi: http://dx.doi.org/10.1175/JAS-D-15-0223.1

[PK2007]

Petters, M. D. & Kreidenweis, S. M. A single parameter representation of hygroscopic growth and cloud condensation nucleus activity. Atmos. Chem. Phys. 7, 1961–1971 (2007).

[Ghan2011]

Ghan, S. J. et al. Droplet nucleation: Physically-based parameterizations and comparative evaluation. J. Adv. Model. Earth Syst. 3, M10001 (2011).

Installation¶

To grab and build the latest version of the model, you should use pip and point it to the source code repository on github:

$ pip install git+git://github.com/darothen/parcel_model.git

This should automatically build the necessary Cython modules and export the code package to your normal package installation directory. If you wish to simply build the code and run it in place, clone the repository, navigate to it in a terminal, and invoke the build command by hand:

$ python setup.py build_ext --inplace

This should produce the compiled file parcel_aux.so in the model package. You can also install the code from the cloned source directory by invoking pip install from within it; this is useful if you’re updating or modifying the model, since you can install an “editable” package which points directly to the git-monitored code:

$ cd path/to/parcel_model/
$ pip install -e .

Dependencies¶

This code was originally written for Python 2.7, and then futurized to Python 3.3+ with hooks for backwards compatibility. By far, the simplest way to run this code is to grab a scientific python distribution, such as Anaconda. This code should work out-of-the box with almost all dependencies filled (exception being numerical solvers) on a recent version (1.2+) of this distribution. To faciliate this, conda environments for Python versions 2.7 and 3.4+ are provided in the parcel_model/ci directory.

Necessary dependencies¶

Note

As of version 1.2.0, the model integration components are being re-written and only the CVODE interface is exposed. As such, Assimulo is temporarily a core and required dependency; in the future the other solvers will be re-enabled. For best results, you will want to manually install Assimulo, as I’ve encountered issues using the available pip or conda packages.

Numerical solver dependencies¶

LSODA - scipy or odespy
VODE, LSODE - odespy
CVODE - Assimulo

Recommended additional packages¶

Testing¶

A nose test-suite is under construction. To check that your model is configured and running correctly, you copy and run the notebook corresponding to the basic run example, or run the command-line interface version of the model with the pre-packed simple run case:

$ cd path/to/parcel_model/
$ ./run_parcel examples/simple.yml

Bugs / Suggestions¶

The code has an issue tracker on github and I strongly encourage you to note any problems with the model there, such as typos or weird behavior and results. Furthermore, I’m looking for ways to expand and extend the model, so if there is something you might wish to see added, please note it there or send me an e-mail. The code was written in such a way that it should be trivial to add physics in a modular fashion.

Example: Activation¶

In this example, we will study the effect of updraft speed on the activation of a lognormal ammonium sulfate accumulation mode aerosol.

# Suppress warnings
import warnings
warnings.simplefilter('ignore')

import parcel_model as pm
import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

First, we indicate the parcel’s initial thermodynamic conditions.

P0 = 100000. # Pressure, Pa
T0 = 279.    # Temperature, K
S0 = -0.1   # Supersaturation, 1-RH

We next define the aerosol distribution to follow the reference simulation from Ghan et al, 2011

aer =  pm.AerosolSpecies('ammonium sulfate',
                          pm.Lognorm(mu=0.05, sigma=2.0, N=1000.),
                          kappa=0.7, bins=100)

Loop over updraft several velocities in the range 0.1 - 10.0 m/s. We will peform a detailed parcel model calculation, as well as calculations with two activation parameterizations. We will also use an accommodation coefficient of \(\alpha_c = 0.1\), following the recommendations of Raatikainen et al (2013).

First, the parcel model calculations:

from parcel_model import binned_activation

Vs = np.logspace(-1, np.log10(10,), 11.)[::-1] # 0.1 - 5.0 m/s
accom = 0.1

smaxes, act_fracs = [], []
for V in Vs:
    # Initialize the model
    model = pm.ParcelModel([aer,], V, T0, S0, P0, accom=accom, console=False)
    par_out, aer_out = model.run(t_end=2500., dt=1.0, solver='cvode',
                                 output='dataframes', terminate=True)
    print(V, par_out.S.max())

    # Extract the supersaturation/activation details from the model
    # output
    S_max = par_out['S'].max()
    time_at_Smax = par_out['S'].argmax()
    wet_sizes_at_Smax = aer_out['ammonium sulfate'].ix[time_at_Smax].iloc[0]
    wet_sizes_at_Smax = np.array(wet_sizes_at_Smax.tolist())

    frac_eq, _, _, _ = binned_activation(S_max, T0, wet_sizes_at_Smax, aer)

    # Save the output
    smaxes.append(S_max)
    act_fracs.append(frac_eq)

[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
10.0 0.0156189147154
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
6.3095734448 0.0116683910368
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
3.98107170553 0.00878287310116
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
2.51188643151 0.00664901290831
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
1.58489319246 0.00505644091867
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
1.0 0.00385393398982
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
0.63095734448 0.00293957320198
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
0.398107170553 0.00224028774582
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
0.251188643151 0.00170480101361
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
0.158489319246 0.0012955732509
[CVode Warning] b'At the end of the first step, there are still some root functions identically 0. This warning will not be issued again.'
0.1 0.000984803827635

Now the activation parameterizations:

smaxes_arg, act_fracs_arg = [], []
smaxes_mbn, act_fracs_mbn = [], []

for V in Vs:
    smax_arg, _, afs_arg = pm.arg2000(V, T0, P0, [aer], accom=accom)
    smax_mbn, _, afs_mbn = pm.mbn2014(V, T0, P0, [aer], accom=accom)

    smaxes_arg.append(smax_arg)
    act_fracs_arg.append(afs_arg[0])
    smaxes_mbn.append(smax_mbn)
    act_fracs_mbn.append(afs_mbn[0])

Finally, we compile our results into a nice plot for visualization.

sns.set(context="notebook", style='ticks')
sns.set_palette("husl", 3)
fig, [ax_s, ax_a] = plt.subplots(1, 2, sharex=True, figsize=(10,4))

ax_s.plot(Vs, np.array(smaxes)*100., color='k', lw=2, label="Parcel Model")
ax_s.plot(Vs, np.array(smaxes_mbn)*100., linestyle='None',
          marker="o", ms=10, label="MBN2014" )
ax_s.plot(Vs, np.array(smaxes_arg)*100., linestyle='None',
          marker="o", ms=10, label="ARG2000" )
ax_s.semilogx()
ax_s.set_ylabel("Superaturation Max, %")
ax_s.set_ylim(0, 2.)

ax_a.plot(Vs, act_fracs, color='k', lw=2, label="Parcel Model")
ax_a.plot(Vs, act_fracs_mbn, linestyle='None',
          marker="o", ms=10, label="MBN2014" )
ax_a.plot(Vs, act_fracs_arg, linestyle='None',
          marker="o", ms=10, label="ARG2000" )
ax_a.semilogx()
ax_a.set_ylabel("Activated Fraction")
ax_a.set_ylim(0, 1.)

plt.tight_layout()
sns.despine()

for ax in [ax_s, ax_a]:
    ax.legend(loc='upper left')
    ax.xaxis.set_ticks([0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0])
    ax.xaxis.set_ticklabels([0.1, 0.2, 0.5, 1.0, 2.0, 5.0, 10.0])
    ax.set_xlabel("Updraft speed, m/s")

Example: Basic Run¶

In this example, we will setup a simple parcel model simulation containing two aerosol modes. We will then run the model with a 1 m/s updraft, and observe how the aerosol population bifurcates into swelled aerosol and cloud droplets.

# Suppress warnings
import warnings
warnings.simplefilter('ignore')

import parcel_model as pm
import numpy as np

%matplotlib inline
import matplotlib.pyplot as plt

Could not find GLIMDA

First, we indicate the parcel’s initial thermodynamic conditions.

P0 = 77500. # Pressure, Pa
T0 = 274.   # Temperature, K
S0 = -0.02  # Supersaturation, 1-RH (98% here)

Next, we define the aerosols present in the parcel. The model itself is agnostic to how the aerosol are specified; it simply expects lists of the radii of wetted aerosol radii, their number concentration, and their hygroscopicity. We can make container objects (:class:AerosolSpecies) that wrap all of this information so that we never need to worry about it.

Here, let’s construct two aerosol modes:

Mode	\(\kappa\) (hygroscopicity)	Mean size (micron)	Std dev	Number Conc (cm**-3)
sulfate	0.54	0.015	1.6	850
sea salt	1.2	0.85	1.2	10

We’ll define each mode using the :class:Lognorm distribution packaged with the model.

sulfate =  pm.AerosolSpecies('sulfate',
                             pm.Lognorm(mu=0.015, sigma=1.6, N=850.),
                             kappa=0.54, bins=200)
sea_salt = pm.AerosolSpecies('sea salt',
                             pm.Lognorm(mu=0.85, sigma=1.2, N=10.),
                             kappa=1.2, bins=40)

The :class:AerosolSpecies class automatically computes gridded/binned representations of the size distributions. Let’s double check that the aerosol distribution in the model will make sense by plotting the number concentration in each bin.

fig = plt.figure(figsize=(10,5))
ax = fig.add_subplot(111)
ax.grid(False, "minor")

sul_c = "#CC0066"
ax.bar(sulfate.rs[:-1], sulfate.Nis*1e-6, np.diff(sulfate.rs),
        color=sul_c, label="sulfate", edgecolor="#CC0066")
sea_c = "#0099FF"
ax.bar(sea_salt.rs[:-1], sea_salt.Nis*1e-6, np.diff(sea_salt.rs),
        color=sea_c, label="sea salt", edgecolor="#0099FF")
ax.semilogx()

ax.set_xlabel("Aerosol dry radius, micron")
ax.set_ylabel("Aerosl number conc., cm$^{-3}$")
ax.legend(loc='upper right')

<matplotlib.legend.Legend at 0x10f4baeb8>

Actually running the model is very straightforward, and involves just two steps:

Instantiate the model by creating a :class:ParcelModel object.
Call the model’s :method:run method.

For convenience this process is encoded into several routines in the driver file, including both a single-strategy routine and an iterating routine which adjusts the the timestep and numerical tolerances if the model crashes. However, we can illustrate the simple model running process here in case you wish to develop your own scheme for running the model.

initial_aerosols = [sulfate, sea_salt]
V = 1.0 # updraft speed, m/s

dt = 1.0 # timestep, seconds
t_end = 250./V # end time, seconds... 250 meter simulation

model = pm.ParcelModel(initial_aerosols, V, T0, S0, P0, console=False, accom=0.3)
parcel_trace, aerosol_traces = model.run(t_end, dt, solver='cvode')

If console is set to True, then some basic debugging output will be written to the terminal, including the initial equilibrium droplet size distribution and some numerical solver diagnostics. The model output can be customized; by default, we get a DataFrame and a Panel of the parcel state vector and aerosol bin sizes as a function of time (and height). We can use this to visualize the simulation results, like in the package’s README.

fig, [axS, axA] = plt.subplots(1, 2, figsize=(10, 4), sharey=True)

axS.plot(parcel_trace['S']*100., parcel_trace['z'], color='k', lw=2)
axT = axS.twiny()
axT.plot(parcel_trace['T'], parcel_trace['z'], color='r', lw=1.5)

Smax = parcel_trace['S'].max()*100
z_at_smax = parcel_trace['z'].ix[parcel_trace['S'].argmax()]
axS.annotate("max S = %0.2f%%" % Smax,
             xy=(Smax, z_at_smax),
             xytext=(Smax-0.3, z_at_smax+50.),
             arrowprops=dict(arrowstyle="->", color='k',
                             connectionstyle='angle3,angleA=0,angleB=90'),
             zorder=10)

axS.set_xlim(0, 0.7)
axS.set_ylim(0, 250)

axT.set_xticks([270, 271, 272, 273, 274])
axT.xaxis.label.set_color('red')
axT.tick_params(axis='x', colors='red')

axS.set_xlabel("Supersaturation, %")
axT.set_xlabel("Temperature, K")
axS.set_ylabel("Height, m")

sulf_array = aerosol_traces['sulfate'].values
sea_array = aerosol_traces['sea salt'].values

ss = axA.plot(sulf_array[:, ::10]*1e6, parcel_trace['z'], color=sul_c,
         label="sulfate")
sa = axA.plot(sea_array*1e6, parcel_trace['z'], color=sea_c, label="sea salt")
axA.semilogx()
axA.set_xlim(1e-2, 10.)
axA.set_xticks([1e-2, 1e-1, 1e0, 1e1], [0.01, 0.1, 1.0, 10.0])
axA.legend([ss[0], sa[0]], ['sulfate', 'sea salt'], loc='upper right')
axA.set_xlabel("Droplet radius, micron")

for ax in [axS, axA, axT]:
    ax.grid(False, 'both', 'both')

In this simple example, the sulfate aerosol population bifurcated into interstitial aerosol and cloud droplets, while the entire sea salt population activated. A peak supersaturation of about 0.63% was reached a few meters above cloud base, where the ambient relative humidity hit 100%.

How many CDNC does this translate into? We can call upon helper methods from the activation package to perform these calculations for us:

from parcel_model import binned_activation

sulf_trace = aerosol_traces['sulfate']
sea_trace = aerosol_traces['sea salt']

ind_final = int(t_end/dt) - 1

T = parcel_trace['T'].iloc[ind_final]
eq_sulf, kn_sulf, alpha_sulf, phi_sulf = \
    binned_activation(Smax/100, T, sulf_trace.iloc[ind_final],  sulfate)
eq_sulf *= sulfate.total_N

eq_sea, kn_sea, alpha_sea, phi_sea = \
    binned_activation(Smax/100, T, sea_trace.iloc[ind_final], sea_salt)
eq_sea *= sea_salt.total_N

print("  CDNC(sulfate) = {:3.1f}".format(eq_sulf))
print(" CDNC(sea salt) = {:3.1f}".format(eq_sea))
print("------------------------")
print("          total = {:3.1f} / {:3.0f} ~ act frac = {:1.2f}".format(
      eq_sulf+eq_sea,
      sea_salt.total_N+sulfate.total_N,
      (eq_sulf+eq_sea)/(sea_salt.total_N+sulfate.total_N)
))

  CDNC(sulfate) = 146.9
 CDNC(sea salt) = 10.0
------------------------
          total = 156.9 / 860 ~ act frac = 0.18

Parcel Model Details¶

Below is the documentation for the parcel model, which is useful for debugging and development. For a higher-level overview, see the scientific description.

Implementation¶

class parcel_model.ParcelModel(aerosols, V, T0, S0, P0, console=False, accom=1.0, truncate_aerosols=False)¶

Wrapper class for instantiating and running the parcel model.

The parcel model has been implemented in an object-oriented format to facilitate easy extensibility to different aerosol and meteorological conditions. A typical use case would involve specifying the initial conditions such as:

>>> import parcel_model as pm
>>> P0 = 80000.
>>> T0 = 283.15
>>> S0 = 0.0
>>> V = 1.0
>>> aerosol1 = pm.AerosolSpecies('sulfate',
...                              Lognorm(mu=0.025, sigma=1.3, N=2000.),
...                              bins=200, kappa=0.54)
>>> initial_aerosols = [aerosol1, ]
>>> z_top = 50.
>>> dt = 0.01

which initializes the model with typical conditions at the top of the boundary layer (800 hPa, 283.15 K, 100% Relative Humidity, 1 m/s updraft), and a simple sulfate aerosol distribution which will be discretized into 200 size bins to track. Furthermore the model was specified to simulate the updraft for 50 meters (z_top) and use a time-discretization of 0.01 seconds. This timestep is used in the model output – the actual ODE solver will generally calculate the trace of the model at many more times.

Running the model and saving the output can be accomplished by invoking:

>>> model = pm.ParcelModel(initial_aerosols, V, T0, S0, P0)
>>> par_out, aer_out = pm.run(z_top, dt)

This will yield par_out, a pandas.DataFrame containing the meteorological conditions in the parcel, and aerosols, a dictionary of DataFrame objects for each species in initial_aerosols with the appropriately tracked size bins and their evolution over time.

Derivative Equation¶

parcel.der(y, t, nr, r_drys, Nis, V, kappas, accom=1.0)¶

Calculates the instantaneous time-derivate of the parcel model system.

Given a current state vector y of the parcel model, computes the tendency of each term including thermodynamic (pressure, temperature, etc) and aerosol terms. The basic aerosol properties used in the model must be passed along with the state vector (i.e. if being used as the callback function in an ODE solver).

This function is implemented in NumPy and Python, and is likely very slow compared to the available Cython version.

Parameters:

y : array_like

Current state of the parcel model system,

y[0] = altitude, m

y[1] = Pressure, Pa

y[2] = temperature, K

y[3] = water vapor mass mixing ratio, kg/kg

y[4] = cloud liquid water mass mixing ratio, kg/kg

y[5] = cloud ice water mass mixing ratio, kg/kg

y[6] = parcel supersaturation

y[7:] = aerosol bin sizes (radii), m

t : float

Current simulation time, in seconds.

nr : Integer

Number of aerosol radii being tracked.

r_drys : array_like

Array recording original aerosol dry radii, m.

Nis : array_like

Array recording aerosol number concentrations, 1/(m**3).

V : float

Updraft velocity, m/s.

kappas : array_like

Array recording aerosol hygroscopicities.

accom : float, optional (default=:const:constants.ac)

Condensation coefficient.

Returns:

x : array_like

Array of shape (``nr``+7, ) containing the evaluated parcel model instaneous derivative.

Notes

This Python sketch of the derivative function shouldn’t really be used for any computational purposes. Instead, see the cythonized version in the auxiliary file, parcel_aux.pyx. In the default configuration, once the code has been built, you can set the environmental variable OMP_NUM_THREADS to control the parallel for loop which calculates the condensational growth rate for each bin.

Reference¶

Main Parcel Model¶

The core of the model has its own documentation page, which you can access here.

ParcelModel(aerosols, V, T0, S0, P0[, ...]) Wrapper class for instantiating and running the parcel model.

Driver Tools¶

Utilities for driving sets of parcel model integration strategies.

Occasionally, a pathological set of input parameters to the parcel model will really muck up the ODE solver’s ability to integrate the model. In that case, it would be nice to quietly adjust some of the numerical parameters for the ODE solver and re-submit the job. This module includes a workhorse function iterate_runs() which can serve this purpose and can serve as an example for more complex integration strategies. Alternatively, :func:`run_model`is a useful shortcut for building/running a model and snagging its output.

`run_model`(V, initial_aerosols, T, P, dt[, ...])	Setup and run the parcel model with given solver configuration.
`iterate_runs`(V, initial_aerosols, T, P[, ...])	Iterate through several different strategies for integrating the parcel model.

Thermodynamics/Kohler Theory¶

Aerosol/atmospheric thermodynamics functions.

The following sets of functions calculate useful thermodynamic quantities that arise in aerosol-cloud studies. Where possible, the source of the parameterization for each function is documented.

`dv`(T, r, P[, accom])	Diffusivity of water vapor in air, modified for non-continuum effects.
`ka`(T, rho, r)	Thermal conductivity of air, modified for non-continuum effects.
`rho_air`(T, P[, RH])	Density of moist air with a given relative humidity, temperature, and pressure.
`es`(T_c)	Calculates the saturation vapor pressure over water for a given temperature.
`sigma_w`(T)	Surface tension of water for a given temperature.
`Seq`(r, r_dry, T, kappa[, neg, approx])	kappa-Kohler theory equilibrium saturation over aerosol.
`kohler_crit`(T, r_dry, kappa[, approx])	Critical radius and supersaturation of an aerosol particle.
`critical_curve`(T, r_a, r_b, kappa[, approx])	Calculates curves of critical radii and supersaturations for aerosol.

Aerosols¶

Container class for encapsulating data about aerosol size distributions.

AerosolSpecies(species, distribution, kappa) Container class for organizing aerosol metadata.

The following are utility functions which might be useful in studying and manipulating aerosol distributions for use in the model or activation routines.

dist_to_conc(dist, r_min, r_max[, rule]) Converts a swath of a size distribution function to an actual number concentration.

Distributions¶

Collection of classes for representing aerosol size distributions.

Most commonly, one would use the Lognorm distribution. However, for the sake of completeness, other canonical distributions will be included here, with the notion that this package could be extended to describe droplet size distributions or other collections of objects.

`BaseDistribution`	Interface for distributions, to ensure that they contain a pdf method.
`Gamma`	Gamma size distribution
`Lognorm`(mu, sigma[, N, base])	Lognormal size distribution.
`MultiModeLognorm`(mus, sigmas, Ns[, base])	Multimode lognormal distribution class.

The following dictionaries containing (multi) Lognormal aerosol size distributions have also been saved for convenience:

FN2005_single_modes: Fountoukis, C., and A. Nenes (2005), Continued development of a cloud droplet formation parameterization for global climate models, J. Geophys. Res., 110, D11212, doi:10.1029/2004JD005591
NS2003_single_modes: Nenes, A., and J. H. Seinfeld (2003), Parameterization of cloud droplet formation in global climate models, J. Geophys. Res., 108, 4415, doi:10.1029/2002JD002911, D14.
whitby_distributions: Whitby, K. T. (1978), The physical characteristics of sulfur aerosols, Atmos. Environ., 12(1-3), 135–159, doi:10.1016/0004-6981(78)90196-8.
jaenicke_distributions: Jaenicke, R. (1993), Tropospheric Aerosols, in Aerosol-Cloud-Climate Interactions, P. V. Hobbs, ed., Academic Press, San Diego, CA, pp. 1-31.

Activation¶

Collection of activation parameterizations.

`lognormal_activation`(smax, mu, sigma, N, kappa)	Compute the activated number/fraction from a lognormal mode
`binned_activation`(Smax, T, rs, aerosol[, approx])	Compute the activation statistics of a given aerosol, its transient size distribution, and updraft characteristics.
`multi_mode_activation`(Smax, T, aerosols, rss)	Compute the activation statistics of a multi-mode, binned_activation aerosol population.
`arg2000`(V, T, P[, aerosols, accom, mus, ...])	Computes droplet activation using a psuedo-analytical scheme.
`mbn2014`(V, T, P[, aerosols, accom, mus, ...])	Computes droplet activation using an iterative scheme.
`shipwayabel2010`(V, T, P, aerosol)	Activation scheme following Shipway and Abel, 2010 (doi:10.1016/j.atmosres.2009.10.005).
`ming2006`(V, T, P, aerosol)	Ming activation scheme.

Constants¶

Commonly used constants in microphysics and aerosol thermodynamics equations as well as important model parameters.

Symbol	Variable	Value	Units	Description
\(g\)	`g`	9.8	m s**-2	gravitational constant
\(C_p\)	`Cp`	1004.0	J/kg	specific heat of dry air at constant pressure
\(\rho_w\)	`rho_w`	1000.0	kg m**-3	density of water at STP
\(R_d\)	`Rd`	287.0	J/kg/K	gas constant for dry air
\(R_v\)	`Rv`	461.5	J/kg/K	gas constant for water vapor
\(R\)	`R`	8.314	J/mol/K	universal gas constant
\(M_w\)	`Mw`	0.018	kg/mol	molecular weight of water
\(M_a\)	`Ma`	0.0289	kg/mol	molecular weight of dry air
\(D_v\)	`Dv`	3e-5	m**2/s	diffusivity of water vapor in air
\(L_v\)	`L`	2.25e6	J/kg/K	latent heat of vaporization of water
\(\alpha_c\)	`ac`	1.0	unitless	condensation coefficient
\(K_a\)	`Ka`	0.02	J/m/s/K	thermal conductivity of air
\(a_T\)	`at`	0.96	unitless	thermal accommodation coefficient
\(\epsilon\)	`epsilon`	0.622	unitless	ratio of \(M_w/M_a\)

Additionally, a reference table containing the 1976 US Standard Atmosphere is implemented in the constant std_atm, which is a pandas DataFrame with the fields

alt, altitude in km
sigma, ratio of density to sea-level density
delta, ratio of pressure to sea-level pressure
theta, ratio of temperature to sea-level temperature
temp, temperature in K
press, pressure in Pa
dens, air density in kg/m**3
k.visc, air kinematic viscosity
ratio, ratio of speed of sound to kinematic viscosity in m**-1

Using default pandas functons, you can interpolate to any reference pressure or height level.

Current version: 1.2.0

Documentation last compiled: Oct 26, 2016

V, T0, S0, P0, aerosols	(floats) Initial parcel settings (see Parameters).
_r0s	(array_like of floats) Initial equilibrium droplet sizes.
_r_drys	(array_like of floats) Dry radii of aerosol population.
_kappas	(array_like of floats) Hygroscopicity of each aerosol size.
_Nis	(array_like of floats) Number concentration of each aerosol size.
_nr	(int) Number of aerosol sizes tracked in model.
_model_set	(boolean) Flag indicating whether or not at any given time the model initialization/equilibration routine has been run with the current model settings.
_y0	(array_like) Initial state vector.

run(t_end, dt, max_steps=1000, solver=”odeint”, output_fmt=”dataframes”, terminate=False, solver_args={})	Execute model simulation.
set_initial_conditions(V=None, T0=None, S0=None, P0=None, aerosols=None)	Re-initialize a model simulation in order to run it.