diff --git a/docs/src/plots/Frostenberg_fig1.jl b/docs/src/plots/Frostenberg_fig1.jl
index 3c47955f2a..589b6e9887 100644
--- a/docs/src/plots/Frostenberg_fig1.jl
+++ b/docs/src/plots/Frostenberg_fig1.jl
@@ -17,7 +17,7 @@ frequency = [
     IN.INP_concentration_frequency(ip, INPC, T) : missing for
     INPC in INPC_range, T in T_range_kelvin
 ]
-mu = [exp(log(-(ip.b * T)^9 * 10^(-9))) for T in T_range_celsius]
+mu = [exp(log(-(T)^9 * 10^(-9))) for T in T_range_celsius]
 
 
 PL.contourf(
diff --git a/parcel/Example_Frostenberg_Immersion_Freezing.jl b/parcel/Example_Frostenberg_Immersion_Freezing.jl
new file mode 100644
index 0000000000..766eea1679
--- /dev/null
+++ b/parcel/Example_Frostenberg_Immersion_Freezing.jl
@@ -0,0 +1,87 @@
+import OrdinaryDiffEq as ODE
+import CairoMakie as MK
+import Thermodynamics as TD
+import CloudMicrophysics as CM
+import CLIMAParameters as CP
+
+# definition of the ODE problem for parcel model
+include(joinpath(pkgdir(CM), "parcel", "Parcel.jl"))
+FT = Float32
+# get free parameters
+tps = TD.Parameters.ThermodynamicsParameters(FT)
+wps = CMP.WaterProperties(FT)
+
+# Initial conditions
+ρₗ = wps.ρw
+Nₐ = FT(0)
+Nₗ = FT(500 * 1e3)
+Nᵢ = FT(0)
+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ᵢ = FT(0)
+x_sulph = FT(0.01)
+
+# Moisture dependent initial conditions
+q = TD.PhasePartition.(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
+R_v = TD.Parameters.R_v(tps)
+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]
+
+# Simulation parameters passed into ODE solver
+w = FT(0.7)                                # updraft speed
+const_dt = FT(1)                           # model timestep
+t_max = FT(1200)
+aerosol = CMP.Illite(FT)
+heterogeneous = "Frostenberg_random"
+condensation_growth = "Condensation"
+deposition_growth = "Deposition"
+DSD = "Monodisperse"
+
+# Plotting
+fig = MK.Figure(resolution = (800, 600))
+ax1 = MK.Axis(fig[1, 1], ylabel = "Supersaturation [-]")
+ax2 = MK.Axis(fig[1, 2], ylabel = "Temperature [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")
+
+params = parcel_params{FT}(
+    const_dt = const_dt,
+    w = w,
+    aerosol = aerosol,
+    heterogeneous = heterogeneous,
+    condensation_growth = condensation_growth,
+    deposition_growth = deposition_growth,
+    size_distribution = DSD,
+)
+# solve ODE
+sol = run_parcel(IC, FT(0), t_max, params)
+
+# Plot results
+MK.lines!(ax1, sol.t / 60, sol[1, :] .-1, label = "S_liq")
+MK.lines!(ax1, sol.t / 60, S_i.(tps, sol[3, :], sol[1, :]) .- 1, label = "S_ice")
+MK.lines!(ax2, sol.t / 60, sol[3, :])
+MK.lines!(ax3, sol.t / 60, sol[6, :] * 1e3)
+MK.lines!(ax4, sol.t / 60, sol[5, :] * 1e3)
+sol_Nₗ = sol[8, :]
+sol_Nᵢ = sol[9, :]
+MK.lines!(ax5, sol.t / 60, sol_Nₗ)
+MK.lines!(ax6, sol.t / 60, sol_Nᵢ)
+
+MK.axislegend(
+    ax1,
+    framevisible = false,
+    labelsize = 12,
+    orientation = :horizontal,
+    nbanks = 2,
+    position = :rt,
+)
+
+MK.save("frostenberg_immersion_freezing.svg", fig)
diff --git a/parcel/ParcelModel.jl b/parcel/ParcelModel.jl
index 61a1e35226..abcdd27fd7 100644
--- a/parcel/ParcelModel.jl
+++ b/parcel/ParcelModel.jl
@@ -22,6 +22,9 @@ Base.@kwdef struct parcel_params{FT} <: CMP.ParametersType{FT}
     const_dt = 1
     w = FT(1)
     r_nuc = FT(0.5 * 1.e-4 * 1e-6)
+    drawing_interval = FT(1)
+    γ = FT(1)
+    ip = CMP.Frostenberg2023(FT)
 end
 
 """
@@ -51,6 +54,7 @@ function parcel_model(dY, Y, p, t)
         Nₗ = clip!(Y[8]),
         Nᵢ = clip!(Y[9]),
         xS = Y[10],
+        t = t,
     )
 
     # Constants
@@ -60,7 +64,7 @@ function parcel_model(dY, Y, p, t)
     ρₗ = wps.ρw
 
     # Get the state values
-    (; Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₗ, Nᵢ) = state
+    (; Sₗ, p_air, T, qᵥ, qₗ, qᵢ, Nₗ, Nᵢ, t) = state
     # Get thermodynamic parameters, phase partition and create thermo state.
     q = TD.PhasePartition(qᵥ + qₗ + qᵢ, qₗ, qᵢ)
     ts = TD.PhaseNonEquil_pTq(tps, p_air, T, q)
@@ -162,6 +166,7 @@ Parcel state vector (all variables are in base SI units):
  - Nₗ    - number concentration of existing water droplets
  - Nᵢ    - concentration of activated ice crystals
  - xS    - percent mass sulphuric acid
+ - t     - time in seconds since the beginning of the simulation 
 
 The parcel parameters struct comes with default values that can be overwritten:
  - deposition - string with deposition ice nucleation parameterization choice ["None" (default), "MohlerAF", "MohlerRate", "ActivityBased", "P3_dep"]
@@ -175,6 +180,9 @@ The parcel parameters struct comes with default values that can be overwritten:
  - const_dt - model timestep [s]. Parcel model is using a simple Euler timestepper. Default value is 1 s
  - w - parcel vertical velocity [m/s]. Default value is 1 m/s
  - r_nuc - assumed size of nucleating ice crystals. Default value is 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_stichastic parametrization of immerison freezing
 """
 function run_parcel(IC, t_0, t_end, pp)
 
@@ -214,6 +222,12 @@ function run_parcel(IC, t_0, t_end, pp)
         imm_params = ABIFM{FT}(pp.H₂SO₄ps, pp.tps, pp.aerosol)
     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)
