Skip to content

Ice thermal conductivity kinetic correction fix and example revision #1923

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
tluettm wants to merge 7 commits into
base: main
Choose a base branch
from
Open

Ice thermal conductivity kinetic correction fix and example revision #1923

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
a5b0b39
Fixed bug for thermal conductivity kinetic correction
tluettm Aug 7, 2026
2c74a34
name change for example
tluettm Aug 7, 2026
44ce3b0
Updated fig_1_2_3.ipynb and fig_S1_S2.ipynb and plot scripts
tluettm Aug 7, 2026
993ffe5
add a regression test checking for the 4-factor, remove empty lines i…
slayoo Aug 10, 2026
7e07d11
rerun of fig_4_5_6_S3_S4 notebook
tluettm Aug 10, 2026
c6206d9
erge branch 'thermal_conductivity_fix_and_paper_revision' of github.c…
tluettm Aug 10, 2026
deb2e67
avoid proper https link not to include it in the docs
slayoo Aug 10, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion PySDM/physics/diffusion_ice_kinetics/standard.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,7 +31,8 @@ def D(const, D, r, lmbd, T):
def K(const, K, r, lmbd, T, rho): # pylint: disable=too-many-arguments
return K / (
r / (r + lmbd)
+ K
+ 4.0
* K
/ const.HAC_ice
/ np.sqrt(8.0 * const.Rd * T / const.PI)
/ const.c_pd
Expand Down
0 ...Luettmer_homogeneous_freezing/__init__.py → ..._examples/Luettmer_et_al_2026/__init__.py
View file
Open in desktop
File renamed without changes.
4 changes: 2 additions & 2 deletions .../Luettmer_homogeneous_freezing/commons.py → ...M_examples/Luettmer_et_al_2026/commons.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
"""commons for homogeneous freezing notebooks"""

from PySDM_examples.Luettmer_homogeneous_freezing.settings import Settings
from PySDM_examples.Luettmer_homogeneous_freezing.simulation import Simulation
from PySDM_examples.Luettmer_et_al_2026.settings import Settings
from PySDM_examples.Luettmer_et_al_2026.simulation import Simulation
from PySDM import Formulae
from PySDM.physics.constants import si
from PySDM.backends import CPU
Expand Down
643 changes: 643 additions & 0 deletions examples/PySDM_examples/Luettmer_et_al_2026/fig_1_2_3.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

1,114 changes: 1,114 additions & 0 deletions examples/PySDM_examples/Luettmer_et_al_2026/fig_4_5_6_S3_S4.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

315 changes: 315 additions & 0 deletions examples/PySDM_examples/Luettmer_et_al_2026/fig_S1_S2.ipynb
View file
Open in desktop

Large diffs are not rendered by default.

69 changes: 35 additions & 34 deletions ...les/Luettmer_homogeneous_freezing/plot.py → ...ySDM_examples/Luettmer_et_al_2026/plot.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

ax_title_size = 18
ax_lab_fsize = 15
ax_lab_fsize2 = 19
tick_fsize = 15
T_frz_bins = np.linspace(-40, -34, num=60, endpoint=True)
T_frz_bins_kelvin = np.linspace(230, 240, num=100, endpoint=True)
Expand Down Expand Up @@ -92,8 +93,9 @@ def plot_thermodynamics_and_bulk(
ax = axs[iax]
ax.plot(time, RH, color="red", linestyle="dashdot", label=r"$S_\text{w}$")
ax.plot(time, RHi, color="blue", linestyle="--", label=r"$S_\text{i}$")
ax.set_ylabel("saturation ratio", fontsize=ax_lab_fsize)
ax.legend(loc="center left", fontsize=ax_lab_fsize)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
ax.set_ylabel("saturation ratio", fontsize=ax_lab_fsize2)
ax.legend(loc="center left", fontsize=ax_lab_fsize2)
ax.set_xlim(time[0], t_lim)
ax.tick_params(labelsize=tick_fsize)
ax.set_title(title_add[iax] + r"ambient thermodynamics", fontsize=ax_lab_fsize)
Expand All @@ -102,9 +104,9 @@ def plot_thermodynamics_and_bulk(

twin = ax.twinx()
twin.plot(time, T, color="black", linestyle="-", label="T")
twin.set_xlabel("time [s]", fontsize=ax_lab_fsize)
twin.set_ylabel("temperature [K]", fontsize=ax_lab_fsize)
twin.legend(loc="upper left", fontsize=ax_lab_fsize)
twin.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
twin.set_ylabel("temperature [K]", fontsize=ax_lab_fsize2)
twin.legend(loc="upper left", fontsize=ax_lab_fsize2)
twin.tick_params(labelsize=tick_fsize)

# mixing ratio and number concentration
Expand All @@ -115,9 +117,9 @@ def plot_thermodynamics_and_bulk(
ax.plot(time, qv, color="black", linestyle="-", label=r"$q_\text{v}$")
ax.set_yscale("log")
ax.set_ylim(1e-5, 1e-2)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
ax.set_ylabel(r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]", fontsize=ax_lab_fsize)
ax.legend(fontsize=ax_lab_fsize)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
ax.set_ylabel(r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]", fontsize=ax_lab_fsize2)
ax.legend(fontsize=ax_lab_fsize2)
ax.tick_params(labelsize=tick_fsize)
ax.set_xlim(time[0], t_lim)
ax.grid(visible=True)
Expand All @@ -128,9 +130,9 @@ def plot_thermodynamics_and_bulk(
twin.plot(time, nc, color="red", linestyle="densly dashdot", label="water")
twin.plot(time, ni, color="blue", linestyle="densly dashed", label="ice")
twin.set_yscale("log")
twin.set_xlabel("time [s]", fontsize=ax_lab_fsize)
twin.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
twin.set_ylabel(
r"number concentration [$\mathrm{kg^{-1}}$]", fontsize=ax_lab_fsize
r"number concentration [$\mathrm{kg^{-1}}$]", fontsize=ax_lab_fsize2
)
twin.tick_params(labelsize=tick_fsize)

Expand All @@ -147,8 +149,8 @@ def plot_thermodynamics_and_bulk(
linewidth=1.5,
)
ax.set_xlim(left=234, right=239)
ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize)
ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize)
ax.set_xlabel("freezing temperature [K]", fontsize=ax_lab_fsize2)
ax.set_ylabel("frozen fraction", fontsize=ax_lab_fsize2)
ax.tick_params(labelsize=tick_fsize)
ax.grid(visible=True)
ax.axvline(x=first_T_frz, color="black", linestyle=":")
Expand All @@ -158,9 +160,9 @@ def plot_thermodynamics_and_bulk(
ax.plot(time, ri * 1e6, color="blue", linestyle="--", label="ice")
ax.set_yscale("log")
ax.set_ylim(1e-2, 1e2)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
ax.set_ylabel("mean radius [µm]", fontsize=ax_lab_fsize)
ax.legend(fontsize=ax_lab_fsize)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
ax.set_ylabel("mean radius [µm]", fontsize=ax_lab_fsize2)
ax.legend(fontsize=ax_lab_fsize2)
ax.set_xlim(time[0], t_lim)
ax.tick_params(labelsize=tick_fsize)
ax.grid(visible=True)
Expand All @@ -173,7 +175,7 @@ def plot_thermodynamics_and_bulk(
if show_jhom:
lin_s_SP2023 = "--"
lin_s_KM2016 = "-"
if simulation["settings"]["hom_freezing"] == "Spichtinger2023":
if simulation["settings"]["hom_freezing"] == "KoopMurray2016_DWA":
lin_s_SP2023 = "-"
lin_s_KM2016 = "--"

Expand All @@ -192,15 +194,15 @@ def plot_thermodynamics_and_bulk(
label="JHOM-DWA",
)
ax.set_ylabel(
r"nucleation rate [$\mathrm{m^{-3} \, s^{-1}}$]", fontsize=ax_lab_fsize
r"nucleation rate [$\mathrm{m^{-3} \, s^{-1}}$]", fontsize=ax_lab_fsize2
)
ax.set_ylim(1e-30, 1e30)
ax.set_title(title_add[iax] + r"nucleation rates", fontsize=ax_lab_fsize)
ax.legend(loc="upper left", fontsize=ax_lab_fsize)
ax.legend(loc="upper left", fontsize=ax_lab_fsize2)
ax.set_yscale("log")
ax.set_xlim(time[0], t_lim)
ax.axvline(x=first_ice_time, color="black", linestyle=":")
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize)
ax.set_xlabel("time [s]", fontsize=ax_lab_fsize2)
if plot_daw:
twin = ax.twinx()
twin.plot(
Expand All @@ -211,9 +213,9 @@ def plot_thermodynamics_and_bulk(
label=r"$\Delta a_{w}$",
)
twin.set_ylim(0.2, 0.35)
twin.set_ylabel("water activity difference", fontsize=ax_lab_fsize)
twin.set_ylabel("water activity difference", fontsize=ax_lab_fsize2)
twin.tick_params(labelsize=tick_fsize)
twin.legend(loc="lower left", fontsize=ax_lab_fsize)
twin.legend(loc="lower left", fontsize=ax_lab_fsize2)
else:
twin = ax.twinx()
twin.plot(
Expand All @@ -224,24 +226,23 @@ def plot_thermodynamics_and_bulk(
linestyle="dashdot",
)
twin.set_ylim(-10, 10)
twin.set_ylabel("relative error", fontsize=ax_lab_fsize)
twin.set_ylabel("relative error", fontsize=ax_lab_fsize2)
twin.tick_params(labelsize=tick_fsize)
twin.legend(loc="lower left", fontsize=ax_lab_fsize)
twin.legend(loc="lower left", fontsize=ax_lab_fsize2)
else:
ax.scatter(radius * 1e6, multiplicity)
ax.set_yscale("log")
ax.set_xscale("log")
ax.set_xlim(1e-3, 5e-0)
ax.set_xlabel("initial radius [µm]", fontsize=ax_lab_fsize)
ax.set_ylabel("multiplicity", fontsize=ax_lab_fsize)
ax.set_xlabel("initial radius [µm]", fontsize=ax_lab_fsize2)
ax.set_ylabel("multiplicity", fontsize=ax_lab_fsize2)
ax.set_title(title_add[iax] + r"CCN size distribution", fontsize=ax_lab_fsize)

ax.tick_params(labelsize=tick_fsize)
ax.grid(visible=True)


def plot_freezing_temperatures_histogram(ax, simulation, plot_rhi=False):

number_of_ensemble_runs = simulation["settings"]["number_of_ensemble_runs"]

for i in range(number_of_ensemble_runs):
Expand Down Expand Up @@ -333,7 +334,6 @@ def plot_freezing_temperatures_2d_histogram_seaborn(
height=4,
width=5,
):

sns.set_theme(style="ticks")

second_axis = True
Expand Down Expand Up @@ -381,7 +381,7 @@ def plot_freezing_temperatures_2d_histogram_seaborn(
ens_variable_sec = ens_variable_sec * 1e3
binwidth = 0.25
elif ens_variable_name == "n_ccn":
y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]"
y_label = r"ccn concentration at STP [$\mathrm{cm^{-3}}$]"
y_label_sec = r"radius [$\mathrm{\mu m}$]"
ens_variable = ens_variable / 1.0e6
ens_variable_sec = ens_variable_sec * 1.0e6
Expand Down Expand Up @@ -453,7 +453,6 @@ def plot_freezing_temperatures_2d_histogram_seaborn(
def plot_ensemble_bulk(
ax, ensemble_simulations, var_name, title_add=""
): # pylint: disable=too-many-nested-blocks

colors = ["blue", "red", "cyan"]
linestyles = ["-", "--", ":"]
pyplot.rcParams["axes.prop_cycle"] = cycler(color=colors) + cycler(
Expand Down Expand Up @@ -517,10 +516,10 @@ def plot_ensemble_bulk(
title = r"$n_\text{i}$"
ax.set_xlim(1e6, 1e10)
elif var_name == "IWC":
ax.set_xscale("log")
ax.set_xscale("linear")
x_label = r"mixing ratio [$\mathrm{kg \, kg^{-1}}$]"
title = "ice mixing ratio"
ax.set_xlim(1e-4, 1e-3)
ax.set_xlim(3e-4, 1e-3)
elif var_name == "freezing_fraction":
title = r"$n_\text{frz}$"
x_label = r"frozen fraction [$\mathrm{\%}$]"
Expand All @@ -531,7 +530,7 @@ def plot_ensemble_bulk(
y_label = r"ccn concentration [$\mathrm{cm^{-3}}$]"
ens_label = r"$n_\text{ccn}$ ensemble"
elif ens_var_name == "w_updraft":
ax.set_yscale("log")
ax.set_yscale("linear")
y_label = r"vertical updraft [$\mathrm{m \, s^{-1}}$]"
ens_label = "w ensemble"
elif ens_var_name == "sigma_droplet_distribution":
Expand All @@ -541,9 +540,11 @@ def plot_ensemble_bulk(
ax.set_yscale("log")
y_label = "number of super-particles"
ens_label = r"$n_\text{sd}$ ensemble"

ax.set_title(title_add + " " + title + " for " + ens_label, fontsize=ax_lab_fsize)
ax.set_xlabel(x_label, fontsize=ax_lab_fsize)
ax.set_ylabel(y_label, fontsize=ax_lab_fsize)
ax.set_xlabel(x_label, fontsize=ax_lab_fsize2)
ax.set_ylabel(y_label, fontsize=ax_lab_fsize2)
ax.grid(visible=True)
ax.legend(fontsize=ax_lab_fsize)
ax.tick_params(labelsize=tick_fsize + 2)
return ax
0 ...Luettmer_homogeneous_freezing/settings.py → ..._examples/Luettmer_et_al_2026/settings.py
View file
Open in desktop
File renamed without changes.
0 .../simple_homogenous_freezing_example.ipynb → .../simple_homogenous_freezing_example.ipynb
View file
Open in desktop
File renamed without changes.
0 ...ettmer_homogeneous_freezing/simulation.py → ...xamples/Luettmer_et_al_2026/simulation.py
View file
Open in desktop
File renamed without changes.
Loading
Loading