diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.C b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.C
index 53fba052e5..8085b1c43b 100644
--- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.C
+++ b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,14 +57,14 @@ dripping::dripping
         dict.optionalSubDict(typeName + "Coeffs")
        .lookupOrDefault("minParticlesPerParcel", 1)
     ),
-    rndGen_(label(0)),
     parcelDistribution_
     (
         distribution::New
         (
+            dimLength,
             dict.optionalSubDict(typeName + "Coeffs")
            .subDict("parcelDistribution"),
-            rndGen_,
+            0,
             0
         )
     )
diff --git a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.H b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.H
index b5255364e2..24224bb0a0 100644
--- a/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.H
+++ b/applications/modules/isothermalFilm/fvModels/filmCloudTransfer/ejectionModels/dripping/dripping.H
@@ -96,9 +96,6 @@ class dripping
         //- Minimum number of droplets per parcel
         scalar minParticlesPerParcel_;
 
-        //- Random number generator
-        randomGenerator rndGen_;
-
         //- Parcel size PDF model
         const autoPtr<distribution> parcelDistribution_;
 
diff --git a/applications/test/distribution/RosinRammler b/applications/test/distribution/RosinRammler
index 3bba84a91f..f4d76f77db 100644
--- a/applications/test/distribution/RosinRammler
+++ b/applications/test/distribution/RosinRammler
@@ -1,5 +1,5 @@
 type RosinRammler;
-min 0.001;
-max 0.015;
-d 0.014;
-n 2;
+min 1 [mm];
+max 15 [mm];
+d 14 [mm];
+n 2 [];
diff --git a/applications/test/distribution/Test-distribution.C b/applications/test/distribution/Test-distribution.C
index f4d8a30c96..641cb17dc9 100644
--- a/applications/test/distribution/Test-distribution.C
+++ b/applications/test/distribution/Test-distribution.C
@@ -59,15 +59,12 @@ int main(int argc, char *argv[])
     static const label nSampled = 100;
     static const label nSamples = 10000000;
 
-    // Initialise the random number generator
-    randomGenerator rndGen(clock::getTime());
-
     // Create a zero-size-exponent, zero-sampling-size-exponent distribution
     // from which to get X-coordinates
     dict.set("Q", 0);
     autoPtr<distribution> distribution00Ptr
     (
-        distribution::New(dict, rndGen, 0)
+        distribution::New(unitAny, dict, 0, clock::getTime())
     );
 
     // Get the X-coordinates for the plots
@@ -105,7 +102,7 @@ int main(int argc, char *argv[])
         (
             Q == 0
           ? distribution00Ptr->clone(0)
-          : distribution::New(dict, rndGen, 0)
+          : distribution::New(unitAny, dict, 0, clock::getTime(), false, false)
         );
 
         // Resize
