diff --git a/docs/src/IceNucleationParcel0D.md b/docs/src/IceNucleationParcel0D.md
index baaccb2eba..2da76dd83d 100644
--- a/docs/src/IceNucleationParcel0D.md
+++ b/docs/src/IceNucleationParcel0D.md
@@ -296,7 +296,7 @@ include("../../parcel/Example_Liquid_only.jl")
 The plots below are the results of the adiabatic parcel model
   with immersion freezing, condensation growth, and deposition growth for
   both a monodisperse and gamma size distribution. Note that this has not
-  yet been validated against literature. Results are very sensitive to 
+  yet been validated against literature. Results are very sensitive to
   initial conditions.
 
 ```@example
@@ -338,77 +338,59 @@ include("../../parcel/Example_P3_ice_nuc.jl")
 
 ## Immersion Freezing Parametrization based on Frostenberg et al. 2023
 
-The parcel model also includes the parametrization of immersion freezing based on [Frostenberg2023](@cite).
-The concentration of ice nucleating particles (INPs) is randomly selected based on a lognormal relative frequency distribution, depending only on the temperature. If the random INP concentration exceeds the existing concentration of ice crystals, new ice crystals are created, such that their total concentration is equal to the new INP concentration.
+Here we show the parcel model results when using the parametrization of immersion freezing
+  based on [Frostenberg2023](@cite).
+The concentration of ice nucleating particles (INPC) depends only on air temperature,
+  and is based on a lognormal relative frequency distribution.
+New ice crystals are created if the INPC exceeds the existing concentration of ice crystals,
+  provided there are sufficient numbers of cloud liquid droplets to freeze.
 
 Three different implementations of this parametrization are used in the parcel model:
-- **Frostenberg_mean**, in which the concentration of INPs is equal to its mean value defined in [Frostenberg2023](@cite) for a given temperature,
-- **Frostenberg_random**, in which the concentration of INPs is drawn randomly from the distribution defined in [Frostenberg2023](@cite). The time between draws is set by the *drawing_interval* parameter,
-- **Frostenberg_stochastic**, in which the freezing is treated as a stochastic process, with mean defined in [Frostenberg2023](@cite) and Wiener noise. The timescale of the process is set by the `\gamma` parameter. 
-
-
-The stochastic implementation is based on the equation for a stochastic process $x$ with Gaussian random noise:
-
+- `mean` - in which INPC is equal to its mean value defined in [Frostenberg2023](@cite).
+- `random` - in which INPC is sampled randomly from the distribution defined in [Frostenberg2023](@cite).
+  The number of model time steps between sampling is set by `sampling_interval`.
+- `stochastic` - in which INPC is solved for as a stochastic process,
+  with the mean and standard deviation defined in [Frostenberg2023](@cite).
+  The inverse timescale of the process is set by ``\gamma``.
+
+The stochastic implementation is based on the equation for a generic
+  stochastic process ``x`` with Gaussian random noise (i.e. Wiener process):
 ```math
 \begin{equation}
-  dx = - \gamma \; x \; dt + g dW
+  dx = - \gamma \; x \; dt + g dW; \;\;\;\;\;\;\;   dW = {\bf{N}}(0, dt^2),
 \end{equation}
 ```
-
-```math
-\begin{equation}
-  dW = {\bf{N}}(0, dt^2)
-\end{equation}
-```
-where ``\bf{N}`` is the lognormal distribution.
-
+where ``\bf{N}`` is the normal distribution.
+For constant ``\gamma`` and ``g``, ``x`` has solution
 ```math
 \begin{equation}
   x(t) = x_0 \; e^{-\gamma t} + g \; \int_0^t e^{\gamma(s-t)} dW
 \end{equation}
 ```
 where ``1/\gamma`` is the assumed timescale of the process.
+We can calculate the variance ``\bf{V}`` as,
 ```math
 \begin{equation}
