diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
index 185d9c1570..5030c1a3c0 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2020 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,6 +27,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "RosinRammler.H"
+#include "MathFunctions.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -48,11 +50,29 @@ Foam::distributionModels::RosinRammler::RosinRammler
 )
 :
     distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue")),
-    d_(distributionModelDict_.get<scalar>("d")),
+    lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
     n_(distributionModelDict_.get<scalar>("n"))
 {
+    const word parcelBasisType =
+        dict.getOrDefault<word>("parcelBasisType", "none");
+
+    if (parcelBasisType == "mass")
+    {
+        WarningInFunction
+            << "Selected parcel basis type: " << parcelBasisType << nl
+            << "    Please consider to use massRosinRammler distribution model"
+            << endl;
+    }
+
+    if (lambda_ < VSMALL || n_ < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Scale/Shape parameter cannot be equal to or less than zero:"
+            << "    lambda = " << lambda_
+            << "    n = " << n_
+            << exit(FatalError);
+    }
+
     check();
 }
 
@@ -60,45 +80,33 @@ Foam::distributionModels::RosinRammler::RosinRammler
 Foam::distributionModels::RosinRammler::RosinRammler(const RosinRammler& p)
 :
     distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
-    d_(p.d_),
+    lambda_(p.lambda_),
     n_(p.n_)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::RosinRammler::~RosinRammler()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::distributionModels::RosinRammler::sample() const
 {
-    const scalar minValueByDPowN = pow(minValue_/d_, n_);
-    const scalar K = 1 - exp(- pow(maxValue_/d_, n_) + minValueByDPowN);
-    const scalar y = rndGen_.sample01<scalar>();
-    return d_*pow(minValueByDPowN - log(1 - K*y), 1/n_);
-}
-
-
-Foam::scalar Foam::distributionModels::RosinRammler::minValue() const
-{
-    return minValue_;
-}
-
-
-Foam::scalar Foam::distributionModels::RosinRammler::maxValue() const
-{
-    return maxValue_;
+    const scalar u = rndGen_.sample01<scalar>();
+    const scalar qMin = pow(minValue_/lambda_, n_);
+    const scalar qMax = pow(maxValue_/lambda_, n_);
+    const scalar r = scalar(1) - exp(-qMax + qMin);
+    return lambda_*pow(qMin - log(scalar(1) - u*r), scalar(1)/n_);
 }
 
 
 Foam::scalar Foam::distributionModels::RosinRammler::meanValue() const
 {
-    return d_;
+    // (C:Eq. 5)
+    const scalar a = scalar(1)/lambda_ + scalar(1);
+    const scalar qMax = pow(maxValue_/n_, lambda_);
+    const scalar qMin = pow(minValue_/n_, lambda_);
+    const scalar gMax = Math::incGamma_P(a, qMax);
+    const scalar gMin = Math::incGamma_P(a, qMin);
+
+    return n_/(exp(-qMin) - exp(-qMax))*(gMax - gMin);
 }
 
 
diff --git a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
index ec4f94aa0f..1cdfa532fc 100644
--- a/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
+++ b/src/lagrangian/distributionModels/RosinRammler/RosinRammler.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2020 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,13 +28,107 @@ Class
     Foam::distributionModels::RosinRammler
 
 Description
-    Rosin-Rammler distributionModel
+    Particle-size distribution model wherein random samples are drawn
+    from the doubly-truncated two-parameter Rosin-Rammler (Weibull)
+    probability density function:
 
