Building the Cloud Parcel Model on Matlab

📘 What is an ODE?

ODE stands for Ordinary Differential Equation.

An ODE is an equation that relates a function and its derivatives with respect to a single independent variable, usually time.

🔢 General Form of an ODE

A first-order ODE:

\frac{dy}{dt} = f(t, y)

$y(t)$ is the unknown function we want to solve.
$f(t, y)$ describes how $y$ changes over time.

You solve it to find $y(t)$ , given:

A starting point $t_0$ ,
An initial value $y(t_0) = y_0$

🧪 Example: Cooling of a Cup of Coffee

\frac{dT}{dt} = -k(T - T_{\text{room}})

$T(t)$ : temperature of the coffee.
$T_{\text{room}}$ : ambient temperature.
$k$ : cooling constant.
The equation says: the hotter the coffee is compared to the room, the faster it cools.

This is a first-order ODE.

🌦 In Your Cloud Parcel Model

Your ODE system includes:

$\frac{dT}{dt}$ : change in temperature
$\frac{dw}{dt}$ : change in vapor mixing ratio
$\frac{dr_i}{dt}$ : change in radius of each droplet

These are solved simultaneously as a system of ODEs:

\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y})

Where $\mathbf{y}$ is a vector like:

\mathbf{y} = [T, w, r_1, r_2, \dots, r_{100}]

🧮 Solving ODEs

In MATLAB, functions like ode45, ode15s, etc., numerically solve ODEs — that is, they compute approximate values of $y(t)$ over time using time stepping methods.

🔄 What is Considered an ODE System?

An ODE system is when you have multiple ODEs that are coupled together — meaning, they influence each other and must be solved simultaneously.

🧮 General Form of an ODE System

Instead of one equation:

\frac{dy}{dt} = f(t, y)

You have many equations, typically written as a vector system:

\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y})

Where:

$\mathbf{y} \in \mathbb{R}^n$ is a vector of unknowns (e.g., temperature, humidity, radii, etc.).
$\mathbf{f}$ is a vector-valued function returning one derivative per variable.

🌦 Example in Cloud Parcel Model

Your parcel model is a classic example of an ODE system.

Suppose:

$T$ : temperature,
$w$ : vapor mixing ratio,
$r_i$ : radius of droplet $i$ , for $i = 1$ to 100 bins.

Then your system is:

\begin{aligned} \frac{dT}{dt} &= f_1(t, T, w, r_1, ..., r_{100}) \\ \frac{dw}{dt} &= f_2(t, T, w, r_1, ..., r_{100}) \\ \frac{dr_1}{dt} &= f_3(t, T, w, r_1) \\ &\vdots \\ \frac{dr_{100}}{dt} &= f_{102}(t, T, w, r_{100}) \end{aligned}

So, this is a 102-dimensional ODE system, where the change in each variable depends on the others.

✅ When You Know It's a System

You solve more than one differential equation at the same time.
The variables are interdependent.
You pass a vector of initial conditions to the solver.
The solver gives you back a matrix of solutions (each column for one variable).

Would you like an illustrative diagram or simple code example of a small ODE system to compare with your parcel model?

Great! Let's break this down step by step.

🧩 What Does "Interdependent" Mean?

Interdependent means the equations (or variables) affect each other.

Example:

If the growth of variable $y_1$ depends on the current value of $y_2$ ,
and the change in $y_2$ also depends on $y_1$ ,
then they are interdependent and must be solved as a system.

🔁 Simple Example of an Interdependent ODE System

Let’s consider a very basic 2-variable ODE system:

\begin{cases} \frac{dy_1}{dt} = -0.5 y_1 + y_2 \\ \frac{dy_2}{dt} = -y_1 - 0.1 y_2

$y_1$ and $y_2$ influence each other's rate of change.
You cannot solve them independently.

💻 MATLAB Code Example

matlab
% Define the system of ODEs
function dydt = simple_system(t, y)
    dydt = zeros(2,1);
    dydt(1) = -0.5 * y(1) + y(2);   % dy1/dt
    dydt(2) = -y(1) - 0.1 * y(2);   % dy2/dt
end

% Initial condition: y1(0)=1, y2(0)=0
y0 = [1; 0];

% Time span
tspan = [0 20];

% Solve using ode45
[t, y] = ode45(@simple_system, tspan, y0);

% Plot
figure;
plot(t, y(:,1), 'r-', t, y(:,2), 'b--');
legend('y_1(t)', 'y_2(t)');
xlabel('Time'); ylabel('Values');
title('Interdependent ODE System');
grid on;

📉 Diagram of Interdependency