@@ -152,14 +149,19 @@ int main(int argc, char *argv[])
         for (label sampleQ = 0; sampleQ <= nQ; ++ sampleQ)
         {
             // Create a Q-size-exponent, Q-sampling-size-exponent distribution
-            Info<< incrIndent;
             autoPtr<distribution> distributionQSampleQPtr
             (
                 distributionQ0Ptr->clone(sampleQ)
             );
-            distributionQSampleQPtr->mean();
-            distributionQSampleQPtr->report();
-            Info<< decrIndent;
+            {
+                OStringStream oss;
+                distributionQSampleQPtr->write(oss, unitAny);
+                writeEntry(oss, "sampleQ", sampleQ);
+                writeEntry(oss, "min", distributionQSampleQPtr->min());
+                writeEntry(oss, "mean", distributionQSampleQPtr->mean());
+                writeEntry(oss, "max", distributionQSampleQPtr->max());
+                Info<< dictionary(IStringStream(oss.str())());
+            }
 
             // Resize
             const label QSampleQi = yNames.size();
diff --git a/applications/test/distribution/tabulatedCumulative b/applications/test/distribution/tabulatedCumulative
index dd490a8c6c..0b2a610021 100644
--- a/applications/test/distribution/tabulatedCumulative
+++ b/applications/test/distribution/tabulatedCumulative
@@ -1,5 +1,6 @@
 type tabulatedCumulative;
-values
+units [mm];
+distribution
 (
     (1    0)
     (1.73 0.0281384)
diff --git a/applications/test/distribution/tabulatedDensity b/applications/test/distribution/tabulatedDensity
index 7e2795e99b..6130a82ae2 100644
--- a/applications/test/distribution/tabulatedDensity
+++ b/applications/test/distribution/tabulatedDensity
@@ -1,5 +1,6 @@
 type tabulatedDensity;
-values
+units [mm];
+distribution
 (
     (1    0.0025)
     (1.73 0.0528)
diff --git a/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H b/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
index 7126d221be..a8089e7b44 100644
--- a/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
+++ b/applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H
@@ -24,17 +24,7 @@
 
     mkDir(pdfPath);
 
-    randomGenerator rndGen(label(0));
-
-    autoPtr<distribution> p
-    (
-        distribution::New
-        (
-            pdfDictionary,
-            rndGen,
-            0
-        )
-    );
+    autoPtr<distribution> p(distribution::New(unitAny, pdfDictionary, 0, 0));
 
     const scalar xMin = p->min();
     const scalar xMax = p->max();
diff --git a/src/OpenFOAM/distributions/RosinRammler/RosinRammler.C b/src/OpenFOAM/distributions/RosinRammler/RosinRammler.C
index fdd6417a42..8e211469d8 100644
--- a/src/OpenFOAM/distributions/RosinRammler/RosinRammler.C
+++ b/src/OpenFOAM/distributions/RosinRammler/RosinRammler.C
@@ -73,27 +73,28 @@ Foam::tmp<Foam::scalarField> Foam::distributions::RosinRammler::Phi
 
 Foam::distributions::RosinRammler::RosinRammler
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<unintegrableForNonZeroQ, RosinRammler>
     (
         typeName,
+        units,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     ),
-    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"})),
-    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"})),
-    d_(dict.lookup<scalar>("d")),
-    n_(dict.lookup<scalar>("n"))
+    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"}, units)),
+    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"}, units)),
+    d_(dict.lookup<scalar>("d", units)),
+    n_(dict.lookup<scalar>("n", unitless))
 {
     validateBounds(dict);
     validatePositive(dict);
     mean();
-    report();
 }
 
 
@@ -147,6 +148,21 @@ Foam::scalar Foam::distributions::RosinRammler::max() const
 }
 
 
+void Foam::distributions::RosinRammler::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<unintegrableForNonZeroQ, RosinRammler>::write(os, units);
+
+    writeEntry(os, "min", units, min_);
+    writeEntry(os, "max", units, max_);
+    writeEntry(os, "d", units, d_);
+    writeEntry(os, "n", unitless, n_);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::RosinRammler::x(const label n) const
 {
diff --git a/src/OpenFOAM/distributions/RosinRammler/RosinRammler.H b/src/OpenFOAM/distributions/RosinRammler/RosinRammler.H
index 8d30e6838d..9597628e42 100644
--- a/src/OpenFOAM/distributions/RosinRammler/RosinRammler.H
+++ b/src/OpenFOAM/distributions/RosinRammler/RosinRammler.H
@@ -124,9 +124,10 @@ public:
         //- Construct from a dictionary
         RosinRammler
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -157,6 +158,9 @@ public:
         //- Return the maximum value
         virtual scalar max() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 };
diff --git a/src/OpenFOAM/distributions/distribution/distribution.C b/src/OpenFOAM/distributions/distribution/distribution.C
index 9f9adfbf25..b85f888355 100644
--- a/src/OpenFOAM/distributions/distribution/distribution.C
+++ b/src/OpenFOAM/distributions/distribution/distribution.C
@@ -76,14 +76,15 @@ Foam::distribution::clipPDF
 Foam::distribution::distribution
 (
     const word& name,
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
-    rndGen_(rndGen),
     Q_(dict.lookup<scalar>("Q")),
-    sampleQ_(sampleQ)
+    sampleQ_(sampleQ),
+    rndGen_("rndGen", dict, std::move(rndGen))
 {
     if (Q_ < 0)
     {
@@ -103,22 +104,22 @@ Foam::distribution::distribution
 
 Foam::distribution::distribution
 (
-    randomGenerator& rndGen,
     const label Q,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
-    rndGen_(rndGen),
     Q_(Q),
-    sampleQ_(sampleQ)
+    sampleQ_(sampleQ),
+    rndGen_(rndGen)
 {}
 
 
 Foam::distribution::distribution(const distribution& d, const label sampleQ)
 :
-    rndGen_(d.rndGen_),
     Q_(d.Q_),
-    sampleQ_(sampleQ)
+    sampleQ_(sampleQ),
+    rndGen_(d.rndGen_)
 {}
 
 
@@ -130,10 +131,16 @@ Foam::distribution::~distribution()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::distribution::report() const
+void Foam::distribution::write(Ostream& os, const unitConversion& units) const
 {
-    Info<< indent << "min/average/max value = " << min() << '/' << mean()
-        << '/' << max() << endl;
+    writeEntry(os, "type", type());
+    writeEntry(os, "Q", Q_);
+}
+
+
+void Foam::distribution::writeState(Ostream& os) const
+{
+    writeEntry(os, "rndGen", rndGen_);
 }
 
 
@@ -163,4 +170,24 @@ Foam::tmp<Foam::scalarField> Foam::distribution::x(const label n) const
 }
 
 
+// * * * * * * * * * * * * * * * Global Functions  * * * * * * * * * * * * * //
+
+void Foam::writeEntry
+(
+    Ostream& os,
+    const unitConversion& units,
+    const distribution& d,
+    const bool write,
+    const bool writeState
+)
+{
+    os  << nl << indent << token::BEGIN_BLOCK << nl << incrIndent;
+
+    if (write) d.write(os, units);
+    if (writeState) d.writeState(os);
+
+    os  << decrIndent << indent << token::END_BLOCK << endl;
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/distributions/distribution/distribution.H b/src/OpenFOAM/distributions/distribution/distribution.H
index 99a580256c..6fb946719c 100644
--- a/src/OpenFOAM/distributions/distribution/distribution.H
+++ b/src/OpenFOAM/distributions/distribution/distribution.H
@@ -79,15 +79,15 @@ protected:
 
     // Protected data
 
-        //- Reference to a random number generator
-        randomGenerator& rndGen_;
-
         //- Distribution size exponent
         const label Q_;
 
         //- Sample size exponent
         const label sampleQ_;
 
+        //- Random number generator
+        mutable randomGenerator rndGen_;
+
 
     // Protected Member Functions
 
@@ -130,11 +130,12 @@ public:
         distribution,
         dictionary,
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         ),
-        (dict, rndGen, sampleQ)
+        (units, dict, sampleQ, std::move(rndGen))
     );
 
 
@@ -144,17 +145,18 @@ public:
         distribution
         (
             const word& name,
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct from components
         distribution
         (
-            randomGenerator& rndGen,
             const label Q,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -170,13 +172,35 @@ public:
         }
 
 
-    //- Selector
-    static autoPtr<distribution> New
-    (
-        const dictionary& dict,
-        randomGenerator& rndGen,
-        const label sampleQ
-    );
+    // Selectors
+
+        //- Select from dictionary and a random generator
+        static autoPtr<distribution> New
+        (
+            const unitConversion& units,
+            const dictionary& dict,
+            const label sampleQ,
+            randomGenerator&& rndGen,
+            const bool report = true
+        );
+
+        //- Select from a dictionary and a random generator seed and global flag
+        static autoPtr<distribution> New
+        (
+            const unitConversion& units,
+            const dictionary& dict,
+            const label sampleQ,
+            const randomGenerator::seed& s,
+            const bool global = false,
+            const bool report = true
+        );
+
+        //- Re-select with a different sample size exponent
+        static autoPtr<distribution> New
+        (
+            autoPtr<distribution>& dPtr,
+            const label sampleQ
+        );
 
 
     //- Destructor
@@ -204,8 +228,11 @@ public:
         //- Return the mean value
         virtual scalar mean() const = 0;
 
-        //- Report
-        void report() const;
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
+        //- Write the state to a stream
+        virtual void writeState(Ostream& os) const;
 
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
@@ -225,6 +252,16 @@ FOR_ALL_FIELD_TYPES(DISTRIBUTION_TEMPLATED_SAMPLE_TYPE);
 #undef DISTRIBUTION_TEMPLATED_SAMPLE_TYPE
 
 
+void writeEntry
+(
+    Ostream& os,
+    const unitConversion&,
+    const distribution& d,
+    const bool write = true,
+    const bool writeState = true
+);
+
+
 /*---------------------------------------------------------------------------*\
                       Class FieldDistribution Declaration
 \*---------------------------------------------------------------------------*/
@@ -263,31 +300,10 @@ public:
 
         //- Sample the distribution into components of a primitive type
         template<class Type>
-        Type sample() const
-        {
-            Type value;
-            for (direction i = 0; i < pTraits<Type>::nComponents; ++ i)
-            {
-                setComponent(value, i) =
-                    static_cast<const Derived&>(*this).sample();
-            }
-            return value;
-        }
+        Type sample() const;
 
         //- Sample the distribution into a field
-        virtual tmp<scalarField> sample(const label n) const
-        {
-            tmp<scalarField> tResult(new scalarField(n));
-            scalarField& result = tResult.ref();
-
-            forAll(result, i)
-            {
-                result[i] =
-                    static_cast<const Derived&>(*this).sample();
-            }
-
-            return tResult;
-        }
+        virtual tmp<scalarField> sample(const label n) const;
 };
 
 
@@ -297,6 +313,12 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#ifdef NoRepository
+    #include "distributionTemplates.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/OpenFOAM/distributions/distribution/distributionNew.C b/src/OpenFOAM/distributions/distribution/distributionNew.C
index c0d3fa81a4..596a76bd10 100644
--- a/src/OpenFOAM/distributions/distribution/distributionNew.C
+++ b/src/OpenFOAM/distributions/distribution/distributionNew.C
@@ -29,14 +29,19 @@ License
 
 Foam::autoPtr<Foam::distribution> Foam::distribution::New
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen,
+    const bool report
 )
 {
     const word distributionType(dict.lookup("type"));
 
-    Info<< "Selecting " << typeName << " type " << distributionType << endl;
+    if (report)
+    {
+        Info<< "Selecting " << typeName << " type " << distributionType << endl;
+    }
 
     dictionaryConstructorTable::iterator cstrIter =
         dictionaryConstructorTablePtr_->find(distributionType);
@@ -50,22 +55,62 @@ Foam::autoPtr<Foam::distribution> Foam::distribution::New
             << exit(FatalError);
     }
 
-    Info<< incrIndent;
-
     autoPtr<distribution> distributionPtr
     (
         cstrIter()
         (
+            units,
             dict.optionalSubDict(distributionType + "Distribution"),
-            rndGen,
-            sampleQ
+            sampleQ,
+            std::move(rndGen)
         )
     );
 
-    Info<< decrIndent;
+    if (report)
+    {
+        Info<< incrIndent << indent
+            << "min/average/max value = "
+            << distributionPtr->min() << '/'
+            << distributionPtr->mean() << '/'
+            << distributionPtr->max()
+            << decrIndent << endl;
+    }
 
     return distributionPtr;
 }
 
 
+Foam::autoPtr<Foam::distribution> Foam::distribution::New
+(
+    const unitConversion& units,
+    const dictionary& dict,
+    const label sampleQ,
+    const randomGenerator::seed& s,
+    const bool global,
+    const bool report
+)
+{
+    return New(units, dict, sampleQ, randomGenerator(s, global), report);
+}
+
+
+Foam::autoPtr<Foam::distribution> Foam::distribution::New
+(
+    autoPtr<distribution>& dPtr,
+    const label sampleQ
+)
+{
+    if (dPtr->sampleQ_ == sampleQ)
+    {
+        return autoPtr<distribution>(dPtr.ptr());
+    }
+    else
+    {
+        autoPtr<distribution> result(dPtr->clone(sampleQ));
+        dPtr.clear();
+        return result;
+    }
+}
+
+
 // ************************************************************************* //
diff --git a/src/OpenFOAM/distributions/distribution/distributionTemplates.C b/src/OpenFOAM/distributions/distribution/distributionTemplates.C
new file mode 100644
index 0000000000..14d00ba773
--- /dev/null
+++ b/src/OpenFOAM/distributions/distribution/distributionTemplates.C
@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2024 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "distribution.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Base, class Derived>
+template<class Type>
+Type Foam::FieldDistribution<Base, Derived>::sample() const
+{
+    Type value;
+    for (direction i = 0; i < pTraits<Type>::nComponents; ++ i)
+    {
+        setComponent(value, i) =
+            static_cast<const Derived&>(*this).sample();
+    }
+    return value;
+}
+
+
+template<class Base, class Derived>
+Foam::tmp<Foam::scalarField> Foam::FieldDistribution<Base, Derived>::sample
+(
+    const label n
+) const
+{
+    tmp<scalarField> tResult(new scalarField(n));
+    scalarField& result = tResult.ref();
+
+    forAll(result, i)
+    {
+        result[i] =
+            static_cast<const Derived&>(*this).sample();
+    }
+
+    return tResult;
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/distributions/exponential/exponential.C b/src/OpenFOAM/distributions/exponential/exponential.C
index 9cb6ca51cb..cc2ab20d70 100644
--- a/src/OpenFOAM/distributions/exponential/exponential.C
+++ b/src/OpenFOAM/distributions/exponential/exponential.C
@@ -71,26 +71,27 @@ Foam::tmp<Foam::scalarField> Foam::distributions::exponential::Phi
 
 Foam::distributions::exponential::exponential
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<unintegrableForNonZeroQ, exponential>
     (
         typeName,
+        units,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     ),
-    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"})),
-    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"})),
-    lambda_(dict.lookup<scalar>("lambda"))
+    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"}, units)),
+    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"}, units)),
+    lambda_(dict.lookup<scalar>("lambda", unitless))
 {
     validateBounds(dict);
     validatePositive(dict);
     mean();
-    report();
 }
 
 
@@ -143,6 +144,20 @@ Foam::scalar Foam::distributions::exponential::max() const
 }
 
 
+void Foam::distributions::exponential::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<unintegrableForNonZeroQ, exponential>::write(os, units);
+
+    writeEntry(os, "min", units, min_);
+    writeEntry(os, "max", units, max_);
+    writeEntry(os, "lambda", unitless, lambda_);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::exponential::x(const label n) const
 {
diff --git a/src/OpenFOAM/distributions/exponential/exponential.H b/src/OpenFOAM/distributions/exponential/exponential.H
index 569063a3fe..c7305c7b19 100644
--- a/src/OpenFOAM/distributions/exponential/exponential.H
+++ b/src/OpenFOAM/distributions/exponential/exponential.H
@@ -114,9 +114,10 @@ public:
         //- Construct from a dictionary
         exponential
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -147,6 +148,9 @@ public:
         //- Return the maximum value
         virtual scalar max() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 };
diff --git a/src/OpenFOAM/distributions/fixedValue/fixedValue.C b/src/OpenFOAM/distributions/fixedValue/fixedValue.C
index f500427414..4b9543b7ba 100644
--- a/src/OpenFOAM/distributions/fixedValue/fixedValue.C
+++ b/src/OpenFOAM/distributions/fixedValue/fixedValue.C
@@ -42,16 +42,20 @@ namespace distributions
 
 Foam::distributions::fixedValue::fixedValue
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label
+    const label,
+    randomGenerator&& rndGen
 )
 :
-    FieldDistribution<distribution, fixedValue>(rndGen, -labelMax, -labelMax),
-    value_(dict.lookup<scalar>("value"))
-{
-    report();
-}
+    FieldDistribution<distribution, fixedValue>
+    (
+        -labelMax,
+        -labelMax,
+        std::move(rndGen)
+    ),
+    value_(dict.lookup<scalar>("value", units))
+{}
 
 
 Foam::distributions::fixedValue::fixedValue
@@ -97,6 +101,18 @@ Foam::scalar Foam::distributions::fixedValue::mean() const
 }
 
 
