diff --git a/cara/monte_carlo/data.py b/cara/monte_carlo/data.py
index bf9c148c..e70a65de 100644
--- a/cara/monte_carlo/data.py
+++ b/cara/monte_carlo/data.py
@@ -3,6 +3,7 @@ import typing
 
 import numpy as np
 from scipy import special as sp
+from scipy.stats import weibull_min
 
 import cara.monte_carlo as mc
 from cara.monte_carlo.sampleable import LogNormal,LogCustomKernel,CustomKernel,Uniform, Custom
@@ -101,13 +102,25 @@ symptomatic_vl_frequencies = LogCustomKernel(
     kernel_bandwidth=0.1
 )
 
+
+# Weibull distribution with a shape factor of 3.47 and a scale factor of 7.01.
+# From https://elifesciences.org/articles/65774 and first line of the figure in
+# https://iiif.elifesciences.org/lax:65774%2Felife-65774-fig4-figsupp3-v2.tif/full/1500,/0/default.jpg
+viral_load = np.linspace(weibull_min.ppf(0.01, c=3.47, scale=7.01),
+                weibull_min.ppf(0.99, c=3.47, scale=7.01), 30)
+frequencies = weibull_min.pdf(viral_load, c=3.47, scale=7.01)
+covid_overal_vl_data = Custom(bounds=(2, 10), function=lambda d: np.interp(d, viral_load, frequencies, right=0., left=0.), max_function=0.16)
+
+
 # Derived from data in doi.org/10.1016/j.ijid.2020.09.025 and
 # https://iosh.com/media/8432/aerosol-infection-risk-hospital-patient-care-full-report.pdf (page 60)
 viable_to_RNA_ratio_distribution = Uniform(0.01, 0.6)
 
+
 # From discussion with virologists
 infectious_dose_distribution = Uniform(10., 100.)
 
+
 # From https://doi.org/10.1101/2021.10.14.21264988 and refererences therein
 virus_distributions = {
     'SARS_CoV_2': mc.SARSCoV2(