diff --git a/cara/models.py b/cara/models.py
index 88e5493a..cf2ffb8f 100644
--- a/cara/models.py
+++ b/cara/models.py
@@ -37,6 +37,7 @@ import typing
 
 import numpy as np
 from scipy.interpolate import interp1d
+import scipy.integrate
 
 if not typing.TYPE_CHECKING:
     from memoization import cached
@@ -477,7 +478,7 @@ class Mask:
     η_inhale: _VectorisedFloat
 
     #: Global factor applied to filtration efficiency of masks when exhaling.
-    factor_exhale: _VectorisedFloat = 1.
+    factor_exhale: float = 1.
 
     #: Pre-populated examples of Masks.
     types: typing.ClassVar[typing.Dict[str, "Mask"]]
@@ -529,31 +530,55 @@ class _ExpirationBase:
     types: typing.ClassVar[typing.Dict[str, "_ExpirationBase"]]
 
     def aerosols(self, mask: Mask):
-        # total volume of aerosols expired (cm^3).
+        """
+        total volume of aerosols expired per volume of air (mL/cm^3).
+        """
         raise NotImplementedError("Subclass must implement")
 
 
 @dataclass(frozen=True)
 class Expiration(_ExpirationBase):
     """
-    Simple model based on four different sizes of particles emitted,
-    with different ejection factors. See Fig. 4 in L. Morawska et al,
-    Size distribution and sites of origin of droplets expelled from the
-    human respiratory tract during expiratory activities,
-    Aerosol Science 40 (2009) pp. 256 - 269.
+    BLO model for the expiration (G. Johnson et al., Modality of human
+    expired aerosol size distributions, Journal of Aerosol Science,
+    vol. 42, no. 12, pp. 839 – 851, 2011,
+    https://doi.org/10.1016/j.jaerosci.2011.07.009).
+    Here all diameters (d) are in microns.
     """
-    ejection_factor: typing.Tuple[float, ...]
-    particle_sizes: typing.Tuple[float, ...] = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4)  # In cm.
+    #: factors assigned to resp. the B, L and O modes. They are
+    # charateristics of the kind of expiratory activity (e.g. breathing,
+    # speaking, singing, or shouting).
+    BLO_factors: typing.Tuple[float, float, float]
 
+    @cached()
     def aerosols(self, mask: Mask):
-        def volume(diameter):
-            return (4 * np.pi * (diameter/2)**3) / 3
-        total = 0
-        for diameter, factor in zip(self.particle_sizes, self.ejection_factor):
-            contribution = (volume(diameter) * factor *
-                            (1 - mask.exhale_efficiency(diameter*1e4)))
-            total += contribution
-        return total
+        """ Result is in mL.cm^-3 """
+        def volume(d):
+            return (np.pi * d**3) / 6.
+
+        def _Bmode(d: float) -> float:
+            # B-mode (see ref. above).
+            return ( (1 / d) * (0.1 / (np.sqrt(2 * np.pi) * 0.262364)) *
+                    np.exp(-1 * (np.log(d) - 0.989541) ** 2 / (2 * 0.262364 ** 2)))
+
+        def _Lmode(d: float) -> float:
+            # L-mode (see ref. above).
+            return ( (1 / d) * (1.0 / (np.sqrt(2 * np.pi) * 0.506818)) *
+                    np.exp(-1 * (np.log(d) - 1.38629) ** 2 / (2 * 0.506818 ** 2)))
+
+        def _Omode(d: float) -> float:
+            # O-mode (see ref. above).
+            return ( (1 / d) * (0.0010008 / (np.sqrt(2 * np.pi) * 0.585005)) *
+                    np.exp(-1 * (np.log(d) - 4.97673) ** 2 / (2 * 0.585005 ** 2)))
+
+        def integrand(d: float) -> float:
+            return (self.BLO_factors[0] * _Bmode(d) +
+                    self.BLO_factors[1] * _Lmode(d) +
+                    self.BLO_factors[2] * _Omode(d)
+                    ) * volume(d) * (1 - mask.exhale_efficiency(d))
+
+        # final result converted from microns^3/cm3 to mL/cm^3
+        return scipy.integrate.quad(integrand, 0.1, 30.)[0]*1e-12
 
 
 @dataclass(frozen=True)