+void Foam::distributions::fixedValue::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<distribution, fixedValue>::write(os, units);
+
+    writeEntry(os, "value", units, value_);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::fixedValue::x(const label n) const
 {
diff --git a/src/OpenFOAM/distributions/fixedValue/fixedValue.H b/src/OpenFOAM/distributions/fixedValue/fixedValue.H
index 2deb70425f..f0f73dd76c 100644
--- a/src/OpenFOAM/distributions/fixedValue/fixedValue.H
+++ b/src/OpenFOAM/distributions/fixedValue/fixedValue.H
@@ -81,9 +81,10 @@ public:
         //- Construct from a dictionary
         fixedValue
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -117,6 +118,9 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 
diff --git a/src/OpenFOAM/distributions/multiNormal/multiNormal.C b/src/OpenFOAM/distributions/multiNormal/multiNormal.C
index 971bf96365..d146bb297e 100644
--- a/src/OpenFOAM/distributions/multiNormal/multiNormal.C
+++ b/src/OpenFOAM/distributions/multiNormal/multiNormal.C
@@ -112,35 +112,37 @@ Foam::tmp<Foam::scalarField> Foam::distributions::multiNormal::Phi
 
 Foam::distributions::multiNormal::multiNormal
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<unintegrableForNonZeroQ, multiNormal>
     (
         typeName,
+        units,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     ),
     cumulativeStrengths_(readCumulativeStrengths(dict))
 {
     const scalar min
     (
-        dict.lookupBackwardsCompatible<scalar>({"min", "minValue"})
+        dict.lookupBackwardsCompatible<scalar>({"min", "minValue"}, units)
     );
     const scalar max
     (
-        dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"})
+        dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"}, units)
     );
     const scalarList mu
     (
-        dict.lookupBackwardsCompatible<scalarList>({"mu", "expectation"})
+        dict.lookupBackwardsCompatible<scalarList>({"mu", "expectation"}, units)
     );
     const scalarList sigma
     (
-        dict.lookup<scalarList>("sigma")
+        dict.lookup<scalarList>("sigma", units)
     );
 
     if
@@ -164,9 +166,9 @@ Foam::distributions::multiNormal::multiNormal
             i,
             new normal
             (
-                rndGen_,
                 0,
                 0,
+                rndGen_.generator(),
                 -1,
                 min,
                 max,
@@ -179,7 +181,6 @@ Foam::distributions::multiNormal::multiNormal
     validateBounds(dict);
     if (q() != 0) validatePositive(dict);
     mean();
-    report();
 }
 
 
@@ -200,9 +201,9 @@ Foam::distributions::multiNormal::multiNormal
             i,
             new normal
             (
-                rndGen_,
                 0,
                 0,
+                randomGenerator(d.distributions_[i].rndGen_),
                 -1,
                 d.distributions_[i].min_,
                 d.distributions_[i].max_,
@@ -257,6 +258,29 @@ Foam::scalar Foam::distributions::multiNormal::max() const
 }
 
 
+void Foam::distributions::multiNormal::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<unintegrableForNonZeroQ, multiNormal>::write(os, units);
+
+    writeEntry(os, "min", units, distributions_[0].min());
+    writeEntry(os, "max", units, distributions_[0].max());
+
+    scalarList mu(distributions_.size()), sigma(distributions_.size());
+    forAll(distributions_, i)
+    {
+        mu[i] = distributions_[i].mu();
+        sigma[i] = distributions_[i].sigma();
+    }
+
+    writeEntry(os, "mu", units, mu);
+    writeEntry(os, "sigma", units, sigma);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::multiNormal::x(const label n) const
 {
diff --git a/src/OpenFOAM/distributions/multiNormal/multiNormal.H b/src/OpenFOAM/distributions/multiNormal/multiNormal.H
index 8536992d3f..814dfe9351 100644
--- a/src/OpenFOAM/distributions/multiNormal/multiNormal.H
+++ b/src/OpenFOAM/distributions/multiNormal/multiNormal.H
@@ -118,9 +118,10 @@ public:
         //- Construct from a dictionary
         multiNormal
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -151,6 +152,9 @@ public:
         //- Return the maximum value
         virtual scalar max() const;
 
+        //- Write stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 };
diff --git a/src/OpenFOAM/distributions/normal/normal.C b/src/OpenFOAM/distributions/normal/normal.C
index bbadcc29f4..852b880194 100644
--- a/src/OpenFOAM/distributions/normal/normal.C
+++ b/src/OpenFOAM/distributions/normal/normal.C
@@ -87,35 +87,36 @@ Foam::scalar Foam::distributions::normal::sampleForZeroQ(const scalar s) const
 
 Foam::distributions::normal::normal
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<unintegrableForNonZeroQ, normal>
     (
         typeName,
+        units,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     ),
-    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"})),
-    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"})),
-    mu_(dict.lookupBackwardsCompatible<scalar>({"mu", "expectation"})),
-    sigma_(dict.lookup<scalar>("sigma"))
+    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"}, units)),
+    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"}, units)),
+    mu_(dict.lookupBackwardsCompatible<scalar>({"mu", "expectation"}, units)),
+    sigma_(dict.lookup<scalar>("sigma", units))
 {
     validateBounds(dict);
     if (q() != 0) validatePositive(dict);
     mean();
-    report();
 }
 
 
 Foam::distributions::normal::normal
 (
-    randomGenerator& rndGen,
     const label Q,
     const label sampleQ,
+    randomGenerator&& rndGen,
     const label n,
     const scalar min,
     const scalar max,
@@ -123,7 +124,13 @@ Foam::distributions::normal::normal
     const scalar sigma
 )
 :
-    FieldDistribution<unintegrableForNonZeroQ, normal>(rndGen, Q, sampleQ, n),
+    FieldDistribution<unintegrableForNonZeroQ, normal>
+    (
+        Q,
+        sampleQ,
+        std::move(rndGen),
+        n
+    ),
     min_(min),
     max_(max),
     mu_(mu),
@@ -195,6 +202,21 @@ Foam::scalar Foam::distributions::normal::mean() const
 }
 
 
+void Foam::distributions::normal::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<unintegrableForNonZeroQ, normal>::write(os, units);
+
+    writeEntry(os, "min", units, min_);
+    writeEntry(os, "max", units, max_);
+    writeEntry(os, "mu", units, mu_);
+    writeEntry(os, "sigma", units, sigma_);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::normal::x(const label n) const
 {
diff --git a/src/OpenFOAM/distributions/normal/normal.H b/src/OpenFOAM/distributions/normal/normal.H
index 7f9c9b5aee..380b94b15f 100644
--- a/src/OpenFOAM/distributions/normal/normal.H
+++ b/src/OpenFOAM/distributions/normal/normal.H
@@ -126,17 +126,18 @@ public:
         //- Construct from a dictionary
         normal
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct from components
         normal
         (
-            randomGenerator& rndGen,
             const label Q,
             const label sampleQ,
+            randomGenerator&& rndGen,
             const label n,
             const scalar min,
             const scalar max,
@@ -175,6 +176,21 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Return the mean value
+        inline scalar mu() const
+        {
+            return mu_;
+        }
+
+        //- Return the standard deviation
+        inline scalar sigma() const
+        {
+            return sigma_;
+        }
+
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 };
diff --git a/src/OpenFOAM/distributions/standardNormal/standardNormal.C b/src/OpenFOAM/distributions/standardNormal/standardNormal.C
index e0e0bd8734..d2f2e611ea 100644
--- a/src/OpenFOAM/distributions/standardNormal/standardNormal.C
+++ b/src/OpenFOAM/distributions/standardNormal/standardNormal.C
@@ -64,15 +64,25 @@ Foam::scalar Foam::distributions::standardNormal::approxErfInv(const scalar y)
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::distributions::standardNormal::standardNormal(randomGenerator& rndGen)
+Foam::distributions::standardNormal::standardNormal(randomGenerator&& rndGen)
 :
-    FieldDistribution<distribution, standardNormal>(rndGen, 0, 0)
+    FieldDistribution<distribution, standardNormal>(0, 0, std::move(rndGen))
+{}
+
+
+Foam::distributions::standardNormal::standardNormal
+(
+    const randomGenerator::seed& s,
+    const bool global
+)
+:
+    standardNormal(randomGenerator(s, global))
 {}
 
 
 Foam::distributions::standardNormal::standardNormal(const standardNormal& d)
 :
-    FieldDistribution<distribution, standardNormal>(d.rndGen_, 0, 0)
+    FieldDistribution<distribution, standardNormal>(d, 0)
 {}
 
 
diff --git a/src/OpenFOAM/distributions/standardNormal/standardNormal.H b/src/OpenFOAM/distributions/standardNormal/standardNormal.H
index 822c6d2411..553e321889 100644
--- a/src/OpenFOAM/distributions/standardNormal/standardNormal.H
+++ b/src/OpenFOAM/distributions/standardNormal/standardNormal.H
@@ -87,7 +87,14 @@ public:
     // Constructors
 
         //- Construct from a random generator
-        standardNormal(randomGenerator& rndGen);
+        standardNormal(randomGenerator&& rndGen);
+
+        //- Construct from a seed
+        standardNormal
+        (
+            const randomGenerator::seed& s,
+            const bool global = false
+        );
 
         //- Construct copy
         standardNormal(const standardNormal& d);
diff --git a/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.C b/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.C
index e578088ad3..bc68de0809 100644
--- a/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.C
+++ b/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.C
@@ -44,20 +44,22 @@ namespace distributions
 
 Foam::distributions::tabulatedCumulative::tabulatedCumulative
 (
+    const unitConversion& defaultUnits,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<distribution, tabulatedCumulative>
     (
         typeName,
+        defaultUnits,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     )
 {
-    List<Pair<scalar>> values(dict.lookup("distribution"));
+    List<Tuple2<scalar, scalar>> values(dict.lookup("distribution"));
 
     // Checks
     if (values.first().second() != 0)
@@ -84,11 +86,15 @@ Foam::distributions::tabulatedCumulative::tabulatedCumulative
         }
     }
 
+    // Optionally read units
+    unitConversion units(defaultUnits);
+    units.readIfPresent("units", dict);
+
     // Copy the coordinates
     x_.resize(values.size());
     forAll(values, i)
     {
-        x_[i] = values[i].first();
+        x_[i] = units.toStandard(values[i].first());
     }
 
     // Set the CDF. Copy if q == 0. Re-integrated if q != 0.
@@ -115,7 +121,6 @@ Foam::distributions::tabulatedCumulative::tabulatedCumulative
     // More checks
     validateBounds(dict);
     if (q() != 0) validatePositive(dict);
-    report();
 }
 
 
@@ -185,6 +190,39 @@ Foam::scalar Foam::distributions::tabulatedCumulative::mean() const
 }
 
 
+void Foam::distributions::tabulatedCumulative::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<distribution, tabulatedCumulative>::write(os, units);
+
+    // Recover the CDF
+    scalarField CDF(CDF_.size());
+    CDF[0] = 0;
+    for (label i = 1; i < CDF_.size(); ++ i)
+    {
+        CDF[i] =
+            CDF[i - 1]
+          + integerPow((x_[i] + x_[i - 1])/2, -q())
+           *(CDF_[i] - CDF_[i - 1]);
+    }
+
+    // Normalise the CDF
+    CDF /= CDF.last();
+
+    // Construct and write the values
+    List<Tuple2<scalar, scalar>> values(CDF_.size());
+    forAll(values, i)
+    {
+        values[i].first() = units.toUser(x_[i]);
+        values[i].second() = CDF[i];
+    }
+    writeEntry(os, "distribution", values);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::distributions::tabulatedCumulative::x
 (
     const label
diff --git a/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.H b/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.H
index 2074a40954..295b8426e0 100644
--- a/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.H
+++ b/src/OpenFOAM/distributions/tabulatedCumulative/tabulatedCumulative.H
@@ -102,9 +102,10 @@ public:
         //- Construct from a dictionary
         tabulatedCumulative
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -141,6 +142,9 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 
diff --git a/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.C b/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.C
index 688573cce9..cc7aad08a5 100644
--- a/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.C
+++ b/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.C
@@ -44,20 +44,22 @@ namespace distributions
 
 Foam::distributions::tabulatedDensity::tabulatedDensity
 (
+    const unitConversion& defaultUnits,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
     FieldDistribution<distribution, tabulatedDensity>
     (
         typeName,
+        defaultUnits,
         dict,
-        rndGen,
-        sampleQ
+        sampleQ,
+        std::move(rndGen)
     )
 {
-    List<Pair<scalar>> values(dict.lookup("distribution"));
+    List<Tuple2<scalar, scalar>> values(dict.lookup("distribution"));
 
     // Checks
     forAll(values, i)
@@ -77,11 +79,15 @@ Foam::distributions::tabulatedDensity::tabulatedDensity
         }
     }
 
+    // Optionally read units
+    unitConversion units(defaultUnits);
+    units.readIfPresent("units", dict);
+
     // Copy the coordinates
     x_.resize(values.size());
     forAll(values, i)
     {
-        x_[i] = values[i].first();
+        x_[i] = units.toStandard(values[i].first());
     }
 
     // Copy the PDF. Scale if q != 0.
@@ -102,7 +108,6 @@ Foam::distributions::tabulatedDensity::tabulatedDensity
     // More checks
     validateBounds(dict);
     if (q() != 0) validatePositive(dict);
-    report();
 }
 
 
@@ -163,6 +168,31 @@ Foam::scalar Foam::distributions::tabulatedDensity::mean() const
 }
 
 
