PySDM_examples.Luettmer_et_al_2026
Homogeneous freezing example
{
"cells": [
{
"cell_type": "markdown",
"id": "75997570-1e10-497f-b966-db30eab8b700",
"metadata": {},
"source": [
" {
"cells": [
{
"cell_type": "markdown",
"id": "1f58ba8b-3ba1-4f69-8057-4bc85e1f5d18",
"metadata": {},
"source": [
" {
"cells": [
{
"cell_type": "markdown",
"id": "5ce42476-e0f9-4b2a-8fe4-6fb69f332da3",
"metadata": {},
"source": [
" {
"cells": [
{
"cell_type": "markdown",
"id": "cb2f8cf9-cbb8-4df3-818c-4e04e5957af2",
"metadata": {},
"source": [
"
\n",
"
"
]
},
{
"cell_type": "markdown",
"id": "ba7ba36d3363edbb",
"metadata": {},
"source": [
"Fig. 1, 2, 3"
]
},
{
"cell_type": "code",
"id": "059aa4b6",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:47:16.324870Z",
"start_time": "2026-08-14T12:47:16.317929Z"
}
},
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
"if 'google.colab' in sys.modules:\n",
" !pip --quiet install open-atmos-jupyter-utils\n",
" from open_atmos_jupyter_utils import pip_install_on_colab\n",
" pip_install_on_colab('PySDM-examples', 'PySDM')"
],
"outputs": [],
"execution_count": 1
},
{
"cell_type": "code",
"id": "789deaf73c1061e4",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:47:28.445396Z",
"start_time": "2026-08-14T12:47:17.377461Z"
}
},
"source": [
"import numpy as np\n",
"from PySDM.physics.constants import si\n",
"from PySDM_examples.Luettmer_et_al_2026.commons import run_simulations, hom_pure_droplet_freezing_backend, hom_pure_droplet_freezing_standard_setup\n",
"from PySDM_examples.Luettmer_et_al_2026 import plot\n",
"from matplotlib import pyplot\n",
"from open_atmos_jupyter_utils import show_plot"
],
"outputs": [],
"execution_count": 2
},
{
"cell_type": "code",
"id": "9bb8100b503bb58d",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:47:29.323161Z",
"start_time": "2026-08-14T12:47:28.451689Z"
}
},
"source": [
"# \"\"\" General settings \"\"\"\n",
"hom_freezing_types = [ \"KoopMurray2016\", \"KoopMurray2016_DWA\" ]\n",
"hom_freezing_labels = [ \"JHOM-T\", \"JHOM-DWA\"]\n",
"number_of_nsd = (1e1, 1e2, 1e3) if 'CI' not in os.environ else (1e2,)\n",
"vertical_updrafts = np.geomspace(0.2,10,num=9 if 'CI' not in os.environ else 2) * si.meter / si.second\n",
"number_concentrations = np.geomspace(100, 10000, num=10 if 'CI' not in os.environ else 2) / si.cm ** 3\n",
"\n",
"backends = hom_pure_droplet_freezing_backend()\n",
"standard = hom_pure_droplet_freezing_standard_setup()"
],
"outputs": [],
"execution_count": 3
},
{
"cell_type": "code",
"id": "58a5c533e06a93d2",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:06:04.007148Z",
"start_time": "2026-08-12T13:05:45.527298Z"
}
},
"source": [
"# Fig. 1 - homogeneous nucleation rates\n",
"def d_aw_ice(S_i,pv_sat_w,pv_sat_i):\n",
" return (S_i - 1) * pv_sat_i / pv_sat_w\n",
"\n",
"# backends setup\n",
"formulae = backends[\"KoopMurray2016\"].formulae\n",
"KM16_rate_formulae = backends[\"KoopMurray2016\"].formulae.homogeneous_ice_nucleation_rate.j_hom\n",
"KM16DWA_rate_formulae = backends[\"KoopMurray2016_DWA\"].formulae.homogeneous_ice_nucleation_rate.j_hom\n",
"SP23_rate_formulae = backends[\"Spichtinger2023\"].formulae.homogeneous_ice_nucleation_rate.j_hom\n",
"KP00_rate_formulae = backends[\"Koop2000\"].formulae.homogeneous_ice_nucleation_rate.j_hom\n",
"\n",
"# Set temperature and saturation ratio wtr water\n",
"fig1_T = np.arange(230,240,0.5) * si.K\n",
"fig1_s_w = np.arange(0.95,1.05,0.01)\n",
"T_mesh, s_w_mesh = np.meshgrid( fig1_T, fig1_s_w)\n",
"e_sat_w = formulae.saturation_vapour_pressure.pvs_water(T_mesh) * si.Pa\n",
"e_sat_i = formulae.saturation_vapour_pressure.pvs_ice(T_mesh) * si.Pa\n",
"s_i_mesh = s_w_mesh * e_sat_w / e_sat_i\n",
"d_aw_i_mesh = d_aw_ice(s_i_mesh,e_sat_w,e_sat_i)\n",
"\n",
"# Calculate rates\n",
"rates_KM16 = np.log10(KM16_rate_formulae( T_mesh, d_aw_i_mesh ))\n",
"rates_KM16DWA = np.log10(KM16DWA_rate_formulae( T_mesh, d_aw_i_mesh ))\n",
"rates_SP23 = np.log10(SP23_rate_formulae( T_mesh, d_aw_i_mesh ))\n",
"rates_KP00 = np.log10(KP00_rate_formulae( T_mesh, d_aw_i_mesh ))\n",
"\n",
"# plot settings\n",
"ax_title_size = 18\n",
"ax_lab_fsize = 15\n",
"tick_fsize = 15\n",
"\n",
"fig, axs = pyplot.subplots(1, 2, figsize=(12, 4), constrained_layout=True)\n",
"fig.suptitle(\"Homogeneous nucleation rates\",fontsize=20)\n",
"\n",
"axs = axs.ravel()\n",
"\n",
"idx = np.where(fig1_s_w == 1)[0]\n",
"axs[0]PySDM_examples.Luettmer_et_al_2026.plot(T_mesh[idx,:][0], rates_KM16[idx,:][0], color=\"black\", label=\"JHOM-T & JHOM-DWA $S_{w}=1$\")\n",
"axs[1]PySDM_examples.Luettmer_et_al_2026.plot(T_mesh[idx,:][0], rates_KM16[idx,:][0], color=\"black\", label=\"JHOM-T\")\n",
"axs[1]PySDM_examples.Luettmer_et_al_2026.plot(np.nan, np.nan, color=\"red\", label=\"JHOM-DWA\")\n",
"axs[0]PySDM_examples.Luettmer_et_al_2026.plot(T_mesh[idx,:][0], rates_KP00[idx,:][0], color=\"blue\", label=r\"KP00 $S_{w}=1$\")\n",
"\n",
"CS = axs[1].contour(\n",
" T_mesh,\n",
" rates_KM16DWA,\n",
" s_w_mesh,\n",
" levels=fig1_s_w[:-1],\n",
" colors=\"red\",\n",
" linewidths=0.9\n",
")\n",
"idx = 18\n",
"manual_locations = [(fig1_T[idx], val) for val in rates_KM16DWA[:,idx]]\n",
"axs[1].clabel(CS, fmt=lambda v: f\"{v:.2f}\", inline=True, fontsize=10, manual=manual_locations[:-1])\n",
"\n",
"axs[0].set_title( r\"(a) at water saturation\", fontsize=ax_lab_fsize )\n",
"axs[1].set_title( r\"(b) for various water saturation ratios\", fontsize=ax_lab_fsize )\n",
"axs[0].legend(loc=\"lower left\",fontsize=ax_lab_fsize)\n",
"axs[0].set_ylim(5,25)\n",
"axs[1].set_ylim(-20,20)\n",
"\n",
"for ax in axs:\n",
" ax.set_xlim(230,240)\n",
" ax.set_xlabel(r\"temperature [K]\",fontsize=ax_lab_fsize)\n",
" xticks = np.arange(230, 241, 1)\n",
" ax.set_xticks(xticks)\n",
" ax.tick_params(axis='both', labelsize=tick_fsize)\n",
" ax.set_yticks(np.arange(ax.get_ylim()[0], ax.get_ylim()[1] + 1, 5))\n",
" ax.set_ylabel(r\"$\log_{10} J_{\mathrm{hom}}\,[\mathrm{m^{-3}\,s^{-1}}]$\", fontsize=ax_lab_fsize)\n",
" ax.grid(visible=True)\n",
" ax.vlines(235, -20, 50, color=\"gray\")\n",
"axs[1].legend(loc=\"lower left\",fontsize=ax_lab_fsize)\n",
"show_plot(\"fig_1\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\"standard,\n",
" \"backend\": backends[\"KoopMurray2016_DWA\"],\n",
" \"hom_freezing\": \"KoopMurray2016_DWA\",\n",
" \"w_updraft\": 2.5 * si.meter / si.second,\n",
" \"deposition_enable\": False,\n",
" \"n_output\": 1,\n",
" \"dz\": 2.0 * si.meter,\n",
" \"n_sd\": int(1e3 if 'CI' not in os.environ else 1e1),\n",
" }\n",
"KoopMurray2016_DWA_reference_high_w_simulation = run_simulations(setting_dict)\n",
"\n"
],
"outputs": [],
"execution_count": 4
},
{
"cell_type": "code",
"id": "dcfcec994864fd60",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T13:07:50.204768Z",
"start_time": "2026-08-14T13:07:42.611450Z"
}
},
"source": [
"panel_labels = [\"(a)\", \"(b)\", \"(c)\", \"(d)\"]\n",
"frz_label = \" JHOM-T \"\n",
"panel_labels = [label + frz_label for label in panel_labels]\n",
"plot.plot_thermodynamics_and_bulk(KoopMurray2016_reference_high_w_simulation, title_add=panel_labels, t_lim=450.)\n",
"show_plot(\"fig_21\")\n",
"\n",
"panel_labels = [\"(e)\", \"(f)\", \"(g)\", \"(h)\"]\n",
"frz_label = \" JHOM-DWA \"\n",
"panel_labels = [label + frz_label for label in panel_labels]\n",
"plot.plot_thermodynamics_and_bulk(KoopMurray2016_DWA_reference_high_w_simulation, title_add=panel_labels, t_lim=450.)\n",
"show_plot(\"fig_22\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig_22.pdf
\"), HTML(value=\"standard,\n",
" \"n_ccn\": n_ccn,\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" \"deposition_enable\": False,\n",
" \"dz\": 2.0 * si.meter,\n",
" \"n_sd\": int(1e3 if 'CI' not in os.environ else 1e1),\n",
" }\n",
" ccn_ensemble_no_deposition[hom_freezing_type].append(run_simulations(setting_dict))"
],
"outputs": [],
"execution_count": 9
},
{
"cell_type": "code",
"id": "406b69b1022934eb",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:22:41.116261Z",
"start_time": "2026-08-12T13:22:35.887834Z"
}
},
"source": [
"# Figure 3\n",
"fig_id=\"fig3_\"\n",
"fig, axs = pyplot.subplots(1, 2, figsize=(12, 4), constrained_layout=True)\n",
"fig.suptitle(\"Ensemble simulations without DEP-ICE\", fontsize=20)\n",
"\n",
"axs = axs.ravel()\n",
"\n",
"axs[0] = plot.plot_freezing_temperatures_histogram_allinone(axs[0], nsd_ensemble_no_deposition[\"KoopMurray2016\"], title = r\"(a) $n_\text{sd}$ ensemble: JHOM-T\")\n",
"axs[1] = plot.plot_freezing_temperatures_histogram_allinone(axs[1], nsd_ensemble_no_deposition[\"KoopMurray2016_DWA\"], title = r\"(b) $n_\text{sd}$ ensemble: JHOM-DWA\")\n",
"show_plot(fig_id+\"ab\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_no_deposition, \"KoopMurray2016\", title = \"(c) w ensemble: JHOM-T\")\n",
"show_plot(fig_id+\"c\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_no_deposition, \"KoopMurray2016_DWA\", title = \"(d) w ensemble: JHOM-DWA\")\n",
"show_plot(fig_id+\"d\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_no_deposition, \"KoopMurray2016\", title = r\"(e) $n_\text{ccn}$ ensemble: JHOM-T\")\n",
"show_plot(fig_id+\"e\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_no_deposition, \"KoopMurray2016_DWA\", title = r\"(f) $n_\text{ccn}$ ensemble: JHOM-DWA\")\n",
"show_plot(fig_id+\"f\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig3_c.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig3_d.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig3_e.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig3_f.pdf
\"), HTML(value=\"
\n",
""
]
},
{
"cell_type": "markdown",
"id": "d31081c6317556e1",
"metadata": {},
"source": [
"Fig. 4, 5, 6 and Supplement Fig. 3 and 4"
]
},
{
"cell_type": "code",
"id": "f82b5b42",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:31:45.068601Z",
"start_time": "2026-08-12T13:31:45.062888Z"
}
},
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
"if 'google.colab' in sys.modules:\n",
" !pip --quiet install open-atmos-jupyter-utils\n",
" from open_atmos_jupyter_utils import pip_install_on_colab\n",
" pip_install_on_colab('PySDM-examples', 'PySDM')"
],
"outputs": [],
"execution_count": 1
},
{
"cell_type": "code",
"id": "e2b74470f7577bba",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"ExecuteTime": {
"end_time": "2026-08-14T11:48:08.412718Z",
"start_time": "2026-08-14T11:48:08.400694Z"
}
},
"source": [
"import numpy as np\n",
"from PySDM.physics.constants import si\n",
"from PySDM_examples.Luettmer_et_al_2026.commons import run_simulations, hom_pure_droplet_freezing_backend, hom_pure_droplet_freezing_standard_setup\n",
"from PySDM_examples.Luettmer_et_al_2026 import plot\n",
"from open_atmos_jupyter_utils import show_plot\n",
"from matplotlib import pyplot"
],
"outputs": [],
"execution_count": 39
},
{
"cell_type": "code",
"id": "7fb198e98a9e71c9",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:31:56.035800Z",
"start_time": "2026-08-12T13:31:55.051752Z"
}
},
"source": [
"# General settings\n",
"hom_freezing_types = [ \"KoopMurray2016\", \"KoopMurray2016_DWA\" ]\n",
"hom_freezing_labels = [ \"JHOM-T\", \"JHOM-DWA\"]\n",
"number_of_nsd = (1e2, 1e3, 1e4) if 'CI' not in os.environ else (1e1, 2e1)\n",
"vertical_updrafts = np.geomspace(0.2,10,num=9 if 'CI' not in os.environ else 2) * si.meter / si.second\n",
"number_concentrations = np.geomspace(100, 10000, num=10 if 'CI' not in os.environ else 2) / si.cm ** 3\n",
"\n",
"backends = hom_pure_droplet_freezing_backend()\n",
"standard = hom_pure_droplet_freezing_standard_setup()\n"
],
"outputs": [],
"execution_count": 3
},
{
"cell_type": "code",
"id": "dd8e02fbe466a63",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:39:17.753326Z",
"start_time": "2026-08-12T13:31:56.087484Z"
}
},
"source": [
"# Simulations for n_sd ensemble with deposition\n",
"nsd_ensemble_deposition = {}\n",
"nsd_ensemble_deposition[\"ens_variable\"] = number_of_nsd\n",
"nsd_ensemble_deposition[\"ens_variable_name\"] = \"n_sd\"\n",
"nsd_ensemble_deposition[\"hom_freezing_types\"] = hom_freezing_types\n",
"nsd_ensemble_deposition[\"hom_freezing_labels\"] = hom_freezing_labels\n",
"\n",
"for hom_freezing_type in hom_freezing_types:\n",
" nsd_ensemble_deposition[hom_freezing_type] = []\n",
" for n_sd in number_of_nsd:\n",
" setting_dict = {\n",
" **standard,\n",
" \"n_sd\": int(n_sd),\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" \"number_of_ensemble_runs\": 5 if 'CI' not in os.environ else 1,\n",
" \"deposition_enable\": True,\n",
" }\n",
" nsd_ensemble_deposition[hom_freezing_type].append( run_simulations(setting_dict) )"
],
"outputs": [],
"execution_count": 4
},
{
"cell_type": "code",
"id": "4113037920c9097f",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T13:53:56.281340Z",
"start_time": "2026-08-12T13:39:17.811196Z"
}
},
"source": [
"# Simulations for updraft histogram with deposition\n",
"updraft_ensemble_deposition = {}\n",
"updraft_ensemble_deposition[\"ens_variable\"] = vertical_updrafts\n",
"updraft_ensemble_deposition[\"ens_variable_name\"] = \"w_updraft\"\n",
"updraft_ensemble_deposition[\"hom_freezing_types\"] = hom_freezing_types\n",
"updraft_ensemble_deposition[\"hom_freezing_labels\"] = hom_freezing_labels\n",
"\n",
"for hom_freezing_type in hom_freezing_types:\n",
" updraft_ensemble_deposition[hom_freezing_type] = []\n",
" for updraft in vertical_updrafts:\n",
" setting_dict = {\n",
" **standard,\n",
" \"w_updraft\": updraft * si.meter / si.second,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" \"number_of_ensemble_runs\": 1,\n",
" }\n",
"\n",
" updraft_ensemble_deposition[hom_freezing_type].append( run_simulations(setting_dict) )"
],
"outputs": [],
"execution_count": 5
},
{
"cell_type": "code",
"id": "a4f8263d88be6374",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T14:04:24.497429Z",
"start_time": "2026-08-12T13:53:56.343164Z"
}
},
"source": [
"# Simulations for CCN concentration histogram\n",
"ccn_ensemble_deposition = {}\n",
"ccn_ensemble_deposition[\"ens_variable\"] = number_concentrations\n",
"ccn_ensemble_deposition[\"ens_variable_name\"] = \"n_ccn\"\n",
"ccn_ensemble_deposition[\"hom_freezing_types\"] = hom_freezing_types\n",
"ccn_ensemble_deposition[\"hom_freezing_labels\"] = hom_freezing_labels\n",
"\n",
"for hom_freezing_type in hom_freezing_types:\n",
" ccn_ensemble_deposition[hom_freezing_type] = []\n",
" for n_dv in number_concentrations:\n",
" setting_dict = {\n",
" *standard,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" \"n_ccn\": n_dv,\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" \"dz\": 0.5 * si.meter,\n",
" }\n",
" ccn_ensemble_deposition[hom_freezing_type].append(run_simulations(setting_dict))"
],
"outputs": [],
"execution_count": 6
},
{
"cell_type": "code",
"id": "223a5e5c3ae6641f",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T14:04:31.861198Z",
"start_time": "2026-08-12T14:04:24.563667Z"
}
},
"source": [
"# Plot Figure 4\n",
"fig_name = \"fig4_\"\n",
"fig, axs = pyplot.subplots(1, 2, figsize=(10, 4), constrained_layout=True)\n",
"fig.suptitle(\"Ensemble simulations with DEP-ICE\", fontsize=20)\n",
"\n",
"axs = axs.ravel()\n",
"\n",
"axs[0] = plot.plot_freezing_temperatures_histogram_allinone(axs[0], nsd_ensemble_deposition[\"KoopMurray2016\"], title = r\"(a) $n_\text{sd}$ ensemble: JHOM-T\", lloc=\"lower left\")\n",
"axs[1] = plot.plot_freezing_temperatures_histogram_allinone(axs[1], nsd_ensemble_deposition[\"KoopMurray2016_DWA\"], title = r\"(b) $n_\text{sd}$ ensemble: JHOM-DWA\", lloc=\"lower left\")\n",
"show_plot(fig_name+\"ab\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_deposition, \"KoopMurray2016\", title = \"(c) w ensemble: JHOM-T\")\n",
"show_plot(fig_name+\"c\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_deposition, \"KoopMurray2016_DWA\", title = \"(d) w ensemble: JHOM-DWA\")\n",
"show_plot(fig_name+\"d\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_deposition, \"KoopMurray2016\", title = r\"(e) $n_\text{ccn}$ ensemble: JHOM-T\")\n",
"show_plot(fig_name+\"e\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_deposition, \"KoopMurray2016_DWA\", title = r\"(f) $n_\text{ccn}$ ensemble: JHOM-DWA\")\n",
"show_plot(fig_name+\"f\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig4_c.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig4_d.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig4_e.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig4_f.pdf
\"), HTML(value=\"standard_sig,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" \"sigma_droplet_distribution\": sig,\n",
" \"deposition_enable\": False,\n",
" }\n",
" sig_ensemble_no_deposition_sig[hom_freezing_type].append(run_simulations(setting_dict))"
],
"outputs": [],
"execution_count": 16
},
{
"cell_type": "code",
"id": "d42dadd688c8c021",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:02:24.117077Z",
"start_time": "2026-08-14T12:02:14.854561Z"
}
},
"source": [
"# Plot Figure 6 and Sup 3\n",
"fig_name = \"fig_sup_3\"\n",
"fig, axs = pyplot.subplots(1, 2, figsize=(10, 4), constrained_layout=True)\n",
"fig.suptitle(\"Ensemble simulations for lognormal CCN size spectra\", fontsize=20)\n",
"\n",
"axs = axs.ravel()\n",
"axs[0] = plot.plot_freezing_temperatures_histogram_allinone(axs[0], nsd_ensemble_no_deposition_sig[\"KoopMurray2016\"], title = r\"(a) $n_\text{sd}$ ensemble: JHOM-T no DEP\", lloc=\"lower left\")\n",
"axs[1] = plot.plot_freezing_temperatures_histogram_allinone(axs[1], nsd_ensemble_deposition_sig[\"KoopMurray2016\"], title = r\"(b) $n_\text{sd}$ ensemble: JHOM-T with DEP\", lloc=\"lower left\")\n",
"show_plot(fig_name)\n",
"\n",
"\n",
"fig_name = \"fig6_\"\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(sig_ensemble_no_deposition_sig, \"KoopMurray2016\", title = r\"(a) $\mathrm{\sigma}$ ensemble: JHOM-T no DEP\")\n",
"show_plot(fig_name+\"a\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(sig_ensemble_deposition_sig, \"KoopMurray2016\", title = r\"(b) $\mathrm{\sigma}$ ensemble: JHOM-T with DEP\")\n",
"show_plot(fig_name+\"b\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_no_deposition_sig, \"KoopMurray2016\", title = \"(c) w ensemble: JHOM-T no DEP\")\n",
"show_plot(fig_name+\"c\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(updraft_ensemble_deposition_sig, \"KoopMurray2016\", title = \"(d) w ensemble: JHOM-T with DEP\")\n",
"show_plot(fig_name+\"d\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_no_deposition_sig, \"KoopMurray2016\", title = r\"(e) $n_\text{ccn}$ ensemble: JHOM-T no DEP\")\n",
"show_plot(fig_name+\"e\")\n",
"\n",
"plot.plot_freezing_temperatures_2d_histogram_seaborn(ccn_ensemble_deposition_sig, \"KoopMurray2016\", title = r\"(f) $n_\text{ccn}$ ensemble: JHOM-T with DEP\")\n",
"show_plot(fig_name+\"f\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\"
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig6_b.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig6_c.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig6_d.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig6_e.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig6_f.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig5.pdf
\"), HTML(value=\"./fig_sup5.pdf
\"), HTML(value=\""
],
"image/svg+xml": "\n\n\n"
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"HBox(children=(HTML(value=\"./fig_sup_4.pdf
\"), HTML(value=\"\n",
""
]
},
{
"cell_type": "markdown",
"id": "516eb145effbc34e",
"metadata": {},
"source": [
"Supplement Fig. 1 and 2"
]
},
{
"cell_type": "code",
"id": "14b48fe6",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:31:45.319943Z",
"start_time": "2026-08-14T12:31:45.306075Z"
}
},
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
"if 'google.colab' in sys.modules:\n",
" !pip --quiet install open-atmos-jupyter-utils\n",
" from open_atmos_jupyter_utils import pip_install_on_colab\n",
" pip_install_on_colab('PySDM-examples', 'PySDM')"
],
"outputs": [],
"execution_count": 1
},
{
"cell_type": "code",
"id": "f17f75bc6ed1d7d",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"ExecuteTime": {
"end_time": "2026-08-14T12:32:07.372590Z",
"start_time": "2026-08-14T12:31:45.353073Z"
}
},
"source": [
"from PySDM_examples.Luettmer_et_al_2026.commons import run_simulations, hom_pure_droplet_freezing_backend, hom_pure_droplet_freezing_standard_setup\n",
"from PySDM_examples.Luettmer_et_al_2026 import plot\n",
"from open_atmos_jupyter_utils import show_plot"
],
"outputs": [],
"execution_count": 2
},
{
"cell_type": "code",
"id": "1905ef6351c078c3",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:32:08.986536Z",
"start_time": "2026-08-14T12:32:07.724127Z"
}
},
"source": [
"backends = hom_pure_droplet_freezing_backend()\n",
"standard = hom_pure_droplet_freezing_standard_setup()\n",
"\n",
"standard_sig = hom_pure_droplet_freezing_standard_setup()\n",
"standard_sig[\"type_droplet_distribution\"] = \"lognormal\"\n",
"standard_sig[\"sigma_droplet_distribution\"] = 2.0\n",
"standard_sig[\"r_ccn\"] = 0.05e-6\n",
"standard_sig[\"n_sd\"] = int(1e4)"
],
"outputs": [],
"execution_count": 3
},
{
"cell_type": "code",
"id": "d8af8a8a03687ada",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:35:52.565719Z",
"start_time": "2026-08-14T12:32:09.048834Z"
}
},
"source": [
"# High output step reference simulations\n",
"setting_dict = {\n",
" **standard,\n",
" \"backend\": backends[\"KoopMurray2016\"],\n",
" \"hom_freezing\": \"KoopMurray2016\",\n",
" \"w_updraft\": 5.,\n",
" \"deposition_enable\": True,\n",
" \"n_output\": 1,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" }\n",
"reference_high_dep_KoopMurray2016_simulation = run_simulations(setting_dict)\n",
"setting_dict = {\n",
" **standard,\n",
" \"backend\": backends[\"KoopMurray2016_DWA\"],\n",
" \"hom_freezing\": \"KoopMurray2016_DWA\",\n",
" \"w_updraft\": 5.,\n",
" \"deposition_enable\": True,\n",
" \"n_output\": 1,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" }\n",
"reference_high_dep_Spichtinger2023_simulation = run_simulations(setting_dict)"
],
"outputs": [],
"execution_count": 4
},
{
"cell_type": "code",
"id": "5fb02478c183102b",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T13:08:29.092087Z",
"start_time": "2026-08-14T13:08:21.627333Z"
}
},
"source": [
"panel_labels = [\"(a)\", \"(b)\", \"(c)\", \"(d)\"]\n",
"frz_label = \" JHOM-T \"\n",
"panel_labels = [label + frz_label for label in panel_labels]\n",
"plot.plot_thermodynamics_and_bulk(reference_high_dep_KoopMurray2016_simulation, title_add=panel_labels, t_lim=300.)\n",
"show_plot(\"fig_sup_11\")\n",
"\n",
"panel_labels = [\"(e)\", \"(f)\", \"(g)\", \"(h)\"]\n",
"frz_label = \" JHOM-DWA \"\n",
"panel_labels = [label + frz_label for label in panel_labels]\n",
"plot.plot_thermodynamics_and_bulk(reference_high_dep_Spichtinger2023_simulation, title_add=panel_labels, t_lim=300.)\n",
"show_plot(\"fig_sup_12\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\"…"
],
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "25012e62506d4fad925933fc4caf9f1b"
}
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
},
{
"data": {
"text/plain": [
"
\"), HTML(value=\"…"
],
"application/vnd.jupyter.widget-view+json": {
"version_major": 2,
"version_minor": 0,
"model_id": "ba29f192ef1747838a626890d8954f37"
}
},
"metadata": {},
"output_type": "display_data",
"jetTransient": {
"display_id": null
}
}
],
"execution_count": 10
},
{
"cell_type": "code",
"id": "e149c6c838849dab",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:43:04.951533Z",
"start_time": "2026-08-14T12:42:12.854649Z"
}
},
"source": [
"# High output step reference simulations\n",
"setting_dict = {\n",
" *standard_sig,\n",
" \"backend\": backends[\"KoopMurray2016\"],\n",
" \"hom_freezing\": \"KoopMurray2016\",\n",
" \"w_updraft\": 2.5,\n",
" \"n_output\": 1,\n",
" \"silent\": False,\n",
" \"deposition_enable\": False,\n",
" \"n_sd\": int(1e4 if 'CI' not in os.environ else 1e1),\n",
" \"sigma_droplet_distribution\": 2.0,\n",
" \"r_ccn\": 50.e-9,\n",
" }\n",
"reference_high_sig20_simulation = run_simulations(setting_dict)"
],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Setting up simulation for KoopMurray2016 with wpdraft=2.5 and N_sd=10000 and n_ccn=749999999.9999999\n",
"Starting simulation...\n",
"all particles frozen or evaporated\n"
]
}
],
"execution_count": 7
},
{
"cell_type": "code",
"id": "37587d30ed4ddc1a",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-14T12:46:52.855958Z",
"start_time": "2026-08-14T12:46:48.623206Z"
}
},
"source": [
"panel_labels = [\"(a)\", \"(b)\", \"(c)\", \"(d)\"]\n",
"frz_label = \" \"\n",
"panel_labels = [label + frz_label for label in panel_labels]\n",
"plot.plot_thermodynamics_and_bulk(reference_high_sig20_simulation, title_add=panel_labels, show_jhom=False, show_tf=True)\n",
"show_plot(\"fig_sup_2\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
\"), HTML(value=\"\n",
""
]
},
{
"cell_type": "markdown",
"id": "940be08d99e094e",
"metadata": {},
"source": [
"Simple homogeneous freezing example for testing"
]
},
{
"cell_type": "code",
"id": "4a2317c7",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T12:04:36.267918Z",
"start_time": "2026-08-12T12:04:36.261795Z"
}
},
"source": [
"import os, sys\n",
"os.environ['NUMBA_THREADING_LAYER'] = 'workqueue' # PySDM & PyMPDATA don't work with TBB; OpenMP has extra dependencies on macOS\n",
"if 'google.colab' in sys.modules:\n",
" !pip --quiet install open-atmos-jupyter-utils\n",
" from open_atmos_jupyter_utils import pip_install_on_colab\n",
" pip_install_on_colab('PySDM-examples', 'PySDM')"
],
"outputs": [],
"execution_count": 1
},
{
"cell_type": "code",
"id": "9df5e9883a926bd4",
"metadata": {
"collapsed": true,
"jupyter": {
"outputs_hidden": true
},
"ExecuteTime": {
"end_time": "2026-08-12T12:04:46.532551Z",
"start_time": "2026-08-12T12:04:36.324129Z"
}
},
"source": [
"from PySDM_examples.Luettmer_et_al_2026.commons import run_simulations, hom_pure_droplet_freezing_backend, \\n",
" hom_pure_droplet_freezing_standard_setup\n",
"from PySDM_examples.Luettmer_et_al_2026 import plot\n",
"from matplotlib import pyplot\n",
"from open_atmos_jupyter_utils import show_plot"
],
"outputs": [],
"execution_count": 2
},
{
"cell_type": "code",
"id": "f4c34819839791db",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T12:04:47.556060Z",
"start_time": "2026-08-12T12:04:46.681120Z"
}
},
"source": [
"hom_freezing_types_all = [ \"KoopMurray2016\", \"Spichtinger2023\", \"KoopMurray2016_DWA\", \"threshold\" ]\n",
"backends = hom_pure_droplet_freezing_backend()\n",
"standard = hom_pure_droplet_freezing_standard_setup()\n",
"standard[\"n_sd\"] = 100\n",
"standard[\"dz\"] = standard[\"dz\"] if 'CI' not in os.environ else 25\n",
"standard[\"deposition_enable\"] = False"
],
"outputs": [],
"execution_count": 3
},
{
"cell_type": "code",
"id": "f2ce64e861da85f1",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T12:08:21.466882Z",
"start_time": "2026-08-12T12:04:47.617663Z"
}
},
"source": [
"simulations = []\n",
"for hom_freezing_type in hom_freezing_types_all:\n",
" setting_dict = {\n",
" **standard,\n",
" \"backend\": backends[hom_freezing_type],\n",
" \"hom_freezing\": hom_freezing_type,\n",
" }\n",
" simulations.append(run_simulations(setting_dict))"
],
"outputs": [],
"execution_count": 4
},
{
"cell_type": "code",
"id": "e94f833a36fe611d",
"metadata": {
"ExecuteTime": {
"end_time": "2026-08-12T12:08:23.426727Z",
"start_time": "2026-08-12T12:08:21.541711Z"
}
},
"source": [
"fig, axs = pyplot.subplots(2, 2, figsize=(8, 8), constrained_layout=True)\n",
"axs = axs.ravel()\n",
"for i,simulation in enumerate(simulations):\n",
" axs[i] = plot.plot_freezing_temperatures_histogram(axs[i], simulation,plot_rhi=False)\n",
"show_plot(\"fig_tfrz_simple_hom_freezing_histogram\")"
],
"outputs": [
{
"data": {
"text/plain": [
"
1""" 2Homogeneous freezing example 3 4.. include:: ./fig_1_2_3.ipynb 5.. include:: ./fig_4_5_6_S3_S4.ipynb 6.. include:: ./fig_S1_S2.ipynb 7.. include:: ./simple_homogenous_freezing_example.ipynb 8""" 9 10from .settings import Settings 11from .simulation import Simulation 12from .plot import ( 13 plot_thermodynamics_and_bulk, 14 plot_freezing_temperatures_histogram, 15 plot_freezing_temperatures_histogram_allinone, 16 plot_freezing_temperatures_2d_histogram_seaborn, 17)