Skip to content

Implementation of Twomey (1959) eq'n for CCN activation spectra (journal club) #1508

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.

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
wants to merge 4 commits into
base: main
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
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
4 changes: 3 additions & 1 deletion PySDM/formulae.py
Original file line number Diff line number Diff line change
Expand Up @@ -57,6 +57,7 @@ def __init__( # pylint: disable=too-many-locals
terminal_velocity: str = "GunnKinzer1949",
air_dynamic_viscosity: str = "ZografosEtAl1987",
bulk_phase_partitioning: str = "Null",
ccn_activation_spectrum: str = "Null",
handle_all_breakups: bool = False,
):
# initialisation of the fields below is just to silence pylint and to enable code hints
Expand Down Expand Up @@ -90,7 +91,8 @@ def __init__( # pylint: disable=too-many-locals
self.air_dynamic_viscosity = air_dynamic_viscosity
self.terminal_velocity = terminal_velocity
self.bulk_phase_partitioning = bulk_phase_partitioning

self.ccn_activation_spectrum = ccn_activation_spectrum

self._components = tuple(
i
for i in dir(self)
Expand Down
1 change: 1 addition & 0 deletions PySDM/physics/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -47,5 +47,6 @@
air_dynamic_viscosity,
terminal_velocity,
bulk_phase_partitioning,
ccn_activation_spectrum,
)
from .constants import convert_to, in_unit, si
2 changes: 2 additions & 0 deletions PySDM/physics/ccn_activation_spectrum/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
from PySDM.impl.null_physics_class import Null
from .twomey_1959 import Twomey1959
15 changes: 15 additions & 0 deletions PySDM/physics/ccn_activation_spectrum/twomey_1959.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,15 @@
import numpy as np
"""
[Twomey 1959](https://doi.org/10.1007/BF01993560). Note that supersaturation is expressed in temperature unit as the elevation of the dew point or as percent of supersaturation; and concentrations are reported for 10C and 800 mb.
"""

class Twomey1959:
def __init__(self, const):
assert np.isfinite(const.TWOMEY_K)
assert np.isfinite(const.TWOMEY_N0)

@staticmethod
def ccn_concentration(const, saturation_ratio):
return const.TWOMEY_N0*np.power(saturation_ratio-1, const.TWOMEY_K)


31 changes: 31 additions & 0 deletions PySDM/physics/constants_defaults.py
Original file line number Diff line number Diff line change
Expand Up @@ -684,3 +684,34 @@ def compute_derived_values(c: dict):
c["water_molar_volume"] = c["Mv"] / c["rho_w"]
c["rho_STP"] = c["p_STP"] / c["Rd"] / c["T_STP"]
c["H_u"] = c["M"] / c["p_STP"]

W76W_G0 = -2.9912729e3 * si.K**2
W76W_G1 = -6.0170128e3 * si.K
W76W_G2 = 1.887643854e1
W76W_G3 = -2.8354721e-2 * si.K**-1
W76W_G4 = 1.7838301e-5 * si.K**-2
W76W_G5 = -8.4150417e-10 * si.K**-3
W76W_G6 = 4.4412543e-13 * si.K**-4
W76W_G7 = 2.858487
W76W_G8 = 1 * si.Pa

B80W_G0 = 6.112 * si.hPa
B80W_G1 = 17.67 * si.dimensionless
B80W_G2 = 243.5 * si.K

one_kelvin = 1 * si.K

bulk_phase_partitioning_T_cold = 235 * si.K
bulk_phase_partitioning_T_warm = 273 * si.K
bulk_phase_partitioning_exponent = np.nan


BOLIN_ISOTOPE_TIMESCALE_COEFF_C1 = np.nan * si.dimensionless
"""
Coeffitient c1 used in [Bolin 1958](https://https://digitallibrary.un.org/record/3892725)
for the falling drop evaporation timescale of equilibration with ambient air void of a given
isotopologue; in the paper timescale is calculated for tritium with assumption of no tritium
in the environment around the drop (Table 1).
"""
TWOMEY_K = np.nan
TWOMEY_N0 = np.nan
30 changes: 30 additions & 0 deletions tests/unit_tests/physics/test_ccn_activation_spectrum.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,30 @@
import numpy as np
from matplotlib import pyplot

from PySDM import Formulae
from PySDM.physics.constants import PER_CENT, si, in_unit

def test_twomey_and_wojciechowski_1969_fig1(plot=True):
"""[Twomey, Wojciechowski 1969](https://doi.org/10.1175/1520-0469(1969)26%3C648:OOTGVO%3E2.0.CO;2)
"""
#arange
for k, N in zip([0.5,0.7], [100,500]):
formulae = Formulae(ccn_activation_spectrum="Twomey1959",constants={"TWOMEY_K": k, "TWOMEY_N0": N/si.cm**3})
supersaturation = np.logspace(np.log10(.2), np.log10(9))*PER_CENT
#act
activated_nuclei_concentration = formulae.ccn_activation_spectrum.ccn_concentration(saturation_ratio=supersaturation+1)
#plot
pyplot.plot(in_unit(supersaturation,PER_CENT),
in_unit(activated_nuclei_concentration,si.cm**-3),
label=f"{k=}"
)
pyplot.xlim(0.1,10)
pyplot.ylim(1,1000)
pyplot.xscale("log")
pyplot.yscale("log")
pyplot.xlabel("Percent supersaturation")
pyplot.ylabel("Nuclei [cm$^{-3}$]")
pyplot.grid()
pyplot.legend()
pyplot.show()
#assert
Loading