diff --git a/PySDM/initialisation/spectra/lognormal.py b/PySDM/initialisation/spectra/lognormal.py
index 1a318dee90..3d799a7a83 100644
--- a/PySDM/initialisation/spectra/lognormal.py
+++ b/PySDM/initialisation/spectra/lognormal.py
@@ -11,23 +11,21 @@
 
 
 class Lognormal(Spectrum):
-    def __init__(self, norm_factor: float, m_mode: float, s_geom: float):
+    def __init__(self, *, norm_factor: float, median=None, s_geom: float, mode=None):
         """`norm_factor=1` corresponds to standard normalised probability density,
         other settings allow to express, e.g., size or mass distributions;
         `m_mode` is the median value, `s_geom` is the geometric standard deviation"""
-        super().__init__(lognorm, (math.log(s_geom), 0, m_mode), norm_factor)
+        assert (median is None) ^ (mode is None)
+        self.median = median or mode * math.exp(s_geom**2)
+        super().__init__(lognorm, (math.log(s_geom), 0, self.median), norm_factor)
 
     @property
     def s_geom(self):
         return math.exp(self.distribution_params[0])
 
     @property
-    def m_mode(self):
-        return self.distribution_params[2]
-
-    @property
-    def median(self):
-        return self.m_mode
+    def mode(self):
+        return self.median / math.exp(self.s_geom**2)
 
     @property
     def geometric_mean(self):
@@ -37,6 +35,6 @@ def __str__(self):
         return (
             f"{self.__class__.__name__}:"
             f" (N={self.norm_factor:.3g},"
-            f" m_mode={self.m_mode:.3g},"
+            f" median={self.median:.3g},"
             f" s_geom={self.s_geom:.3g})"
         )
diff --git a/docs/markdown/pysdm_landing.md b/docs/markdown/pysdm_landing.md
index 350cf70971..6969120e68 100644
--- a/docs/markdown/pysdm_landing.md
+++ b/docs/markdown/pysdm_landing.md
@@ -402,7 +402,7 @@ env = Parcel(
     T0=300 * si.K,
     w= 2.5 * si.m / si.s
 )
-spectrum = Lognormal(norm_factor=1e4/si.mg, m_mode=50*si.nm, s_geom=1.4)
+spectrum = Lognormal(norm_factor=1e4/si.mg, median=50*si.nm, s_geom=1.4)
 kappa = .5 * si.dimensionless
 cloud_range = (.5 * si.um, 25 * si.um)
 output_interval = 4
@@ -483,7 +483,7 @@ env = Parcel(pyargs( ...
     'T0', 300 * si.K, ...
     'w', 2.5 * si.m / si.s ...
 ));
-spectrum = Lognormal(pyargs('norm_factor', 1e4/si.mg, 'm_mode', 50 * si.nm, 's_geom', 1.4));
+spectrum = Lognormal(pyargs('norm_factor', 1e4/si.mg, 'median', 50 * si.nm, 's_geom', 1.4));
 kappa = .5;
 cloud_range = py.tuple({.5 * si.um, 25 * si.um});
 output_interval = 4;
@@ -584,7 +584,7 @@ env = Parcel(
   T0=300 * si.K,
   w=2.5 * si.m / si.s
 )
-spectrum = Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)
+spectrum = Lognormal(norm_factor=1e4 / si.mg, median=50 * si.nm, s_geom=1.5)
 kappa = .5 * si.dimensionless
 cloud_range = (.5 * si.um, 25 * si.um)
 output_interval = 4