-  {\bf{V}}(x) = g^2 \int_0^t e^{2\gamma(s-t)} ds = \frac{g^2}{2\gamma}(1 - e^{-2\gamma t})
+  {\bf{V}}(t) = g^2 \int_0^t e^{2\gamma(s-t)} ds = \frac{g^2}{2\gamma}(1 - e^{-2\gamma t}).
 \end{equation}
 ```
-where ``\bf{V}`` is the variance.
-We map this notation onto our problem:
 
+We use this process to model ``x=\log(\text{INPC})``,
+  which tends toward the mean ``\mu(T)``.
+If we denote ``\tau = 1 / \gamma`` as the process timescale and
+  ``\sigma = g / \sqrt{2 / \gamma}`` as the process uncertainty,
+  then
 ```math
 \begin{equation}
-  x = ln(INPC)
+  d\log(\text{INPC}) = - \frac{\log(\text{INPC}) - μ}{\tau} + \sigma \sqrt{\frac{2}{\tau \; dt}} \; {\bf{N}}(0, 1))
 \end{equation}
 ```
 
-```math
-\begin{equation}
-  dx = -\gamma (x - \mu(T)) dt + g dW
-\end{equation}
-```
-
-```math
-\begin{equation}
-  x(t) = x e^{-\gamma t} + \mu(T) (1 - e^{-\gamma t}) + g \int_0^t e^{-\gamma (t-s)} dW
-\end{equation}
-```
-
-```math
-\begin{equation}
-  x(t) = {\bf{N}}(x_0 e^{-\gamma t} + (1-e^{-\gamma t}) \mu(T), \frac{g^2}{2\gamma} (1-e^{-2\gamma t}))
-\end{equation}
-```
-If ``\gamma >> 1``, then
-```math
-\begin{equation}
-  x(t) = {\bf{N}}(\mu(T), \frac{g^2}{2\gamma})
-\end{equation}
-```
-Hence ``g = \sigma * \sqrt{2 \gamma}``.
-
-
-The following plot shows resuls of the parcel model with Frostenberg_mean parametrization, Frostenberg_random parametrization with different drawing intervals `t_d` and Frostenberg_stochastic parametrization with different time scales `\gamma`. The timestep of the model is 1 s.
+The following plot shows resuls of the parcel model with the `mean` (black line),
+  `random` (dotted lines) and `stochastic` (solid lines) parameterization options.
+We show results for two sampling intervals ``\Delta t`` (random),
+  two process time scales ``\tau`` (stochastic),
+  and two model time steps `dt`.
 
 ```@example
 include("../../parcel/Example_Frostenberg_Immersion_Freezing.jl")
diff --git a/parcel/Example_Deposition_Nucleation.jl b/parcel/Example_Deposition_Nucleation.jl
index 64aaa97d26..d086792cd5 100644
--- a/parcel/Example_Deposition_Nucleation.jl
+++ b/parcel/Example_Deposition_Nucleation.jl
@@ -21,6 +21,7 @@ qᵥ = FT(3.3e-4)
 qₗ = FT(0)
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -32,7 +33,7 @@ eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 e = eᵥ(qᵥ, p₀, Rₐ, Rᵥ)
 
 Sₗ = FT(e / eₛ)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles
diff --git a/parcel/Example_Frostenberg_Immersion_Freezing.jl b/parcel/Example_Frostenberg_Immersion_Freezing.jl
index a71b08d7c2..a13babc47a 100644
--- a/parcel/Example_Frostenberg_Immersion_Freezing.jl
+++ b/parcel/Example_Frostenberg_Immersion_Freezing.jl
@@ -2,9 +2,12 @@ import OrdinaryDiffEq as ODE
 import CairoMakie as MK
 import Thermodynamics as TD
 import CloudMicrophysics as CM
-import CLIMAParameters as CP
+import CloudMicrophysics.HetIceNucleation as CMI
+import ClimaParams as CP
+
 import Random as RAND
-random_seeds = [0, 1234, 5678, 3443]
+RAND.seed!(44)
+N_ensemble = 32
 
 # definition of the ODE problem for parcel model
 include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
@@ -22,9 +25,10 @@ r₀ = FT(1e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(251)
 qᵥ = FT(8.1e-4)
-qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
+qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # ρₐ ~ 1.2
 qᵢ = FT(0)
 x_sulph = FT(0.01)
+ln_INPC₀ = FT(CMI_het.INP_concentration_mean(T₀))
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -33,162 +37,163 @@ Rₐ = TD.gas_constant_air(tps, q)
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 e = eᵥ(qᵥ, p₀, Rₐ, R_v)
 Sₗ = FT(e / eₛ)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC₀]
 
 # Simulation parameters passed into ODE solver
