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.


class pyrcel.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 pyrcel 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 =, 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.

See also

companion routine which computes equilibrium droplet sizes and sets the model’s state vectors.


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.
run(t_end, output_dt=1.0, solver_dt=None, max_steps=1000, solver='odeint', output_fmt='dataframes', terminate=False, terminate_depth=100.0, **solver_args)

Run the parcel model simulation.

Once the model has been instantiated, a simulation can immediately be performed by invoking this method. The numerical details underlying the simulation and the times over which to integrate can be flexibly set here.

Time – The user must specify two timesteps: output_dt, which is the timestep between output snapshots of the state of the parcel model, and solver_dt, which is the the interval of time before the ODE integrator is paused and re-started. It’s usually okay to use a very large solver_dt, as output_dt can be interpolated from the simulation. In some cases though a small solver_dt could be useful to force the solver to use smaller internal timesteps.

Numerical Solver – By default, the model will use the odeint wrapper of LSODA shipped by default with scipy. Some fine-tuning of the solver tolerances is afforded here through the max_steps. For other solvers, a set of optional arguments solver_args can be passed.

Solution Output – Several different output formats are available by default. Additionally, the output arrays are saved with the ParcelModel instance so they can be used later.


t_end : float

Total time over interval over which the model should be integrated

output_dt : float

Timestep intervals to report model output.

solver_dt : float

Timestep interval for calling solver integration routine.

max_steps : int

Maximum number of steps allowed by solver to satisfy error tolerances per timestep.

solver : {‘odeint’, ‘lsoda’, ‘lsode’, ‘vode’, cvode’}

Choose which numerical solver to use: * ‘odeint’: LSODA implementation from ODEPACK via

SciPy’s integrate module

  • ‘lsoda’: LSODA implementation from ODEPACK via odespy
  • ‘lsode’: LSODE implementation from ODEPACK via odespy
  • ‘vode’ : VODE implementation from ODEPACK via odespy
  • ‘cvode’ : CVODE implementation from Sundials via Assimulo
  • ‘lsodar’ : LSODAR implementation from Sundials via Assimulo

output_fmt : str, one of {‘dataframes’, ‘arrays’, ‘smax’}

Choose format of solution output.

terminate : boolean

End simulation at or shortly after a maximum supersaturation has been achieved

terminate_depth : float, optional (default=100.)

Additional depth (in meters) to integrate after termination criterion eached.


DataFrames, array, or float

Depending on what was passed to the output argument, different types of data might be returned:

  • `dataframes’: (default) will process the output into
    two pandas DataFrames - the first one containing profiles of the meteorological quantities tracked in the model, and the second a dictionary of DataFrames with one for each AerosolSpecies, tracking the growth in each bin for those species.
  • ‘arrays’: will return the raw output from the solver
    used internally by the parcel model - the state vector y and the evaluated timesteps converted into height coordinates.
  • ‘smax’: will only return the maximum supersaturation
    value achieved in the simulation.


The parcel model failed to complete successfully or failed to initialize.

See also

right-hand side derivative evaluated during model integration.
set_initial_conditions(V=None, T0=None, S0=None, P0=None, aerosols=None)

Set the initial conditions and parameters for a new parcel model run without having to create a new ParcelModel instance.

Based on the aerosol population which has been stored in the model, this method will finish initializing the model. This has three major parts:

  1. concatenate the aerosol population information (their dry radii, hygroscopicities, etc) into single arrays which can be placed into the state vector for forward integration.
  2. Given the initial ambient water vapor concentration (computed from the temperature, pressure, and supersaturation), determine how much water must already be coated on the aerosol particles in order for their size to be in equilibrium.
  3. Set-up the state vector with these initial conditions.

Once the state vector has been set up, the setup routine will record attributes in the parent instance of the ParcelModel.


V, T0, S0, P0 : floats

The updraft speed and initial temperature (K), pressure (Pa), supersaturation (percent, with 0.0 = 100% RH).

aerosols : array_like sequence of AerosolSpecies

The aerosols contained in the parcel.



If an equilibrium droplet size distribution could not be calculated.


The actual setup occurs in the private method _setup_run(); this method is simply an interface that can be used to modify an existing ParcelModel.

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.


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 (

Condensation coefficient.


x : array_like

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


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.