+    elseif pp.heterogeneous == "Frostenberg_mean"
+        imm_params = Frostenberg_mean{FT}(pp.ip)
+    elseif pp.heterogeneous == "Frostenberg_stochastic"
+        imm_params = Frostenberg_stochastic{FT}(pp.ip, pp.γ)
     else
         throw("Unrecognized heterogeneous mode")
     end
diff --git a/parcel/ParcelParameters.jl b/parcel/ParcelParameters.jl
index ce3abf8792..0644c19356 100644
--- a/parcel/ParcelParameters.jl
+++ b/parcel/ParcelParameters.jl
@@ -38,6 +38,20 @@ struct P3_het{FT} <: CMP.ParametersType{FT}
     const_dt::FT
 end
 
+struct Frostenberg_random{FT} <: CMP.ParametersType{FT}
+    ip::CMP.ParametersType{FT}
+    drawing_interval::FT
+end
+
+struct Frostenberg_stochastic{FT} <: CMP.ParametersType{FT}
+    ip::CMP.ParametersType{FT}
+    γ::FT
+end
+
+struct Frostenberg_mean{FT} <: CMP.ParametersType{FT} 
+    ip::CMP.ParametersType{FT}
+end
+
 struct ABHOM{FT} <: CMP.ParametersType{FT}
     tps::TDP.ThermodynamicsParameters{FT}
     ips::CMP.ParametersType{FT}
