Building the Cloud Parcel Model on Matlab - Part 2: Create the core ODE system (B. use a realistic particle size distribution (PSD) instead of a single size)

INTRODUCTION

✅ What does generate the aerosol size distribution?

They generate the initial aerosol size distribution using a sum of lognormal modes, defined by:

Number concentration (N),
Geometric mean diameter (rg),
Geometric standard deviation (σ).

This gives you dN/dlogDp as a function of Dp (dry diameter), independently of composition.

Black carbon doesn't influence:

Mode centers (rg),
Widths (σ),
Total number per mode (N).

🟢 These are purely geometrical/statistical parameters for the aerosol distribution.

This code block sets up a realistic aerosol size distribution based on a tri-modal lognormal fit commonly used in cloud activation studies, including in Nenes et al. (2001, 2002). Let's break it down and explain each part:

🔧 Function Overview

function [Dp, Nd, kappa, bc_frac] = setup_aerosol_distribution(N_bins)

This function generates:

Dp: Particle diameters [m]
Nd: Number distribution dN/dlogDp [#/m³]
kappa: Hygroscopicity for each size bin
bc_frac: Black carbon mass fraction

🟦 1. Define Mode Parameters (based on Nenes 2002)

rg = [0.007, 0.027, 0.43]*1e-6;        % Geometric mean radii [m]
sigma = [1.8, 2.16, 2.21];            % Geometric standard deviations
N_cm3 = [106000, 32000, 5.4];         % Number concentrations [#/cm³]
N = N_cm3 * 1e6;                      % Convert to [#/m³]

These values correspond to:

Aitken, accumulation, and coarse modes.
These are typical for urban/continental air (based on observational fits).

🟦 2. Create Log-Spaced Particle Diameter Grid


Dp = logspace(log10(0.003e-6), log10(3e-6), N_bins);

Creates a logarithmic array of particle diameters between 0.003 µm and 3 µm. This covers most CCN-active particle sizes.

🟦 3. Build Multimodal dN/dlogDp Size Distribution

dNdlogDp = zeros(size(Dp));
for i = 1:3
    dNdlogDp = dNdlogDp + lognormal_dist(Dp, rg(i)*2, sigma(i), N(i));
end

For each mode (Aitken, accumulation, coarse), we:

Generate a lognormal distribution using lognormal_dist (a separate function - check in the upcoming sections).
Sum them together to get total aerosol size distribution.
Note: rg(i)*2 gives modal diameter, not radius.

🟦 4. Assign κ to Each Size Bin


kappa = assign_kappa(Dp);

This function should return a kappa value for each diameter bin. It may return:

A fixed value (e.g., 0.3 for ammonium sulfate)
Or vary with size (e.g., accounting for internal mixing or BC content)

🟦 5. Set Black Carbon Fraction

bc_frac = 0.1;

Sets 10% of aerosol mass to black carbon (BC), a common value for polluted urban air.
This will be used in heating scenarios (e.g. Scenario 3, 4 in Nenes 2002) to add feedback from BC absorption.

✅ Output Explanation

Output	Description
`Dp`	Diameter grid [m]
`Nd`	Size-resolved number concentration: dN/dlogDp [#/m³]
`kappa`	Array of κ values for each size
`bc_frac`	Fraction of aerosol mass that is BC

📌 Summary of Physical Meaning

This function enables the cloud parcel model to:

Use a realistic aerosol population (not just mono-disperse).
Simulate how particles of different sizes activate differently.
Model feedbacks from black carbon heating, if included.

Here are the supporting functions for lognormal_dist and assign_kappa used in the setup_aerosol_distribution function. These match the formulations typically used in atmospheric modeling studies like Nenes et al. (2001, 2002):

🔹 `lognormal_dist` — Generates a lognormal number distribution (dN/dlogDp)

To use the aerosol size di

This defines how particles are distributed across diameters in each mode.


function dNdlogDp = lognormal_dist(Dp, Dg, sigma, N)
    % Dp: particle diameter array [m]
    % Dg: geometric mean diameter [m]
    % sigma: geometric standard deviation
    % N: total number concentration for the mode [#/m^3]

    dNdlogDp = (N ./ (log(sigma) * sqrt(2 * pi))) .* ...
               exp(- (log(Dp ./ Dg)).^2 ./ (2 * log(sigma)^2));
end

📝 Explanation:

dNdlogDp gives particle number per unit logarithmic diameter.
This is the lognormal form typically used in aerosol modeling.
The integral of dNdlogDp over logDp gives total particle number N.

🔹 `assign_kappa` — Assign hygroscopicity (κ) for each bin

You can use a simple constant, or make it size-dependent to simulate mixing.

Option A: Constant κ (e.g., all ammonium sulfate)


function kappa = assign_kappa(Dp)
    kappa = 0.3 * ones(size(Dp)); % uniform κ value
end

Option B: Size-dependent κ (simple BC-core / sulfate-shell model)

function kappa = assign_kappa(Dp)
    % Assume small particles more BC-like, large ones more sulfate-like
    kappa = 0.1 + 0.2 * (Dp / max(Dp)); % κ increases with size
end

This kappa = 0.1 + 0.2 * (Dp / max(Dp)) line means:

κ ranges from 0.1 (smallest particles) to 0.3 (largest particles).
Assumes a continuous variation of hygroscopicity with size.
Implies internal mixing, where BC and sulfate are mixed in varying proportions as size increases.

This is heuristic, not physically tied to the three lognormal modes defined earlier.

Option C: different k for different modes of the full PSD

Nenes et al. (2002) used three lognormal modes, each likely representing different aerosol types:

Mode	Mode Size (rg)	Likely Composition	κ Value
1	~0.007 µm	Fresh BC (hydrophobic)	κ ≈ 0.0–0.05
2	~0.027 µm	Aged organics/BC mix	κ ≈ 0.1–0.2
3	~0.43 µm	Sulfate-rich (soluble)	κ ≈ 0.3

So to assign a different κ per mode, here's how you do it.

✅ Option C1: Assign κ by bin membership per mode

Assuming you've already summed the three lognormal modes into a single size distribution:

Step-by-step:


function kappa = assign_kappa_modewise(Dp, rg, sigma, N)
    % Define κ per mode (matching rg)
    kappa_modes = [0.05, 0.15, 0.3];

    % Preallocate numerator and denominator for weighted κ
    kappa = zeros(size(Dp));
    total_contrib = zeros(size(Dp));

    for i = 1:length(rg)
        % Get number concentration from each mode at all Dp
        dNdlogDp_mode = lognormal_dist(Dp, rg(i), sigma(i), N(i));

        % Add weighted contribution to kappa
        kappa = kappa + dNdlogDp_mode .* kappa_modes(i);
        total_contrib = total_contrib + dNdlogDp_mode;
    end

    % Normalize to get mean κ at each Dp
    kappa = kappa ./ total_contrib;
end

This gives you a size-resolved κ that:

Reflects the contribution of each mode at each size bin.
Still assumes internal mixing (each bin can have mixed modes).

✅ Option C2: External mixture (Mode-split κ assignment)

If instead you treat the modes separately (external mixture), you could:

Compute each mode's dN/dlogDp separately.
Assign a single κ value to all particles from that mode.
Concatenate into one vector.


Dp = logspace(log10(0.003e-6), log10(3e-6), N_bins);
kappa_all = [];
Nd_all = [];

for i = 1:3
    dNdlogDp = lognormal_dist(Dp, rg(i), sigma(i), N(i));
    kappa_mode = kappa_modes(i) * ones(size(Dp));
    
    Nd_all = [Nd_all, dNdlogDp];
    kappa_all = [kappa_all, kappa_mode];
end

Now you have:

Nd_all: number distribution from all modes (separately).
kappa_all: matched κ values.
And this now represents external mixing.

🔬 Conclusion

If you want internal mixing → assign κ based on weighted mode contributions at each size.
If you want external mixing → treat each mode as distinct and assign one κ per mode.
Both are valid — your choice depends on the physical assumptions you're modeling (e.g. BC aging, aerosol processing, etc.).

✅ Summary Table

Function	Purpose	Output
`lognormal_dist`	Generate tri-modal aerosol distribution	`dN/dlogDp` [#/m³]
`assign_kappa`	Assign hygroscopicity for each size bin	`κ` array

Recap of the differences between Option C1 and C2:

Let’s walk through the logic step by step:

🧠 What Does `assign_kappa_modewise` (Option C1) Do?

This function does not assign a κ value per mode to the particles. Instead, it assigns a κ value per particle size, by combining the contribution from all modes that overlap at that size.

How?

It loops over all three modes (with different κ_mode(i) values).
For each diameter bin Dp(j):
- It computes the number of particles in that bin from each mode.
- It multiplies each mode's contribution by its respective κ value.
- It sums these and divides by the total number in that bin.
The result is:
$\kappa(D_p) = \frac{\sum_i \kappa_i \cdot \frac{dN_i}{d\log D_p}}{\sum_i \frac{dN_i}{d\log D_p}}$
where $\kappa_i$ is the κ for mode $i$ .

🎯 What Does This Mean Physically?

You’re not assuming each particle is of one type (pure BC, pure sulfate, etc).

Instead, you’re saying:

“At each size $D_p$ , particles come from multiple modes, and therefore are a mixture of components.”

So at Dp = 0.2 µm, for example:

Mode 1 contributes a bit (κ = 0.05),
Mode 2 contributes more (κ = 0.15),
Mode 3 also contributes (κ = 0.3),
→ the effective κ is a weighted blend.

🔁 This is Internal Mixture:

Each particle of size $D_p$ is assumed to contain material from multiple sources.
You assign one κ per particle, but that κ is the average composition of that size, based on all contributing modes.

❌ What Would External Mixture (Option C2) Be Instead?

For an external mixture, you would:

Treat each mode completely separately in the simulation.
Each mode has its own size distribution and fixed κ.
You would simulate growth per mode, not mix them.
There would be no averaging of κ across modes.

So:

Mode 1 → simulate growth with κ = 0.05
Mode 2 → simulate growth with κ = 0.15
Mode 3 → simulate growth with κ = 0.3
Each produces its own droplet growth trajectory.

✅ Summary Table

Approach	κ assignment	Particle composition	Mixture type
`assign_kappa_modewise`	Weighted κ per Dp based on all modes	Each size contains mixed material	Internal
Fixed κ per mode, no mixing	Separate growth per mode	Each particle is pure material	External

🧪 Analogy

Imagine again a fruit salad:

You measure how much apple, orange, banana exist at each spoonful size.
You average the composition by weight for each spoon size.
Now each spoon has a different mix — this is internal mixture.

If instead:

You say “this spoon is only apple slices,” and another is “only banana,”
That's an external mixture.

🔹 Step-by-Step Conversion: From dN/dlogDp to Initial Parcel Inputs

To use the aerosol size distribution (`dN/dlogDp`) as input to the parcel model, we need to convert it into initial radius and number of droplets per bin. These will initialize the growth of droplets in the ODE system.
Here is how to set it up:

🔹 Step 1: Setup aerosol distribution

Using the function you've already defined:


[Dp, Nd, kappa, bc_frac] = setup_aerosol_distribution(50);

Dp is the particle diameter in [m] (center of each bin)
Nd is dN/dlogDp [#/m³]
kappa is assigned κ for each bin

🔹 Step 2: Convert to droplet number per bin

To use in the parcel model, we need total number per bin N_bin, not just dN/dlogDp.

Use numerical integration over logDp to get per-bin values:


% Compute d(logDp)
logDp = log10(Dp);
dlogDp = diff(logDp);
dlogDp = [dlogDp; dlogDp(end)]; % repeat last for equal length

% Approximate bin number: Nd * dlogDp
N_bin = Nd .* dlogDp; % units [#/m³]

🔹 Step 3: Define initial dry radius

Initial radius is dry particle radius before water uptake:


ri0 = Dp(:) / 2;  % radius in meters

🔹 Step 4: Bundle into initial condition `y0`

Parcel model needs:


T0 = 273.65;     % Initial temperature [K]
w0 = 0.010;      % Initial water vapor mixing ratio [kg/kg]
y0 = [T0; w0; ri0];  % Initial state vector

Then pass N_bin and kappa into params:


params = struct( ...
    'P', 85000, ...
    'Nd', N_bin, ...      % now size-resolved
    'kappa', kappa, ...
    'rho_air', 1.2 ...
);

🔹 Step 5: ODE solver


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

🧠 Why This Matters

Quantity	Purpose
`N_bin`	Number of droplets per size bin
`ri0`	Dry radius per bin
`kappa`	Growth sensitivity to water

This setup mimics how lognormal aerosol distributions are used in cloud parcel models like Nenes et al. (2002), allowing multimodal, size-resolved activation and latent heating feedback (if enabled).

NOTE:

in Nenes et al. (2002), the black carbon (BC) fraction does not affect the initial size distribution generation. Here's how it works, clearly:

🚫 What BC does not do:

It does not create a new mode in the size distribution.
It does not change the total number concentration or the shape of the distribution.

✅ What BC does affect:

Activation behavior:
BC-containing particles have a lower hygroscopicity (κ) → they require higher supersaturation to activate into cloud droplets.
Thermal feedback:
BC absorbs radiation, locally heats the parcel → changes temperature (T) → affects vapor pressure and supersaturation.
Heterogeneous internal mixing:
You can tag a fraction of particles (e.g. 10%) in each size bin as BC-containing and assign:
- Lower κ values to them,
- Heating terms in the energy equation.

But again: this tagging happens after generating the initial dN/dlogDp.

🧪 Optional: Making a BC-sensitive size distribution

If you want BC to affect the distribution itself, you would need:

Source-specific emissions (e.g. biomass vs fossil),
Fresh vs aged BC (coated with organics or sulfate),
A mixing model that couples chemistry and microphysics.

That goes beyond Nenes (2002), into chemical transport models or aerosol aging modules.

So, how do they differentiate internal and external mixed BC? Both have the same size distribution?

This is great and an important one for understanding how black carbon (BC) influences cloud activation. Let’s break it down using the framework in Nenes et al. (2002) and general aerosol-cloud interaction modeling:

🔹 Internal vs External Mixture — What’s the Difference?

Type of Mixing	Definition	Effect on κ	Effect on Heating	Impact on Activation
Internal mixing	Each particle contains both BC and soluble material (e.g. sulfate)	Intermediate κ (e.g. 0.1–0.3) depending on the fraction of soluble	All particles contribute to absorption (if BC is well-distributed)	Lower activation fraction (suppressed growth due to BC core)
External mixing	BC and soluble particles are separate populations	Two κ values: low for BC (κ ≈ 0.0), higher for soluble (κ ≈ 0.3–0.6)	Only BC particles heat parcel	Only soluble particles activate (unless high S)

🔸 Do They Share the Same Size Distribution?

Yes and no.

1. Geometric Distribution (dN/dlogDp):

Both are often initialized using the same size distribution shape (same bins, same rg and σ) — especially in simplified models like the one in Nenes et al. (2002).

You might assign, say, 90% of the particles as soluble and 10% as BC, within each bin → same shape, different κ or heating effect.
This is common for representing internally mixed aerosols.

2. Differentiated Sub-populations (External):

In external mixing:

You could split the total distribution into two separate number distributions.
For example:
- One lognormal mode for sulfate (e.g. κ = 0.3),
- Another for BC (e.g. κ = 0.0).

Each may have:

Different number concentration (N),
Different rg and σ.

So, in external mixing, BC and soluble aerosol may have distinct size distributions.

🔹 How Did Nenes et al. (2002) Handle It?

In their paper:

They mostly assumed internal mixing for simplicity in the parcel model (more consistent with atmospheric aging),
Used effective κ values to represent mixtures,
Explored heating scenarios separately, allowing them to simulate BC absorption effects even if not fully coupled to aerosol composition.

🔸 Implementation Example in Code

Let’s say you want 10% of particles in each bin to be BC internally mixed:


% In each bin:
kappa = (1 - bc_frac) * kappa_soluble + bc_frac * kappa_BC;

Or, for external mixing (two separate populations):

% Define two Nd vectors:
Nd_soluble = Nd * (1 - bc_frac);
Nd_bc = Nd * bc_frac;

% Assign separate κ values:
kappa_soluble = 0.3;  % sulfate
kappa_bc = 0.0;       % pure BC

🔍 Summary

Internal mixing: One size distribution, adjusted κ and possibly heating.
External mixing: Separate size distributions for BC and soluble particles.
Nenes (2002) mostly used internal mixing, with heating effects added separately.
The choice impacts droplet activation, heating, and ultimately cloud microphysics and radiative effects.

SCRIPT....

%%% INTERNAL vs EXTERNAL MIXTURE COMPARISON

%%% (1) Generate aerosol distribution using 3 lognormal modes

function [Dp, dNdlogDp, kappa_internal, kappa_external] = aerosol_mixture_comparison(N_bins)

% Define mode parameters (same as Nenes 2001/2002)

rg = [0.007, 0.027, 0.43] * 1e-6; % geometric mean radius [m]

sigma = [1.8, 2.16, 2.21]; % geometric standard deviation

N_cm3 = [106000, 32000, 5.4]; % number concentration [cm^-3]

N = N_cm3 * 1e6; % convert to #/m^3

% Kappa values per mode (e.g. BC, organics, sulfate)

kappa_mode = [0.05, 0.15, 0.3];

% Generate diameter bins

Dp = logspace(log10(0.003e-6), log10(3e-6), N_bins); % particle diameter [m]

dlogDp = log10(Dp(2)) - log10(Dp(1));

dNdlogDp = zeros(N_bins, 3);

% Fill in size distribution for each mode

for i = 1:3

dNdlogDp(:, i) = lognormal_dist(Dp, rg(i)*2, sigma(i), N(i));

end

% Total size distribution (sum of all modes)

dNdlogDp_total = sum(dNdlogDp, 2);

% ----- INTERNAL MIXTURE CASE -----

kappa_internal = zeros(size(Dp));

for j = 1:N_bins

if dNdlogDp_total(j) > 0

kappa_internal(j) = sum(kappa_mode .* dNdlogDp(j, :) ) / dNdlogDp_total(j);

else

kappa_internal(j) = NaN;

end

% ----- EXTERNAL MIXTURE CASE -----

% Assign fixed kappa to each mode; each mode handled separately

% This just returns 3 separate kappa arrays

kappa_external = repmat(kappa_mode, N_bins, 1);

% Optional plot

figure;

subplot(2,1,1);

loglog(Dp*1e6, dNdlogDp_total, 'k-', 'LineWidth', 2); hold on;

loglog(Dp*1e6, dNdlogDp, '--');

xlabel('Dp [μm]'); ylabel('dN/dlogDp [#/m^3]');

legend('Total', 'Mode 1', 'Mode 2', 'Mode 3');

title('Size Distribution');

subplot(2,1,2);

semilogx(Dp*1e6, kappa_internal, 'r-', 'LineWidth', 2); hold on;

semilogx(Dp*1e6, kappa_external, '--');

xlabel('Dp [μm]'); ylabel('κ');

legend('Internal mixture', 'Mode 1 (ext)', 'Mode 2 (ext)', 'Mode 3 (ext)');

title('Kappa Values');

end

%%% Helper: lognormal distribution

function dNdlogDp = lognormal_dist(Dp, rg, sigma, N)

% dN/dlogDp = N / (sqrt(2pi) * log(sigma)) * exp(-(ln(Dp/rg))^2 / (2 log^2(sigma)))

lnDp = log(Dp);

lnsigma = log(sigma);

lnrg = log(rg);

dNdlogDp = N ./ (sqrt(2*pi) * lnsigma) .* exp( - (lnDp - lnrg).^2 ./ (2 * lnsigma^2));

end

Explanations..

Here's the full MATLAB script that shows the difference between internal and external mixtures of aerosol particles based on the κ-Köhler theory:

Key Concepts:

Internal mixture: All modes are treated as mixed together. A size bin has contributions from all modes, so its κ is the weighted average of those contributions.
External mixture: Each mode retains its identity. The κ for particles in each mode is fixed and distinct. No averaging across modes.

What the script does:

It defines a trimodal aerosol population based on Nenes (2002).
Assigns different κ values to each mode (e.g., 0.05 for BC-like, 0.3 for sulfate).
Plots both total number size distribution and the κ value as a function of diameter for both internal and external mixtures.

Let me know if you want to simulate activation, hygroscopic growth, or CCN activity from this basis.