y₁ ———► y₂
↑        ↓
└────────┘

y₁ affects y₂, and y₂ affects y₁.
This loop creates feedback → requires solving together.

🌩 Comparison to Your Parcel Model

In a parcel model:

$T$ (temperature) affects supersaturation.
Supersaturation affects droplet growth $r_i$ .
Growing droplets consume vapor, which reduces $w$ .
Change in $w$ changes $T$ via latent heat release.

So all variables — T, w, r₁, r₂, ..., rₙ — are interdependent.

If you're simulating a process where some variables change very fast while others change slowly, ode15s is more efficient and stable than standard solvers like ode45.

🛠️ Why We Use `ode15s` Instead of `ode45` for the Cloud Parcel ODE System in MATLAB

ode15s in MATLAB is a numerical solver used to integrate stiff systems of ordinary differential equations (ODEs).

📌 In simple terms:

If you're simulating a process where some variables change very fast while others change slowly, ode15s is more efficient and stable than standard solvers like ode45.

🔧 Technical Description:

matlab
[t, y] = ode15s(odefun, tspan, y0)

odefun — function handle to your system of ODEs (e.g., @(t,y) f(t,y)).
tspan — time interval to simulate (e.g., [0 1800] seconds).
y0 — initial condition vector for the state variables.
t — output time points.
y — solution matrix: each row corresponds to t, each column to a state variable.

🧠 When to Use `ode15s` (instead of `ode45`):

Use ode15s when:

You have stiff equations — common in chemical kinetics, cloud microphysics, combustion, etc.
The solver with ode45 becomes very slow or unstable.
There are widely varying timescales in your model (e.g., rapid condensation with slow temperature rise).

🔍 Why It's Used in Cloud Parcel Models:

In cloud microphysics (like your script simulating droplet growth and supersaturation evolution), the ODEs become stiff because:

Droplet radii can grow fast under supersaturation.
Supersaturation itself responds quickly to droplet growth.
Heat and mass transfer processes evolve at different speeds.

Using ode15s helps avoid tiny time steps and instabilities, making the simulation accurate and efficient.

CLARITY....
Again, what is stiffness?

What is a stiff system?

A stiff system of ODEs has:

Rapidly changing components alongside slow ones.
The system requires extremely small time steps if solved with explicit solvers like ode45, even though the solution itself might be smooth.

In this example:

ode45 struggles to maintain accuracy due to the stiffness, resulting in very small steps and sometimes instability.
ode15s (stiff solver) handles it efficiently, using implicit methods that remain stable even with larger steps.

Summary:

Feature	`ode45` (RK45)	`ode15s` (BDF)
Method Type	Explicit Runge-Kutta	Implicit multistep
Best for	Non-stiff problems	Stiff problems
Speed on stiff	Very slow or unstable	Fast and stable
Step size	Tiny (in stiff regions)	Can be large

In differential equations, a stiff system is one where some components evolve very quickly while others change slowly. This causes trouble for standard solvers like ode45, which must take very small time steps to accurately capture the fast changes — even if you're more interested in the slow dynamics.

🔹 Real-World Analogy:

Imagine you're filming a scene with both a snail and a hummingbird. If you use a regular camera speed (like ode45), the hummingbird's wings become a blur — you'd need a high-speed camera (like ode15s) to capture both properly without wasting effort.

🔸 Example: Chemical Reaction Kinetics

Let’s say you have a simple system:

\begin{aligned} \frac{d y_{1}}{d t} & = - 1000 y_{1} + y_{2} \\ \frac{d y_{2}}{d t} & = y_{1} - y_{2​} \end{aligned}

This models two reacting species:

y₁ changes very fast due to the -1000 y₁ term
y₂ changes more slowly

🔸 MATLAB Demo Code

matlab
function stiff_demo
    tspan = [0 1];
    y0 = [1; 0];

    % Solve with stiff solver
    [t1, y1] = ode15s(@stiff_system, tspan, y0);

    % Solve with non-stiff solver
    [t2, y2] = ode45(@stiff_system, tspan, y0);

    % Plot comparison
    figure; hold on;
    plot(t1, y1(:,1), 'r-', 'DisplayName','y1 - ode15s');
    plot(t1, y1(:,2), 'r--', 'DisplayName','y2 - ode15s');
    plot(t2, y2(:,1), 'b-', 'DisplayName','y1 - ode45');
    plot(t2, y2(:,2), 'b--', 'DisplayName','y2 - ode45');
    xlabel('Time'); ylabel('Solution');
    title('Stiff System: ode15s vs ode45');
    legend; grid on;
end