diff --git a/parcel/ParcelTendencies.jl b/parcel/ParcelTendencies.jl
index e215da7dbe..192812a48d 100644
--- a/parcel/ParcelTendencies.jl
+++ b/parcel/ParcelTendencies.jl
@@ -3,6 +3,7 @@ import CloudMicrophysics.Common as CMO
 import CloudMicrophysics.HetIceNucleation as CMI_het
 import CloudMicrophysics.HomIceNucleation as CMI_hom
 import CloudMicrophysics.Parameters as CMP
+using Distributions
 
 function deposition_nucleation(::Empty, state, dY)
     FT = eltype(state)
@@ -106,6 +107,37 @@ function immersion_freezing(params::P3_het, PSD, state)
     return max(FT(0), (Nᵢ_het - Nᵢ) / const_dt)
 end
 
+function immersion_freezing(params::Frostenberg_random, PSD, state)
+    FT = eltype(state)
+    (; ip, drawing_interval) = params
+    (; T, Nₗ, Nᵢ, t) = state
+    if mod(t, drawing_interval) == 0
+        μ = CMI_het.INP_concentration_mean(T)
+        INPC = exp(rand(Normal(μ, ip.σ)))
+    end
+    return min(Nₗ, max(FT(0), INPC - Nᵢ))
+end
+
+function immersion_freezing(params::Frostenberg_mean, PSD, state)
+    FT = eltype(state)
+    (; ip) = params
+    (; T, Nₗ, Nᵢ) = state
+    INPC = exp(CMI_het.INP_concentration_mean(T))
+    return min(Nₗ, max(FT(0), INPC - Nᵢ))
+end
+
+function immersion_freezing(params::Frostenberg_stochastic, PSD, state)
+    FT = eltype(state)
+    (; ip, γ) = params
+    (; T, Nₗ, Nᵢ, 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(Normal(mean, st_dev)))
+    return min(Nₗ, max(FT(0), INPC - Nᵢ))
+end
+
 function homogeneous_freezing(::Empty, PSD, state)
     FT = eltype(state)
     return FT(0)
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
index f702515b14..6db272c84f 100644
--- a/src/IceNucleation.jl
+++ b/src/IceNucleation.jl
@@ -13,6 +13,7 @@ export ABIFM_J
 export P3_deposition_N_i
 export P3_het_N_i
 export INP_concentration_frequency
+export INP_concentration_mean
 
 """
     dust_activated_number_fraction(dust, ip, Si, T)
@@ -165,7 +166,7 @@ function P3_het_N_i(
 end
 
 """
-    INP_concentration_frequency(INPC,T)
+    INP_concentration_frequency(params,INPC,T)
 
  - `params` - a struct with INPC(T) distribution parameters
  - `INPC` - concentration of ice nucleating particles [m^-3]
@@ -181,12 +182,27 @@ function INP_concentration_frequency(
     T::FT,
 ) where {FT}
 
-    T_celsius = T - 273.15
-    (; σ, a, b) = params
+    μ = INP_concentration_mean(T)
 
-    μ = log(-(b * T_celsius)^9 * 10^(-9))
+    return 1 / (sqrt(2 * FT(π)) * params.σ) * exp(-(log(INPC) - μ)^2 / (2 * params.σ^2))
+end
+
+
+"""
+    INP_concentration_mean(T)
+
+ - `T` - air temperature [K]
+
+Returns the logarithm of mean INP concentration (in m^-3), depending on the temperature.
+Based on the function μ(T) in Frostenberg et al., 2023. See DOI: 10.5194/acp-23-10883-2023
+"""
+function INP_concentration_mean(
+    T::FT,
+) where {FT}
 
-    return 1 / (sqrt(2 * FT(π)) * σ) * exp(-(log(a * INPC) - μ)^2 / (2 * σ^2))
+    T_celsius = T - FT(273.15)
+    
+    return log(-(T_celsius)^9 * 10^(-9))
 end
 
 end # end module