-w = FT(0.7)                                # updraft speed
-const_dt = FT(1)                           # model timestep
+w = FT(0.7)
 t_max = FT(1200)
-aerosol = CMP.Illite(FT)
 condensation_growth = "Condensation"
 deposition_growth = "Deposition"
 DSD = "Monodisperse"
 
+# Model time step
+const_dt_range = [FT(1), FT(0.25)]
+# Every how many time steps will I sample for the "random" option
+n_dt_range_12 = [[1, 60], [4, 240]]
+# Assumed process timescale for the stochastic option
+τ_range = [1, 100]
+
+results_mean = []
+labels_mean = []
+time_mean = []
+for const_dt in const_dt_range
+    # Frostenberg_mean
+    local params = parcel_params{FT}(
+        const_dt = const_dt,
+        w = w,
+        heterogeneous = "Frostenberg_mean",
+        condensation_growth = condensation_growth,
+        deposition_growth = deposition_growth,
+        size_distribution = DSD,
+    )
+    # solve ODE
+    sol = run_parcel(IC, FT(0), t_max, params)
+    push!(time_mean, sol.t)
+    push!(results_mean, sol)
+    push!(labels_mean, "dt = " * string(const_dt))
+end
+
+results_random = []
+labels_random = []
+time_random = []
+for (const_dt, n_dt_range) in zip(const_dt_range, n_dt_range_12)
+    # Frostenberg_random with different sampling frequencies
+    for n_dt in n_dt_range
+        # creating an ensamble of solutions
+        sampling_interval = FT(n_dt * const_dt)
+        solutions_rd = []
+        for ensemble_member in range(1, length = N_ensemble)
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                w = w,
+                heterogeneous = "Frostenberg_random",
+                condensation_growth = condensation_growth,
+                deposition_growth = deposition_growth,
+                size_distribution = DSD,
+                sampling_interval = sampling_interval,
+            )
+            # solve ODE
+            sol_rd = run_parcel(IC, FT(0), t_max, params)
+            push!(solutions_rd, sol_rd)
+        end
+        mean_sol = sum(solutions_rd) / length(solutions_rd)
+        push!(time_random, solutions_rd[1].t)
+        push!(results_random, mean_sol)
+        push!(
+            labels_random,
+            "dt = " * string(const_dt) * " Δt = " * string(sampling_interval),
+        )
+    end
+end
+
+results_stoch = []
+labels_stoch = []
+time_stoch = []
+for const_dt in const_dt_range
+    # Frostenberg_stochastic with different timescales γ
+    for τ in τ_range
+        # creating an ensamble of solutions
+        γ = FT(1) / τ
+        solutions_st = []
+        for ensemble_member in range(1, length = N_ensemble)
+            local params = parcel_params{FT}(
+                const_dt = const_dt,
+                w = w,
+                heterogeneous = "Frostenberg_stochastic",
+                condensation_growth = condensation_growth,
+                deposition_growth = deposition_growth,
+                size_distribution = DSD,
+                γ = γ,
+            )
+            # solve ODE
+            sol_st = run_parcel(IC, FT(0), t_max, params)
+            push!(solutions_st, sol_st)
+        end
+        mean_sol = sum(solutions_st) / length(solutions_st)
+        push!(time_stoch, solutions_st[1].t)
+        push!(results_stoch, mean_sol)
+        push!(labels_stoch, "dt = " * string(const_dt) * " τ = " * string(τ))
+    end
+end
+
 # Plotting
-fig = MK.Figure(resolution = (900, 700))
+fig = MK.Figure(size = (1200, 900))
 ax1 = MK.Axis(fig[1, 1], ylabel = "Ice Supersaturation [-]")
-ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [K]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "T [K]")
 ax3 = MK.Axis(fig[2, 1], ylabel = "q_ice [g/kg]")
 ax4 = MK.Axis(fig[2, 2], ylabel = "q_liq [g/kg]")
-ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_liq")
-ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_ice")
-
-colors = [:blue, :green, :darkorange]
-
-function plot_results(
-    sol_t,
-    sol,
-    variable,
-    label = "",
-    unit = "",
-    linestyle = :solid,
-    color = :black,
+ax5 = MK.Axis(fig[3, 1], xlabel = "Time [min]", ylabel = "N_ice [1/cm3]")
+ax6 = MK.Axis(fig[3, 2], xlabel = "Time [min]", ylabel = "N_liq [1/cm3]")
+ax7 = MK.Axis(
+    fig[4, 1:2],
+    xlabel = "T [K]",
+    ylabel = "INPC [1/m3]",
+    yscale = log10,
 )
 
+colors = [:red, :green, :orange, :limegreen]
+colors_mean = [:black, :gray]
+
+function plot_results(sol, t, ll, cl, ls)
     MK.lines!(
         ax1,
-        sol_t / 60,
+        t / 60,
         S_i.(tps, sol[3, :], sol[1, :]) .- 1,
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
-    )
-    MK.lines!(
-        ax2,
-        sol_t / 60,
-        sol[3, :],
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
-    )
-    MK.lines!(
-        ax3,
-        sol_t / 60,
-        sol[6, :] * 1e3,
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
-    )
-    MK.lines!(
-        ax4,
-        sol_t / 60,
-        sol[5, :] * 1e3,
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
-    )
-    MK.lines!(
-        ax5,
-        sol_t / 60,
-        sol[8, :],
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
+        linestyle = ls,
+        color = cl,
     )
+    MK.lines!(ax2, t / 60, sol[3, :], linestyle = ls, color = cl)
+    MK.lines!(ax3, t / 60, sol[6, :] * 1e3, linestyle = ls, color = cl)
+    MK.lines!(ax4, t / 60, sol[5, :] * 1e3, linestyle = ls, color = cl)
+    MK.lines!(ax5, t / 60, sol[9, :] * 1e-6, linestyle = ls, color = cl)
     MK.lines!(
         ax6,
-        sol_t / 60,
-        sol[9, :],
-        label = label * string(variable) * unit,
-        linestyle = linestyle,
-        color = color,
+        t / 60,
+        sol[8, :] * 1e-6,
+        linestyle = ls,
+        color = cl,
+        label = ll,
     )
-
-    MK.axislegend(ax2, position = :lt)
 end
-
-#Frostenberg_mean
-params = parcel_params{FT}(
-    const_dt = const_dt,
-    w = w,
-    aerosol = aerosol,
-    heterogeneous = "Frostenberg_mean",
-    condensation_growth = condensation_growth,
-    deposition_growth = deposition_growth,
-    size_distribution = DSD,
-)
-# solve ODE
-sol = run_parcel(IC, FT(0), t_max, params)
-plot_results(sol.t, sol, "mean")
-
-
-# Frostenberg_random with different drawing frequencies
-drawing_interval_range = range(FT(1), stop = FT(5), length = 3)
-
-for (drawing_interval, color) in zip(drawing_interval_range, colors)
-
-    # creating an ensamble of solutions
-    solutions = []
-    for random_seed in random_seeds
-
-        RAND.seed!(trunc(Int, random_seed)) #set the random seed
-
-        params = parcel_params{FT}(
-            const_dt = const_dt,
-            w = w,
-            aerosol = aerosol,
-            heterogeneous = "Frostenberg_random",
-            condensation_growth = condensation_growth,
-            deposition_growth = deposition_growth,
-            size_distribution = DSD,
-            drawing_interval = drawing_interval,
-        )
-        # solve ODE
-        sol = run_parcel(IC, FT(0), t_max, params)
-        push!(solutions, sol)
-    end
-    mean_sol = sum(solutions) / length(solutions)
-    plot_results(sol.t, mean_sol, drawing_interval, "t_d=", " s", :dot, color)
+for (sol, t, ll, cl) in zip(results_stoch, time_stoch, labels_stoch, colors)
+    ls = :solid
+    plot_results(sol, t, ll, cl, ls)
+    MK.lines!(ax7, sol[3, :], exp.(sol[11, :]), linestyle = ls, color = cl)
 end
-
-
-# Frostenberg_stochastic with different timescales γ
-γ_range = range(FT(1), stop = FT(5), length = 3)
-
-for (γ, color) in zip(γ_range, colors)
-
-    # creating an ensamble of solutions
-    solutions = []
-    for random_seed in random_seeds
-
-        RAND.seed!(trunc(Int, random_seed)) #set the random seed
-
-        params = parcel_params{FT}(
-            const_dt = const_dt,
-            w = w,
-            aerosol = aerosol,
-            heterogeneous = "Frostenberg_stochastic",
-            condensation_growth = condensation_growth,
-            deposition_growth = deposition_growth,
-            size_distribution = DSD,
-            γ = γ,
-        )
-        # solve ODE
-        sol = run_parcel(IC, FT(0), t_max, params)
-        push!(solutions, sol)
-    end
-    mean_sol = sum(solutions) / length(solutions)
-    plot_results(sol.t, mean_sol, γ, "γ=", " s", :solid, color)
+for (sol, t, ll, cl) in zip(results_random, time_random, labels_random, colors)
+    ls = :dot
+    plot_results(sol, t, ll, cl, ls)
 end
-
+for (sol, t, ll, cl) in zip(results_mean, time_mean, labels_mean, colors_mean)
+    ls = :solid
+    plot_results(sol, t, ll, cl, ls)
+    MK.lines!(
+        ax7,
+        sol[3, :],
+        exp.(CMI_het.INP_concentration_mean.(sol[3, :])),
+        linestyle = ls,
+        color = cl,
+    )
+end
+fig[1:2, 3] = MK.Legend(fig, ax6, framevisible = false)
 MK.save("frostenberg_immersion_freezing.svg", fig)
diff --git a/parcel/Example_Immersion_Freezing.jl b/parcel/Example_Immersion_Freezing.jl
index a0e67f2ffc..6ac7897a10 100644
--- a/parcel/Example_Immersion_Freezing.jl
+++ b/parcel/Example_Immersion_Freezing.jl
@@ -24,6 +24,7 @@ qᵥ = FT(8.1e-4)
 qₗ = Nₗ * 4 / 3 * FT(π) * r₀^3 * ρₗ / FT(1.2) # 1.2 should be ρₐ
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -32,7 +33,7 @@ Rₐ = TD.gas_constant_air(tps, q)
 eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 e = eᵥ(qᵥ, p₀, Rₐ, R_v)
 Sₗ = FT(e / eₛ)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(0.4)                                # updraft speed
diff --git a/parcel/Example_Jensen_et_al_2022.jl b/parcel/Example_Jensen_et_al_2022.jl
index 19f1d0fd39..25c6add7f0 100644
--- a/parcel/Example_Jensen_et_al_2022.jl
+++ b/parcel/Example_Jensen_et_al_2022.jl
@@ -22,6 +22,7 @@ Nᵢ = FT(0)
 T₀ = FT(190)
 cᵥ₀ = FT(5 * 1e-6)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Constants
 ρₗ = wps.ρw
@@ -41,7 +42,7 @@ Sᵢ = FT(1.55)
 Sₗ = Sᵢ / ξ(tps, T₀)
 e = Sₗ * eₛ
 p₀ = e / cᵥ₀
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(1)           # updraft speed
diff --git a/parcel/Example_Liquid_only.jl b/parcel/Example_Liquid_only.jl
index e6d6556305..941277e75a 100644
--- a/parcel/Example_Liquid_only.jl
+++ b/parcel/Example_Liquid_only.jl
@@ -26,6 +26,7 @@ r₀ = FT(8e-6)
 p₀ = FT(800 * 1e2)
 T₀ = FT(273.15 + 7.0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 e = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
 Sₗ = FT(1)
 md_v = (p₀ - e) / R_d / T₀
@@ -34,7 +35,7 @@ ml_v = Nₗ * 4 / 3 * FT(π) * ρₗ * r₀^3
 qᵥ = mv_v / (md_v + mv_v + ml_v)
 qₗ = ml_v / (md_v + mv_v + ml_v)
 qᵢ = FT(0)
-IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 w = FT(10)                                 # updraft speed
diff --git a/parcel/Example_P3_ice_nuc.jl b/parcel/Example_P3_ice_nuc.jl
index 4034244726..f8c904fb72 100644
--- a/parcel/Example_P3_ice_nuc.jl
+++ b/parcel/Example_P3_ice_nuc.jl
@@ -21,6 +21,7 @@ qᵥ = FT(8.3e-4)
 qₗ = FT(Nₗ * 4 / 3 * π * rₗ^3 * wps.ρw / 1.2)
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -60,7 +61,7 @@ for it in [1, 2, 3]
     local ρₐ = TD.air_density(tps, ts)
     local eₛ = TD.saturation_vapor_pressure(tps, T₀, TD.Liquid())
     local Sₗ = FT(e / eₛ)
-    local IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+    local IC = [Sₗ, p₀, T₀, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
     #! format: off
     if mode == "P3_dep"
diff --git a/parcel/Example_P3_vs_activitybased.jl b/parcel/Example_P3_vs_activitybased.jl
index 2a281c9c18..9948d3966b 100644
--- a/parcel/Example_P3_vs_activitybased.jl
+++ b/parcel/Example_P3_vs_activitybased.jl
@@ -31,6 +31,7 @@ qᵥ = FT(8.3e-4)
 qₗ = FT(Nₗ * 4 / 3 * π * rₗ^3 * ρₗ / 1.2)
 qᵢ = FT(0)
 x_sulph = FT(0)
+ln_INPC = FT(0)
 
 # Moisture dependent initial conditions
 q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
@@ -49,9 +50,9 @@ e = eᵥ(qᵥ, p₀, Rₐ, R_v)
 Sₗ_dep = FT(e / eₛ_dep)
 Sₗ_het = FT(e / eₛ_het)
 Sₗ_hom = FT(e / eₛ_hom)
-IC_dep = [Sₗ_dep, p₀, T₀_dep, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
-IC_het = [Sₗ_het, p₀, T₀_het, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
-IC_hom = [Sₗ_hom, p₀, T₀_hom, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph]
+IC_dep = [Sₗ_dep, p₀, T₀_dep, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
+IC_het = [Sₗ_het, p₀, T₀_het, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
+IC_hom = [Sₗ_hom, p₀, T₀_hom, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 
 # Simulation parameters passed into ODE solver
 r_nuc = FT(1.25e-6)                     # assumed size of nucleated particles
@@ -164,7 +165,7 @@ for heterogeneous in heterogeneous_modes
             if aerosol == CMP.DesertDust(FT)
                 line_color = :green
                 MK.lines!(ax4, sol.t, S_i.(tps, sol[3, :], (sol[1, :])), label = heterogeneous, color = line_color) # saturation
-                MK.lines!(ax5, sol.t, sol[9, :] .* 1e-6, color = line_color)                                        # ICNC 
+                MK.lines!(ax5, sol.t, sol[9, :] .* 1e-6, color = line_color)                                        # ICNC
                 MK.lines!(ax6, sol.t, sol[3, :], color = line_color)                                                # Temperature
             elseif aerosol == CMP.Illite(FT)
                 line_color = :orange
diff --git a/parcel/Example_Tully_et_al_2023.jl b/parcel/Example_Tully_et_al_2023.jl
index 7915c5ef65..8720790440 100644
--- a/parcel/Example_Tully_et_al_2023.jl
+++ b/parcel/Example_Tully_et_al_2023.jl
@@ -13,24 +13,25 @@ include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
 """
 function get_initial_condition(
     tps,
-    p_a,
+    p_air,
     T,
-    q_vap,
-    q_liq,
-    q_ice,
-    N_aer,
-    N_liq,
-    N_ice,
+    qᵥ,
+    qₗ,
+    qᵢ,
+    Nₐ,
+    Nₗ,
+    Nᵢ,
     x_sulph,
+    ln_INPC,
 )
-    q = TD.PhasePartition(q_vap + q_liq + q_ice, q_liq, q_ice)
+    q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
     R_a = TD.gas_constant_air(tps, q)
     R_v = TD.Parameters.R_v(tps)
     e_sl = TD.saturation_vapor_pressure(tps, T, TD.Liquid())
-    e = eᵥ(q_vap, p_a, R_a, R_v)
-    S_liq = e / e_sl
+    e = eᵥ(qᵥ, p_air, R_a, R_v)
+    Sₗ = e / e_sl
 
-    return [S_liq, p_a, T, q_vap, q_liq, q_ice, N_aer, N_liq, N_ice, x_sulph]
+    return [Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₐ, Nₗ, Nᵢ, x_sulph, ln_INPC]
 end
 
 """
@@ -55,6 +56,7 @@ function Tully_et_al_2023(FT)
     q_liq_0 = FT(0)
     q_ice_0 = FT(0)
     x_sulph = FT(0)
+    ln_INPC = FT(0)
     # Initial conditions for the 2nd period
     T2 = FT(229.25)
     q_vap2 = FT(3.3e-4)
@@ -106,6 +108,7 @@ function Tully_et_al_2023(FT)
             N_droplets,
             N_0,
             x_sulph,
+            ln_INPC,
         )
         sol1 = run_parcel(IC1, 0, t_max, params)
         if mode == "MohlerAF"
@@ -124,6 +127,7 @@ function Tully_et_al_2023(FT)
                 sol1[8, end],
                 sol1[9, end],
                 x_sulph,
+                ln_INPC,
             )
             sol2 = run_parcel(IC2, sol1.t[end], sol1.t[end] + t_max, params)
 
@@ -142,6 +146,7 @@ function Tully_et_al_2023(FT)
                 sol2[8, end],
                 sol2[9, end],
                 x_sulph,
+                ln_INPC,
             )
             sol3 = run_parcel(IC3, sol2.t[end], sol2.t[end] + t_max, params)
 
diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl
index 81bc2797f5..5b369fcff9 100644
--- a/parcel/ParcelModel.jl
+++ b/parcel/ParcelModel.jl
@@ -23,7 +23,7 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     w = FT(1)
     r_nuc = FT(0.5 * 1.e-4 * 1e-6)
     A_aer = FT(1e-9)
-    drawing_interval = FT(1)
+    sampling_interval = FT(1)
     γ = FT(1)
     ip = CMP.Frostenberg2023(FT)
 end
@@ -54,7 +54,8 @@ function parcel_model(dY, Y, p, t)
         Nₐ = clip!(Y[7]),
         Nₗ = clip!(Y[8]),
         Nᵢ = clip!(Y[9]),
-        xS = Y[10],
+        xS = Y[10],        # TODO - to be deleted
+        ln_INPC = Y[11],   # needed only in stochastic Frostenberg
         t = t,
     )
 
@@ -94,6 +95,8 @@ function parcel_model(dY, Y, p, t)
     dqᵢ_dt_dep = dNᵢ_dt_dep * 4 / 3 * FT(π) * r_nuc^3 * ρᵢ / ρ_air
 
     # Heterogeneous ice nucleation
+    #@info(imm_params)
+    dln_INPC_imm = INPC_model(imm_params, state)
     dNᵢ_dt_imm = immersion_freezing(imm_params, PSD, state)
     dqᵢ_dt_imm = dNᵢ_dt_imm * PSD.Vₗ * ρᵢ / ρ_air
 
@@ -143,6 +146,7 @@ function parcel_model(dY, Y, p, t)
     dY[8] = dNₗ_dt      # mumber concentration of droplets
     dY[9] = dNᵢ_dt      # number concentration of activated particles
     dY[10] = FT(0)      # sulphuric acid concentration
+    dY[11] = dln_INPC_imm
 end
 
 """
@@ -183,8 +187,8 @@ The parcel parameters struct comes with default values that can be overwritten:
  - r_nuc - assumed size of nucleating ice crystals. Default value is 5e-11 [m]
  - A_aer - assumed surface area of ice nucleating aerosol. Default value assumes radius of 5e-11 [m]
  - ip - parameters of INPC frequency distribution for Frostenberg parametrization of immerison freezing
- - drawing_interval - time in seconds between random draws in Frostenberg_random parametrization of immerison freezing
- - γ - timescale for the stochastic process in Frostenberg_stochastic parametrization of immerison freezing
+ - sampling_interval - number of time steps between random draws in Frostenberg_random parametrization of immerison freezing
+ - γ - inverse timescale for the stochastic process in Frostenberg_stochastic parametrization of immerison freezing
 """
 function run_parcel(IC, t_0, t_end, pp)
 
@@ -225,11 +229,12 @@ function run_parcel(IC, t_0, t_end, pp)
     elseif pp.heterogeneous == "P3_het"
         imm_params = P3_het{FT}(pp.ips, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_random"
-        imm_params = Frostenberg_random{FT}(pp.ip, pp.drawing_interval)
+        imm_params =
+            Frostenberg_random{FT}(pp.ip, pp.sampling_interval, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_mean"
-        imm_params = Frostenberg_mean{FT}(pp.ip)
+        imm_params = Frostenberg_mean{FT}(pp.ip, pp.const_dt)
     elseif pp.heterogeneous == "Frostenberg_stochastic"
-        imm_params = Frostenberg_stochastic{FT}(pp.ip, pp.γ)
+        imm_params = Frostenberg_stochastic{FT}(pp.ip, pp.γ, pp.const_dt)
     else
         throw("Unrecognized heterogeneous mode")
     end
diff --git a/parcel/ParcelParameters.jl b/parcel/ParcelParameters.jl
index ecb80dae16..8875f25405 100644
--- a/parcel/ParcelParameters.jl
+++ b/parcel/ParcelParameters.jl
@@ -41,16 +41,19 @@ end
 
 struct Frostenberg_random{FT} <: CMP.ParametersType{FT}
     ip::CMP.ParametersType{FT}
-    drawing_interval::FT
+    sampling_interval::FT
+    const_dt::FT
 end
 
 struct Frostenberg_stochastic{FT} <: CMP.ParametersType{FT}
     ip::CMP.ParametersType{FT}
     γ::FT
+    const_dt::FT
 end
 
 struct Frostenberg_mean{FT} <: CMP.ParametersType{FT}
     ip::CMP.ParametersType{FT}
+    const_dt::FT
 end
 
 struct ABHOM{FT} <: CMP.ParametersType{FT}
diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl
index 06463d6a6c..21d147df12 100644
--- a/parcel/ParcelTendencies.jl
+++ b/parcel/ParcelTendencies.jl
@@ -109,9 +109,9 @@ end
 
 function immersion_freezing(params::Frostenberg_random, PSD, state)
     FT = eltype(state)
-    (; ip, drawing_interval) = params
+    (; ip, sampling_interval, const_dt) = params
     (; T, Nₗ, Nᵢ, t) = state
-    if mod(t, drawing_interval) == 0
+    if mod(t, sampling_interval) == 0
         μ = CMI_het.INP_concentration_mean(T)
         INPC = exp(rand(DS.Normal(μ, ip.σ)))
     else
@@ -122,22 +122,37 @@ end
 
 function immersion_freezing(params::Frostenberg_mean, PSD, state)
     FT = eltype(state)
-    (; ip) = params
+    (; ip, const_dt) = params
     (; T, Nₗ, Nᵢ) = state
     INPC = exp(CMI_het.INP_concentration_mean(T))
     return min(Nₗ, max(FT(0), INPC - Nᵢ)) / const_dt
 end
 
-function immersion_freezing(params::Frostenberg_stochastic, PSD, state)
+function INPC_model(params, state)
     FT = eltype(state)
-    (; ip, γ) = params
-    (; T, Nₗ, Nᵢ, t) = state
+    return FT(0)
+end
+
+function INPC_model(params::Frostenberg_stochastic, state)
+    FT = eltype(state)
+    (; ip, γ, const_dt) = params
+    (; T, ln_INPC, t) = state
+
     μ = CMI_het.INP_concentration_mean(T)
     g = ip.σ * sqrt(2 * γ)
-    mean = (1 - exp(-γ * t)) * μ
-    st_dev = sqrt(g^2 / (2 * γ) * (1 - exp(-2 * γ * t)))
-    INPC = exp(rand(DS.Normal(mean, st_dev)))
-    return min(Nₗ, max(FT(0), INPC - Nᵢ)) / const_dt
+
+    dln_INPC =
+        -γ * (ln_INPC - μ) * const_dt +
+        g * sqrt(const_dt) * rand(DS.Normal(0, 1))
+
+    return dln_INPC / const_dt
+end
+
+function immersion_freezing(params::Frostenberg_stochastic, PSD, state)
+    FT = eltype(state)
+    (; ip, γ, const_dt) = params
+    (; T, Nₗ, Nᵢ, ln_INPC, t) = state
+    return min(Nₗ, max(FT(0), exp(ln_INPC) - Nᵢ)) / const_dt
 end
 
 function homogeneous_freezing(::Empty, PSD, state)