diff --git a/tests/smoke_tests/parcel_a/pyrcel/test_parcel_example.py b/tests/smoke_tests/parcel_a/pyrcel/test_parcel_example.py
index 0f2f9398c9..3157013cf1 100644
--- a/tests/smoke_tests/parcel_a/pyrcel/test_parcel_example.py
+++ b/tests/smoke_tests/parcel_a/pyrcel/test_parcel_example.py
@@ -25,12 +25,13 @@ def test_humidity_and_temperature_profile(s_max, s_250m, T_250m, scipy_solver):
         settings = Settings(
             dz=1 * si.m,
             n_sd_per_mode=(5, 5),
+            # see also https://github.com/darothen/pyrcel/pull/26
             aerosol_modes_by_kappa={
                 0.54: Lognormal(
-                    norm_factor=850 / si.cm**3, m_mode=15 * si.nm, s_geom=1.6
+                    norm_factor=850 / si.cm**3, median=15 * si.nm, s_geom=1.6
                 ),
                 1.2: Lognormal(
-                    norm_factor=10 / si.cm**3, m_mode=850 * si.nm, s_geom=1.2
+                    norm_factor=10 / si.cm**3, median=850 * si.nm, s_geom=1.2
                 ),
             },
             vertical_velocity=1.0 * si.m / si.s,
diff --git a/tests/smoke_tests/parcel_c/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py b/tests/smoke_tests/parcel_c/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py
index 9f2b64bcbf..1aa970478c 100644
--- a/tests/smoke_tests/parcel_c/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py
+++ b/tests/smoke_tests/parcel_c/abdul_razzak_ghan_2000/test_single_supersaturation_peak.py
@@ -42,7 +42,7 @@ def test_single_saturation_peak(
     n_steps = 70
     n_sd = 2
     kappa = 0.4
-    spectrum = Lognormal(norm_factor=5000 / si.cm**3, m_mode=50.0 * si.nm, s_geom=2.0)
+    spectrum = Lognormal(norm_factor=5000 / si.cm**3, median=50.0 * si.nm, s_geom=2.0)
     formulae = Formulae(constants=CONSTANTS_ARG)
     builder = Builder(
         backend=CPU(formulae),
diff --git a/tests/unit_tests/dynamics/collisions/test_croupiers.py b/tests/unit_tests/dynamics/collisions/test_croupiers.py
index 421ada1eae..b461746d3a 100644
--- a/tests/unit_tests/dynamics/collisions/test_croupiers.py
+++ b/tests/unit_tests/dynamics/collisions/test_croupiers.py
@@ -16,14 +16,14 @@ def test_final_state(croupier, backend_class):
 
     # Arrange
     n_part = 100000
-    v_mean = 2e-6
+    v_median = 2e-6
     d = 1.2
     n_sd = 32
     x = 4
     y = 4
 
     attributes = {}
-    spectrum = Lognormal(n_part, v_mean, d)
+    spectrum = Lognormal(norm_factor=n_part, median=v_median, s_geom=d)
     attributes["volume"], attributes["multiplicity"] = Linear(
         spectrum
     ).sample_deterministic(n_sd)
diff --git a/tests/unit_tests/dynamics/condensation/test_parcel_sanity_checks.py b/tests/unit_tests/dynamics/condensation/test_parcel_sanity_checks.py
index ce99960e77..a13cf4fefc 100644
--- a/tests/unit_tests/dynamics/condensation/test_parcel_sanity_checks.py
+++ b/tests/unit_tests/dynamics/condensation/test_parcel_sanity_checks.py
@@ -17,7 +17,7 @@
 from PySDM.physics import si
 
 FORMULAE = Formulae()
-SPECTRUM = Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)
+SPECTRUM = Lognormal(norm_factor=1e4 / si.mg, median=50 * si.nm, s_geom=1.5)
 N_SD = 64
 R_DRY, specific_concentration = spectral_sampling.Logarithmic(
     SPECTRUM
diff --git a/tests/unit_tests/impl/test_moments.py b/tests/unit_tests/impl/test_moments.py
index b6ad286f99..a08cf148b2 100644
--- a/tests/unit_tests/impl/test_moments.py
+++ b/tests/unit_tests/impl/test_moments.py
@@ -14,11 +14,11 @@ class TestMaths:
     def test_moment_0d(backend_class):
         # Arrange
         n_part = 100000
-        v_mean = 2e-6
+        v_median = 2e-6
         d = 1.2
         n_sd = 32
 
-        spectrum = Lognormal(n_part, v_mean, d)
+        spectrum = Lognormal(norm_factor=n_part, median=v_median, s_geom=d)
         v, n = Linear(spectrum).sample_deterministic(n_sd)
         T = 300.0
         n = discretise_multiplicities(n)
@@ -75,11 +75,11 @@ def test_moment_0d(backend_class):
     def test_spectrum_moment_0d(backend_class):
         # Arrange
         n_part = 100000
-        v_mean = 2e-6
+        v_median = 2e-6
         d = 1.2
         n_sd = 32
 
-        spectrum = Lognormal(n_part, v_mean, d)
+        spectrum = Lognormal(norm_factor=n_part, median=v_median, s_geom=d)
         v, n = Linear(spectrum).sample_deterministic(n_sd)
         T = 300.0
         n = discretise_multiplicities(n)
diff --git a/tests/unit_tests/initialisation/test_spectra_lognormal.py b/tests/unit_tests/initialisation/test_spectra_lognormal.py
index 7ca5a1d64e..cd3568f616 100644
--- a/tests/unit_tests/initialisation/test_spectra_lognormal.py
+++ b/tests/unit_tests/initialisation/test_spectra_lognormal.py
@@ -10,43 +10,59 @@ class TestSpectraLognormal:
 
     @staticmethod
     @pytest.mark.parametrize(
-        "m_mode, s_geom",
+        "median, s_geom",
         (
             (0.01, 1.5),
             (0.1, 1.1),
             (1, 1.01),
         ),
     )
-    def test_median(m_mode, s_geom):
+    def test_median(median, s_geom):
         # arrange
-        sut = Lognormal(m_mode=m_mode, norm_factor=1, s_geom=s_geom)
+        sut = Lognormal(median=median, norm_factor=1, s_geom=s_geom)
 
         # act
-        median = sut.percentiles(0.5)
+        actual = sut.percentiles(0.5)
 
         # assert
-        assert median == m_mode
+        assert median == actual
 
     @staticmethod
     @pytest.mark.parametrize(
-        "m_mode, s_geom",
+        "median, s_geom",
         (
             (0.01, 3.5),
             (0.1, 2.1),
             (1, 1.5),
         ),
     )
-    def test_mean(m_mode, s_geom):
+    def test_mean(median, s_geom):
         # arrange
-        sut = Lognormal(m_mode=m_mode, norm_factor=1, s_geom=s_geom)
-        x = np.linspace(m_mode / 1000, m_mode * 1000, num=10000)
+        sut = Lognormal(median=median, norm_factor=1, s_geom=s_geom)
+        x = np.linspace(median / 1000, median * 1000, num=10000)
 
         # act
         mean = np.sum(sut.pdf(x) * x) / np.sum(sut.pdf(x))
 
         # assert
         np.testing.assert_approx_equal(
-            actual=np.log(m_mode) + 0.5 * np.log(s_geom) ** 2,
+            actual=np.log(median) + 0.5 * np.log(s_geom) ** 2,
             desired=np.log(mean),
             significant=3,
         )
+
+    @staticmethod
+    def test_mode():
+        # arrange
+        spectrum = Lognormal(mode=44, s_geom=1.2, norm_factor=666)
+
+        # assert
+        assert spectrum.mode == 44
+
+    @staticmethod
+    def test_instantiation_passing_mode():
+        # arrange & act
+        spectrum = Lognormal(mode=1, s_geom=1, norm_factor=1)
+
+        # assert
+        np.testing.assert_approx_equal(actual=spectrum.median, desired=np.e)
diff --git a/tests/unit_tests/initialisation/test_spectral_discretisation.py b/tests/unit_tests/initialisation/test_spectral_discretisation.py
index 91b98d8935..398cc48917 100644
--- a/tests/unit_tests/initialisation/test_spectral_discretisation.py
+++ b/tests/unit_tests/initialisation/test_spectral_discretisation.py
@@ -7,10 +7,10 @@
 from PySDM.initialisation import spectra
 from PySDM.physics import si
 
-m_mode = 0.5e-5
+median = 0.5e-5
 n_part = 256 * 16
 s_geom = 1.5
-spectrum = spectra.Lognormal(n_part, m_mode, s_geom)
+spectrum = spectra.Lognormal(norm_factor=n_part, median=median, s_geom=s_geom)
 m_range = (0.1 * si.um, 100 * si.um)
 formulae = Formulae()
 
@@ -91,5 +91,6 @@ def test_error_threshold_with_deterministic_sampling(
         sampling_class, error_threshold
     ):
         sampling_class(
-            spectra.Lognormal(n_part, m_mode, s_geom), error_threshold=error_threshold
+            spectra.Lognormal(norm_factor=n_part, median=median, s_geom=s_geom),
+            error_threshold=error_threshold,
         ).sample_deterministic(n_sd=10)
diff --git a/tests/unit_tests/initialisation/test_spectro_glacial_discretisation.py b/tests/unit_tests/initialisation/test_spectro_glacial_discretisation.py
index 00a72cc85f..0f6080e7ad 100644
--- a/tests/unit_tests/initialisation/test_spectro_glacial_discretisation.py
+++ b/tests/unit_tests/initialisation/test_spectro_glacial_discretisation.py
@@ -13,7 +13,7 @@
 )
 spectrum = Lognormal(
     norm_factor=1,
-    m_mode=formulae.trivia.sphere_surface(diameter=0.74 * si.um),
+    median=formulae.trivia.sphere_surface(diameter=0.74 * si.um),
     s_geom=np.exp(0.25),
 )
 m_range = (
diff --git a/tests/unit_tests/physics/test_spectra.py b/tests/unit_tests/physics/test_spectra.py
index 2321f55406..4cc3c0d30d 100644
--- a/tests/unit_tests/physics/test_spectra.py
+++ b/tests/unit_tests/physics/test_spectra.py
@@ -13,7 +13,7 @@ def test_size_distribution_n_part():
         # Arrange
         s = 1.5
         n_part = 256
-        sut = Lognormal(n_part, 0.5e-5, s)
+        sut = Lognormal(norm_factor=n_part, median=0.5e-5, s_geom=s)
 
         # Act
         m, dm = np.linspace(0.1e-6, 100e-6, 100, retstep=True)
@@ -27,7 +27,7 @@ def test_size_distribution_r_mode():
         # Arrange
         s = 1.001
         r_mode = 1e-6
-        sut = Lognormal(norm_factor=1, m_mode=r_mode, s_geom=s)
+        sut = Lognormal(norm_factor=1, mode=r_mode, s_geom=s)
 
         # Act
         m = np.linspace(start=0.01e-6, stop=100e-6, num=1000)
@@ -69,7 +69,7 @@ class TestSum:
 
     s = 1.001
     r_mode = 1e-6
-    lognormal = Lognormal(1, r_mode, s)
+    lognormal = Lognormal(norm_factor=1, mode=r_mode, s_geom=s)
 
     @staticmethod
     def test_size_distribution():
diff --git a/tests/unit_tests/products/test_activation_criteria.py b/tests/unit_tests/products/test_activation_criteria.py
index 37b4dde0f9..970bdb4214 100644
--- a/tests/unit_tests/products/test_activation_criteria.py
+++ b/tests/unit_tests/products/test_activation_criteria.py
@@ -44,7 +44,7 @@ def test_activation_criteria(backend, plot=False):
     builder.add_dynamic(Condensation())
 
     r_dry, specific_concentration = spectral_sampling.ConstantMultiplicity(
-        Lognormal(norm_factor=1e4 / si.mg, m_mode=50 * si.nm, s_geom=1.5)
+        Lognormal(norm_factor=1e4 / si.mg, median=50 * si.nm, s_geom=1.5)
     ).sample_deterministic(builder.particulator.n_sd)
 
     particulator = builder.build(