+void Foam::distributions::tabulatedDensity::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    FieldDistribution<distribution, tabulatedDensity>::write(os, units);
+
+    // Recover the PDF
+    scalarField PDF(integerPow(x_, -q())*PDF_);
+
+    // Normalise the PDF
+    PDF /= unintegrable::integrate(x_, PDF)->last();
+
+    // Construct and write the values
+    List<Tuple2<scalar, scalar>> values(PDF_.size());
+    forAll(values, i)
+    {
+        values[i].first() = units.toUser(x_[i]);
+        values[i].second() = PDF[i];
+    }
+    writeEntry(os, "distribution", values);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::distributions::tabulatedDensity::x
 (
     const label
diff --git a/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.H b/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.H
index 9513a16d8c..3ebc0aeef9 100644
--- a/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.H
+++ b/src/OpenFOAM/distributions/tabulatedDensity/tabulatedDensity.H
@@ -103,9 +103,10 @@ public:
         //- Construct from a dictionary
         tabulatedDensity
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -139,6 +140,9 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         virtual tmp<scalarField> x(const label n) const;
 
diff --git a/src/OpenFOAM/distributions/uniform/uniform.C b/src/OpenFOAM/distributions/uniform/uniform.C
index fa7501c775..359caadb4c 100644
--- a/src/OpenFOAM/distributions/uniform/uniform.C
+++ b/src/OpenFOAM/distributions/uniform/uniform.C
@@ -57,20 +57,27 @@ Foam::scalar Foam::distributions::uniform::Phi(const scalar x) const
 
 Foam::distributions::uniform::uniform
 (
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
-    FieldDistribution<distribution, uniform>(typeName, dict, rndGen, sampleQ),
-    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"})),
-    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"})),
+    FieldDistribution<distribution, uniform>
+    (
+        typeName,
+        units,
+        dict,
+        sampleQ,
+        std::move(rndGen)
+    ),
+    min_(dict.lookupBackwardsCompatible<scalar>({"min", "minValue"}, units)),
+    max_(dict.lookupBackwardsCompatible<scalar>({"max", "maxValue"}, units)),
     Phi0_(Phi(min_)),
     Phi1_(Phi(max_))
 {
     validateBounds(dict);
     if (q() != 0) validatePositive(dict);
-    report();
 }
 
 
@@ -137,6 +144,19 @@ Foam::scalar Foam::distributions::uniform::mean() const
 }
 
 
+void Foam::distributions::uniform::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    distribution::write(os, units);
+
+    writeEntry(os, "min", units, min_);
+    writeEntry(os, "max", units, max_);
+}
+
+
 Foam::tmp<Foam::scalarField>
 Foam::distributions::uniform::PDF(const scalarField& x) const
 {
diff --git a/src/OpenFOAM/distributions/uniform/uniform.H b/src/OpenFOAM/distributions/uniform/uniform.H
index fdc0d0e32a..2e21d2d18c 100644
--- a/src/OpenFOAM/distributions/uniform/uniform.H
+++ b/src/OpenFOAM/distributions/uniform/uniform.H
@@ -99,9 +99,10 @@ public:
         //- Construct from a dictionary
         uniform
         (
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct copy
@@ -135,6 +136,9 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return the distribution probability density function
         virtual tmp<scalarField> PDF(const scalarField& x) const;
 };
diff --git a/src/OpenFOAM/distributions/unintegrable/unintegrable.C b/src/OpenFOAM/distributions/unintegrable/unintegrable.C
index 4575ab59d6..4ec642cea5 100644
--- a/src/OpenFOAM/distributions/unintegrable/unintegrable.C
+++ b/src/OpenFOAM/distributions/unintegrable/unintegrable.C
@@ -305,25 +305,26 @@ Foam::distributions::unintegrable::Phi01() const
 Foam::distributions::unintegrable::unintegrable
 (
     const word& name,
+    const unitConversion& units,
     const dictionary& dict,
-    randomGenerator& rndGen,
-    const label sampleQ
+    const label sampleQ,
+    randomGenerator&& rndGen
 )
 :
-    distribution(name, dict, rndGen, sampleQ),
+    distribution(name, units, dict, sampleQ, std::move(rndGen)),
     n_((1<<dict.lookupOrDefault<label>("level", 16)) + 1)
 {}
 
 
 Foam::distributions::unintegrable::unintegrable
 (
-    randomGenerator& rndGen,
     const label Q,
     const label sampleQ,
+    randomGenerator&& rndGen,
     const label n
 )
 :
-    distribution(rndGen, Q, sampleQ),
+    distribution(Q, sampleQ, std::move(rndGen)),
     n_(n)
 {}
 
@@ -373,6 +374,22 @@ Foam::scalar Foam::distributions::unintegrable::mean() const
 }
 
 
