Building the Cloud Parcel Model on Matlab - Part 2: Create the core ODE system (A. basic)

Here’s how you can construct and describe the ODE system of a cloud parcel model—based on Nenes et al. (2001, 2002) and examples like Pyrcel—for full implementation in MATLAB. This clarifies the variables, equations, and overall structure you'll include in your solver (ode15s).

🌫️ 1. State Vector (`y`)

Your ODE state vector y should stack the following:

Temperature: $T(t)$ — air parcel temperature [K]
Water vapor: $w(t)$ — mass mixing ratio [kg/kg]
Droplet radii: $r_i(t)$ , for $i = 1,\dots,N$ — each aerosol bin radius [m]

So in MATLAB:


y = [ T; w; r1; r2; …; rN ];

🔢 2. ODEs Formulation

You’ll define the right-hand side function parcel_ode(t, y, params), which returns:

$\frac{d y}{d t} = [d T / d t; d w / d t; d r_{1} / d t; \dots; d r_{N} / d t]$

2.1 Droplet Growth (`dr_i/dt`)

Using κ-Köhler theory (as in Nenes 2001, Pyrcel implementation

\frac{d r_{i}}{d t} = \frac{G (T, p, \dots) S}{r_{i​}}

$G$ : depends on diffusivity, latent heat, etc.
$S = \frac{e}{e_s} - 1$ : supersaturation from temperature and vapor.

2.2 Vapor Change (`dw/dt`)

Mass conservation: condensation on all droplets reduces vapor:

\frac{d w}{d t} = - \frac{1}{ρ_{air}} \sum_{i} 4 π r_{i}^{2} N_{i} \frac{d r_{i}}{d t​}

Where $N_i$ is the aerosol number concentration.

2.3 Temperature Change (`dT/dt`)

Includes:

Adiabatic cooling from ascent.
Latent heat from condensation: $-\frac{L_v}{c_p} dw/dt$
Potential heating from BC scenarios: integrated $Q_{\text{ext}}$ and $g_a$ terms.

Example:

dTdt = -g*V/cp + Qext_term - (Lv/cp)*(dw_dt);

Depending on your scenario, Qext_term toggles on or off.

🔄 3. Putting It All Together in MATLAB

matlab

function dydt = parcel_ode(t, y, p)
    T = y(1); w = y(2); ri = y(3:end);
    S = compute_supersat(T, w, p);
    dri = droplet_growth_rate(ri, S, T, p);
    dw = -sum(4*pi*ri.^2 .* p.Nd .* dri) / p.rho_air;
    dT = -p.g * p.V / p.Cp - p.Lv/p.Cp * dw + heating_term(ri, p);
    dydt = [dT; dw; dri];
end

You'd define heating_term(ri, p) differently depending on No Heat, Droplet, External, or Internal heating.

🗺️ 4. Why Is This a "System"?

Because the derivatives are coupled:

$dr_i/dt$ depends on $T$ and $w$ via supersaturation.
$dw/dt$ depends on all $dr_i/dt$ via water mass uptake.
$dT/dt$ depends on $dw/dt$ (latent heating) and sometimes $ri$ (BC heating).

You must solve these equations together, which is why you pass them as a single function to ode15s.

✅ Summary

The ODE system includes variables: T, w, {r₁...rₙ}
Derivative functions respond to each other (coupling)
Implemented as one MATLAB function returning dydt
Solved together with ode15s, a solver adept at stiff, coupled systems

LOOK CLOSER...

"Where is the function definition of each component in the core ODE system? Why are some directly defined, some not?"

Let's break this down logically and transparently.

🔹 The Core ODE System

The heart of the cloud parcel model is an ODE system that describes how three things evolve with time:

Temperature (T)
Water vapor mixing ratio (w)
Droplet radius (ri) for each size bin

This is often written in the form:

\frac{d}{dt} \begin{bmatrix} T \\ w \\ r_1 \\ r_2 \\ \dots \\ r_N \end{bmatrix} = \text{f}(t, y, \text{params})

In MATLAB, this is what ode15s integrates numerically.

🔹 Why Some Components Are Defined Inline

In your script, the ODE function parcel_ode_kappa looks like this:

matlab
function dydt = parcel_ode_kappa(t, y, p)
    T = y(1);
    w = y(2);
    ri = y(3:end);

    S = compute_supersat(T, w, p);          % <--- called from another function
    dri_dt = droplet_growth_rate(ri, S, T, p);  % <--- another function

    dT_dt = 0;                              % <--- defined inline
    dw_dt = -sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) / p.rho_air;

    dydt = [dT_dt; dw_dt; dri_dt];
end

📌 Explained:

✅ compute_supersat(...) → defined as a separate function
Reason: It's a reusable physics-based function, involving vapor pressure and Clausius-Clapeyron relation.
✅ droplet_growth_rate(...) → also a separate function
Reason: Includes G(T) physics term, reused across variants and better organized when modular.
🔁 dT_dt = 0 → defined inline
Reason: In the no-heating case, temperature change is neglected or simplified. It's a fixed condition, not a formula needing modularization.