-   \f[
-       CDF(x) =
-            (1 - exp(-(x/d)^n + (d_0/d)^n)
-           /(1 - exp(-(d_1/d)^n + (d_0/d)^n)
-   \f]
+    \f[
+        f(x; \lambda, n, A, B) =
+            \frac{
+                \frac{n}{\lambda}
+                \left( \frac{x}{\lambda} \right)^{n-1}
+                \exp\{ -(\frac{x}{\lambda})^n \}
+            }{
+                \exp\{- (\frac{A}{\lambda})^n \}
+              - \exp\{- (\frac{B}{\lambda})^n \}
+            }
+    \f]
+    where
+
+    \vartable
+        f(x; \lambda, n, A, B) | Rosin-Rammler probability density function
+        \lambda   | Scale parameter
+        n         | Shape parameter
+        x         | Sample
+        A         | Minimum of the distribution
+        B         | Maximum of the distribution
+    \endvartable
+
+    Constraints:
+    - \f$ \infty > B > A > 0\f$
+    - \f$ x \in [B,A] \f$
+    - \f$ \lambda > 0 \f$
+    - \f$ n > 0 \f$
+
+    Random samples are generated by the inverse transform sampling technique
+    by using the quantile function of the doubly-truncated two-parameter
+    Rosin-Rammler (Weibull) probability density function:
+
+    \f[
+        x = \lambda \left( q_{min} - \ln (1 - u r) \right)^{1/n}
+    \f]
+    with
+
+    \f[
+        r = 1 - \exp(-q_{max} + q_{min})
+    \f]
+
+    \f[
+        q_{min} = \left(\frac{A}{\lambda}\right)^n
+    \f]
+
+    \f[
+        q_{max} = \left(\frac{B}{\lambda}\right)^n
+    \f]
+    where \f$ u \f$ is sample drawn from the uniform probability
+    density function on the unit interval \f$ (0, 1) \f$.
+
+    Reference:
+    \verbatim
+    Doubly-truncated two-parameter Weibull distribution and its moments (tag:C):
+        Crénin, F. (2015).
+        Truncated Weibull Distribution Functions and Moments.
+        SSRN 2690255.
+        DOI:10.2139/ssrn.2690255
+    \endverbatim
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        RosinRammler;
+                    RosinRammlerDistribution
+                    {
+                        lambda      <scaleParameterValue>;
+                        n           <shapeParameterValue>;
+                        minValue    <minValue>;
+                        maxValue    <maxValue>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: RosinRammler           | word | yes  | -
+      RosinRammlerDistribution | Distribution settings | dict | yes  | -
+      lambda       | Scale parameter                   | scalar | yes  | -
+      n            | Shape parameter                   | scalar | yes  | -
+      minValue     | Minimum of the distribution       | scalar | yes  | -
+      maxValue     | Maximum of the distribution       | scalar | yes  | -
+    \endtable
 
 SourceFiles
     RosinRammler.C
@@ -62,19 +157,11 @@ class RosinRammler
 {
     // Private Data
 
-        //- Distribution minimum
-        scalar minValue_;
+        //- Scale parameter
+        scalar lambda_;
 
-        //- Distribution maximum
-        scalar maxValue_;
-
-        // Model coefficients
-
-            //- Scale parameter
-            scalar d_;
-
-            //- Shape parameter
-            scalar n_;
+        //- Shape parameter
+        scalar n_;
 
 
 public:
@@ -88,7 +175,7 @@ public:
         //- Construct from components
         RosinRammler(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         RosinRammler(const RosinRammler& p);
 
         //- Construct and return a clone
@@ -97,23 +184,20 @@ public:
             return autoPtr<distributionModel>(new RosinRammler(*this));
         }
 
+        //- No copy assignment
+        void operator=(const RosinRammler&) = delete;
+
 
     //- Destructor
-    virtual ~RosinRammler();
+    virtual ~RosinRammler() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/binned/binned.C b/src/lagrangian/distributionModels/binned/binned.C
index 6393ed3846..36d7370387 100644
--- a/src/lagrangian/distributionModels/binned/binned.C
+++ b/src/lagrangian/distributionModels/binned/binned.C
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2016 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,11 +32,11 @@ License
 
 namespace Foam
 {
-    namespace distributionModels
-    {
-        defineTypeNameAndDebug(binned, 0);
-        addToRunTimeSelectionTable(distributionModel, binned, dictionary);
-    }
+namespace distributionModels
+{
+    defineTypeNameAndDebug(binned, 0);
+    addToRunTimeSelectionTable(distributionModel, binned, dictionary);
+}
 }
 
 
@@ -58,6 +58,16 @@ void Foam::distributionModels::binned::initialise()
 
     // Normalise
     scalar sumProb = xy_.last()[1];
+
+    if (sumProb < VSMALL)
+    {
+        FatalErrorInFunction
+            << type() << "distribution: "
+            << "The sum of elements in the second column cannot be zero." << nl
+            << "sum = " << sumProb
+            << exit(FatalError);
+    }
+
     forAll(xy_, bini)
     {
         xy_[bini][1] /= sumProb;
@@ -90,14 +100,10 @@ Foam::distributionModels::binned::binned
     xy_(distributionModelDict_.lookup("distribution")),
     meanValue_(0)
 {
-    if (maxValue() < minValue())
-    {
-        FatalErrorInFunction
-            << "Maximum value is smaller than the minimum value:"
-            << "    maxValue = " << maxValue()
-            << ", minValue = " << minValue()
-            << exit(FatalError);
-    }
+    minValue_ = xy_[0][0];
+    maxValue_ = xy_[xy_.size()-1][0];
+
+    check();
 
     initialise();
 }
@@ -114,16 +120,16 @@ Foam::distributionModels::binned::binned
     xy_(),
     meanValue_(0)
 {
-    scalar minValue = GREAT;
-    scalar maxValue = -GREAT;
+    minValue_ = GREAT;
+    maxValue_ = -GREAT;
     forAll(sampleData, i)
     {
-        minValue = min(minValue, sampleData[i]);
-        maxValue = max(maxValue, sampleData[i]);
+        minValue_ = min(minValue_, sampleData[i]);
+        maxValue_ = max(maxValue_, sampleData[i]);
     }
 
-    const label bin0 = floor(minValue/binWidth);
-    const label bin1 = ceil(maxValue/binWidth);
+    const label bin0 = floor(minValue_/binWidth);
+    const label bin1 = ceil(maxValue_/binWidth);
     const label nBin = bin1 - bin0;
 
     if (nBin == 0)
@@ -177,39 +183,21 @@ Foam::distributionModels::binned::binned(const binned& p)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::binned::~binned()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::distributionModels::binned::sample() const
 {
-    scalar y = rndGen_.sample01<scalar>();
+    const scalar u = rndGen_.sample01<scalar>();
 
     for (label i = 0; i < xy_.size() - 1; ++i)
     {
-        if (xy_[i][1] > y)
+        if (xy_[i][1] > u)
         {
             return xy_[i][0];
         }
     }
 
-    return xy_.last()[0];
-}
-
-
-Foam::scalar Foam::distributionModels::binned::minValue() const
-{
-    return xy_.first()[0];
-}
-
-
-Foam::scalar Foam::distributionModels::binned::maxValue() const
-{
-    return xy_.last()[0];
+    return maxValue_;
 }
 
 
diff --git a/src/lagrangian/distributionModels/binned/binned.H b/src/lagrangian/distributionModels/binned/binned.H
index 63d88dc9f1..20e36899ab 100644
--- a/src/lagrangian/distributionModels/binned/binned.H
+++ b/src/lagrangian/distributionModels/binned/binned.H
@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2015-2016 OpenCFD Ltd.
+    Copyright (C) 2015-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +27,57 @@ Class
     Foam::distributionModels::binned
 
 Description
-    A binned distribution model where the random sample is taken from a
-    discrete set of bins
+    Particle-size distribution model wherein random samples are
+    drawn from a given discrete set of (\c bin, \c probability) pairs,
+    where in terms of its meaning, \c bins correspond to particle sizes
+    and \c probabilities correspond to (relative) probability of occurrences.
+
+    The second column (i.e. \c probability) are normalised by the sum of all
+    its values, resulting in a normalised column in the domain of [0,1].
+    To generate a sample, first a sample drawn from the uniform probability
+    density function on the unit interval (i.e. \c u), and then, the \c bin
+    corresponding to the first \c probability larger than \c u is fetched
+    as the particle size to be further processed.
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        binned;
+                    binnedDistribution
+                    {
+                        distribution
+                        (
+                            (<bin1>    <probability1>)
+                            (<bin2>    <probability2>)
+                            ...
+                            (<binN>    <probabilityN>)
+                        );
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: binned                 | word | yes  | -
+      binnedDistribution | Distribution settings       | dict | yes  | -
+      distribution | \<bin\>-\<probability\> pairs     | dict | yes  | -
+      \<bin\>      | Particle size                     | scalar | yes | -
+      \<probability\> | Probability of occurrence      | scalar | yes | -
+    \endtable
 
 SourceFiles
     binned.C
@@ -70,18 +119,18 @@ class binned
 :
     public distributionModel
 {
-    // Private data
+    typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
 
-        typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
+    // Private Data
 
-        // List of (bin probability)
+        // List of (bin probability) pairs
         List<pair> xy_;
 
-        //- Distribution mean value
+        //- Mean of the distribution
         scalar meanValue_;
 
 
-    // Private member functions
+    // Private Member Functions
 
         //- Initialise the distribution parameters
         void initialise();
@@ -101,6 +150,7 @@ public:
         binned(const dictionary& dict, Random& rndGen);
 
         //- Construct from components
+        //  Allows negative entries
         binned
         (
             const UList<scalar>& sampleData,
@@ -108,7 +158,7 @@ public:
             Random& rndGen
         );
 
-        //- Construct copy
+        //- Copy construct
         binned(const binned& p);
 
         //- Construct and return a clone
@@ -117,23 +167,20 @@ public:
             return autoPtr<distributionModel>(new binned(*this));
         }
 
+        //- No copy assignment
+        void operator=(const binned&) = delete;
+
 
     //- Destructor
-    virtual ~binned();
+    virtual ~binned() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the arithmetic mean of the distribution data
         virtual scalar meanValue() const;
 
         //- Write data to stream
diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.C b/src/lagrangian/distributionModels/distributionModel/distributionModel.C
index 4817bd026c..93d13e802a 100644
--- a/src/lagrangian/distributionModels/distributionModel/distributionModel.C
+++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -43,19 +44,31 @@ void Foam::distributionModel::check() const
     if (minValue() < 0)
     {
         FatalErrorInFunction
-            << type() << "distribution: Minimum value must be greater than "
-            << "zero." << nl << "Supplied minValue = " << minValue()
+            << type() << "Distribution: "
+            << "Minimum value must be greater than zero." << nl
+            << "Supplied minValue = " << minValue()
             << abort(FatalError);
     }
 
     if (maxValue() < minValue())
     {
         FatalErrorInFunction
-            << type() << "distribution: Maximum value is smaller than the "
-            << "minimum value:" << nl << "    maxValue = " << maxValue()
+            << type() << "Distribution: "
+            << "Maximum value cannot be smaller than minimum value" << nl
+            << "    maxValue = " << maxValue()
             << ", minValue = " << minValue()
             << abort(FatalError);
     }
+
+    if (maxValue() == minValue())
+    {
+        WarningInFunction
+            << type() << "Distribution: "
+            << "Maximum and minimum values are equal to each other" << nl
+            << "    maxValue = " << maxValue()
+            << ", minValue = " << minValue()
+            << endl;
+    }
 }
 
 
@@ -69,7 +82,9 @@ Foam::distributionModel::distributionModel
 )
 :
     distributionModelDict_(dict),
-    rndGen_(rndGen)
+    rndGen_(rndGen),
+    minValue_(distributionModelDict_.getOrDefault<scalar>("minValue", GREAT)),
+    maxValue_(distributionModelDict_.getOrDefault<scalar>("maxValue", -GREAT))
 {}
 
 
@@ -79,14 +94,24 @@ Foam::distributionModel::distributionModel
 )
 :
     distributionModelDict_(p.distributionModelDict_),
-    rndGen_(p.rndGen_)
+    rndGen_(p.rndGen_),
+    minValue_(p.minValue_),
+    maxValue_(p.maxValue_)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::distributionModel::~distributionModel()
-{}
+Foam::scalar Foam::distributionModel::minValue() const
+{
+    return minValue_;
+}
+
+
+Foam::scalar Foam::distributionModel::maxValue() const
+{
+    return maxValue_;
+}
 
 
 // ************************************************************************* //
diff --git a/src/lagrangian/distributionModels/distributionModel/distributionModel.H b/src/lagrangian/distributionModels/distributionModel/distributionModel.H
index c866d972fe..dc4215faf5 100644
--- a/src/lagrangian/distributionModels/distributionModel/distributionModel.H
+++ b/src/lagrangian/distributionModels/distributionModel/distributionModel.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,15 +24,21 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
+Namespace
+    Foam::distributionModels
+
+Description
+    A namespace for various probability distribution model implementations.
+
 Class
     Foam::distributionModel
 
 Description
-    A library of runtime-selectable distribution models.
+    A library of runtime-selectable doubly-truncated probability distribution
+    models. Returns random samples based on given distribution parameters.
 
-    Returns a sampled value given the expectation (nu) and variance (sigma^2)
-
-    Current distribution models include:
+    Available distribution models include:
+    - binned
     - exponential
     - fixedValue
     - general
@@ -41,11 +48,6 @@ Description
     - Mass-based Rosin-Rammler
     - uniform
 
-    The distributionModel is tabulated in equidistant nPoints, in an interval.
-    These values are integrated to obtain the cumulated distribution model,
-    which is then used to change the distribution from unifrom to
-    the actual distributionModel.
-
 SourceFiles
     distributionModel.C
     distributionModelNew.C
@@ -70,10 +72,9 @@ namespace Foam
 
 class distributionModel
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Coefficients dictionary
         const dictionary distributionModelDict_;
@@ -81,6 +82,12 @@ protected:
         //- Reference to the random number generator
         Random& rndGen_;
 
+        //- Minimum of the distribution
+        scalar minValue_;
+
+        //- Maximum of the distribution
+        scalar maxValue_;
+
 
     // Protected Member Functions
 
@@ -118,7 +125,7 @@ public:
             Random& rndGen
         );
 
-        //- Construct copy
+        //- Copy construct
         distributionModel(const distributionModel& p);
 
         //- Construct and return a clone
@@ -134,21 +141,22 @@ public:
 
 
     //- Destructor
-    virtual ~distributionModel();
+    virtual ~distributionModel() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const = 0;
 
-        //- Return the minimum value
-        virtual scalar minValue() const = 0;
+        //- Return the minimum of the distribution
+        virtual scalar minValue() const;
 
-        //- Return the maximum value
-        virtual scalar maxValue() const = 0;
+        //- Return the maximum of the distribution
+        virtual scalar maxValue() const;
 
-        //- Return the maximum value
+        //- Return the theoretical mean of the distribution, or
+        //- the arithmetic mean of the distribution data
         virtual scalar meanValue() const = 0;
 };
 
diff --git a/src/lagrangian/distributionModels/exponential/exponential.C b/src/lagrangian/distributionModels/exponential/exponential.C
index 7586507f86..96d3df2c15 100644
--- a/src/lagrangian/distributionModels/exponential/exponential.C
+++ b/src/lagrangian/distributionModels/exponential/exponential.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,10 +49,16 @@ Foam::distributionModels::exponential::exponential
 )
 :
     distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue")),
     lambda_(distributionModelDict_.get<scalar>("lambda"))
 {
+    if (lambda_ < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Rate parameter cannot be equal to or less than zero:" << nl
+            << "    lambda = " << lambda_
+            << exit(FatalError);
+    }
+
     check();
 }
 
@@ -59,43 +66,24 @@ Foam::distributionModels::exponential::exponential
 Foam::distributionModels::exponential::exponential(const exponential& p)
 :
     distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
     lambda_(p.lambda_)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::exponential::~exponential()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::distributionModels::exponential::sample() const
 {
-    scalar y = rndGen_.sample01<scalar>();
-    scalar K = exp(-lambda_*maxValue_) - exp(-lambda_*minValue_);
-    return -(1.0/lambda_)*log(exp(-lambda_*minValue_) + y*K);
-}
-
-
-Foam::scalar Foam::distributionModels::exponential::minValue() const
-{
-    return minValue_;
-}
-
-
-Foam::scalar Foam::distributionModels::exponential::maxValue() const
-{
-    return maxValue_;
+    const scalar u = rndGen_.sample01<scalar>();
+    const scalar qMin = exp(-lambda_*minValue_);
+    const scalar qMax = exp(-lambda_*maxValue_);
+    return -(scalar(1)/lambda_)*log(qMin + u*(qMax - qMin));
 }
 
 
 Foam::scalar Foam::distributionModels::exponential::meanValue() const
 {
-    return 1.0/lambda_;
+    return scalar(1)/lambda_;
 }
 
 
diff --git a/src/lagrangian/distributionModels/exponential/exponential.H b/src/lagrangian/distributionModels/exponential/exponential.H
index a9777a65c3..ca35bebe66 100644
--- a/src/lagrangian/distributionModels/exponential/exponential.H
+++ b/src/lagrangian/distributionModels/exponential/exponential.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,87 @@ Class
     Foam::distributionModels::exponential
 
 Description
-    exponential distribution model
+    Particle-size distribution model wherein random samples are drawn
+    from the doubly-truncated exponential probability density function:
+
+    \f[
+        f(x; \lambda, A, B) =
+            \lambda \frac{\exp(-\lambda (x - A))}{1 - \exp(-\lambda(B-A))}
+    \f]
+    where
+
+    \vartable
+        f(x; \lambda, A, B) | Exponential probability density function
+        \lambda   | Rate parameter
+        x         | Sample
+        A         | Minimum of the distribution
+        B         | Maximum of the distribution
+    \endvartable
+
+    Constraints:
+    - \f$ \infty > B > A > 0 \f$
+    - \f$ x \in [B,A] \f$
+    - \f$ \lambda > 0 \f$
+
+    Random samples are generated by the inverse transform sampling technique
+    by using the quantile function of the doubly-truncated exponential
+    probability density function:
+
+    \f[
+        x = - \frac{1}{\lambda} \ln (r)
+    \f]
+
+    with
+
+    \f[
+        r = q_{min} + u (q_{max} - q_{min})
+    \f]
+
+    \f[
+        q_{min} = \exp(-\lambda A)
+    \f]
+
+    \f[
+        q_{max} = \exp(-\lambda B)
+    \f]
+    where \f$ u \f$ is a sample drawn from the uniform probability
+    density function on the unit interval \f$ (0, 1) \f$.
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        exponential;
+                    exponentialDistribution
+                    {
+                        lambda      <lambdaValue>;
+                        minValue    <minValue>;
+                        maxValue    <maxValue>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: exponential            | word | yes  | -
+      exponentialDistribution | Distribution settings  | dict | yes  | -
+      lambda       | Rate parameter                    | scalar | yes  | -
+      minValue     | Minimum of the distribution       | scalar | yes  | -
+      maxValue     | Maximum of the distribution       | scalar | yes  | -
+    \endtable
 
 SourceFiles
     exponential.C
@@ -54,18 +135,10 @@ class exponential
 :
     public distributionModel
 {
-    // Private data
+    // Private Data
 
-        //- Distribution minimum
-        scalar minValue_;
-
-        //- Distribution maximum
-        scalar maxValue_;
-
-
-        // Model coefficients
-
-            scalar lambda_;
+        //- Rate parameter
+        scalar lambda_;
 
 
 public:
@@ -79,7 +152,7 @@ public:
         //- Construct from components
         exponential(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         exponential(const exponential& p);
 
         //- Construct and return a clone
@@ -88,23 +161,20 @@ public:
             return autoPtr<distributionModel>(new exponential(*this));
         }
 
+        //- No copy assignment
+        void operator=(const exponential&) = delete;
+
 
     //- Destructor
-    virtual ~exponential();
+    virtual ~exponential() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/fixedValue/fixedValue.C b/src/lagrangian/distributionModels/fixedValue/fixedValue.C
index ab7df3a950..9cc2f8dbb3 100644
--- a/src/lagrangian/distributionModels/fixedValue/fixedValue.C
+++ b/src/lagrangian/distributionModels/fixedValue/fixedValue.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -49,7 +50,15 @@ Foam::distributionModels::fixedValue::fixedValue
 :
     distributionModel(typeName, dict, rndGen),
     value_(distributionModelDict_.get<scalar>("value"))
-{}
+{
+    if (value_ < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Fixed value cannot be equal to or less than zero:" << nl
+            << "    value = " << value_
+            << exit(FatalError);
+    }
+}
 
 
 Foam::distributionModels::fixedValue::fixedValue(const fixedValue& p)
@@ -59,12 +68,6 @@ Foam::distributionModels::fixedValue::fixedValue(const fixedValue& p)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::fixedValue::~fixedValue()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::distributionModels::fixedValue::fixedValue::sample() const
diff --git a/src/lagrangian/distributionModels/fixedValue/fixedValue.H b/src/lagrangian/distributionModels/fixedValue/fixedValue.H
index 9b07b2670e..d8437df290 100644
--- a/src/lagrangian/distributionModels/fixedValue/fixedValue.H
+++ b/src/lagrangian/distributionModels/fixedValue/fixedValue.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,39 @@ Class
     Foam::distributionModels::fixedValue
 
 Description
-    Returns a fixed value
+    Particle-size distribution model wherein samples are given fixed values.
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        fixedValue;
+                    fixedValueDistribution
+                    {
+                        value   <value>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: fixedValue             | word | yes  | -
+      fixedValueDistribution | Distribution settings   | dict | yes  | -
+      value        | Fixed value for size              | scalar | yes | -
+    \endtable
 
 SourceFiles
     fixedValue.C
@@ -53,9 +86,9 @@ class fixedValue
 :
     public distributionModel
 {
-    // Private data
+    // Private Data
 
-        //- Fixed value
+        //- Fixed value for size
         scalar value_;
 
 
@@ -70,7 +103,7 @@ public:
         //- Construct from components
         fixedValue(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         fixedValue(const fixedValue& p);
 
         //- Construct and return a clone
@@ -79,23 +112,26 @@ public:
             return autoPtr<distributionModel>(new fixedValue(*this));
         }
 
+        //- No copy assignment
+        void operator=(const fixedValue&) = delete;
+
 
     //- Destructor
-    virtual ~fixedValue();
+    virtual ~fixedValue() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
+        //- Return the minimum of the distribution
         virtual scalar minValue() const;
 
-        //- Return the maximum value
+        //- Return the maximum of the distribution
         virtual scalar maxValue() const;
 
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/general/general.C b/src/lagrangian/distributionModels/general/general.C
index 08f2daf01d..04ce24fefa 100644
--- a/src/lagrangian/distributionModels/general/general.C
+++ b/src/lagrangian/distributionModels/general/general.C
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2011-2019 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,13 +47,14 @@ void Foam::distributionModels::general::initialise()
 {
     static scalar eps = ROOTVSMALL;
 
-    const label nEntries = xy_.size();
+    integral_.setSize(nEntries_);
 
-    integral_.setSize(nEntries);
+    // Fill out the integral table (x, P(x<=0)) and calculate mean
+    // For density function: P(x<=0) = int f(x) and mean = int x*f(x)
+    // For cumulative function: mean = int 1-P(x<=0) = maxValue_ - int P(x<=0)
+    integral_[0] = 0;
 
-    // Normalize the cumulative distribution
-    integral_[0] = 0.0;
-    for (label i = 1; i < nEntries; i++)
+    for (label i = 1; i < nEntries_; ++i)
     {
         scalar k = (xy_[i][1] - xy_[i-1][1])/(xy_[i][0] - xy_[i-1][0] + eps);
         scalar d = xy_[i-1][1] - k*xy_[i-1][0];
@@ -61,18 +62,32 @@ void Foam::distributionModels::general::initialise()
         scalar y0 = xy_[i-1][0]*(0.5*k*xy_[i-1][0] + d);
         scalar area = y1 - y0;
 
-        integral_[i] = area + integral_[i-1];
+        if (cumulative_)
+        {
+            integral_[i] = xy_[i][1];
+            meanValue_ += area;
+        }
+        else
+        {
+            integral_[i] = area + integral_[i-1];
+
+            y1 = sqr(xy_[i][0])*(1.0/3.0*k*xy_[i][0] + 0.5*d);
+            y0 = sqr(xy_[i-1][0])*(1.0/3.0*k*xy_[i-1][0] + 0.5*d);
+            meanValue_ += y1 - y0;
+        }
     }
 
+    // normalize the distribution
     scalar sumArea = integral_.last();
 
-    meanValue_ = sumArea/(maxValue() - minValue() + eps);
-
-    for (label i=0; i < nEntries; i++)
+    for (label i = 0; i < nEntries_; ++i)
     {
         xy_[i][1] /= sumArea + eps;
         integral_[i] /= sumArea + eps;
     }
+
+    meanValue_ /= sumArea;
+    meanValue_ = cumulative_ ? (maxValue_ - meanValue_ + eps) : meanValue_;
 }
 
 
@@ -86,11 +101,64 @@ Foam::distributionModels::general::general
 :
     distributionModel(typeName, dict, rndGen),
     xy_(distributionModelDict_.lookup("distribution")),
-    meanValue_(0.0),
-    integral_()
+    nEntries_(xy_.size()),
+    meanValue_(0),
+    integral_(nEntries_),
+    cumulative_(distributionModelDict_.getOrDefault("cumulative", false))
 {
+    minValue_ = xy_[0][0];
+    maxValue_ = xy_[nEntries_-1][0];
+
     check();
 
+    // Additional sanity checks
+    if (cumulative_ && xy_[0][1] != 0)
+    {
+        FatalErrorInFunction
+            << type() << "distribution: "
+            << "Elements in the second column for cumulative "
+            << "distribution functions must start from zero." << nl
+            << "First element = " << xy_[0][1]
+            << exit(FatalError);
+    }
+
+    for (label i = 0; i < nEntries_; ++i)
+    {
+        if (i > 0 && xy_[i][0] <= xy_[i-1][0])
+        {
+            FatalErrorInFunction
+                << type() << "distribution: "
+                << "Elements in the first column must "
+                << "be specified in an ascending order." << nl
+                << "Please see the row i = " << i << nl
+                << "xy[i-1] = " << xy_[i-1] << nl
+                << "xy[i] = " << xy_[i]
+                << exit(FatalError);
+        }
+
+        if (xy_[i][0] < 0 || xy_[i][1] < 0)
+        {
+            FatalErrorInFunction
+                << type() << "distribution: "
+                << "Input pairs cannot contain any negative element." << nl
+                << "Please see the row i = " << i << nl
+                << "xy[i] = " << xy_[i]
+                << exit(FatalError);
+        }
+
+        if (cumulative_ && i > 0 && xy_[i][1] < xy_[i-1][1])
+        {
+            FatalErrorInFunction
+                << type() << "distribution: "
+                << "Elements in the second column for cumulative "
+                << "distribution functions must be non-decreasing." << nl
+                << "Please see the row i = " << i << nl
+                << "xy[i-1] = " << xy_[i-1] << nl
+                << "xy[i] = " << xy_[i]
+                << exit(FatalError);
+        }
+    }
+
     initialise();
 }
 
@@ -104,22 +172,24 @@ Foam::distributionModels::general::general
 :
     distributionModel(typeName, dictionary::null, rndGen),
     xy_(),
-    meanValue_(0.0),
-    integral_()
+    nEntries_(0),
+    meanValue_(0),
+    integral_(),
+    cumulative_(false)
 {
-    scalar minValue = GREAT;
-    scalar maxValue = -GREAT;
+    minValue_ = GREAT;
+    maxValue_ = -GREAT;
     forAll(sampleData, i)
     {
-        minValue = min(minValue, sampleData[i]);
-        maxValue = max(maxValue, sampleData[i]);
+        minValue_ = min(minValue_, sampleData[i]);
+        maxValue_ = max(maxValue_, sampleData[i]);
     }
 
-    label bin0 = floor(minValue/binWidth);
-    label bin1 = ceil(maxValue/binWidth);
-    label nEntries = bin1 - bin0;
+    label bin0 = floor(minValue_/binWidth);
+    label bin1 = ceil(maxValue_/binWidth);
+    nEntries_ = bin1 - bin0;
 
-    if (nEntries == 0)
+    if (nEntries_ == 0)
     {
         WarningInFunction
             << "Data cannot be binned - zero bins generated" << nl
@@ -130,10 +200,10 @@ Foam::distributionModels::general::general
         return;
     }
 
-    xy_.setSize(nEntries);
+    xy_.setSize(nEntries_);
 
     // Populate bin boundaries and initialise occurrences
-    for (label bini = 0; bini < nEntries; ++bini)
+    for (label bini = 0; bini < nEntries_; ++bini)
     {
         xy_[bini][0] = (bin0 + bini)*binWidth;
         xy_[bini][1] = 0;
@@ -154,14 +224,10 @@ Foam::distributionModels::general::general(const general& p)
 :
     distributionModel(p),
     xy_(p.xy_),
+    nEntries_(p.nEntries_),
     meanValue_(p.meanValue_),
-    integral_(p.integral_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::general::~general()
+    integral_(p.integral_),
+    cumulative_(p.cumulative_)
 {}
 
 
@@ -169,57 +235,46 @@ Foam::distributionModels::general::~general()
 
 Foam::scalar Foam::distributionModels::general::sample() const
 {
-    scalar y = rndGen_.sample01<scalar>();
+    const scalar u = rndGen_.sample01<scalar>();
 
-    // Find the interval where y is in the table
+    // Find the interval where u is in the table
     label n = 1;
-    while (integral_[n] <= y)
+    while (integral_[n] <= u)
     {
         n++;
     }
 
-    scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
-    scalar d = xy_[n-1][1] - k*xy_[n-1][0];
+    const scalar k = (xy_[n][1] - xy_[n-1][1])/(xy_[n][0] - xy_[n-1][0]);
+    const scalar d = xy_[n-1][1] - k*xy_[n-1][0];
 
-    scalar alpha = y + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
-    scalar x = 0.0;
+    if (cumulative_)
+    {
+        return (u - d)/k;
+    }
+
+    const scalar alpha =
+        u + xy_[n-1][0]*(0.5*k*xy_[n-1][0] + d) - integral_[n-1];
 
     // If k is small it is a linear equation, otherwise it is of second order
     if (mag(k) > SMALL)
     {
-        scalar p = 2.0*d/k;
-        scalar q = -2.0*alpha/k;
-        scalar sqrtEr = sqrt(0.25*p*p - q);
+        const scalar p = 2.0*d/k;
+        const scalar q = -2.0*alpha/k;
+        const scalar sqrtEr = sqrt(0.25*p*p - q);
 
-        scalar x1 = -0.5*p + sqrtEr;
-        scalar x2 = -0.5*p - sqrtEr;
+        const scalar x1 = -0.5*p + sqrtEr;
+        const scalar x2 = -0.5*p - sqrtEr;
         if ((x1 >= xy_[n-1][0]) && (x1 <= xy_[n][0]))
         {
-            x = x1;
+            return x1;
         }
         else
         {
-            x = x2;
+            return x2;
         }
     }
-    else
-    {
-        x = alpha/d;
-    }
 
-    return x;
-}
-
-
-Foam::scalar Foam::distributionModels::general::minValue() const
-{
-    return xy_.first()[0];
-}
-
-
-Foam::scalar Foam::distributionModels::general::maxValue() const
-{
-    return xy_.last()[0];
+    return alpha/d;
 }
 
 
@@ -278,11 +333,11 @@ void Foam::distributionModels::general::readDict(const dictionary& dict)
 Foam::tmp<Foam::Field<Foam::scalar>>
 Foam::distributionModels::general::x() const
 {
-    tmp<Field<scalar>> tx(new Field<scalar>(xy_.size()));
-    scalarField& xi = tx.ref();
+    auto tx = tmp<scalarField>::New(xy_.size());
+    auto& x = tx.ref();
     forAll(xy_, i)
     {
-        xi[i] = xy_[i][0];
+        x[i] = xy_[i][0];
     }
 
     return tx;
@@ -292,11 +347,11 @@ Foam::distributionModels::general::x() const
 Foam::tmp<Foam::Field<Foam::scalar>>
 Foam::distributionModels::general::y() const
 {
-    tmp<Field<scalar>> ty(new Field<scalar>(xy_.size()));
-    scalarField& yi = ty.ref();
+    auto ty = tmp<scalarField>::New(xy_.size());
+    auto& y = ty.ref();
     forAll(xy_, i)
     {
-        yi[i] = xy_[i][1];
+        y[i] = xy_[i][1];
     }
 
     return ty;
diff --git a/src/lagrangian/distributionModels/general/general.H b/src/lagrangian/distributionModels/general/general.H
index 39d974a7d1..e3f41e8b16 100644
--- a/src/lagrangian/distributionModels/general/general.H
+++ b/src/lagrangian/distributionModels/general/general.H
@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2011-2019 OpenFOAM Foundation
+    Copyright (C) 2016-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +28,66 @@ Class
     Foam::distributionModels::general
 
 Description
-    general distribution model
+    Particle-size distribution model wherein random samples are
+    drawn from a given arbitrary probability density function
+    or cumulative distribution function. Input distributions are
+    specified as pairs of (size, probability) (i.e. (point, value) ).
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        general;
+                    generalDistribution
+                    {
+                        cumulative      false;
+
+                        distribution
+                        (
+                            (<size1>    <probability1>)
+                            (<size2>    <probability2>)
+                            ...
+                            (<sizeN>    <probabilityN>)
+                        );
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: general                | word | yes  | -
+      generalDistribution | Distribution settings      | dict | yes  | -
+      distribution | \<size\>-\<probability\> pairs    | dict | yes  | -
+      \<size\>     | Particle size                     | scalar | yes | -
+      \<probability\> | Volume fraction/probability    | scalar | yes | -
+      cumulative   | Flag to determine if input distribution is specified <!--
+                   --> as cumulative or as density     | bool | no   | false
+    \endtable
+
+Note
+  - An example for a pair within \c distribution subdictionary
+  can be given as follows: "(1e-07 1.3)" means 1.3\% of particles
+  are modelled with a particle diameter of "1e-7" [m], and so on.
+  - Variation between input pairs is assumed to be linear.
+  - Elements in the second column (i.e. probability) are normalised.
+  - Elements in the second column for cumulative distribution
+  functions must start from zero and must be non-decreasing (i.e. monotonic).
+  - Elements in the first column (i.e. size) must be specified
+  in an ascending order.
+  - Input pairs cannot contain any negative element.
 
 SourceFiles
     general.C
@@ -65,26 +124,34 @@ namespace distributionModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class general Declaration
+                          Class general Declaration
 \*---------------------------------------------------------------------------*/
 
 class general
 :
     public distributionModel
 {
-    // Private data
+    typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
 
-        typedef VectorSpace<Vector<scalar>, scalar, 2> pair;
+    // Private Data
 
-        // List of (bin probability)
+        //- List of (x, y=f(x)) pairs
         List<pair> xy_;
 
+        //- Number of entries in the xy_ list
+        label nEntries_;
+
+        //- Mean of the distribution
         scalar meanValue_;
 
+        //- Values of cumulative distribution function
         List<scalar> integral_;
 
+        //- Flag to determine if input is specified as cumulative or as density
+        bool cumulative_;
 
-    // Private member functions
+
+    // Private Member Functions
 
         //- Initialise the distribution parameters
         void initialise();
@@ -102,6 +169,7 @@ public:
         general(const dictionary& dict, Random& rndGen);
 
         //- Construct from components
+        //  Allows negative entries
         general
         (
             const UList<scalar>& sampleData,
@@ -109,7 +177,7 @@ public:
             Random& rndGen
         );
 
-        //- Construct copy
+        //- Copy construct
         general(const general& p);
 
         //- Construct and return a clone
@@ -118,29 +186,26 @@ public:
             return autoPtr<distributionModel>(new general(*this));
         }
 
+        //- No copy assignment
+        void operator=(const general&) = delete;
+
 
     //- Destructor
-    virtual ~general();
+    virtual ~general() = default;
 
 
     // Member Functions
 
         //- Bin boundaries
-        virtual tmp<Field<scalar>> x() const;
+        virtual tmp<scalarField> x() const;
 
         //- Probabilities
-        virtual tmp<Field<scalar>> y() const;
+        virtual tmp<scalarField> y() const;
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the arithmetic mean of the distribution data
         virtual scalar meanValue() const;
 
         //- Write data to stream
diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
index 130e4a9759..0342f1b75d 100644
--- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
+++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.C
@@ -50,11 +50,18 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
 )
 :
     distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue")),
-    d_(distributionModelDict_.get<scalar>("d")),
+    lambda_(distributionModelDict_.getCompat<scalar>("lambda", {{"d", 2112}})),
     n_(distributionModelDict_.get<scalar>("n"))
 {
+    if (lambda_ < VSMALL || n_ < VSMALL)
+    {
+        FatalErrorInFunction
+            << "Scale/Shape parameter cannot be equal to or less than zero:"
+            << "    lambda = " << lambda_
+            << "    n = " << n_
+            << exit(FatalError);
+    }
+
     check();
 }
 
@@ -65,53 +72,37 @@ Foam::distributionModels::massRosinRammler::massRosinRammler
 )
 :
     distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
-    d_(p.d_),
+    lambda_(p.lambda_),
     n_(p.n_)
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::massRosinRammler::~massRosinRammler()
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::scalar Foam::distributionModels::massRosinRammler::sample() const
 {
-    scalar d;
-
-    // Re-sample if the calculated d is out of the physical range
+    scalar d = 0;
     do
     {
-        const scalar a = 3/n_ + 1;
-        const scalar P = rndGen_.sample01<scalar>();
-        const scalar x = Math::invIncGamma(a, P);
-        d = d_*pow(x, 1/n_);
-    } while (d < minValue_ || d > maxValue_);
+        // (YHD:Inverse of Eq. 10)
+        const scalar a = scalar(3)/n_ + scalar(1);
+        const scalar cdfA = Math::incGamma_P(a, pow(minValue_/lambda_, n_) );
+        const scalar cdfB = Math::incGamma_P(a, pow(maxValue_/lambda_, n_) );
+
+        const scalar u = rndGen_.position<scalar>(cdfA, cdfB);
+        const scalar x = Math::invIncGamma(a, u);
+        d = lambda_*pow(x, scalar(1)/n_);
+
+    } while (std::isnan(d));
 
     return d;
 }
 
 
-Foam::scalar Foam::distributionModels::massRosinRammler::minValue() const
-{
-    return minValue_;
-}
-
-
-Foam::scalar Foam::distributionModels::massRosinRammler::maxValue() const
-{
-    return maxValue_;
-}
-
-
 Foam::scalar Foam::distributionModels::massRosinRammler::meanValue() const
 {
-    return d_;
+    // (YHD:Eqs. 11-12)
+    return lambda_*tgamma(scalar(1)/n_ + scalar(1));
 }
 
 
diff --git a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
index 57bb9b2d84..2699203246 100644
--- a/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
+++ b/src/lagrangian/distributionModels/massRosinRammler/massRosinRammler.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,27 +25,141 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::massRosinRammler
+    Foam::distributionModels::massRosinRammler
 
 Description
-    Mass-based Rosin-Rammler distributionModel.
+    Particle-size distribution model wherein random samples are drawn
+    from the two-parameter Rosin-Rammler (Weibull) probability density
+    function corrected to take into account varying number of particles
+    per parcels for fixed-mass parcels.
 
-    Corrected form of the Rosin-Rammler distribution taking into account the
-    varying number of particels per parces for for fixed-mass parcels.  This
-    distribution should be used when
+    "There is a factor of \f$x^3\f$ difference between the size distribution
+    probability density function of individual particles and modeled parcels
+    of particles, \f$ f_{parcels}(D) = x^3 f_{particles}(D) \f$, because the
+    submodel parameter, \f$ PPP \f$ (number of particles per parcel) is based
+    on a fixed mass per parcel which weights the droplet distribution by
+    a factor proportional to \f$ 1/x^3 \f$ (i.e. \f$ PPP = \dot{m}_{total}
+    \Delta_{t_{inj(per-parcel)}} / {m}_{particle}  \f$)." (YHD:p. 1374-1375).
+
+    \c massRosinRammler should be preferred over \c RosinRammler
+    when \c parcelBasisType is based on \c mass:
     \verbatim
         parcelBasisType mass;
     \endverbatim
 
-    See equation 10 in reference:
+    The doubly-truncated mass-corrected
+    Rosin-Rammler probability density function:
+
+    \f[
+        f(x; \lambda, n, A, B) =
+                x^3
+                \frac{n}{\lambda}
+                \left( \frac{x}{\lambda} \right)^{n-1}
+                \frac{
+                    \exp\{ -(\frac{x}{\lambda})^n \}
+                }
+                {
+                    \exp\{ -(\frac{A}{\lambda})^n \}
+                  - \exp\{ -(\frac{B}{\lambda})^n \}
+                }
+    \f]
+    where
+
+    \vartable
+        f(x; \lambda, n, A, B) | Rosin-Rammler probability density function
+        \lambda   | Scale parameter
+        n         | Shape parameter
+        x         | Sample
+        A         | Minimum of the distribution
+        B         | Maximum of the distribution
+    \endvartable
+
+    Constraints:
+    - \f$ \lambda > 0 \f$
+    - \f$ n > 0 \f$
+
+    Random samples are generated by the inverse transform sampling technique
+    by using the quantile function of the doubly-truncated two-parameter
+    mass-corrected Rosin-Rammler (Weibull) probability density function:
+
+    \f[
+        d = \lambda \, Q(a, u)^{1/n}
+    \f]
+    with
+
+    \f[
+        a = \frac{3}{n} + 1
+    \f]
+    where \f$ Q(z_1, z_2) \f$ is the inverse of regularised lower incomplete
+    gamma function, and \f$ u \f$ is a sample drawn from the uniform
+    probability density function on the interval \f$ (x, y) \f$:
+
+    \f[
+        x = \gamma \left( a, \frac{A}{\lambda}^n \right)
+    \f]
+
+    \f[
+        y = \gamma \left( a, \frac{B}{\lambda}^n \right)
+    \f]
+    where \f$ \gamma(z_1, z_2) \f$ is the lower incomplete gamma function.
+
+    Reference:
     \verbatim
-        Yoon, S. S., Hewson, J. C., DesJardin, P. E., Glaze, D. J.,
-        Black, A. R., & Skaggs, R. R. (2004).
-        Numerical modeling and experimental measurements of a high speed
-        solid-cone water spray for use in fire suppression applications.
-        International Journal of Multiphase Flow, 30(11), 1369-1388.
+        Standard model (tag:YHD):
+            Yoon, S. S., Hewson, J. C., DesJardin, P. E.,
+            Glaze, D. J., Black, A. R., & Skaggs, R. R. (2004).
+            Numerical modeling and experimental measurements of a high speed
+            solid-cone water spray for use in fire suppression applications.
+            International Journal of Multiphase Flow, 30(11), 1369-1388.
+            DOI:10.1016/j.ijmultiphaseflow.2004.07.006
     \endverbatim
 
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+                parcelBasisType  mass;
+                ...
+
+                sizeDistribution
+                {
+                    type        massRosinRammler;
+                    massRosinRammlerDistribution
+                    {
+                        lambda      <scaleParameterValue>;
+                        n           <shapeParameterValue>;
+                        minValue    <minValue>;
+                        maxValue    <maxValue>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: massRosinRammler       | word | yes  | -
+      massRosinRammlerDistribution | Distribution settings | dict | yes  | -
+      lambda       | Scale parameter                   | scalar | yes  | -
+      n            | Shape parameter                   | scalar | yes  | -
+      minValue     | Minimum of the distribution       | scalar | yes  | -
+      maxValue     | Maximum of the distribution       | scalar | yes  | -
+    \endtable
+
+Note
+  - In the current context, \c lambda (or \c d) is called "characteristic
+    droplet size", and \c n "empirical dimensionless constant to specify
+    the distribution width, sometimes referred to as the dispersion
+    coefficient." (YHD:p. 1374).
+
 SourceFiles
     massRosinRammler.C
 
@@ -70,19 +185,12 @@ class massRosinRammler
 :
     public distributionModel
 {
-    // Private data
+    // Private Data
 
-        //- Distribution minimum
-        scalar minValue_;
+        //- Scale parameter
+        scalar lambda_;
 
-        //- Distribution maximum
-        scalar maxValue_;
-
-        //- Characteristic droplet size
-        scalar d_;
-
-        //- Empirical dimensionless constant to specify the distribution width,
-        //  sometimes referred to as the dispersion coefficient
+        //- Shape parameter
         scalar n_;
 
 
@@ -97,7 +205,7 @@ public:
         //- Construct from components
         massRosinRammler(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         massRosinRammler(const massRosinRammler& p);
 
         //- Construct and return a clone
@@ -106,23 +214,20 @@ public:
             return autoPtr<distributionModel>(new massRosinRammler(*this));
         }
 
+        //- No copy assignment
+        void operator=(const massRosinRammler&) = delete;
+
 
     //- Destructor
-    virtual ~massRosinRammler();
+    virtual ~massRosinRammler() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/multiNormal/multiNormal.C b/src/lagrangian/distributionModels/multiNormal/multiNormal.C
index d026c6b7dd..43761f9326 100644
--- a/src/lagrangian/distributionModels/multiNormal/multiNormal.C
+++ b/src/lagrangian/distributionModels/multiNormal/multiNormal.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,6 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "multiNormal.H"
+#include "mathematicalConstants.H"
+#include "MathFunctions.H"
+#include "ListOps.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -48,37 +52,67 @@ Foam::distributionModels::multiNormal::multiNormal
 )
 :
     distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue")),
-    range_(maxValue_ - minValue_),
-    expectation_(distributionModelDict_.lookup("expectation")),
-    variance_(distributionModelDict_.lookup("variance")),
-    strength_(distributionModelDict_.lookup("strength"))
+    mu_
+    (
+        distributionModelDict_.lookupCompat
+        (
+            "mu",
+            {{"expectation", 2112}}
+        )
+    ),
+    sigma_
+    (
+        distributionModelDict_.lookupCompat
+        (
+            "sigma",
+            {{"variance", 2112}}
+        )
+    ),
+    weight_
+    (
+        distributionModelDict_.lookupCompat
+        (
+            "weight",
+            {{"strength", 2112}}
+        )
+    )
 {
     check();
 
-    scalar sMax = 0;
-    label n = strength_.size();
-    for (label i=0; i<n; i++)
+    scalar sum = 0;
+    for (label i = 0; i < weight_.size(); ++i)
     {
-        scalar x = expectation_[i];
-        scalar s = strength_[i];
-        for (label j=0; j<n; j++)
+        if (i > 0 && weight_[i] < weight_[i-1])
         {
-            if (i!=j)
-            {
-                scalar x2 = (x - expectation_[j])/variance_[j];
-                scalar y = exp(-0.5*x2*x2);
-                s += strength_[j]*y;
-            }
+            FatalErrorInFunction
+                << type() << "distribution: "
+                << "Weights must be specified in a monotonic order." << nl
+                << "Please see the row i = " << i << nl
+                << "weight[i-1] = " << weight_[i-1] << nl
+                << "weight[i] = " << weight_[i]
+                << exit(FatalError);
         }
 
-        sMax = max(sMax, s);
+        sum += weight_[i];
     }
 
-    for (label i=0; i<n; i++)
+    if (sum < VSMALL)
     {
-        strength_[i] /= sMax;
+        FatalErrorInFunction
+            << type() << "distribution: "
+            << "The sum of weights cannot be zero." << nl
+            << "weight = " << weight_
+            << exit(FatalError);
+    }
+
+    for (label i = 1; i < weight_.size(); ++i)
+    {
+        weight_[i] += weight_[i-1];
+    }
+
+    for (auto& w : weight_)
+    {
+        w /= sum;
     }
 }
 
@@ -86,18 +120,9 @@ Foam::distributionModels::multiNormal::multiNormal
 Foam::distributionModels::multiNormal::multiNormal(const multiNormal& p)
 :
     distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
-    range_(p.range_),
-    expectation_(p.expectation_),
-    variance_(p.variance_),
-    strength_(p.strength_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::multiNormal::~multiNormal()
+    mu_(p.mu_),
+    sigma_(p.sigma_),
+    weight_(p.weight_)
 {}
 
 
@@ -105,54 +130,54 @@ Foam::distributionModels::multiNormal::~multiNormal()
 
 Foam::scalar Foam::distributionModels::multiNormal::sample() const
 {
-    scalar y = 0;
-    scalar x = 0;
-    label n = expectation_.size();
-    bool success = false;
+    const scalar u = rndGen_.sample01<scalar>();
 
-    while (!success)
+    for (label i = 0; i < weight_.size(); ++i)
     {
-        x = minValue_ + range_*rndGen_.sample01<scalar>();
-        y = rndGen_.sample01<scalar>();
-        scalar p = 0.0;
-
-        for (label i=0; i<n; i++)
+        if (weight_[i] > u)
         {
-            scalar nu = expectation_[i];
-            scalar sigma = variance_[i];
-            scalar s = strength_[i];
-            scalar v = (x - nu)/sigma;
-            p += s*exp(-0.5*v*v);
-        }
-
-        if (y<p)
-        {
-            success = true;
+            return sample(mu_[i], sigma_[i]);
         }
     }
 
-    return x;
+    const label last = weight_.size() - 1;
+
+    return sample(mu_[last], sigma_[last]);
 }
 
 
-Foam::scalar Foam::distributionModels::multiNormal::minValue() const
+Foam::scalar Foam::distributionModels::multiNormal::sample
+(
+    const scalar mu,
+    const scalar sigma
+) const
 {
-    return minValue_;
-}
+    const scalar a = (minValue_ - mu)/sigma;
+    const scalar b = (maxValue_ - mu)/sigma;
 
+    const scalar aPhi = 0.5*(1.0 + erf(a/Foam::sqrt(2.0)));
+    const scalar bPhi = 0.5*(1.0 + erf(b/Foam::sqrt(2.0)));
 
-Foam::scalar Foam::distributionModels::multiNormal::maxValue() const
-{
-    return maxValue_;
+    const scalar u = rndGen_.sample01<scalar>();
+    const scalar p = u*(bPhi - aPhi) + aPhi;
+
+    // (B:p. 20-24)
+    const scalar x =
+        mu + sigma*Foam::sqrt(2.0)*Math::erfInv(2.0*p - 1.0);
+
+    // Note: numerical approximation of the inverse function yields slight
+    //       inaccuracies
+
+    return min(max(x, minValue_), maxValue_);
 }
 
 
 Foam::scalar Foam::distributionModels::multiNormal::meanValue() const
 {
-    scalar mean = 0.0;
-    forAll(strength_, i)
+    scalar mean = 0;
+    forAll(weight_, i)
     {
-        mean += strength_[i]*expectation_[i];
+        mean += weight_[i]*mu_[i];
     }
 
     return mean;
diff --git a/src/lagrangian/distributionModels/multiNormal/multiNormal.H b/src/lagrangian/distributionModels/multiNormal/multiNormal.H
index b544a746e1..8639cbf123 100644
--- a/src/lagrangian/distributionModels/multiNormal/multiNormal.H
+++ b/src/lagrangian/distributionModels/multiNormal/multiNormal.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,12 +28,149 @@ Class
     Foam::distributionModels::multiNormal
 
 Description
-    A multiNormal distribution model
+    Particle-size distribution model wherein random samples are drawn
+    from a mixture of a finite set of doubly-truncated univariate normal
+    probability density functions:
 
+    \f[
+        g (\mathbf{x}; \mathbf{\mu}, \mathbf{\sigma}, A, B) =
+            \sum_i w_i f(x_i; \mu_i, \sigma_i, A, B)
+    \f]
+    with for any distribution:
+
+    \f[
+        f(x; \mu, \sigma, A, B) =
+            \frac{1}{\sigma}
+            \frac{
+                \phi \left( \frac{x - \mu}{\sigma} \right)
+            }{
+                \Phi \left( \frac{B - \mu}{\sigma} \right)
+              - \Phi \left( \frac{A - \mu}{\sigma} \right)
+            }
+    \f]
+    where
+
+    \vartable
+      f(x; \mu, \sigma, A, B) | Doubly-truncated univariate normal distribution
+      \mu         | Mean of the parent general normal distribution
+      \sigma      | Standard deviation of the parent general normal distribution
+      \phi(\cdot) | General normal probability density function
+      \Phi(\cdot) | General normal cumulative distribution function
+      x           | Sample
+      A           | Minimum of the distribution (the same for each distribution)
+      B           | Maximum of the distribution (the same for each distribution)
+      w_i         | Weighting factor
+    \endvartable
+
+    Constraints:
+    - \f$ \infty > B > A > 0 \f$
+    - \f$ x \in [B,A] \f$
+    - \f$ \sigma^2 > 0 \f$
+    - \f$ w_i >= 0 \f$
+
+    Random samples are generated by a combination of the inverse transform
+    sampling technique and categorical sampling in three steps:
+    - Draw a sample from the uniform probability density function
+    on the unit interval \f$u = (0, 1)\f$
+    - Find the interval among normalised cumulative weight intervals
+    wherein \f$ u \f$ resides
+    - Draw a sample from the distribution corresponding to the interval by
+    using the quantile function of the doubly-truncated univariate normal
+    probability density function by the following expressions (similar to
+    \c distributionModels::normal):
+
+    \f[
+        x = \mu + \sigma \sqrt{2} \, {erf}^{-1} \left( 2 p - 1 \right)
+    \f]
+    with
+
+    \f[
+        p = u \,
+            \left(
+                \Phi\left(
+                    \frac{B - \mu}{\sigma}
+                    \right)
+              - \Phi\left(
+                    \frac{A - \mu}{\sigma}
+                    \right)
+            \right)
+          + \Phi\left( \frac{A - \mu}{\sigma} \right)
+    \f]
+
+    \f[
+        \Phi(\xi) =
+            \frac{1}{2}
+            \left(
+                1 + {erf}\left(\frac{\xi - \mu}{\sigma \sqrt{2} }\right)
+            \right)
+    \f]
+    where \f$ u \f$ is another sample drawn from the uniform probability
+    density function on the unit interval \f$ (0, 1) \f$.
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
     \verbatim
-        model = sum_i strength_i * exp(-0.5*((x - expectation_i)/variance_i)^2 )
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        multiNormal;
+                    multiNormalDistribution
+                    {
+                        minValue  <min>;
+                        maxValue  <max>;
+                        mu
+                        (
+                            <mean1>
+                            <mean2>
+                            ...
+                        );
+                        sigma
+                        (
+                            <standard deviation1>
+                            <standard deviation2>
+                            ...
+                        );
+                        weight
+                        (
+                            <weight1>
+                            <weight2>
+                            ...
+                        );
+                    }
+                }
+            }
+        }
+    }
     \endverbatim
 
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd   | Deflt
+      type       | Type name: multiNormal            | word | yes    | -
+      multiNormalDistribution | Distribution settings  | dict | yes  | -
+      minValue   | Minimum of the distribution       | scalar | yes  | -
+      maxValue   | Maximum of the distribution       | scalar | yes  | -
+      mu         | List of means of parent general normal <!--
+                 --> distributions                | scalarList | yes | -
+      sigma      | List of standard deviations of parent  <!--
+                 --> general normal distributions | scalarList | yes | -
+      weight   | List of weights of a given distribution in <!--
+                 --> the distribution mixture     | scalarList | yes | -
+    \endtable
+
+Notes
+  - The sum of normal distributions (i.e. a mixture distribution) should not
+  be confused with the sum of normally-distributed random variables.
+  - \c minValue and \c maxValue are the same for all distributions
+  in the distribution mixture.
+  - \c weight should always be input in a non-decreasing (i.e. monotonic) order.
 
 SourceFiles
     multiNormal.C
@@ -59,23 +197,17 @@ class multiNormal
 :
     public distributionModel
 {
-    // Private data
+    // Private Data
 
-        //- Distribution minimum
-        scalar minValue_;
+        //- List of means of the parent general normal distributions
+        List<scalar> mu_;
 
-        //- Distribution maximum
-        scalar maxValue_;
+        //- List of standard deviations of
+        //- the parent general normal distributions
+        List<scalar> sigma_;
 
-        //- Distribution range
-        scalar range_;
-
-
-        // Model coefficients
-
-            List<scalar> expectation_;
-            List<scalar> variance_;
-            List<scalar> strength_;
+        //- List of weights of a given distribution in the mixture
+        List<scalar> weight_;
 
 
 public:
@@ -89,7 +221,7 @@ public:
         //- Construct from components
         multiNormal(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         multiNormal(const multiNormal& p);
 
         //- Construct and return a clone
@@ -98,23 +230,23 @@ public:
             return autoPtr<distributionModel>(new multiNormal(*this));
         }
 
+        //- No copy assignment
+        void operator=(const multiNormal&) = delete;
+
 
     //- Destructor
-    virtual ~multiNormal();
+    virtual ~multiNormal() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
+        //- Sample the normal distribution
+        scalar sample(const scalar mu, const scalar sigma) const;
 
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/normal/normal.C b/src/lagrangian/distributionModels/normal/normal.C
index 803c555b2e..e953ec6b05 100644
--- a/src/lagrangian/distributionModels/normal/normal.C
+++ b/src/lagrangian/distributionModels/normal/normal.C
@@ -27,8 +27,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "normal.H"
-#include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
 #include "MathFunctions.H"
+#include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -50,44 +51,40 @@ Foam::distributionModels::normal::normal
 )
 :
     distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue")),
-    expectation_(distributionModelDict_.get<scalar>("expectation")),
-    variance_(distributionModelDict_.get<scalar>("variance")),
-    a_(0.147)
+    mu_
+    (
+        distributionModelDict_.getCompat<scalar>
+        (
+            "mu",
+            {{"expectation", 2112}}
+        )
+    ),
+    sigma_
+    (
+        distributionModelDict_.getCompat<scalar>
+        (
+            "sigma",
+            {{"variance", 2112}}
+        )
+    )
 {
-    if (minValue_ < 0)
+    if (mag(sigma_) == 0)
     {
         FatalErrorInFunction
-            << "Minimum value must be greater than zero. "
-            << "Supplied minValue = " << minValue_
-            << abort(FatalError);
+            << "Standard deviation cannot be zero." << nl
+            << "    sigma = " << sigma_ << nl
+            << exit(FatalError);
     }
 
-    if (maxValue_ < minValue_)
-    {
-        FatalErrorInFunction
-            << "Maximum value is smaller than the minimum value:"
-            << "    maxValue = " << maxValue_ << ", minValue = " << minValue_
-            << abort(FatalError);
-    }
+    check();
 }
 
 
 Foam::distributionModels::normal::normal(const normal& p)
 :
     distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_),
-    expectation_(p.expectation_),
-    variance_(p.variance_),
-    a_(p.a_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::normal::~normal()
+    mu_(p.mu_),
+    sigma_(p.sigma_)
 {}
 
 
@@ -95,37 +92,45 @@ Foam::distributionModels::normal::~normal()
 
 Foam::scalar Foam::distributionModels::normal::sample() const
 {
+    const scalar a = (minValue_ - mu_)/sigma_;
+    const scalar b = (maxValue_ - mu_)/sigma_;
 
-    scalar a = erf((minValue_ - expectation_)/variance_);
-    scalar b = erf((maxValue_ - expectation_)/variance_);
+    const scalar aPhi = 0.5*(scalar(1) + erf(a/Foam::sqrt(scalar(2))));
+    const scalar bPhi = 0.5*(scalar(1) + erf(b/Foam::sqrt(scalar(2))));
 
-    scalar y = rndGen_.sample01<scalar>();
-    scalar x = Math::erfInv(y*(b - a) + a)*variance_ + expectation_;
+    const scalar u = rndGen_.sample01<scalar>();
+    const scalar p = u*(bPhi - aPhi) + aPhi;
+
+    // (B:p. 20-24)
+    const scalar x =
+        mu_ + sigma_*Foam::sqrt(scalar(2))*Math::erfInv(scalar(2)*p - scalar(1));
 
     // Note: numerical approximation of the inverse function yields slight
     //       inaccuracies
 
-    x = min(max(x, minValue_), maxValue_);
-
-    return x;
-}
-
-
-Foam::scalar Foam::distributionModels::normal::minValue() const
-{
-    return minValue_;
-}
-
-
-Foam::scalar Foam::distributionModels::normal::maxValue() const
-{
-    return maxValue_;
+    return min(max(x, minValue_), maxValue_);
 }
 
 
 Foam::scalar Foam::distributionModels::normal::meanValue() const
 {
-    return expectation_;
+    const scalar a = (minValue_ - mu_)/sigma_;
+    const scalar b = (maxValue_ - mu_)/sigma_;
+
+    // (B:p. 2)
+    const scalar aphi =
+        scalar(1)/Foam::sqrt(scalar(2)*constant::mathematical::pi)
+       *exp(-0.5*sqr(a));
+    const scalar bphi =
+        scalar(1)/Foam::sqrt(scalar(2)*constant::mathematical::pi)
+       *exp(-0.5*sqr(b));
+
+    // (B:p. 4)
+    const scalar aPhi = 0.5*(scalar(1) + erf(a/Foam::sqrt(scalar(2))));
+    const scalar bPhi = 0.5*(scalar(1) + erf(b/Foam::sqrt(scalar(2))));
+
+    // (B:p. 25)
+    return mu_ - sigma_*(bphi - aphi)/(bPhi - aPhi + VSMALL);
 }
 
 
diff --git a/src/lagrangian/distributionModels/normal/normal.H b/src/lagrangian/distributionModels/normal/normal.H
index f98206b25b..dfadd2d513 100644
--- a/src/lagrangian/distributionModels/normal/normal.H
+++ b/src/lagrangian/distributionModels/normal/normal.H
@@ -28,13 +28,124 @@ Class
     Foam::distributionModels::normal
 
 Description
-    A normal distribution model
+    Particle-size distribution model wherein random samples are drawn
+    from the doubly-truncated univariate normal probability density function:
 
+    \f[
+        f(x; \mu, \sigma, A, B) =
+            \frac{1}{\sigma}
+            \frac{
+                \phi \left( \frac{x - \mu}{\sigma} \right)
+            }{
+                \Phi \left( \frac{B - \mu}{\sigma} \right)
+              - \Phi \left( \frac{A - \mu}{\sigma} \right)
+            }
+    \f]
+    where
+
+    \vartable
+      f(x; \mu, \sigma, A, B) | Doubly-truncated univariate normal distribution
+      \mu         | Mean of the parent general normal distribution
+      \sigma      | Standard deviation of the parent general normal distribution
+      \phi(\cdot) | General normal probability density function
+      \Phi(\cdot) | General normal cumulative distribution function
+      x           | Sample
+      A           | Minimum of the distribution
+      B           | Maximum of the distribution
+    \endvartable
+
+    Constraints:
+    - \f$ \infty > B > A > 0 \f$
+    - \f$ x \in [B,A] \f$
+    - \f$ \sigma^2 > 0 \f$
+
+    Random samples are generated by the inverse transform sampling technique
+    by using the quantile function of the doubly-truncated univariate normal
+    probability density function:
+
+    \f[
+        x = \mu + \sigma \sqrt{2} \, {erf}^{-1} \left( 2 p - 1 \right)
+    \f]
+    with
+
+    \f[
+        p = u \,
+            \left(
+                \Phi\left(
+                    \frac{B - \mu}{\sigma}
+                    \right)
+              - \Phi\left(
+                    \frac{A - \mu}{\sigma}
+                    \right)
+            \right)
+          + \Phi\left( \frac{A - \mu}{\sigma} \right)
+    \f]
+
+    \f[
+        \Phi(\xi) =
+            \frac{1}{2}
+            \left(
+                1 + {erf}\left(\frac{\xi - \mu}{\sigma \sqrt{2} }\right)
+            \right)
+    \f]
+    where \f$ u \f$ is sample drawn from the uniform probability
+    density function on the unit interval \f$ (0, 1) \f$.
+
+    Reference:
     \verbatim
-        model = strength * exp(-0.5*((x - expectation)/variance)^2 )
+        Governing expressions (tag:B):
+            Burkardt, J. (2014).
+            The truncated normal distribution.
+            Department of Scientific Computing Website,
+            Florida State University, 1-35.
+            URL:people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
+            (Retrieved on: 19 Feb 2021)
     \endverbatim
 
-    strength only has meaning if there's more than one distribution model
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        normal;
+                    normalDistribution
+                    {
+                        mu        <mean>;
+                        sigma     <stardard deviation>;
+                        minValue  <min>;
+                        maxValue  <max>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd   | Deflt
+      type       | Type name: normal                 | word | yes    | -
+      normalDistribution | Distribution settings     | dict | yes    | -
+      mu         | Mean of the parent general normal distribution <!--
+                 -->                                 | scalar | yes  | -
+      sigma      | Standard deviation of the parent general normal <!--
+                 --> distribution                    | scalar | yes  | -
+      minValue   | Minimum of the distribution       | scalar | yes  | -
+      maxValue   | Maximum of the distribution       | scalar | yes  | -
+    \endtable
+
+Note
+  - The mean and standard deviation of the parent general normal probability
+  distribution function (i.e. input) are not the same with those of the
+  truncated probability distribution function.
 
 SourceFiles
     normal.C
@@ -63,19 +174,11 @@ class normal
 {
     // Private Data
 
-        //- Distribution minimum
-        scalar minValue_;
+        //- Mean of the parent general normal distribution
+        scalar mu_;
 
-        //- Distribution maximum
-        scalar maxValue_;
-
-
-        // Model coefficients
-
-            scalar expectation_;
-            scalar variance_;
-
-            scalar a_;
+        //- Standard deviation of the parent general normal distribution
+        scalar sigma_;
 
 
 public:
@@ -89,7 +192,7 @@ public:
         //- Construct from components
         normal(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         normal(const normal& p);
 
         //- Construct and return a clone
@@ -98,23 +201,20 @@ public:
             return autoPtr<distributionModel>(new normal(*this));
         }
 
+        //- No copy assignment
+        void operator=(const normal&) = delete;
+
 
     //- Destructor
-    virtual ~normal();
+    virtual ~normal() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };
 
diff --git a/src/lagrangian/distributionModels/uniform/uniform.C b/src/lagrangian/distributionModels/uniform/uniform.C
index ef665a2858..690ac67dd4 100644
--- a/src/lagrangian/distributionModels/uniform/uniform.C
+++ b/src/lagrangian/distributionModels/uniform/uniform.C
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,9 +48,7 @@ Foam::distributionModels::uniform::uniform
     Random& rndGen
 )
 :
-    distributionModel(typeName, dict, rndGen),
-    minValue_(distributionModelDict_.get<scalar>("minValue")),
-    maxValue_(distributionModelDict_.get<scalar>("maxValue"))
+    distributionModel(typeName, dict, rndGen)
 {
     check();
 }
@@ -57,15 +56,7 @@ Foam::distributionModels::uniform::uniform
 
 Foam::distributionModels::uniform::uniform(const uniform& p)
 :
-    distributionModel(p),
-    minValue_(p.minValue_),
-    maxValue_(p.maxValue_)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::distributionModels::uniform::~uniform()
+    distributionModel(p)
 {}
 
 
@@ -77,18 +68,6 @@ Foam::scalar Foam::distributionModels::uniform::sample() const
 }
 
 
-Foam::scalar Foam::distributionModels::uniform::minValue() const
-{
-    return minValue_;
-}
-
-
-Foam::scalar Foam::distributionModels::uniform::maxValue() const
-{
-    return maxValue_;
-}
-
-
 Foam::scalar Foam::distributionModels::uniform::meanValue() const
 {
     return 0.5*(minValue_ + maxValue_);
diff --git a/src/lagrangian/distributionModels/uniform/uniform.H b/src/lagrangian/distributionModels/uniform/uniform.H
index 1b99c6c9c1..55b35b595c 100644
--- a/src/lagrangian/distributionModels/uniform/uniform.H
+++ b/src/lagrangian/distributionModels/uniform/uniform.H
@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,7 +28,69 @@ Class
     Foam::distributionModels::uniform
 
 Description
-    Uniform/equally-weighted distribution model
+    Particle-size distribution model wherein random samples are drawn
+    from the doubly-truncated uniform probability density function:
+
+    \f[
+        f(x; A, B) = \frac{1}{B - A}
+    \f]
+    where
+
+    \vartable
+      f(x; A, B)  | Doubly-truncated uniform distribution
+      x           | Sample
+      A           | Minimum of the distribution
+      B           | Maximum of the distribution
+    \endvartable
+
+    Constraints:
+    - \f$ \infty > B > A > 0 \f$
+    - \f$ x \in [B,A] \f$
+    - \f$ \sigma^2 > 0 \f$
+
+    Random samples are generated by the inverse transform sampling technique
+    by using the quantile function of the uniform probability density function:
+
+    \f[
+        x = u \, (B - A) + A
+    \f]
+
+    where \f$ u \f$ is sample drawn from the uniform probability
+    density function on the unit interval \f$ (0, 1) \f$.
+
+Usage
+    Minimal example by using \c constant/\<CloudProperties\>:
+    \verbatim
+    subModels
+    {
+        injectionModels
+        {
+            <name>
+            {
+                ...
+
+                sizeDistribution
+                {
+                    type        uniform;
+                    uniformDistribution
+                    {
+                        minValue  <min>;
+                        maxValue  <max>;
+                    }
+                }
+            }
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                       | Type | Reqd   | Deflt
+      type       | Type name: uniform                | word | yes    | -
+      uniformDistribution | Distribution settings    | dict | yes    | -
+      minValue   | Minimum of the distribution       | scalar | yes  | -
+      maxValue   | Maximum of the distribution       | scalar | yes  | -
+    \endtable
 
 SourceFiles
     uniform.C
@@ -54,15 +117,6 @@ class uniform
 :
     public distributionModel
 {
-    // Private data
-
-        //- Distribution minimum
-        scalar minValue_;
-
-        //- Distribution maximum
-        scalar maxValue_;
-
-
 public:
 
     //- Runtime type information
@@ -74,7 +128,7 @@ public:
         //- Construct from components
         uniform(const dictionary& dict, Random& rndGen);
 
-        //- Construct copy
+        //- Copy construct
         uniform(const uniform& p);
 
         //- Construct and return a clone
@@ -83,23 +137,20 @@ public:
             return autoPtr<distributionModel>(new uniform(*this));
         }
 
+        //- No copy assignment
+        void operator=(const uniform&) = delete;
+
 
     //- Destructor
-    virtual ~uniform();
+    virtual ~uniform() = default;
 
 
     // Member Functions
 
-        //- Sample the distributionModel
+        //- Sample the distribution
         virtual scalar sample() const;
 
-        //- Return the minimum value
-        virtual scalar minValue() const;
-
-        //- Return the maximum value
-        virtual scalar maxValue() const;
-
-        //- Return the mean value
+        //- Return the theoretical mean of the distribution
         virtual scalar meanValue() const;
 };