+void Foam::distributions::unintegrable::write
+(
+    Ostream& os,
+    const unitConversion& units
+) const
+{
+    distribution::write(os, units);
+
+    // Recover the level
+    label n = n_ - 1, level = 0;
+    while (n >>= 1) ++ level;
+
+    writeEntryIfDifferent(os, "level", 16, level);
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::distributions::unintegrable::PDF
 (
     const scalarField& x
diff --git a/src/OpenFOAM/distributions/unintegrable/unintegrable.H b/src/OpenFOAM/distributions/unintegrable/unintegrable.H
index fc09f61e41..03a386e9d3 100644
--- a/src/OpenFOAM/distributions/unintegrable/unintegrable.H
+++ b/src/OpenFOAM/distributions/unintegrable/unintegrable.H
@@ -187,17 +187,18 @@ public:
         unintegrable
         (
             const word& name,
+            const unitConversion& units,
             const dictionary& dict,
-            randomGenerator& rndGen,
-            const label sampleQ
+            const label sampleQ,
+            randomGenerator&& rndGen
         );
 
         //- Construct from components
         unintegrable
         (
-            randomGenerator& rndGen,
             const label Q,
             const label sampleQ,
+            randomGenerator&& rndGen,
             const label n
         );
 
@@ -220,6 +221,9 @@ public:
         //- Return the mean value
         virtual scalar mean() const;
 
+        //- Write to a stream
+        virtual void write(Ostream& os, const unitConversion& units) const;
+
         //- Return coordinates to plot across the range of the distribution
         using distribution::x;
 
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
index db61f055c0..c06df4ebbd 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.C
@@ -30,7 +30,6 @@ License
 #include "constants.H"
 #include "zeroGradientFvPatchFields.H"
 #include "polyMeshTetDecomposition.H"
-#include "standardNormal.H"
 
 using namespace Foam::constant;
 
@@ -651,6 +650,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
     ),
     constProps_(),
     rndGen_(label(149382906)),
+    stdNormal_(rndGen_.generator()),
     boundaryT_
     (
         volScalarField
@@ -905,6 +905,7 @@ Foam::DSMCCloud<ParcelType>::DSMCCloud
     ),
     constProps_(),
     rndGen_(label(971501)),
+    stdNormal_(rndGen_.generator()),
     boundaryT_
     (
         volScalarField
@@ -1043,7 +1044,7 @@ Foam::vector Foam::DSMCCloud<ParcelType>::equipartitionLinearVelocity
 {
     return
         sqrt(physicoChemical::k.value()*temperature/mass)
-       *distributions::standardNormal(rndGen_).sample<vector>();
+       *stdNormal_.sample<vector>();
 }
 
 
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
index 06240813fe..b67153e598 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloud.H
@@ -40,6 +40,7 @@ SourceFiles
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "randomGenerator.H"
+#include "standardNormal.H"
 #include "fvMesh.H"
 #include "volFields.H"
 #include "scalarIOField.H"
@@ -143,6 +144,9 @@ class DSMCCloud
         //- Random number generator
         randomGenerator rndGen_;
 
+        //- Standard normal distribution
+        distributions::standardNormal stdNormal_;
+
 
         // boundary value fields
 
@@ -267,9 +271,12 @@ public:
                 inline const typename ParcelType::constantProperties&
                     constProps(label typeId) const;
 
-                //- Return references to the random object
+                //- Return reference to the random generator
                 inline randomGenerator& rndGen();
 
+                //- Return reference to the standard normal distribution
+                inline distributions::standardNormal& stdNormal();
+
 
             // References to the boundary fields for surface data collection
 
diff --git a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloudI.H b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloudI.H
index 9a66d94106..3f2ca7c6b2 100644
--- a/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloudI.H
+++ b/src/lagrangian/DSMC/clouds/Templates/DSMCCloud/DSMCCloudI.H
@@ -124,6 +124,14 @@ inline Foam::randomGenerator& Foam::DSMCCloud<ParcelType>::rndGen()
 }
 
 
+template<class ParcelType>
+inline Foam::distributions::standardNormal&
+Foam::DSMCCloud<ParcelType>::stdNormal()
+{
+    return stdNormal_;
+}
+
+
 template<class ParcelType>
 inline Foam::volScalarField::Boundary&
 Foam::DSMCCloud<ParcelType>::qBF()
diff --git a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
index 7a2e4cff85..748fcd678e 100644
--- a/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
+++ b/src/lagrangian/DSMC/submodels/InflowBoundaryModel/FreeStream/FreeStream.C
@@ -146,8 +146,8 @@ void Foam::FreeStream<CloudType>::inflow()
 
     const scalar deltaT = mesh.time().deltaTValue();
 
-    randomGenerator& rndGen(cloud.rndGen());
-    distributions::standardNormal stdNormal(rndGen);
+    randomGenerator& rndGen = cloud.rndGen();
+    distributions::standardNormal& stdNormal = cloud.stdNormal();
 
     scalar sqrtPi = sqrt(pi);
 
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
index 4faa887396..2ec16fa23e 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MaxwellianThermal/MaxwellianThermal.C
@@ -81,8 +81,8 @@ void Foam::MaxwellianThermal<CloudType>::correct
 
     CloudType& cloud(this->owner());
 
-    randomGenerator& rndGen(cloud.rndGen());
-    distributions::standardNormal stdNormal(rndGen);
+    randomGenerator& rndGen = cloud.rndGen();
+    distributions::standardNormal& stdNormal = cloud.stdNormal();
 
     while (mag(Ut) < small)
     {
diff --git a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
index 00292f46c5..439978362d 100644
--- a/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
+++ b/src/lagrangian/DSMC/submodels/WallInteractionModel/MixedDiffuseSpecular/MixedDiffuseSpecular.C
@@ -79,8 +79,8 @@ void Foam::MixedDiffuseSpecular<CloudType>::correct
 
     CloudType& cloud(this->owner());
 
-    randomGenerator& rndGen(cloud.rndGen());
-    distributions::standardNormal stdNormal(rndGen);
+    randomGenerator& rndGen = cloud.rndGen();
+    distributions::standardNormal& stdNormal = cloud.stdNormal();
 
     if (diffuseFraction_ > rndGen.scalar01())
     {
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
index a8f90a5041..c02b215dd5 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C
@@ -1080,7 +1080,8 @@ Foam::moleculeCloud::moleculeCloud
     cellOccupancy_(mesh_.nCells()),
     il_(mesh_, pot_.pairPotentials().rCutMax(), false),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(clock::getTime()),
+    stdNormal_(rndGen_.generator())
 {
     if (readFields)
     {
@@ -1110,7 +1111,8 @@ Foam::moleculeCloud::moleculeCloud
     pot_(pot),
     il_(mesh_, 0.0, false),
     constPropList_(),
-    rndGen_(clock::getTime())
+    rndGen_(clock::getTime()),
+    stdNormal_(rndGen_.generator())
 {
     if (readFields)
     {
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
index 5df439ebef..fa6bdb5ce8 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H
@@ -43,6 +43,7 @@ SourceFiles
 #include "InteractionLists.H"
 #include "labelVector.H"
 #include "randomGenerator.H"
+#include "standardNormal.H"
 #include "fileName.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -72,6 +73,8 @@ class moleculeCloud
 
         randomGenerator rndGen_;
 
+        distributions::standardNormal stdNormal_;
+
 
     // Private Member Functions
 
@@ -193,6 +196,8 @@ public:
 
             inline randomGenerator& rndGen();
 
+            inline distributions::standardNormal& stdNormal();
+
 
     // Member Functions
 
diff --git a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
index 5725aecdb2..3477ca9a16 100644
--- a/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
+++ b/src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H
@@ -24,7 +24,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "constants.H"
-#include "standardNormal.H"
 
 using namespace Foam::constant;
 
@@ -304,7 +303,7 @@ inline Foam::vector Foam::moleculeCloud::equipartitionLinearVelocity
 {
     return
         sqrt(physicoChemical::k.value()*temperature/mass)
-       *distributions::standardNormal(rndGen_).sample<vector>();
+       *stdNormal_.sample<vector>();
 }
 
 
@@ -316,24 +315,22 @@ inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
 {
     scalar sqrtKbT = sqrt(physicoChemical::k.value()*temperature);
 
-    distributions::standardNormal stdNormal(rndGen_);
-
     if (cP.linearMolecule())
     {
         return sqrtKbT*vector
         (
             0.0,
-            sqrt(cP.momentOfInertia().yy())*stdNormal.sample(),
-            sqrt(cP.momentOfInertia().zz())*stdNormal.sample()
+            sqrt(cP.momentOfInertia().yy())*stdNormal_.sample(),
+            sqrt(cP.momentOfInertia().zz())*stdNormal_.sample()
         );
     }
     else
     {
         return sqrtKbT*vector
         (
-            sqrt(cP.momentOfInertia().xx())*stdNormal.sample(),
-            sqrt(cP.momentOfInertia().yy())*stdNormal.sample(),
-            sqrt(cP.momentOfInertia().zz())*stdNormal.sample()
+            sqrt(cP.momentOfInertia().xx())*stdNormal_.sample(),
+            sqrt(cP.momentOfInertia().yy())*stdNormal_.sample(),
+            sqrt(cP.momentOfInertia().zz())*stdNormal_.sample()
         );
     }
 }
@@ -387,4 +384,10 @@ inline Foam::randomGenerator& Foam::moleculeCloud::rndGen()
 }
 
 
+inline Foam::distributions::standardNormal& Foam::moleculeCloud::stdNormal()
+{
+    return stdNormal_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.C b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.C
index 9b241304ea..3c86941965 100644
--- a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.C
+++ b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.C
@@ -239,8 +239,6 @@ void Foam::MomentumCloud<CloudType>::cloudReset(MomentumCloud<CloudType>& c)
 {
     CloudType::cloudReset(c);
 
-    rndGen_ = c.rndGen_;
-
     forces_.transfer(c.forces_);
 
     functions_.transfer(c.functions_);
@@ -302,6 +300,7 @@ Foam::MomentumCloud<CloudType>::MomentumCloud
     ),
     cpuLoad_(particleProperties_.lookupOrDefault("cpuLoad", false)),
     rndGen_(0),
+    stdNormal_(rndGen_.generator()),
     cellOccupancyPtr_(),
     cellLengthScale_(mag(cbrt(this->mesh().V()))),
     rho_(rho),
@@ -414,6 +413,7 @@ Foam::MomentumCloud<CloudType>::MomentumCloud
     subModelProperties_(c.subModelProperties_),
     cpuLoad_(c.cpuLoad_),
     rndGen_(c.rndGen_),
+    stdNormal_(c.stdNormal_),
     cellOccupancyPtr_(nullptr),
     cellLengthScale_(c.cellLengthScale_),
     rho_(c.rho_),
@@ -504,6 +504,7 @@ Foam::MomentumCloud<CloudType>::MomentumCloud
     subModelProperties_(dictionary::null),
     cpuLoad_(c.cpuLoad_),
     rndGen_(0),
+    stdNormal_(rndGen_.generator()),
     cellOccupancyPtr_(nullptr),
     cellLengthScale_(c.cellLengthScale_),
     rho_(c.rho_),
diff --git a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.H b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.H
index 3a6b2e8080..b6d23c8e70 100644
--- a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.H
+++ b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloud.H
@@ -56,6 +56,7 @@ SourceFiles
 #include "timeIOdictionary.H"
 #include "autoPtr.H"
 #include "randomGenerator.H"
+#include "standardNormal.H"
 #include "fvMesh.H"
 #include "volFields.H"
 #include "fvMatrices.H"
@@ -161,6 +162,9 @@ protected:
         //- Random number generator - used by some injection routines
         mutable randomGenerator rndGen_;
 
+        //- Standard normal distribution
+        mutable distributions::standardNormal stdNormal_;
+
         //- Cell occupancy information for each parcel, (demand driven)
         autoPtr<List<DynamicList<parcelType*>>> cellOccupancyPtr_;
 
@@ -367,9 +371,12 @@ public:
 
             // Cloud data
 
-                //- Return reference to the random object
+                //- Return reference to the random generator
                 inline randomGenerator& rndGen() const;
 
+                //- Return reference to the standard normal distribution
+                inline distributions::standardNormal& stdNormal() const;
+
                 //- Return the cell occupancy information for each
                 //  parcel, non-const access, the caller is
                 //  responsible for updating it for its own purposes
diff --git a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloudI.H b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloudI.H
index 596a6a0860..2f9a31198a 100644
--- a/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloudI.H
+++ b/src/lagrangian/parcel/clouds/Templates/MomentumCloud/MomentumCloudI.H
@@ -366,6 +366,14 @@ inline Foam::randomGenerator& Foam::MomentumCloud<CloudType>::rndGen() const
 }
 
 
+template<class CloudType>
+inline Foam::distributions::standardNormal&
+Foam::MomentumCloud<CloudType>::stdNormal() const
+{
+    return stdNormal_;
+}
+
+
 template<class CloudType>
 inline Foam::List<Foam::DynamicList<typename CloudType::particleType*>>&
 Foam::MomentumCloud<CloudType>::cellOccupancy()
diff --git a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
index 94d7456800..293234ebd5 100644
--- a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
+++ b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.C
@@ -101,7 +101,7 @@ Foam::vector Foam::GradientDispersionRAS<CloudType>::update
     scalar& tTurb
 )
 {
-    distributions::standardNormal stdNormal(this->owner().rndGen());
+    distributions::standardNormal& stdNormal = this->owner().stdNormal();
 
     const scalar cps = 0.16432;
 
diff --git a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
index 42e24bd7e4..6fb246b06c 100644
--- a/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
+++ b/src/lagrangian/parcel/submodels/Momentum/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C
@@ -73,7 +73,7 @@ Foam::vector Foam::StochasticDispersionRAS<CloudType>::update
 )
 {
     randomGenerator& rndGen = this->owner().rndGen();
-    distributions::standardNormal stdNormal(rndGen);
+    distributions::standardNormal& stdNormal = this->owner().stdNormal();
 
     const scalar cps = 0.16432;
 
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.C
index d15df7ba46..c464e1ecdf 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/CellZoneInjection/CellZoneInjection.C
@@ -173,9 +173,10 @@ Foam::CellZoneInjection<CloudType>::CellZoneInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     )
 {
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.C
index dd9b360b24..0e705d376c 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ConeInjection/ConeInjection.C
@@ -199,9 +199,10 @@ Foam::ConeInjection<CloudType>::ConeInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     ),
     dInner_(vGreat),
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
index 9df0a07b4c..513bca28f5 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/FieldActivatedInjection/FieldActivatedInjection.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,9 +83,10 @@ Foam::FieldActivatedInjection<CloudType>::FieldActivatedInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     )
 {
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.C
index 31e4b37d3c..424cbfd294 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/ManualInjection/ManualInjection.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2024 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,9 +63,10 @@ Foam::ManualInjection<CloudType>::ManualInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     ),
     ignoreOutOfBounds_
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
index cf61fc8057..215792deda 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchFlowRateInjection/PatchFlowRateInjection.C
@@ -61,9 +61,10 @@ Foam::PatchFlowRateInjection<CloudType>::PatchFlowRateInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     )
 {}
diff --git a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.C b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.C
index 218e687704..8beafc3cf8 100644
--- a/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.C
+++ b/src/lagrangian/parcel/submodels/Momentum/InjectionModel/PatchInjection/PatchInjection.C
@@ -98,9 +98,10 @@ Foam::PatchInjection<CloudType>::PatchInjection
     (
         distribution::New
         (
+            dimLength,
             this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen(),
-            this->sizeSampleQ()
+            this->sizeSampleQ(),
+            owner.rndGen().generator()
         )
     )
 {
diff --git a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
index 2869f45072..35ee6d8907 100644
--- a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
+++ b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.C
@@ -70,7 +70,6 @@ Foam::BrownianMotionForce<CloudType>::BrownianMotionForce
 )
 :
     ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
-    rndGen_(owner.rndGen()),
     lambda_(this->coeffs().template lookup<scalar>("lambda")),
     turbulence_(readBool(this->coeffs().lookup("turbulence"))),
     kPtr_(nullptr),
@@ -85,7 +84,6 @@ Foam::BrownianMotionForce<CloudType>::BrownianMotionForce
 )
 :
     ParticleForce<CloudType>(bmf),
-    rndGen_(bmf.rndGen_),
     lambda_(bmf.lambda_),
     turbulence_(bmf.turbulence_),
     kPtr_(nullptr),
@@ -172,7 +170,7 @@ Foam::forceSuSp Foam::BrownianMotionForce<CloudType>::calcCoupled
     }
 
     randomGenerator& rndGen = this->owner().rndGen();
-    distributions::standardNormal stdNormal(rndGen);
+    distributions::standardNormal& stdNormal = this->owner().stdNormal();
 
     // To generate a cubic distribution (i.e., 3 independent directions):
     // value.Su() = f*stdNormal.sample<vector>();
diff --git a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
index 0ee5f1562f..edc3d50700 100644
--- a/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
+++ b/src/lagrangian/parcel/submodels/Thermodynamic/ParticleForces/BrownianMotion/BrownianMotionForce.H
@@ -64,9 +64,6 @@ class BrownianMotionForce
 {
     // Private Data
 
-        //- Reference to the cloud random number generator
-        randomGenerator& rndGen_;
-
         //- Molecular free path length [m]
         const scalar lambda_;
 
diff --git a/src/randomProcesses/processes/OUprocess/OUprocess.C b/src/randomProcesses/processes/OUprocess/OUprocess.C
index 68870835a3..78b990d8a3 100644
--- a/src/randomProcesses/processes/OUprocess/OUprocess.C
+++ b/src/randomProcesses/processes/OUprocess/OUprocess.C
@@ -25,7 +25,6 @@ License
 
 #include "OUprocess.H"
 #include "Kmesh.H"
-#include "standardNormal.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -36,13 +35,11 @@ namespace Foam
 
 complexVector OUprocess::WeinerProcess(const scalar deltaT) const
 {
-    distributions::standardNormal stdNormal(rndGen_);
-
     return sqrt(deltaT)*complexVector
     (
-        complex(stdNormal.sample(), stdNormal.sample()),
-        complex(stdNormal.sample(), stdNormal.sample()),
-        complex(stdNormal.sample(), stdNormal.sample())
+        complex(stdNormal_.sample(), stdNormal_.sample()),
+        complex(stdNormal_.sample(), stdNormal_.sample()),
+        complex(stdNormal_.sample(), stdNormal_.sample())
     );
 }
 
@@ -56,7 +53,7 @@ OUprocess::OUprocess
     const dictionary& OUdict
 )
 :
-    rndGen_(label(0)),
+    stdNormal_(0),
     Kmesh_(kmesh),
     OUfield_(Kmesh_.size()),
 
diff --git a/src/randomProcesses/processes/OUprocess/OUprocess.H b/src/randomProcesses/processes/OUprocess/OUprocess.H
index ae63830535..4b8022f8d3 100644
--- a/src/randomProcesses/processes/OUprocess/OUprocess.H
+++ b/src/randomProcesses/processes/OUprocess/OUprocess.H
@@ -37,7 +37,7 @@ SourceFiles
 
 #include "complexFields.H"
 #include "scalar.H"
-#include "randomGenerator.H"
+#include "standardNormal.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -54,7 +54,7 @@ class OUprocess
 {
     // Private Data
 
-        mutable randomGenerator rndGen_;
+        const distributions::standardNormal stdNormal_;
 
         const Kmesh& Kmesh_;
         mutable complexVectorField OUfield_;