🔹 When Are Components Defined Inline vs Separately?

Component	Defined Separately?	Why?
`compute_supersat(...)`	✅ Yes	Nontrivial physics; reused across scenarios
`droplet_growth_rate(...)`	✅ Yes	Uses G(T); modularity & clarity
`compute_G(...)`	✅ Yes	Fundamental thermodynamic function
`dT_dt` in base case	❌ No	Trivial or scenario-specific
Heating terms	❌ or ✅	Depends on scenario (defined inline in scenario ODEs)

🔹 Takeaway

Rule of thumb in your script:

🔹 If the component is physics-based, reused, or depends on parameters, it gets its own function.
🔸 If it’s scenario-specific, trivial, or static (like dT_dt = 0 in the base case), it stays inline.

✅ Example: In Heating Scenarios

For heating cases (e.g. parcel_ode_internal_kappa), you'll see:

matlab
dT_dt = ... + p.Qext .* p.ga / p.Cp;

This may be left inline because the expression is still simple and doesn't benefit from being extracted unless used multiple times.

The script below contains all of the core component functions of the ODE system (compute_supersat, droplet_growth_rate, compute_G, latent_heating_rate, and vapor_consumption_rate) that are now explicitly defined.

The script contains two main parts: 2nd part -> script defining the ODE system function itself along with the constructing functions; 1st part-> the driver code that call the ODE functions.

This gives clarity and full control over how each physical process is modeled in both heating and non-heating scenarios. You can now modify or extend any part as needed.

%%%Part 1: The Driver Code That Calls the ODE Function

% Define parameters

params_base = struct( ...

'P', 85000, ... % Pressure [Pa]

'Nd', 100e6, ... % Droplet number concentration [#/m^3]

'kappa', 0.3, ... % Hygroscopicity

'rho_air', 1.2 ... % Air density [kg/m^3]

);

% Initial conditions: T0 [K], w0 [kg/kg], ri0 [m] for each droplet

T0 = 273.15 + 0.5; % Slightly above freezing

w0 = 0.010; % Initial water vapor mixing ratio

ri0 = 0.05e-6 * ones(50,1); % 50 small droplets, radius = 0.05 micron

y0 = [T0; w0; ri0]; % Concatenate all initial conditions

% Call ODE solver

[t, y] = ode15s(@(t,y) parcel_ode_kappa(t, y, params_base), [0 1800], y0);

%%%Part 2: The core ODE system function

function dydt = parcel_ode_kappa(t, y, p)

% Extract state variables

T = y(1);

w = y(2);

ri = y(3:end);

% Compute supersaturation

S = compute_supersat(T, w, p);

% Compute growth rate of droplets

dri_dt = droplet_growth_rate(ri, S, T, p);

% Compute latent heating (set to zero in no-heating scenario)

dT_dt = latent_heating_rate(ri, dri_dt, p, T, S);

% Compute change in water vapor

dw_dt = vapor_consumption_rate(ri, dri_dt, p);

% Combine derivatives

dydt = [dT_dt; dw_dt; dri_dt];

end

%%%%Each of the individual composing functions definition

%%supersaturation

function S = compute_supersat(T, w, p)

es = 611 * exp((17.27 * (T - 273.15)) / (T - 35.85));

ws = 0.622 * es / (p.P - es);

S = (w - ws) / ws;

end

%%growth rate of droplet_main

function dri_dt = droplet_growth_rate(ri, S, T, p)

G = compute_G(T, p);

ak = kappa_kohler(ri, p);

Seff = (1 + ak) .* S;

dri_dt = G .* Seff ./ ri;

end

%%growth rate cofficient-part of the growth rate of droplet function

function G = compute_G(T, p)

L = 2.5e6; % latent heat [J/kg]

Rv = 461; % gas constant for water vapor [J/kg/K]

Dv = 2.5e-5; % diffusivity of water vapor [m^2/s]

ka = 0.024; % thermal conductivity [W/m/K]

rho_w = 1000; % density of water [kg/m^3]

es = 611 * exp((17.27 * (T - 273.15)) / (T - 35.85));

G = (rho_w ./ (Rv * T)) .* (Dv + ka * L^2 / (Rv * T^2));

end

%%latent heating

function dT_dt = latent_heating_rate(ri, dri_dt, p, T, S)

% No heating by default

dT_dt = 0;

% Can override in scenarios with heating by implementing additional terms

end

%%change in water vapor

function dw_dt = vapor_consumption_rate(ri, dri_dt, p)

dw_dt = -sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) / p.rho_air;

end

%%kappa value

function ak = kappa_kohler(ri, p)

ak = p.kappa ./ ri; % simple kappa-Köhler expression

end

Important Note!!!

First....

