diff --git a/docs/bibliography.bib b/docs/bibliography.bib
index 28220fe244..fb88ca5ebc 100644
--- a/docs/bibliography.bib
+++ b/docs/bibliography.bib
@@ -143,6 +143,18 @@ @article{Dunne2016
   doi = {10.1126/science.aaf2649}
 }
 
+@Article{Frostenberg2023,
+author = {Frostenberg, H. C. and Welti, A. and Luhr, M. and Savre, J. and Thomson, E. S. and Ickes, L.},
+title = {The chance of freezing -- a conceptional study to parameterize temperature-dependent freezing by including randomness of ice-nucleating particle concentrations},
+journal = {Atmospheric Chemistry and Physics},
+volume = {23},
+year = {2023},
+number = {19},
+pages = {10883--10900},
+url = {https://acp.copernicus.org/articles/23/10883/2023/},
+doi = {10.5194/acp-23-10883-2023}
+}
+
 @article{Glassmeier2016,
   author = {Franziska Glassmeier and Ulrike Lohmann},
   title = {Constraining Precipitation Susceptibility of Warm-, Ice-, and Mixed-Phase Clouds with Microphysical Equations},
diff --git a/docs/src/API.md b/docs/src/API.md
index 62cba66ec1..6b7e209f92 100644
--- a/docs/src/API.md
+++ b/docs/src/API.md
@@ -92,6 +92,7 @@ HetIceNucleation.deposition_J
 HetIceNucleation.ABIFM_J
 HetIceNucleation.P3_deposition_N_i
 HetIceNucleation.P3_het_N_i
+HetIceNucleation.INP_concentration_frequency
 ```
 
 # Homogeneous ice nucleation
@@ -138,6 +139,7 @@ Parameters.Seasalt
 Parameters.Sulfate
 Parameters.AerosolActivationParameters
 Parameters.IceNucleationParameters
+Parameters.Frostenberg2023
 Parameters.H2SO4SolutionParameters
 Parameters.Mohler2006
 Parameters.Koop2000
diff --git a/docs/src/IceNucleation.md b/docs/src/IceNucleation.md
index 9315e8de34..65eee37649 100644
--- a/docs/src/IceNucleation.md
+++ b/docs/src/IceNucleation.md
@@ -6,7 +6,9 @@ The `IceNucleation.jl` module includes:
   - the parameterization of activation of dust aerosol particles into ice crystals
     via deposition of water vapor,
   - water activity based parameterization of immersion freezing,
-  - water activity based parameterization of homogeneous freezing.
+  - water activity based parameterization of homogeneous freezing,
+  - parametrization of temperature-dependent aerosol-independent ice nucleating particles concentration for immersion freezing.
+
 The parameterization for deposition on dust particles is an implementation of
   the empirical formulae from [Mohler2006](@cite)
   and is valid for two types of dust particles:
@@ -14,6 +16,7 @@ The parameterization for deposition on dust particles is an implementation of
   The parameterization for immersion freezing is an implementation of [KnopfAlpert2013](@cite)
   and is valid for droplets containing sulphuric acid.
   The parameterization for homogeneous freezing is an implementation of [Koop2000](@cite).
+  The parametrization for ice nucleating particles concentration for immersion freezing is an implementation of [Frostenberg2023](@cite).
 
 !!! note
 
@@ -239,3 +242,37 @@ Multiple sulphuric acid concentrations, ``x``,
     The current parameterization in CloudMicrophysics.jl is not valid for \Delta a_w
     values that are obtained from pure water droplets. Though CliMA lines look far
     from the Spichtinger 2023 line, the lines seem to move closer as x approaches 0.
+
+
+## INP Concentration Frequency
+
+With this parametrization, the concentration of ice nucleating particles (INPs) is found based on the relative frequency distribution, which depends on the temperature, but does not depend on aerosol types. It is based on [Frostenberg2023](@cite) and is derived from measurements, for marine
+data sets. It is a lognormal distribution, described by
+
+```math
+\begin{equation}
+  D(\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi} \sigma} exp \left(- \frac{[ln(a \cdot INPC) - \mu(T)]^2}{2 \sigma^2} \right),
+\end{equation}
+```
+where ``T`` is the temperature in degrees Celsius, ``INPC`` is the INP concentration in m``^{-3}`` and the temperature-dependent mean ``\mu(T)`` is given by
+```math
+\begin{equation}
+  \mu(T) = ln(-(b \cdot T)^9 \times 10^{-9})
+\end{equation}
+```
+ ``\sigma^2`` is the variance, ``a`` and ``b`` are coefficients. The parameters defined in [Frostenberg2023](@cite) for marine data sets are ``\sigma=1.37``, ``a=1`` m``^3``, ``b=1``/C.
+
+!!! note
+
+    Our implementation uses base SI units and takes ``T`` in Kelvin.
+
+The following plot shows the relative frequency distribution for INPCs, as a function of temperature (the same as figure 1 in [Frostenberg2023](@cite)). 
+
+```@example
+include("plots/Frostenberg_fig1.jl")
+```
+![](Frostenberg_fig1.svg)
+
+The following plot shows the relative frequency distribution for INPCs at the temperature ``T=-16 C``.
+
+![](Frostenberg_fig1_T16.svg)
\ No newline at end of file
diff --git a/docs/src/plots/Frostenberg_fig1.jl b/docs/src/plots/Frostenberg_fig1.jl
new file mode 100644
index 0000000000..3c47955f2a
--- /dev/null
+++ b/docs/src/plots/Frostenberg_fig1.jl
@@ -0,0 +1,72 @@
+import Plots as PL
+
+import CloudMicrophysics as CM
+import CloudMicrophysics.HetIceNucleation as IN
+import CLIMAParameters as CP
+import CloudMicrophysics.Parameters as CMP
+
+FT = Float32
+ip = CMP.Frostenberg2023(FT)
+
+T_range_celsius = range(-40, stop = -2, length = 500)  # air temperature [°C]
+T_range_kelvin = T_range_celsius .+ 273.15  # air temperature [K]
+INPC_range = 10.0 .^ (range(-5, stop = 7, length = 500)) #ice nucleating particle concentration
+
+frequency = [
+    IN.INP_concentration_frequency(ip, INPC, T) > 0.015 ?
+    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]
+
+
+PL.contourf(
+    T_range_kelvin,
+    INPC_range,
+    frequency,
+    xlabel = "T [K]",
+    ylabel = "INPC [m⁻³]",
+    colorbar_title = "frequency",
+    yaxis = :log,
+    color = :lighttest,
+    gridlinewidth = 3,
+    ylims = (1e-5, 1e7),
+    lw = 0,
+)
+
+PL.plot!(
+    T_range_kelvin,
+    mu,
+    label = "median INPC",
+    legend = :topright,
+    color = :darkred,
+)
+
+PL.plot!(
+    repeat([257], 50),
+    10.0 .^ (range(-2, stop = 6, length = 50)),
+    label = "T = 257 K",
+    color = :darkgreen,
+    linestyle = :dash,
+)
+PL.savefig("Frostenberg_fig1.svg")
+
+
+#plotting the distribution for T=257 K
+
+T = FT(257.0) # [K]
+INPC_range = FT(10) .^ range(FT(-1), 4, length = 100)
+frequency = [IN.INP_concentration_frequency(ip, INPC, T) for INPC in INPC_range]
+
+PL.plot(
+    INPC_range,
+    frequency,
+    xlabel = "INPC [m⁻³]",
+    ylabel = "frequency",
+    xaxis = :log,
+    color = :darkgreen,
+    label = "T = 257 K",
+    legend = :topright,
+)
+
+PL.savefig("Frostenberg_fig1_T16.svg")
diff --git a/src/IceNucleation.jl b/src/IceNucleation.jl
index d70b262344..d896c1d1a7 100644
--- a/src/IceNucleation.jl
+++ b/src/IceNucleation.jl
@@ -12,6 +12,8 @@ export deposition_J
 export ABIFM_J
 export P3_deposition_N_i
 export P3_het_N_i
+export INP_concentration_frequency
+
 
 """
     dust_activated_number_fraction(dust, ip, Si, T)
@@ -155,6 +157,31 @@ function P3_het_N_i(T::FT, N_l::FT, V_l::FT, Δt::FT) where {FT}
     return N_l * (1 - exp(-B * V_l_converted * Δt * exp(a * Tₛ)))
 end
 
+"""
+    INP_concentration_frequency(INPC,T)
+
+ - `params` - a struct with INPC(T) distribution parameters
+ - `INPC` - concentration of ice nucleating particles [m^-3]
+ - `T` - air temperature [K]
+
+Returns the relative frequency of a given INP concentration,
+depending on the temperature.
+Based on Frostenberg et al., 2023. See DOI: 10.5194/acp-23-10883-2023
+"""
+function INP_concentration_frequency(
+    params::CMP.Frostenberg2023,
+    INPC::FT,
+    T::FT,
+) where {FT}
+
+    T_celsius = T - 273.15
+    (; σ, a, b) = params
+
+    μ = log(-(b * T_celsius)^9 * 10^(-9))
+
+    return 1 / (sqrt(2 * FT(π)) * σ) * exp(-(log(a * INPC) - μ)^2 / (2 * σ^2))
+end
+
 end # end module
 
 """
diff --git a/src/parameters/IceNucleation.jl b/src/parameters/IceNucleation.jl
index 9f313738ff..1923957afb 100644
--- a/src/parameters/IceNucleation.jl
+++ b/src/parameters/IceNucleation.jl
@@ -1,4 +1,5 @@
 export IceNucleationParameters
+export Frostenberg2023
 
 """
     Mohler2006{FT}
@@ -87,3 +88,36 @@ function IceNucleationParameters(
     HOM = typeof(homogeneous)
     return IceNucleationParameters{FT, DEP, HOM}(deposition, homogeneous)
 end
+
+
+"""
+    Frostenberg2023{FT}
+
+Parameters for frequency distribution of INP concentration
+DOI: 10.5194/acp-23-10883-2023
+
+# Fields
+$(DocStringExtensions.FIELDS)
+"""
+Base.@kwdef struct Frostenberg2023{FT} <: ParametersType{FT}
+    "standard deviation"
+    σ::FT
+    "coefficient"
+    a::FT
+    "coefficient"
+    b::FT
+end
+
+Frostenberg2023(::Type{FT}) where {FT <: AbstractFloat} =
+    Frostenberg2023(CP.create_toml_dict(FT))
+
+function Frostenberg2023(td::CP.AbstractTOMLDict)
+    name_map = (;
+        :Frostenberg2023_standard_deviation => :σ,
+        :Frostenberg2023_a_coefficient => :a,
+        :Frostenberg2023_b_coefficient => :b,
+    )
+    parameters = CP.get_parameter_values(td, name_map, "CloudMicrophysics")
+    FT = CP.float_type(td)
+    return Frostenberg2023{FT}(; parameters...)
+end
diff --git a/test/gpu_tests.jl b/test/gpu_tests.jl
index a46fde12b0..65b9ff65a9 100644
--- a/test/gpu_tests.jl
+++ b/test/gpu_tests.jl
@@ -595,6 +595,20 @@ end
     end
 end
 
+@kernel function IceNucleation_INPC_frequency_kernel!(
+    output::AbstractArray{FT},
+    ip_frostenberg,
+    INPC,
+    T,
+) where {FT}
+
+    i = @index(Group, Linear)
+
+    @inbounds begin
+        output[1] = CMI_het.INP_concentration_frequency(ip_frostenberg, INPC, T)
+    end
+end
+
 @kernel function IceNucleation_homogeneous_J_kernel!(
     ip,
     output::AbstractArray{FT},
@@ -681,6 +695,7 @@ function test_gpu(FT)
     feldspar = CMP.Feldspar(FT)
     ferrihydrite = CMP.Ferrihydrite(FT)
     ip = CMP.IceNucleationParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
 
     @testset "Aerosol activation kernels" begin
         dims = (6, 2)
@@ -1026,6 +1041,18 @@ function test_gpu(FT)
         dims = (1, 1)
         (; output, ndrange) = setup_output(dims, FT)
 
+        T = FT(233)
+        INPC = FT(220000)
+
+        kernel! = IceNucleation_INPC_frequency_kernel!(backend, work_groups)
+        kernel!(output, ip_frostenberg, INPC, T; ndrange)
+
+        # test INPC_frequency is callable and returns a reasonable value
+        @test Array(output)[1] ≈ FT(0.26) rtol = 0.1
+
+        dims = (1, 1)
+        (; output, ndrange) = setup_output(dims, FT)
+
         T = ArrayType([FT(220)])
         x_sulph = ArrayType([FT(0.15)])
         Delta_a_w = ArrayType([FT(0.2907389666103033)])
diff --git a/test/heterogeneous_ice_nucleation_tests.jl b/test/heterogeneous_ice_nucleation_tests.jl
index 84aab00405..f09de630c8 100644
--- a/test/heterogeneous_ice_nucleation_tests.jl
+++ b/test/heterogeneous_ice_nucleation_tests.jl
@@ -16,6 +16,7 @@ function test_heterogeneous_ice_nucleation(FT)
     tps = TD.Parameters.ThermodynamicsParameters(FT)
     H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
     ip = CMP.IceNucleationParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
     # more parameters for aerosol properties
     ATD = CMP.ArizonaTestDust(FT)
     desert_dust = CMP.DesertDust(FT)
@@ -207,6 +208,21 @@ function test_heterogeneous_ice_nucleation(FT)
         TT.@test CMI_het.P3_het_N_i(T_cold, N_liq, V_l, Δt) >
                  CMI_het.P3_het_N_i(T_warm, N_liq, V_l, Δt)
     end
+
+    TT.@testset "Frostenberg" begin
+
+        temperatures = FT.([233, 257])
+        INPCs = FT.([220000, 9])
+        frequencies = FT.([0.26, 0.08])
+
+        for (T, INPC, frequency) in zip(temperatures, INPCs, frequencies)
+            TT.@test CMI_het.INP_concentration_frequency(
+                ip_frostenberg,
+                INPC,
+                T,
+            ) ≈ frequency rtol = 0.1
+        end
+    end
 end
 
 println("Testing Float64")
diff --git a/test/performance_tests.jl b/test/performance_tests.jl
index b18c5b8a6d..b677164236 100644
--- a/test/performance_tests.jl
+++ b/test/performance_tests.jl
@@ -63,6 +63,7 @@ function benchmark_test(FT)
     kaolinite = CMP.Kaolinite(FT)
     ip = CMP.IceNucleationParameters(FT)
     H2SO4_prs = CMP.H2SO4SolutionParameters(FT)
+    ip_frostenberg = CMP.Frostenberg2023(FT)
     # aerosol nucleation parameters
     mixed_nuc = CMP.MixedNucleationParameters(FT)
     h2so4_nuc = CMP.H2S04NucleationParameters(FT)
@@ -116,6 +117,8 @@ function benchmark_test(FT)
     V_liq = FT(4 / 3 * π * r_liq^3)
     Δt = FT(25)
 
+    INPC = FT(1e5)
+
     # P3 scheme
     bench_press(P3.thresholds, (p3, ρ_r, F_r), 12e6, 2048, 80)
 
@@ -151,6 +154,11 @@ function benchmark_test(FT)
     bench_press(CMI_het.P3_deposition_N_i, (T_air_cold), 230)
     bench_press(CMI_het.ABIFM_J, (desert_dust, Delta_a_w), 230)
     bench_press(CMI_het.P3_het_N_i, (T_air_cold, N_liq, V_liq, Δt), 230)
+    bench_press(
+        CMI_het.INP_concentration_frequency,
+        (ip_frostenberg, INPC, T_air_cold),
+        110,
+    )
     bench_press(CMI_hom.homogeneous_J, (ip.homogeneous, Delta_a_w), 230)
 
     # non-equilibrium