function dydt = stiff_system(t, y)
    dydt = [-1000*y(1) + y(2);
             y(1) - y(2)];
end

🔸 Observation:

ode15s takes fewer, well-spaced steps and remains stable.
ode45 takes many small steps, may be slow or unstable due to the stiff term -1000*y₁.

✅ Conclusion:

Stiffness = Fast + Slow dynamics mixed together.
Use ode15s, ode23s, or ode23t for such problems.

In a cloud parcel model, stiffness arises naturally because different physical processes occur on vastly different time scales. Here's a real example of stiffness in the cloud microphysics context:

🌧️ Example: Aerosol Activation and Droplet Growth

In the parcel model, you typically integrate variables like:

T(t) = temperature
w(t) = water vapor mixing ratio
rᵢ(t) = droplet radius for each aerosol bin (many sizes)
S(t) = supersaturation (depends on T and w)

⚠️ Where Stiffness Appears

Droplet Growth Equation (Fast dynamics):
$\frac{dr_i}{dt} = \frac{G S}{r_i}$
- Early in the simulation, when droplets are tiny (small rᵢ), this derivative becomes very large.
- The droplets grow extremely rapidly over milliseconds, especially for the most activated particles.
Thermodynamics & Condensation (Slow dynamics):
- T(t) and w(t) change more gradually as the air parcel ascends.
- These change over seconds to minutes — much slower than droplet growth in initial phases.

🧠 Summary of Time Scales:

Process	Time Scale	Comment
Droplet activation	milliseconds	Very rapid radius change
Supersaturation response	milliseconds–sec	Tight feedback with condensation
Temperature evolution	seconds	Slow due to latent heat balance
Aerosol processing	seconds–minutes	Gradual

💻 Stiffness in Practice:

Using ode45 on this system might:
- Take tiny time steps initially for stability.
- Fail or take excessive time if radii grow fast.
ode15s can adapt by treating fast processes implicitly, skipping over the small time scales where needed.

🔬 Real Implication in Research:

In simulations like the one you're doing (e.g., with BC heating), early droplet growth can dominate dynamics briefly, but you're mostly interested in CDNC and effective radius at 1800s. You don’t need millisecond resolution — hence, stiff solvers are ideal.

Here’s a breakdown of which terms in your cloud parcel model script are stiff, and how ode15s helps:

⚙️ Core ODE System in Your Script

In the function parcel_ode_kappa and its variants, you solve the following:

matlab
dydt = [dT_dt; dw_dt; dri_dt];

Where:

T(t) is temperature
w(t) is water vapor mixing ratio
rᵢ(t) is the radius of the i-th droplet bin (for i = 1 to 100 bins)

🔍 Identify Stiff Components

1. Droplet Growth (`dri_dt`) — STIFF

matlab
dri_dt = G * S ./ ri;

This term is very large when ri is small (i.e., just activated or pre-activated aerosol).
Early in the simulation, ri ≈ 0 ⇒ rapid radius increase.
Small ri → division by small number → large derivatives → numerical stiffness.

👉 This makes rᵢ(t) components evolve on sub-second or millisecond time scales.

2. Water Vapor (`dw_dt`) — MODERATELY FAST

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

This is a cumulative sink of water vapor due to condensation.
Strongly coupled to dri_dt: as droplets grow fast, w must drop quickly.
Can be somewhat stiff depending on how many droplets activate at once.

3. Temperature (`dT_dt`) — SLOW (in your base model)

matlab
dT_dt = 0;

In your current no-heating scenario, dT_dt is explicitly set to zero.
In the heating scenarios (droplet, external, internal), it will be non-zero due to latent heat or radiation.

Still, temperature changes slowly, unless external energy sources (e.g. BC absorption) dominate.

🧮 Why You Use `ode15s`

The fast radius growth (dri_dt) and its interaction with dw_dt cause stiffness.
If you use ode45:
- MATLAB would need tiny time steps to handle early rapid growth.
- Slow-changing variables like T would be updated unnecessarily often.
ode15s recognizes that ri grows fast, but T and w don’t need to be resolved at such high resolution.
It uses implicit numerical techniques that are stable for stiff systems, allowing:
- Larger time steps when possible.
- Efficient long-term simulation, even if early parts are rapid.

📌 Where to Annotate in Your Code

Here’s how you can annotate your script (you can add this comment block above your ODE function definitions):

matlab

%% Note on Stiffness
% This system is stiff due to rapid growth of droplet radius (ri)
% immediately after activation. Small ri values cause large dri_dt,
% creating rapidly changing dynamics. Meanwhile, w and T evolve more slowly.
% To handle this multiscale behavior, we use the stiff solver ode15s.