In the above script of the cloud parcel model, only the total number concentration (Nd) is specified, rather than a full aerosol size distribution. This is a simplifying assumption, and here’s why and how it's used:

🔹 Why Only Total Number Concentration (`Nd`)?

Simplification for computational efficiency:
- Modeling the full aerosol size distribution requires solving many more ODEs (e.g., using bin-resolved or modal schemes), which makes the system larger and stiffer.
- For many idealized studies (e.g., understanding the effects of heating, supersaturation dynamics), a monodisperse assumption is sufficient.
Each droplet in the ODE is already being tracked individually:
- In your model, you initialize, say, 50 droplets (ri0 = 0.05e-6 * ones(50,1)), all with the same initial size.
- This represents a subset of the entire droplet population, assuming that all particles activate and grow in a similar way.
- You scale up the total water vapor consumption and heating effects by multiplying with Nd, like this:
```
dw_dt = -sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) / p.rho_air;
```

🔹 What If You Want to Use Size Distribution?

If you want to incorporate an aerosol size distribution, you need to:

Define the initial distribution (e.g., lognormal with parameters: N, Dg, σg).
Initialize droplets at multiple sizes representing the distribution.
Weight each droplet’s effect based on its share in the distribution (e.g., dN per bin).
Use either:
- Sectional (bin) method: divide particles into bins of different sizes.
- Modal method: track modes of the lognormal distribution analytically.

This makes the model more physically realistic but also much more computationally intensive.

🔹 Example: Lognormal Distribution with Bins

If you wanted to track size distribution:

matlab
% Generate 10 size bins from lognormal
Dp = logspace(log10(0.01e-6), log10(0.5e-6), 10);  % in meters
dN = lognormal_pdf(Dp, N, Dg, σg);                 % your own function

Then for each ri, you solve a droplet growth equation, and in the dw_dt, you multiply each term by its corresponding dN(i) instead of a uniform Nd.

✅ Summary

Feature	Current Model	Full Size Distribution Model
Aerosol size info	Only total `Nd`, monodisperse `ri0`	Lognormal or measured dN/dlogDp
Complexity	Low, idealized	High, realistic
Used for	Heating sensitivity, conceptual studies	Activation spectra, CCN, entrainment effects
Computation	Fast	Slow

More details: refer to the other posts

2nd...

the script does not currently include latent heat release upon condensation. The function responsible for that effect:

function dT_dt = latent_heating_rate(ri, dri_dt, p, T, S)
    % No heating by default
    dT_dt = 0;

    % Can override in scenarios with heating by implementing additional terms
end

is a placeholder — it explicitly sets dT_dt = 0, meaning temperature is not affected by condensation, which is not physically accurate in real cloud processes.

The basic assumption in the cloud parcel model is that:

Latent heat is released upon condensation on droplets.

Explanation:

The cloud parcel model follows the adiabatic ascent of an air parcel. As it rises, the air parcel expands and cools.
Once the temperature reaches the dew point, water vapor begins to condense onto aerosol particles (cloud condensation nuclei, or CCN).
This condensation releases latent heat, which partially offsets the cooling due to expansion. This is a critical feedback in the thermodynamics of cloud formation.

Therefore:

❌ No latent heat upon condensation → not assumed in realistic cloud parcel models. This would apply only in dry adiabatic or simplified models before condensation starts.
✅ Latent heat release upon condensation → assumed and included in most cloud parcel models, especially after cloud base.

Summary:

Assumption	Used in Cloud Parcel Model?	Notes
No latent heat released	❌	Only before condensation (dry ascent)
Latent heat released upon condensation	✅	Key to cloud formation, affects temperature and buoyancy

✅ How to Include Latent Heating

To incorporate the latent heating effect, you need to compute the rate of temperature change due to condensation:

🔧 Replace `latent_heating_rate` with:

matlab

function dT_dt = latent_heating_rate(ri, dri_dt, p, T, S)
    L = 2.5e6;           % Latent heat of vaporization [J/kg]
    cp = 1005;           % Specific heat of air at constant pressure [J/kg/K]
    rho_w = 1000;        % Density of water [kg/m^3]

    % Total mass flux of condensed water (kg/m^3/s)
    dmdt = sum(4 * pi * ri.^2 .* p.Nd .* dri_dt) * rho_w;

    % Temperature change due to release of latent heat
    dT_dt = (L / (p.rho_air * cp)) * dmdt;
end

This version computes the heat added to the air due to the mass of water condensing, and divides by air's heat capacity to get the temperature tendency.

🧠 Summary:

Component	Current Script	Physical Reality	Fix Required?
Latent heat release	❌ Ignored	✅ Important	✅ Yes
Function `latent_heating_rate`	`dT_dt = 0`	Should depend on `dri_dt` and latent heat	✅ Replace as shown above

Check here!! for the full corrected complete version of the script including the basic assumption of the cloud parcel model - condensation driven latent heat release

Cari Blog Ini

Aerosol-Awan