@@ -581,11 +606,11 @@ class MultipleExpiration(_ExpirationBase):
 
 
 _ExpirationBase.types = {
-    'Breathing': Expiration((0.084, 0.009, 0.003, 0.002)),
-    'Whispering': Expiration((0.11, 0.014, 0.004, 0.002)),
-    'Talking': Expiration((0.236, 0.068, 0.007, 0.011)),
-    'Unmodulated Vocalization': Expiration((0.751, 0.139, 0.0139, 0.059)),
-    'Superspreading event': Expiration((np.inf, np.inf, np.inf, np.inf)),
+    'Breathing': Expiration((1., 0., 0.)),
+    'Talking': Expiration((1., 1., 1.)),
+    'Shouting': Expiration((1., 5., 5.)),
+    'Singing': Expiration((1., 5., 5.)),
+    'Superspreading event': Expiration((np.inf, 0., 0.)),
 }
 
 
diff --git a/cara/tests/conftest.py b/cara/tests/conftest.py
index 24983cf1..44b25423 100644
--- a/cara/tests/conftest.py
+++ b/cara/tests/conftest.py
@@ -7,11 +7,9 @@ import pytest
 def baseline_model():
     model = models.ConcentrationModel(
         room=models.Room(volume=75),
-        ventilation=models.SlidingWindow(
-            active=models.PeriodicInterval(period=120, duration=120),
-            inside_temp=models.PiecewiseConstant((0,24),(293,)),
-            outside_temp=models.PiecewiseConstant((0,24),(283,)),
-            window_height=1.6, opening_length=0.6,
+        ventilation=models.AirChange(
+            active=models.SpecificInterval(((0,24),)),
+            air_exch=30.,
         ),
         infected=models.InfectedPopulation(
             number=1,
@@ -19,7 +17,7 @@ def baseline_model():
             presence=models.SpecificInterval(((0, 4), (5, 8))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -30,7 +28,7 @@ def baseline_exposure_model(baseline_model):
     return models.ExposureModel(
         baseline_model,
         exposed=models.Population(
-            number=10,
+            number=1000,
             presence=baseline_model.infected.presence,
             activity=baseline_model.infected.activity,
             mask=baseline_model.infected.mask,
diff --git a/cara/tests/models/test_concentration_model.py b/cara/tests/models/test_concentration_model.py
index f7230430..26707236 100644
--- a/cara/tests/models/test_concentration_model.py
+++ b/cara/tests/models/test_concentration_model.py
@@ -14,7 +14,6 @@ from cara import models
         {'air_change': np.array([100, 120])},
         {'viral_load_in_sputum': np.array([5e8, 1e9])},
         {'quantum_infectious_dose': np.array([50, 20])},
-        {'factor_exhale': np.array([0.92, 0.95])},
     ]
 )
 def test_concentration_model_vectorisation(override_params):
@@ -24,7 +23,6 @@ def test_concentration_model_vectorisation(override_params):
         'air_change': 100,
         'viral_load_in_sputum': 1e9,
         'quantum_infectious_dose': 50,
-        'factor_exhale': 0.95,
     }
     defaults.update(override_params)
 
@@ -36,7 +34,7 @@ def test_concentration_model_vectorisation(override_params):
             number=1,
             presence=always,
             mask=models.Mask(
-                factor_exhale=defaults['factor_exhale'],
+                factor_exhale=0.95,
                 η_inhale=0.3,
             ),
             activity=models.Activity(
@@ -47,9 +45,7 @@ def test_concentration_model_vectorisation(override_params):
                 viral_load_in_sputum=defaults['viral_load_in_sputum'],
                 quantum_infectious_dose=defaults['quantum_infectious_dose'],
             ),
-            expiration=models.Expiration(
-                ejection_factor=(0.084, 0.009, 0.003, 0.002),
-            ),
+            expiration=models.Expiration((1., 0., 0.)),
         )
     )
     concentrations = c_model.concentration(10)
diff --git a/cara/tests/test_expiration.py b/cara/tests/test_expiration.py
index 7638d613..434caf8c 100644
--- a/cara/tests/test_expiration.py
+++ b/cara/tests/test_expiration.py
@@ -9,7 +9,7 @@ from cara import models
 
 def test_multiple_wrong_weight_size():
     weights = (1., 2., 3.)
-    e_base = models.Expiration((0.084, 0.009, 0.003, 0.002))
+    e_base = models.Expiration((1., 0., 0.))
     with pytest.raises(
             ValueError,
             match=re.escape("expirations and weigths should contain the"
@@ -19,11 +19,25 @@ def test_multiple_wrong_weight_size():
 
 
 def test_multiple():
-    weights = (1., 2.)
-    e1 = models.Expiration((0.03, 0.02, 0.01, 0.005))
-    e2 = models.Expiration((0.05, 0.04, 0.03, 0.01))
+    weights = (1., 1.)
+    e1 = models.Expiration((1., 0., 0.))
+    e2 = models.Expiration((4., 5., 5.))
+    e_expected = models.Expiration((2.5, 2.5, 2.5))
     e = models.MultipleExpiration([e1, e2], weights)
-    assert e.aerosols(models.Mask.types['No mask']) == (
-              e1.aerosols(models.Mask.types['No mask'])/3. +
-            2*e2.aerosols(models.Mask.types['No mask'])/3.
-    )
+    mask = models.Mask.types['Type I']
+    npt.assert_almost_equal(e_expected.aerosols(mask), e.aerosols(mask))
+
+
+# expected values obtained from another code
+@pytest.mark.parametrize(
+    "BLO_weights, expected_aerosols",
+    [
+        [(1.,0.,0.), 1.38924e-12],
+        [(1.,1.,1.), 1.07129e-10],
+        [(1.,5.,5.), 5.30088e-10],
+    ],
+)
+def test_expiration_aerosols(BLO_weights, expected_aerosols):
+    mask = models.Mask.types['No mask']
+    e = models.Expiration(BLO_weights)
+    npt.assert_allclose(e.aerosols(mask), expected_aerosols, rtol=1e-4)
diff --git a/cara/tests/test_infected_population.py b/cara/tests/test_infected_population.py
index 8ef94729..7acef7f7 100644
--- a/cara/tests/test_infected_population.py
+++ b/cara/tests/test_infected_population.py
@@ -8,7 +8,6 @@ import cara.models
     "override_params", [
         {'viral_load_in_sputum': np.array([5e8, 1e9])},
         {'quantum_infectious_dose': np.array([50, 20])},
-        {'factor_exhale': np.array([0.92, 0.95])},
         {'exhalation_rate': np.array([0.75, 0.81])},
     ]
 )
@@ -16,7 +15,6 @@ def test_infected_population_vectorisation(override_params):
     defaults = {
         'viral_load_in_sputum': 1e9,
         'quantum_infectious_dose': 50,
-        'factor_exhale': 0.95,
         'exhalation_rate': 0.75,
     }
     defaults.update(override_params)
@@ -26,7 +24,7 @@ def test_infected_population_vectorisation(override_params):
             number=1,
             presence=office_hours,
             mask=cara.models.Mask(
-                factor_exhale=defaults['factor_exhale'],
+                factor_exhale=0.95,
                 η_inhale=0.3,
             ),
             activity=cara.models.Activity(
@@ -37,9 +35,7 @@ def test_infected_population_vectorisation(override_params):
                 viral_load_in_sputum=defaults['viral_load_in_sputum'],
                 quantum_infectious_dose=defaults['quantum_infectious_dose'],
             ),
-            expiration=cara.models.Expiration(
-                ejection_factor=(0.084, 0.009, 0.003, 0.002),
-        ),
+            expiration=cara.models.Expiration((1., 0., 0.)),
     )
     emission_rate = infected.emission_rate(10)
     assert isinstance(emission_rate, np.ndarray)
diff --git a/cara/tests/test_known_quantities.py b/cara/tests/test_known_quantities.py
index ea8a1ad0..a3b1e5db 100644
--- a/cara/tests/test_known_quantities.py
+++ b/cara/tests/test_known_quantities.py
@@ -6,21 +6,12 @@ import cara.models as models
 import cara.data as data
 
 
-def test_no_mask_aerosols(baseline_model):
-    exp = models.Expiration.types['Unmodulated Vocalization']
-    npt.assert_allclose(
-        exp.aerosols(models.Mask.types['No mask']),
-        6.077541e-12,
-        rtol=1e-5,
-    )
-
-
-def test_no_mask_emission_rate(baseline_model):
-    rate = 151.938514
+def test_no_mask_superspeading_emission_rate(baseline_model):
+    expected_rate = 970.
     npt.assert_allclose(
         [baseline_model.infected.emission_rate(t) for t in [0, 1, 4, 4.5, 5, 8, 9]],
-        [0, rate, rate, 0, 0, rate, 0],
-        rtol=1e-5
+        [0, expected_rate, expected_rate, 0, 0, expected_rate, 0],
+        rtol=1e-12
     )
 
 
@@ -48,23 +39,24 @@ def baseline_periodic_hepa():
 
 
 def test_concentrations(baseline_model):
+    # expected concentrations were computed analytically
     ts = [0, 4, 5, 7, 10]
     concentrations = [baseline_model.concentration(t) for t in ts]
     npt.assert_allclose(
         concentrations,
-        [0.000000e+00, 2.619538e-01, 1.146999e-04, 2.619537e-01, 5.022289e-08],
-        rtol=1e-5
+        [0.000000e+00, 0.4189594, 1.6422648e-14, 0.4189594, 6.4374587e-28],
+        rtol=1e-6
     )
 
 
 def test_smooth_concentrations(baseline_model):
     # We don't care about the actual concentrations in this test, but rather
     # that the curve itself is smooth.
-    dx = 0.1
-    dy_limit = dx * 2  # Anything more than this is a bit steep.
+    dx = 0.002
+    dy_limit = 0.2  # Anything more than this (in relative) is a bit steep.
     ts = np.arange(0, 10, dx)
     concentrations = [baseline_model.concentration(t) for t in ts]
-    assert np.abs(np.diff(concentrations)).max() < dy_limit
+    assert np.abs(np.diff(concentrations)).max()/np.mean(concentrations) < dy_limit
 
 
 def build_model(interval_duration):
@@ -80,7 +72,7 @@ def build_model(interval_duration):
             presence=models.SpecificInterval(((0, 4), (5, 8))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -98,8 +90,8 @@ def test_concentrations_startup(baseline_model):
 def test_r0(baseline_exposure_model):
     # expected r0 was computed with a trapezoidal integration, using
     # a mesh of 100'000 pts per exposed presence interval.
-    p = baseline_exposure_model.infection_probability()
-    npt.assert_allclose(p, 89.001748)
+    r0 = baseline_exposure_model.reproduction_number()
+    npt.assert_allclose(r0, 973.535888)
 
 
 def test_periodic_window(baseline_periodic_window, baseline_room):
@@ -183,15 +175,6 @@ def test_multiple_ventilation_HEPA_HVAC_AirChange(volume, expected_value):
                         expected_value,rtol=1e-5)
 
 
-def test_expiration_aerosols():
-    mask = models.Mask.types['Type I']
-    exp1 = models.Expiration((0.751, 0.139, 0.0139, 0.059),
-                            particle_sizes = (0.8e-4, 1.8e-4, 3.5e-4, 5.5e-4))
-    exp2 = models.Expiration((0.059, 0.0139, 0.751, 0.139),
-                            particle_sizes = (5.5e-4, 3.5e-4, 0.8e-4, 1.8e-4))
-    npt.assert_allclose(exp1.aerosols(mask), exp2.aerosols(mask), rtol=1e-5)
-
-
 @pytest.mark.parametrize(
     "time, expected_value",
     [
@@ -239,7 +222,7 @@ def build_hourly_dependent_model(month, intervals_open=((7.5, 8.5),),
             presence=models.SpecificInterval(intervals_presence_infected),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -260,7 +243,7 @@ def build_constant_temp_model(outside_temp, intervals_open=((7.5, 8.5),)):
             presence=models.SpecificInterval(((0, 4), (5, 7.5))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -287,7 +270,7 @@ def build_hourly_dependent_model_multipleventilation(month, intervals_open=((7.5
             presence=models.SpecificInterval(((0, 4), (5, 7.5))),
             mask=models.Mask.types['No mask'],
             activity=models.Activity.types['Light activity'],
-            expiration=models.Expiration.types['Unmodulated Vocalization'],
+            expiration=models.Expiration.types['Superspreading event'],
         ),
     )
     return model
@@ -370,16 +353,16 @@ def build_exposure_model(concentration_model):
     )
 
 
-# expected r0 were computed with a trapezoidal integration, using
+# expected quanta were computed with a trapezoidal integration, using
 # a mesh of 100'000 pts per exposed presence interval.
 @pytest.mark.parametrize(
-    "month, expected_r0",
+    "month, expected_quanta",
     [
-        ['Jan', 86.258533],
-        ['Jun', 99.972103],
+        ['Jan', 10.136783],
+        ['Jun', 41.800377],
     ],
 )
-def test_r0_hourly_dep(month,expected_r0):
+def test_quanta_hourly_dep(month,expected_quanta):
     m = build_exposure_model(
         build_hourly_dependent_model(
             month,
@@ -387,19 +370,19 @@ def test_r0_hourly_dep(month,expected_r0):
             intervals_presence_infected=((8, 12), (13, 17))
         )
     )
-    p = m.infection_probability()
-    npt.assert_allclose(p, expected_r0)
+    quanta = m.quanta_exposure()
+    npt.assert_allclose(quanta, expected_quanta)
 
-# expected r0 were computed with a trapezoidal integration, using
+# expected quanta were computed with a trapezoidal integration, using
 # a mesh of 100'000 pts per exposed presence interval.
 @pytest.mark.parametrize(
-    "month, expected_r0",
+    "month, expected_quanta",
     [
-        ['Jan', 86.422736],
-        ['Jun', 99.982323],
+        ['Jan', 10.19818],
+        ['Jun', 44.130683],
     ],
 )
-def test_r0_hourly_dep_refined(month,expected_r0):
+def test_quanta_hourly_dep_refined(month,expected_quanta):
     m = build_exposure_model(
         build_hourly_dependent_model(
             month,
@@ -408,5 +391,5 @@ def test_r0_hourly_dep_refined(month,expected_r0):
             temperatures=data.GenevaTemperatures,
         )
     )
-    p = m.infection_probability()
-    npt.assert_allclose(p, expected_r0)
+    quanta = m.quanta_exposure()
+    npt.assert_allclose(quanta, expected_quanta)
diff --git a/cara/tests/test_monte_carlo.py b/cara/tests/test_monte_carlo.py
index f3d9fb08..c3bfc2fe 100644
--- a/cara/tests/test_monte_carlo.py
+++ b/cara/tests/test_monte_carlo.py
@@ -53,7 +53,7 @@ def baseline_mc_model() -> cara.monte_carlo.ConcentrationModel:
             presence=cara.models.SpecificInterval(((0, 4), (5, 8))),
             mask=cara.models.Mask.types['No mask'],
             activity=cara.models.Activity.types['Light activity'],
-            expiration=cara.models.Expiration.types['Unmodulated Vocalization'],
+            expiration=cara.models.Expiration.types['Breathing'],
         ),
     )
     return mc_model