diff --git a/applications/utilities/preProcessing/setWaves/setWaves.C b/applications/utilities/preProcessing/setWaves/setWaves.C
index a87a882b40..85438a2496 100644
--- a/applications/utilities/preProcessing/setWaves/setWaves.C
+++ b/applications/utilities/preProcessing/setWaves/setWaves.C
@@ -30,14 +30,17 @@ Description
 
 \*---------------------------------------------------------------------------*/
 
-#include "fvCFD.H"
+#include "argList.H"
 #include "levelSet.H"
 #include "pointFields.H"
 #include "timeSelector.H"
+#include "uniformDimensionedFields.H"
+#include "volFields.H"
 #include "wallDist.H"
 #include "waveAlphaFvPatchScalarField.H"
 #include "waveVelocityFvPatchVectorField.H"
-#include "waveSuperposition.H"
+
+using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,18 +48,27 @@ int main(int argc, char *argv[])
 {
     timeSelector::addOptions(false, false);
 
-    Foam::argList::addOption
+    #include "addDictOption.H"
+    #include "addRegionOption.H"
+
+    argList::addOption
+    (
+        "alpha",
+        "name",
+        "name of the volume fraction field, default is \"alpha\""
+    );
+
+    argList::addOption
     (
         "U",
         "name",
         "name of the velocity field, default is \"U\""
     );
 
-    Foam::argList::addOption
+    argList::addBoolOption
     (
-        "alpha",
-        "name",
-        "name of the volume fraction field, default is \"alpha\""
+        "gas",
+        "the volume fraction field is that that of the gas above the wave"
     );
 
     #include "setRootCase.H"
@@ -64,11 +76,31 @@ int main(int argc, char *argv[])
 
     instantList timeDirs = timeSelector::selectIfPresent(runTime, args);
 
-    #include "createMesh.H"
+    #include "createNamedMesh.H"
+
+    const word dictName("setWavesDict");
+    #include "setSystemMeshDictionaryIO.H"
+    Info<< "Reading " << dictName << "\n" << endl;
+    IOdictionary setWavesDict(dictIO);
+
     #include "readGravitationalAcceleration.H"
 
     const pointMesh& pMesh = pointMesh::New(mesh);
 
+    // Parse options
+    const word alphaName = setWavesDict.lookupOrDefault<word>("alpha", "alpha");
+    const word UName = setWavesDict.lookupOrDefault<word>("U", "U");
+    const bool liquid = setWavesDict.lookupOrDefault<bool>("liquid", true);
+    const dimensionedVector UMean
+    (
+        "UMean",
+        dimVelocity,
+        setWavesDict.lookup("UMean")
+    );
+
+    // Get the wave models
+    const waveSuperposition& waves = waveSuperposition::New(mesh);
+
     forAll(timeDirs, timeI)
     {
         runTime.setTime(timeDirs[timeI], timeI);
@@ -83,7 +115,7 @@ int main(int argc, char *argv[])
         (
             IOobject
             (
-                args.optionFound("alpha") ? args["alpha"] : "alpha",
+                alphaName,
                 runTime.timeName(),
                 mesh,
                 IOobject::MUST_READ
@@ -94,7 +126,7 @@ int main(int argc, char *argv[])
         (
             IOobject
             (
-                args.optionFound("U") ? args["U"] : "U",
+                UName,
                 runTime.timeName(),
                 mesh,
                 IOobject::MUST_READ
@@ -125,7 +157,7 @@ int main(int argc, char *argv[])
         (
             IOobject("uGasp", runTime.timeName(), mesh),
             pMesh,
-            dimensionedVector("0", dimLength, vector::zero)
+            dimensionedVector("0", dimVelocity, vector::zero)
         );
         volVectorField uLiq
         (
@@ -137,165 +169,31 @@ int main(int argc, char *argv[])
         (
             IOobject("uLiqp", runTime.timeName(), mesh),
             pMesh,
-            dimensionedVector("0", dimLength, vector::zero)
+            dimensionedVector("0", dimVelocity, vector::zero)
         );
 
-        // The number of wave patches
-        label nWaves = 0;
+        // Offset
+        const vector offset = UMean.value()*t;
 
-        // Whether the alpha conditions refer to the liquid phase
-        bool liquid = false;
+        // Cell centres and points
+        const pointField& ccs = mesh.cellCentres();
+        const pointField& pts = mesh.points();
 
-        // Loop the patches, averaging and superimposing wave model data
-        forAll(mesh.boundary(), patchi)
+        // Internal field
+        h.primitiveFieldRef() = waves.height(t, ccs + offset);
+        hp.primitiveFieldRef() = waves.height(t, pts + offset);
+        uGas.primitiveFieldRef() = waves.UGas(t, ccs + offset);
+        uGasp.primitiveFieldRef() = waves.UGas(t, pts + offset);
+        uLiq.primitiveFieldRef() = waves.ULiquid(t, ccs + offset);
+        uLiqp.primitiveFieldRef() = waves.ULiquid(t, pts + offset);
+
+        // Boundary fields
+        forAll(mesh.boundary(), patchj)
         {
-            fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
-            fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
-
-            const bool isWave = isA<waveAlphaFvPatchScalarField>(alphap);
-
-            if (isA<waveVelocityFvPatchVectorField>(Up) != isWave)
-            {
-                FatalErrorInFunction
-                    << "The alpha condition on patch " << Up.patch().name()
-                    << " is " << alphap.type() << " and the velocity condition"
-                    << " is " << Up.type() << ". Wave boundary conditions must"
-                    << " be set in pairs. If the alpha condition is "
-                    << waveAlphaFvPatchScalarField::typeName
-                    << " then the velocity condition must be "
-                    << waveVelocityFvPatchVectorField::typeName
-                    << " and vice-versa." << exit(FatalError);
-            }
-
-            if (!isWave)
-            {
-                continue;
-            }
-
-            Info<< "Adding waves from patch " << Up.patch().name() << endl;
-
-            const waveSuperposition& waves =
-                refCast<waveVelocityFvPatchVectorField>(Up).waves();
-
-            const bool liquidp =
-                refCast<waveAlphaFvPatchScalarField>(alphap).liquid();
-            if (nWaves > 0 && liquidp != liquid)
-            {
-                FatalErrorInFunction
-                    << "All " << waveAlphaFvPatchScalarField::typeName
-                    << "patch fields must be configured for the same phase,"
-                    << " i.e., the liquid switch must have the same value."
-                    << exit(FatalError);
-            }
-            liquid = liquidp;
-
-            const pointField& ccs = mesh.cellCentres();
-            const pointField& pts = mesh.points();
-
-            // Internal field superposition
-            h.primitiveFieldRef() += waves.height(t, ccs);
-            hp.primitiveFieldRef() += waves.height(t, pts);
-            uGas.primitiveFieldRef() += waves.UGas(t, ccs) - waves.UMean(t);
-            uGasp.primitiveFieldRef() += waves.UGas(t, pts) - waves.UMean(t);
-            uLiq.primitiveFieldRef() += waves.ULiquid(t, ccs) - waves.UMean(t);
-            uLiqp.primitiveFieldRef() += waves.ULiquid(t, pts) - waves.UMean(t);
-
-            // Boundary field superposition
-            forAll(mesh.boundary(), patchj)
-            {
-                const pointField& fcs = mesh.boundary()[patchj].Cf();
-                h.boundaryFieldRef()[patchj] += waves.height(t, fcs);
-                uGas.boundaryFieldRef()[patchj] +=
-                    waves.UGas(t, fcs) - waves.UMean(t);
-                uLiq.boundaryFieldRef()[patchj] +=
-                    waves.ULiquid(t, fcs) - waves.UMean(t);
-            }
-
-            ++ nWaves;
-        }
-
-        // Create the mean velocity field
-        volVectorField UMean
-        (
-            IOobject("UMean", runTime.timeName(), mesh),
-            mesh,
-            dimensionedVector("UMean", dimVelocity, Zero)
-        );
-
-        if (nWaves == 0)
-        {
-            // Warn and skip to the next time if there are no wave patches
-            WarningInFunction
-                << "No " << waveAlphaFvPatchScalarField::typeName << " or "
-                << waveVelocityFvPatchVectorField::typeName << " patch fields "
-                << "were found. No waves have been set." << endl;
-
-            continue;
-        }
-        else if (nWaves == 1)
-        {
-            // Set a mean velocity equal to that on the only wave patch
-            forAll(mesh.boundary(), patchi)
-            {
-                const fvPatchVectorField& Up = U.boundaryField()[patchi];
-                if (!isA<waveVelocityFvPatchVectorField>(Up))
-                {
-                    continue;
-                }
-
-                const waveSuperposition& waves =
-                    refCast<const waveVelocityFvPatchVectorField>(Up).waves();
-
-                UMean ==
-                    dimensionedVector("UMean", dimVelocity, waves.UMean(t));
-            }
-        }
-        else if (nWaves > 1)
-        {
-            // Set the mean velocity by distance weighting from the wave patches
-
-            // Create weighted average fields for the mean velocity
-            volScalarField weight
-            (
-                IOobject("weight", runTime.timeName(), mesh),
-                mesh,
-                dimensionedScalar("0", dimless/dimLength, 0)
-            );
-            volVectorField weightUMean
-            (
-                IOobject("weightUMean", runTime.timeName(), mesh),
-                mesh,
-                dimensionedVector("0", dimVelocity/dimLength, vector::zero)
-            );
-
-            // Loop the patches, inverse-distance weighting the mean velocities
-            forAll(mesh.boundary(), patchi)
-            {
-                const fvPatchVectorField& Up = U.boundaryField()[patchi];
-                if (!isA<waveVelocityFvPatchVectorField>(Up))
-                {
-                    continue;
-                }
-
-                const waveSuperposition& waves =
-                    refCast<const waveVelocityFvPatchVectorField>(Up).waves();
-
-                const volScalarField w
-                (
-                    1
-                   /(
-                       wallDist(mesh, labelList(1, patchi)).y()
-                     + dimensionedScalar("ySmall", dimLength, small)
-                    )
-                );
-
-                weight += w;
-                weightUMean +=
-                    w*dimensionedVector("wUMean", dimVelocity, waves.UMean(t));
-            }
-
-            // Complete the average for the mean velocity
-            UMean = weightUMean/weight;
+            const pointField& fcs = mesh.boundary()[patchj].Cf();
+            h.boundaryFieldRef()[patchj] = waves.height(t, fcs + offset);
+            uGas.boundaryFieldRef()[patchj] = waves.UGas(t, fcs + offset);
+            uLiq.boundaryFieldRef()[patchj] = waves.ULiquid(t, fcs + offset);
         }
 
         // Set the fields
@@ -306,10 +204,13 @@ int main(int argc, char *argv[])
         forAll(mesh.boundary(), patchi)
         {
             fvPatchScalarField& alphap = alpha.boundaryFieldRef()[patchi];
-            fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
             if (isA<waveAlphaFvPatchScalarField>(alphap))
             {
                 alphap == refCast<waveAlphaFvPatchScalarField>(alphap).alpha();
+            }
+            fvPatchVectorField& Up = U.boundaryFieldRef()[patchi];
+            if (isA<waveVelocityFvPatchVectorField>(Up))
+            {
                 Up == refCast<waveVelocityFvPatchVectorField>(Up).U();
             }
         }
diff --git a/etc/caseDicts/annotated/setWavesDict b/etc/caseDicts/annotated/setWavesDict
new file mode 100644
index 0000000000..b603c25136
--- /dev/null
+++ b/etc/caseDicts/annotated/setWavesDict
@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setWavesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Optional name of the volume fraction field
+alpha   alpha.water;
+
+// Optional name of the velocity field
+U       U;
+
+// Optional flag which indicates whether the alpha field is that of the liquid
+// under the waves (true) or the gas over the waves (false)
+liquid  true;
+
+// The mean flow velocity over which to superimpose waves
+UMean   (1 0 0);
+
+
+// ************************************************************************* //
diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
index 77ad8d2153..c1ee40ee97 100644
--- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
+++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
@@ -122,20 +122,19 @@ Foam::waveAlphaFvPatchScalarField::waveAlphaFvPatchScalarField
 Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alpha() const
 {
     const scalar t = db().time().timeOutputValue();
-
+    const waveSuperposition& waves = waveSuperposition::New(db());
     const waveVelocityFvPatchVectorField& Up =
         refCast<const waveVelocityFvPatchVectorField>
         (
             patch().lookupPatchField<volVectorField, scalar>(UName_)
         );
-    const waveSuperposition& waves = Up.waves();
 
     return
         levelSetFraction
         (
             patch(),
-            waves.height(t, patch().Cf()),
-            waves.height(t, patch().patch().localPoints()),
+            waves.height(t, patch().Cf() + Up.offset()),
+            waves.height(t, patch().patch().localPoints() + Up.offset()),
             !liquid_
         );
 }
@@ -143,14 +142,14 @@ Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alpha() const
 
 Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alphan() const
 {
+    const scalar t = db().time().timeOutputValue();
+    const waveSuperposition& waves = waveSuperposition::New(db());
     const waveVelocityFvPatchVectorField& Up =
         refCast<const waveVelocityFvPatchVectorField>
         (
             patch().lookupPatchField<volVectorField, scalar>(UName_)
         );
 
-    const scalar t = db().time().timeOutputValue();
-
     const fvMeshSubset& subset = Up.faceCellSubset();
     const fvMesh& meshs = subset.subMesh();
     const label patchis = findIndex(subset.patchMap(), patch().index());
@@ -160,8 +159,8 @@ Foam::tmp<Foam::scalarField> Foam::waveAlphaFvPatchScalarField::alphan() const
         levelSetFraction
         (
             meshs,
-            Up.waves().height(t, meshs.cellCentres())(),
-            Up.waves().height(t, meshs.points())(),
+            waves.height(t, meshs.cellCentres() + Up.offset())(),
+            waves.height(t, meshs.points() + Up.offset())(),
             !liquid_
         )
     );
diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
index 7c2be193c4..64006317e2 100644
--- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
+++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
@@ -27,7 +27,8 @@ Class
 Description
     This boundary condition provides a waveAlpha condition. This sets the phase
     fraction to that specified by a superposition of wave models. All the
-    parameters are looked up from the corresponding velocity condition.
+    wave modelling parameters are obtained from a centrally registered
+    waveSuperposition class.
 
     Flow reversal will occur in the event that the amplitude of the velocity
     oscillation is greater than the mean flow. This triggers special handling,
@@ -58,10 +59,14 @@ Usage
     {
         type        waveAlpha;
         U           U;
+        liquid      true;
         inletOutlet true;
     }
     \endverbatim
 
+See also
+    Foam::waveSuperposition
+
 SourceFiles
     waveAlphaFvPatchScalarField.C
 
diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C
index 5822e273a4..7fed7ae8f9 100644
--- a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C
+++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.C
@@ -115,38 +115,38 @@ Foam::wavePressureFvPatchScalarField::wavePressureFvPatchScalarField
 
 Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::p() const
 {
+    const scalar t = db().time().timeOutputValue();
+    const waveSuperposition& waves = waveSuperposition::New(db());
     const waveVelocityFvPatchVectorField& Up =
         refCast<const waveVelocityFvPatchVectorField>
         (
             patch().lookupPatchField<volVectorField, scalar>(UName_)
         );
 
-    const scalar t = db().time().timeOutputValue();
-
     return
         levelSetAverage
         (
             patch(),
-            Up.waves().height(t, patch().Cf()),
-            Up.waves().height(t, patch().patch().localPoints()),
-            Up.waves().pGas(t, patch().Cf())(),
-            Up.waves().pGas(t, patch().patch().localPoints())(),
-            Up.waves().pLiquid(t, patch().Cf())(),
-            Up.waves().pLiquid(t, patch().patch().localPoints())()
+            waves.height(t, patch().Cf() + Up.offset()),
+            waves.height(t, patch().patch().localPoints() + Up.offset()),
+            waves.pGas(t, patch().Cf() + Up.offset())(),
+            waves.pGas(t, patch().patch().localPoints() + Up.offset())(),
+            waves.pLiquid(t, patch().Cf() + Up.offset())(),
+            waves.pLiquid(t, patch().patch().localPoints() + Up.offset())()
         );
 }
 
 
 Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::pn() const
 {
+    const scalar t = db().time().timeOutputValue();
+    const waveSuperposition& waves = waveSuperposition::New(db());
     const waveVelocityFvPatchVectorField& Up =
         refCast<const waveVelocityFvPatchVectorField>
         (
             patch().lookupPatchField<volVectorField, scalar>(UName_)
         );
 
-    const scalar t = db().time().timeOutputValue();
-
     const fvMeshSubset& subset = Up.faceCellSubset();
     const fvMesh& meshs = subset.subMesh();
     const label patchis = findIndex(subset.patchMap(), patch().index());
@@ -156,12 +156,12 @@ Foam::tmp<Foam::scalarField> Foam::wavePressureFvPatchScalarField::pn() const
         levelSetAverage
         (
             meshs,
-            Up.waves().height(t, meshs.cellCentres())(),
-            Up.waves().height(t, meshs.points())(),
-            Up.waves().pGas(t, meshs.cellCentres())(),
-            Up.waves().pGas(t, meshs.points())(),
-            Up.waves().pLiquid(t, meshs.cellCentres())(),
-            Up.waves().pLiquid(t, meshs.points())()
+            waves.height(t, meshs.cellCentres() + Up.offset())(),
+            waves.height(t, meshs.points() + Up.offset())(),
+            waves.pGas(t, meshs.cellCentres() + Up.offset())(),
+            waves.pGas(t, meshs.points() + Up.offset())(),
+            waves.pLiquid(t, meshs.cellCentres() + Up.offset())(),
+            waves.pLiquid(t, meshs.points() + Up.offset())()
         )
     );
 
diff --git a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H
index 5199c5d35c..ffbfa05148 100644
--- a/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H
+++ b/src/waves/derivedFvPatchFields/wavePressure/wavePressureFvPatchScalarField.H
@@ -27,7 +27,8 @@ Class
 Description
     This boundary condition provides a wavePressure condition. This sets the
     pressure to a value specified by a superposition of wave models. All the
-    parameters are looked up from the corresponding velocity condition.
+    wave modelling parameters are obtained from a centrally registered
+    waveSuperposition class.
 
     This functions like an outletInlet condition. Faces on which the flow is
     leaving the domain have a value set by the wave model. Faces on which the
@@ -55,6 +56,8 @@ Usage
     }
     \endverbatim
 
+See also
+    Foam::waveSuperposition
 
 SourceFiles
     wavePressureFvPatchScalarField.C
diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C
index 76eddfec3f..75a34f4b84 100644
--- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C
+++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C
@@ -42,7 +42,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_("phi"),
     pName_("p"),
     inletOutlet_(true),
-    waves_(db()),
+    UMean_(nullptr),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {
@@ -63,7 +63,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     pName_(dict.lookupOrDefault<word>("p", "p")),
     inletOutlet_(dict.lookupOrDefault<Switch>("inletOutlet", true)),
-    waves_(db(), dict),
+    UMean_(Function1<vector>::New("UMean", dict)),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {
@@ -94,7 +94,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    waves_(ptf.waves_),
+    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -109,7 +109,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    waves_(ptf.waves_),
+    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -125,7 +125,7 @@ Foam::waveVelocityFvPatchVectorField::waveVelocityFvPatchVectorField
     phiName_(ptf.phiName_),
     pName_(ptf.pName_),
     inletOutlet_(ptf.inletOutlet_),
-    waves_(ptf.waves_),
+    UMean_(ptf.UMean_, false),
     faceCellSubset_(nullptr),
     faceCellSubsetTimeIndex_(-1)
 {}
@@ -162,17 +162,19 @@ Foam::waveVelocityFvPatchVectorField::faceCellSubset() const
 Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::U() const
 {
     const scalar t = db().time().timeOutputValue();
+    const waveSuperposition& waves = waveSuperposition::New(db());
 
     return
-        levelSetAverage
+        UMean_->value(t)
+      + levelSetAverage
         (
             patch(),
-            waves_.height(t, patch().Cf()),
-            waves_.height(t, patch().patch().localPoints()),
-            waves_.UGas(t, patch().Cf())(),
-            waves_.UGas(t, patch().patch().localPoints())(),
-            waves_.ULiquid(t, patch().Cf())(),
-            waves_.ULiquid(t, patch().patch().localPoints())()
+            waves.height(t, patch().Cf() + offset()),
+            waves.height(t, patch().patch().localPoints() + offset()),
+            waves.UGas(t, patch().Cf() + offset())(),
+            waves.UGas(t, patch().patch().localPoints() + offset())(),
+            waves.ULiquid(t, patch().Cf() + offset())(),
+            waves.ULiquid(t, patch().patch().localPoints() + offset())()
         );
 }
 
@@ -180,6 +182,7 @@ Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::U() const
 Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::Un() const
 {
     const scalar t = db().time().timeOutputValue();
+    const waveSuperposition& waves = waveSuperposition::New(db());
 
     const fvMeshSubset& subset = faceCellSubset();
     const fvMesh& meshs = subset.subMesh();
@@ -187,15 +190,16 @@ Foam::tmp<Foam::vectorField> Foam::waveVelocityFvPatchVectorField::Un() const
 
     const vectorField Us
     (
-        levelSetAverage
+        UMean_->value(t)
+      + levelSetAverage
         (
             meshs,
-            waves_.height(t, meshs.cellCentres())(),
-            waves_.height(t, meshs.points())(),
-            waves_.UGas(t, meshs.cellCentres())(),
-            waves_.UGas(t, meshs.points())(),
-            waves_.ULiquid(t, meshs.cellCentres())(),
-            waves_.ULiquid(t, meshs.points())()
+            waves.height(t, meshs.cellCentres() + offset())(),
+            waves.height(t, meshs.points() + offset())(),
+            waves.UGas(t, meshs.cellCentres() + offset())(),
+            waves.UGas(t, meshs.points() + offset())(),
+            waves.ULiquid(t, meshs.cellCentres() + offset())(),
+            waves.ULiquid(t, meshs.points() + offset())()
         )
     );
 
@@ -293,7 +297,7 @@ void Foam::waveVelocityFvPatchVectorField::write
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "p", "p", pName_);
     writeEntryIfDifferent<Switch>(os, "inletOutlet", true, inletOutlet_);
-    waves_.write(os);
+    UMean_->writeData(os);
 }
 
 
diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
index dffcab052a..71de15223a 100644
--- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
+++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
@@ -26,9 +26,9 @@ Class
 
 Description
     This boundary condition provides a waveVelocity condition. This sets the
-    velocity to that specified by a superposition of wave models. The
-    corresponding phase fraction and pressure conditions look this condition up
-    and re-use the wave modelling.
+    velocity to that specified by a superposition of wave models. All the
+    wave modelling parameters are obtained from a centrally registered
+    waveSuperposition class.
 
     Flow reversal will occur in the event that the amplitude of the velocity
     oscillation is greater than the mean flow. This triggers special handling,
@@ -62,14 +62,8 @@ Usage
         phi         | Name of the flux field                 | no     | phi
         p           | Name of the pressure field             | no     | p
         inletOutlet | does the condition behave like inletOutlet | no | true
-        origin      | origin of the wave coordinate system   | yes    |
-        direction   | direction of the mean flow             | yes    |
-        speed       | speed of the mean flow                 | yes    |
-        waves       | list of wave models to superimpose     | yes    |
+        UMean       | velocity of the mean flow              | yes    |
         ramp        | ramping function for the mean flow speed | no   | None
-        scale       | scale factor along the mean flow direction | no | None
-        crossScale  | scale factor across the mean flow direction | no | None
-        heightAboveWave | use with the height above the wave | no     | false
     \endtable
 
     Example of the boundary condition specification:
@@ -77,32 +71,17 @@ Usage
     <patchName>
     {
         type        waveVelocity;
-        origin      (0 25 0);
-        direction   (1 0 0);
-        speed       2;
-        waves
-        (
-            Airy
-            {
-                length      40;
-                amplitude   0.5;
-                phase       0;
-                angle       0;
-            }
-            Airy
-            {
-                length      20;
-                amplitude   0.25;
-                phase       1.5708;
-                angle       0;
-            }
-        );
+        phi         phi;
+        p           p;
+        inletOutlet yes;
+        UMean       (2 0 0);
         ramp        constant 1;
-        scale       table ((100 1) (200 0));
-        crossScale  constant 1;
     }
     \endverbatim
 
+See also
+    Foam::waveSuperposition
+
 SourceFiles
     waveVelocityFvPatchVectorField.C
 
@@ -140,8 +119,8 @@ class waveVelocityFvPatchVectorField
         //- Act as an inlet/outlet patch?
         const Switch inletOutlet_;
 
-        //- Wave superposition
-        const waveSuperposition waves_;
+        //- Mean velocity
+        const autoPtr<Function1<vector>> UMean_;
 
         //- Mesh subset corresponding to the patch adjacent cells
         mutable autoPtr<fvMeshSubset> faceCellSubset_;
@@ -227,12 +206,6 @@ public:
                 return pName_;
             }
 
-            //- Access the wave models
-            const waveSuperposition& waves() const
-            {
-                return waves_;
-            }
-
             //- Access the face-cell subset
             const fvMeshSubset& faceCellSubset() const;
 
@@ -245,6 +218,18 @@ public:
             //- Return the current modelled velocity field in the neighbour cell
             tmp<vectorField> Un() const;
 
+            //- Return the mean velocity
+            inline vector UMean() const
+            {
+                return UMean_->value(db().time().timeOutputValue());
+            }
+
+            //- Return the distance offset for the wave models
+            inline vector offset() const
+            {
+                return - UMean_->integrate(0, db().time().timeOutputValue());
+            }
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
diff --git a/src/waves/waveModels/Airy/Airy.C b/src/waves/waveModels/Airy/Airy.C
index 2d7e3e68dc..e1065948af 100644
--- a/src/waves/waveModels/Airy/Airy.C
+++ b/src/waves/waveModels/Airy/Airy.C
@@ -56,11 +56,10 @@ Foam::scalar Foam::waveModels::Airy::celerity() const
 Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::angle
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
-    return phase_ + k()*(x - (u + celerity())*t);
+    return phase_ + k()*(x - celerity()*t);
 }
 
 
@@ -74,14 +73,13 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::vi
 (
     const label i,
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
     const scalarField x(xz.component(0));
     const scalarField z(xz.component(1));
 
-    const scalarField phi(angle(t, u, x));
+    const scalarField phi(angle(t, x));
     const scalarField kz(k()*z);
 
     if (deep())
@@ -132,38 +130,35 @@ Foam::waveModels::Airy::~Airy()
 Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::elevation
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
-    return amplitude(t)*cos(angle(t, u, x));
+    return amplitude(t)*cos(angle(t, x));
 }
 
 
 Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::velocity
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
     const scalar ka = k()*amplitude(t);
 
-    return celerity()*ka*vi(1, t, u, xz);
+    return celerity()*ka*vi(1, t, xz);
 }
 
 
 Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::pressure
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
     // It is a fluke of the formulation that the time derivative of the velocity
     // potential equals the x-derivative multiplied by the celerity. This allows
     // for this shortcut in evaluating the unsteady pressure.
-    return celerity()*velocity(t, u, xz)->component(0);
+    return celerity()*velocity(t, xz)->component(0);
 }
 
 
diff --git a/src/waves/waveModels/Airy/Airy.H b/src/waves/waveModels/Airy/Airy.H
index 3c30065df7..214f9812a6 100644
--- a/src/waves/waveModels/Airy/Airy.H
+++ b/src/waves/waveModels/Airy/Airy.H
@@ -87,7 +87,6 @@ protected:
         tmp<scalarField> angle
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
@@ -99,7 +98,6 @@ protected:
         (
             const label i,
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 
@@ -151,32 +149,29 @@ public:
                 return depth_;
             }
 
-        //- Get the wave elevation at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity.
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
         virtual tmp<scalarField> elevation
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<vector2DField> velocity
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 
-        //- Get the wave pressure at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave pressure at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<scalarField> pressure
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 
diff --git a/src/waves/waveModels/Stokes2/Stokes2.C b/src/waves/waveModels/Stokes2/Stokes2.C
index c5499e34b6..641783cfd5 100644
--- a/src/waves/waveModels/Stokes2/Stokes2.C
+++ b/src/waves/waveModels/Stokes2/Stokes2.C
@@ -61,7 +61,6 @@ Foam::waveModels::Stokes2::~Stokes2()
 Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes2::elevation
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
@@ -77,15 +76,14 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes2::elevation
     }
 
     return
-        Airy::elevation(t, u, x)
-      + (1/k())*sqr(ka)*B22*cos(2*angle(t, u, x));
+        Airy::elevation(t, x)
+      + (1/k())*sqr(ka)*B22*cos(2*angle(t, x));
 }
 
 
 Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes2::velocity
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
@@ -100,8 +98,8 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes2::velocity
     }
 
     return
-        Airy::velocity(t, u, xz)
-      + celerity()*sqr(ka)*A22ByA11*vi(2, t, u, xz);
+        Airy::velocity(t, xz)
+      + celerity()*sqr(ka)*A22ByA11*vi(2, t, xz);
 }
 
 
diff --git a/src/waves/waveModels/Stokes2/Stokes2.H b/src/waves/waveModels/Stokes2/Stokes2.H
index 99817bd5bd..bd059ec195 100644
--- a/src/waves/waveModels/Stokes2/Stokes2.H
+++ b/src/waves/waveModels/Stokes2/Stokes2.H
@@ -85,22 +85,20 @@ public:
 
     // Member Functions
 
-        //- Get the wave elevation at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity.
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
         virtual tmp<scalarField> elevation
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<vector2DField> velocity
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 };
diff --git a/src/waves/waveModels/Stokes5/Stokes5.C b/src/waves/waveModels/Stokes5/Stokes5.C
index 5ec6001bc9..dfafb0273d 100644
--- a/src/waves/waveModels/Stokes5/Stokes5.C
+++ b/src/waves/waveModels/Stokes5/Stokes5.C
@@ -61,7 +61,6 @@ Foam::waveModels::Stokes5::~Stokes5()
 Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
@@ -110,10 +109,10 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
             << "B55 = " << B55 << endl;
     }
 
-    const scalarField phi(angle(t, u, x));
+    const scalarField phi(angle(t, x));
 
     return
-        Stokes2::elevation(t, u, x)
+        Stokes2::elevation(t, x)
       + (1/k())
        *(
            pow3(ka)*B31*(cos(phi) - cos(3*phi))
@@ -126,7 +125,6 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
 Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes5::velocity
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
@@ -192,15 +190,15 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes5::velocity
             << "A55 = " << A55ByA11*A11 << endl;
     }
 
-    const vector2DField v1(vi(1, t, u, xz)), v3(vi(3, t, u, xz));
+    const vector2DField v1(vi(1, t, xz)), v3(vi(3, t, xz));
 
     return
-        Stokes2::velocity(t, u, xz)
+        Stokes2::velocity(t, xz)
       + celerity()
        *(
             pow3(ka)*(A31ByA11*v1 + A33ByA11*v3)
-          + pow4(ka)*(A42ByA11*vi(2, t, u, xz) + A44ByA11*vi(4, t, u, xz))
-          + pow5(ka)*(A51ByA11*v1 + A53ByA11*v3 + A55ByA11*vi(5, t, u, xz))
+          + pow4(ka)*(A42ByA11*vi(2, t, xz) + A44ByA11*vi(4, t, xz))
+          + pow5(ka)*(A51ByA11*v1 + A53ByA11*v3 + A55ByA11*vi(5, t, xz))
         );
 }
 
diff --git a/src/waves/waveModels/Stokes5/Stokes5.H b/src/waves/waveModels/Stokes5/Stokes5.H
index 7194deb8b7..f0c00b3abc 100644
--- a/src/waves/waveModels/Stokes5/Stokes5.H
+++ b/src/waves/waveModels/Stokes5/Stokes5.H
@@ -84,22 +84,20 @@ public:
 
     // Member Functions
 
-        //- Get the wave elevation at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity.
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
         virtual tmp<scalarField> elevation
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<vector2DField> velocity
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 };
diff --git a/src/waves/waveModels/solitary/solitary.C b/src/waves/waveModels/solitary/solitary.C
index 50cd3bb9ea..bc9cd4bd95 100644
--- a/src/waves/waveModels/solitary/solitary.C
+++ b/src/waves/waveModels/solitary/solitary.C
@@ -62,24 +62,22 @@ Foam::scalar Foam::waveModels::solitary::celerity(const scalar t) const
 Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::parameter
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
-    return k(t)*(x - offset_ - (u + celerity(t))*t);
+    return k(t)*(x - offset_ - celerity(t)*t);
 }
 
 
 Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::Pi
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
     const scalar clip = log(great);
 
-    return 1/sqr(cosh(max(- clip, min(clip, parameter(t, u, x)))));
+    return 1/sqr(cosh(max(- clip, min(clip, parameter(t, x)))));
 }
 
 
@@ -116,24 +114,22 @@ Foam::waveModels::solitary::~solitary()
 Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::elevation
 (
     const scalar t,
-    const scalar u,
     const scalarField& x
 ) const
 {
-    return amplitude(t)*Pi(t, u, x);
+    return amplitude(t)*Pi(t, x);
 }
 
 
 Foam::tmp<Foam::vector2DField> Foam::waveModels::solitary::velocity
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
     const scalar A = alpha(t);
     const scalarField Z(max(scalar(0), 1 + xz.component(1)/depth()));
-    const scalarField P(Pi(t, u, xz.component(0)));
+    const scalarField P(Pi(t, xz.component(0)));
 
     return
         celerity(t)
@@ -144,7 +140,7 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::solitary::velocity
                 (4 + 2*A - 6*A*sqr(Z))*P
               + (- 7*A + 9*A*sqr(Z))*sqr(P)
             ),
-            A*Z*depth()*k(t)*tanh(parameter(t, u, xz.component(0)))
+            A*Z*depth()*k(t)*tanh(parameter(t, xz.component(0)))
            *(
                 (2 + A - A*sqr(Z))*P
               + (- 7*A + 3*A*sqr(Z))*sqr(P)
@@ -156,7 +152,6 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::solitary::velocity
 Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::pressure
 (
     const scalar t,
-    const scalar u,
     const vector2DField& xz
 ) const
 {
diff --git a/src/waves/waveModels/solitary/solitary.H b/src/waves/waveModels/solitary/solitary.H
index d598154628..19c18a41d4 100644
--- a/src/waves/waveModels/solitary/solitary.H
+++ b/src/waves/waveModels/solitary/solitary.H
@@ -86,7 +86,6 @@ class solitary
         tmp<scalarField> parameter
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
@@ -94,7 +93,6 @@ class solitary
         tmp<scalarField> Pi
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
@@ -140,32 +138,29 @@ public:
                 return depth_;
             }
 
-        //- Get the wave elevation at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity.
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
         virtual tmp<scalarField> elevation
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<vector2DField> velocity
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 
-        //- Get the wave pressure at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave pressure at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<scalarField> pressure
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const;
 
diff --git a/src/waves/waveModels/waveModel/waveModel.H b/src/waves/waveModels/waveModel/waveModel.H
index 078a25dc7e..42182ccde0 100644
--- a/src/waves/waveModels/waveModel/waveModel.H
+++ b/src/waves/waveModels/waveModel/waveModel.H
@@ -128,32 +128,29 @@ public:
         //- Get the (scalar) value of gravity.
         scalar g() const;
 
-        //- Get the wave elevation at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity.
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
         virtual tmp<scalarField> elevation
         (
             const scalar t,
-            const scalar u,
             const scalarField& x
         ) const = 0;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<vector2DField> velocity
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const = 0;
 
-        //- Get the wave pressure at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, and z with
-        //  negative gravity.
+        //- Get the wave pressure at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation, and z with negative
+        //  gravity.
         virtual tmp<scalarField> pressure
         (
             const scalar t,
-            const scalar u,
             const vector2DField& xz
         ) const = 0;
 
diff --git a/src/waves/waveSuperposition/waveSuperposition.C b/src/waves/waveSuperposition/waveSuperposition.C
index 275e204d7b..2be492b542 100644
--- a/src/waves/waveSuperposition/waveSuperposition.C
+++ b/src/waves/waveSuperposition/waveSuperposition.C
@@ -26,18 +26,42 @@ License
 #include "waveSuperposition.H"
 #include "uniformDimensionedFields.H"
 
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+const Foam::word Foam::waveSuperposition::dictName("waveProperties");
+
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+const Foam::waveSuperposition& Foam::waveSuperposition::New
+(
+    const objectRegistry& db
+)
+{
+    if (db.foundObject<waveSuperposition>(dictName))
+    {
+        return db.lookupObject<waveSuperposition>(dictName);
+    }
+    else
+    {
+        waveSuperposition* ptr = new waveSuperposition(db);
+        ptr->store();
+        return *ptr;
+    }
+}
+
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::waveSuperposition::transformation
 (
     const vectorField& p,
     tensor& axes,
-    scalar& u,
     vectorField& xyz
 ) const
 {
     const uniformDimensionedVectorField& g =
-        db_.lookupObject<uniformDimensionedVectorField>("g");
+        db().lookupObject<uniformDimensionedVectorField>("g");
     const scalar magG = mag(g.value());
     const vector gHat = g.value()/magG;
 
@@ -47,8 +71,6 @@ void Foam::waveSuperposition::transformation
 
     axes = tensor(dSurfHat, - gHat ^ dSurfHat, - gHat);
 
-    u = speed_*magDSurf;
-
     xyz = axes & (p - origin_);
 }
 
@@ -64,7 +86,7 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::elevation
     forAll(waveModels_, wavei)
     {
         const vector2D d(cos(waveAngles_[wavei]), sin(waveAngles_[wavei]));
-        result += waveModels_[wavei].elevation(t, d.x()*speed_, d & xy);
+        result += waveModels_[wavei].elevation(t, d & xy);
     }
 
     return scale(xy)*result;
@@ -92,7 +114,7 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::velocity
         );
         const vector2DField uw
         (
-            waveModels_[wavei].velocity(t, d.x()*speed_, xz)
+            waveModels_[wavei].velocity(t, xz)
         );
         result += zip
         (
@@ -129,9 +151,9 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pressure
         );
         const vector2DField uw
         (
-            waveModels_[wavei].velocity(t, d.x()*speed_, xz)
+            waveModels_[wavei].velocity(t, xz)
         );
-        result += waveModels_[wavei].pressure(t, d.x()*speed_, xz);
+        result += waveModels_[wavei].pressure(t, xz);
     }
 
     tmp<scalarField> s = scale(zip(xyz.component(0), xyz.component(1)));
@@ -174,67 +196,36 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::scale
 
 Foam::waveSuperposition::waveSuperposition(const objectRegistry& db)
 :
-    db_(db),
-    origin_(vector::zero),
-    direction_(vector(1, 0, 0)),
-    speed_(0),
-    waveModels_(),
-    waveAngles_(),
-    ramp_(),
-    scale_(),
-    crossScale_(),
-    heightAboveWave_(false)
-{}
-
-
-Foam::waveSuperposition::waveSuperposition(const waveSuperposition& waves)
-:
-    db_(waves.db_),
-    origin_(waves.origin_),
-    direction_(waves.direction_),
-    speed_(waves.speed_),
-    waveModels_(waves.waveModels_),
-    waveAngles_(waves.waveAngles_),
-    ramp_(waves.ramp_, false),
-    scale_(waves.scale_, false),
-    crossScale_(waves.crossScale_, false),
-    heightAboveWave_(waves.heightAboveWave_)
-{}
-
-
-Foam::waveSuperposition::waveSuperposition
-(
-    const objectRegistry& db,
-    const dictionary& dict
-)
-:
-    db_(db),
-    origin_(dict.lookup("origin")),
-    direction_(dict.lookup("direction")),
-    speed_(readScalar(dict.lookup("speed"))),
-    waveModels_(),
-    waveAngles_(),
-    ramp_
+    IOdictionary
     (
-        dict.found("ramp")
-      ? Function1<scalar>::New("ramp", dict)
-      : autoPtr<Function1<scalar>>()
+        IOobject
+        (
+            dictName,
+            db.time().constant(),
+            db,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
     ),
+    origin_(lookup("origin")),
+    direction_(lookup("direction")),
+    waveModels_(),
+    waveAngles_(),
     scale_
     (
-        dict.found("scale")
-      ? Function1<scalar>::New("scale", dict)
+        found("scale")
+      ? Function1<scalar>::New("scale", *this)
       : autoPtr<Function1<scalar>>()
     ),
     crossScale_
     (
-        dict.found("crossScale")
-      ? Function1<scalar>::New("crossScale", dict)
+        found("crossScale")
+      ? Function1<scalar>::New("crossScale", *this)
       : autoPtr<Function1<scalar>>()
     ),
-    heightAboveWave_(dict.lookupOrDefault<Switch>("heightAboveWave", false))
+    heightAboveWave_(lookupOrDefault<Switch>("heightAboveWave", false))
 {
-    const PtrList<entry> waveEntries(dict.lookup("waves"));
+    const PtrList<entry> waveEntries(lookup("waves"));
 
     waveModels_.setSize(waveEntries.size());
     waveAngles_.setSize(waveEntries.size());
@@ -267,9 +258,8 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::height
 ) const
 {
     tensor axes;
-    scalar u;
     vectorField xyz(p.size());
-    transformation(p, axes, u, xyz);
+    transformation(p, axes, xyz);
 
     return
         xyz.component(2)
@@ -284,16 +274,15 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::ULiquid
 ) const
 {
     tensor axes;
-    scalar u;
     vectorField xyz(p.size());
-    transformation(p, axes, u, xyz);
+    transformation(p, axes, xyz);
 
     if (heightAboveWave_)
     {
         xyz.replace(2, height(t, p));
     }
 
-    return UMean(t) + (velocity(t, xyz) & axes);
+    return velocity(t, xyz) & axes;
 }
 
 
@@ -304,9 +293,8 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::UGas
 ) const
 {
     tensor axes;
-    scalar u;
     vectorField xyz(p.size());
-    transformation(p, axes, u, xyz);
+    transformation(p, axes, xyz);
 
     axes = tensor(- axes.x(), - axes.y(), axes.z());
 
@@ -317,7 +305,7 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::UGas
 
     xyz.replace(2, - xyz.component(2));
 
-    return UMean(t) + (velocity(t, xyz) & axes);
+    return velocity(t, xyz) & axes;
 }
 
 
@@ -328,9 +316,8 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pLiquid
 ) const
 {
     tensor axes;
-    scalar u;
     vectorField xyz(p.size());
-    transformation(p, axes, u, xyz);
+    transformation(p, axes, xyz);
 
     if (heightAboveWave_)
     {
@@ -348,9 +335,8 @@ Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pGas
 ) const
 {
     tensor axes;
-    scalar u;
     vectorField xyz(p.size());
-    transformation(p, axes, u, xyz);
+    transformation(p, axes, xyz);
 
     axes = tensor(- axes.x(), - axes.y(), axes.z());
 
@@ -369,7 +355,6 @@ void Foam::waveSuperposition::write(Ostream& os) const
 {
     os.writeKeyword("origin") << origin_ << token::END_STATEMENT << nl;
     os.writeKeyword("direction") << direction_ << token::END_STATEMENT << nl;
-    os.writeKeyword("speed") << speed_ << token::END_STATEMENT << nl;
     os.writeKeyword("waves") << nl << token::BEGIN_LIST << nl << incrIndent;
     forAll(waveModels_, wavei)
     {
@@ -380,10 +365,6 @@ void Foam::waveSuperposition::write(Ostream& os) const
             << nl << decrIndent << indent << token::END_BLOCK << nl;
     }
     os  << decrIndent << token::END_LIST << token::END_STATEMENT << nl;
-    if (ramp_.valid())
-    {
-        ramp_->writeData(os);
-    }
     if (scale_.valid())
     {
         scale_->writeData(os);
diff --git a/src/waves/waveSuperposition/waveSuperposition.H b/src/waves/waveSuperposition/waveSuperposition.H
index 0d5be062e0..9b968cac48 100644
--- a/src/waves/waveSuperposition/waveSuperposition.H
+++ b/src/waves/waveSuperposition/waveSuperposition.H
@@ -25,8 +25,48 @@ Class
     Foam::waveSuperposition
 
 Description
-    A wrapper around a list of wave models. Superimposes the modelled values
-    of elevation and velocity.
+    A wrapper around a list of wave models. Superimposes the modelled values of
+    elevation and velocity. The New method looks up or or constructs an
+    instance of this class on demand and returns a reference. Properties are
+    read from a dictionary in constant.
+
+Usage
+    \table
+        Property   | Description                             | Req'd? | Default
+        origin     | origin of the wave coordinate system    | yes    |
+        direction  | direction of the wave coordinate system | yes    |
+        waves      | list of wave models to superimpose      | yes    |
+        scale      | scale factor in the direction           | no     | None
+        crossScale | scale factor perpendicular to the direction | no | None
+        heightAboveWave | use the height above the wave as the vertical \\
+        coordinate | no | false
+    \endtable
+
+    Example specification:
+    \verbatim
+    origin      (0 25 0);
+    direction   (1 0 0);
+    waves
+    (
+        Airy
+        {
+            length      40;
+            amplitude   0.5;
+            phase       0;
+            angle       0;
+        }
+        Airy
+        {
+            length      20;
+            amplitude   0.25;
+            phase       1.5708;
+            angle       0;
+        }
+    );
+    scale       table ((100 1) (200 0));
+    crossScale  constant 1;
+    heightAboveWave no;
+    \endverbatim
 
 SourceFiles
     waveSuperposition.C
@@ -37,6 +77,7 @@ SourceFiles
 #define waveSuperposition_H
 
 #include "waveModel.H"
+#include "IOdictionary.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,34 +89,27 @@ namespace Foam
 \*---------------------------------------------------------------------------*/
 
 class waveSuperposition
+:
+    public IOdictionary
 {
-    // Private data
-
-        //- Reference to the database
-        const objectRegistry& db_;
+    // Private Data
 
         //- The origin of the wave coordinate system
         const vector origin_;
 
-        //- The mean flow direction
+        //- The direction of the wave coordinate system
         const vector direction_;
 
-        //- The mean flow speed
-        const scalar speed_;
-
         //- Wave models to superimpose
         PtrList<waveModel> waveModels_;
 
-        //- The angle relative to the mean velocity at which the waves propagate
+        //- The angle relative to the direction at which the waves propagate
         scalarList waveAngles_;
 
-        //- Ramp for the mean flow speed
-        const autoPtr<Function1<scalar>> ramp_;
-
-        //- Scaling in the flow direction
+        //- Scaling in the local x-direction
         const autoPtr<Function1<scalar>> scale_;
 
-        //- Scaling perpendicular to the flow direction
+        //- Scaling perpendicular to the local x-direction
         const autoPtr<Function1<scalar>> crossScale_;
 
         //- Calculate wave properties using the height above the wave (true) or
@@ -90,27 +124,26 @@ class waveSuperposition
         (
             const vectorField& p,
             tensor& axes,
-            scalar& u,
             vectorField& xyz
         ) const;
 
-        //- Get the wave elevation relative to the mean at a given time, mean
-        //  velocity and local coordinates. Local x is aligned with the mean
-        //  velocity, and y is perpendicular to both x and gravity.
+        //- Get the wave elevation relative to the mean at a given time and
+        //  local coordinates. Local x is aligned with the direction, and y is
+        //  perpendicular to both x and gravity.
         tmp<scalarField> elevation
         (
             const scalar t,
             const vector2DField& xy
         ) const;
 
-        //- Get the wave velocity at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, z with
-        //  negative gravity, and y is perpendicular to both.
+        //- Get the wave velocity at a given time and local coordinates. Local
+        //  x is aligned with the direction, z with negative gravity, and y is
+        //  perpendicular to both.
         tmp<vectorField> velocity(const scalar t, const vectorField& xyz) const;
 
-        //- Get the wave pressure at a given time, mean velocity and local
-        //  coordinates. Local x is aligned with the mean velocity, z with
-        //  negative gravity, and y is perpendicular to both.
+        //- Get the wave pressure at a given time and local coordinates. Local
+        //  x is aligned with the direction, z with negative gravity, and y is
+        //  perpendicular to both.
         tmp<scalarField> pressure(const scalar t, const vectorField& xyz) const;
 
         //- Get the scaling factor, calculated from the optional scaling
@@ -120,6 +153,19 @@ class waveSuperposition
 
 public:
 
+    // Static Data
+
+        //- The name of the dictionary
+        static const word dictName;
+
+
+    // Static Member Functions
+
+        //- Return a reference to the wave model on the given database,
+        //  constructing if it doesn't exist
+        static const waveSuperposition& New(const objectRegistry& db);
+
+
     // Constructors
 
         //- Construct from a database
@@ -153,11 +199,8 @@ public:
         //- Get the gas pressure at a given time and global positions
         tmp<scalarField> pGas(const scalar t, const vectorField& p) const;
 
-        //- Get the mean flow velocity
-        inline vector UMean(const scalar t) const
-        {
-            return (ramp_.valid() ? ramp_->value(t) : 1)*direction_*speed_;
-        }
+        //- Inherit write from regIOobject
+        using regIOobject::write;
 
         //- Write
         void write(Ostream&) const;
diff --git a/tutorials/multiphase/interFoam/RAS/DTCHullWave/0/U.orig b/tutorials/multiphase/interFoam/RAS/DTCHullWave/0/U.orig
index 9ed04308ef..f0e3d4678b 100644
--- a/tutorials/multiphase/interFoam/RAS/DTCHullWave/0/U.orig
+++ b/tutorials/multiphase/interFoam/RAS/DTCHullWave/0/U.orig
@@ -29,20 +29,7 @@ boundaryField
     inlet
     {
         type            waveVelocity;
-        origin          (0 0 0.244);
-        direction       (-1 0 0);
-        speed           $UMean;
-        waves
-        (
-            Stokes2
-            {
-                length      3;
-                amplitude   0.04;
-                phase       0;
-                angle       0;
-            }
-        );
-        scale           table ((4 1) (12 0));
+        UMean           (#neg $UMean 0 0);
     }
 
     outlet
diff --git a/tutorials/multiphase/interFoam/RAS/DTCHullWave/constant/waveProperties b/tutorials/multiphase/interFoam/RAS/DTCHullWave/constant/waveProperties
new file mode 100644
index 0000000000..38c13b9aae
--- /dev/null
+++ b/tutorials/multiphase/interFoam/RAS/DTCHullWave/constant/waveProperties
@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      waveProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+origin          (0 0 0.244);
+
+direction       (-1 0 0);
+
+waves
+(
+    Stokes2
+    {
+        length      3;
+        amplitude   0.04;
+        phase       0;
+        angle       0;
+    }
+);
+
+scale           table ((4 1) (12 0));
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/RAS/DTCHullWave/system/setWavesDict b/tutorials/multiphase/interFoam/RAS/DTCHullWave/system/setWavesDict
new file mode 100644
index 0000000000..362f80e95f
--- /dev/null
+++ b/tutorials/multiphase/interFoam/RAS/DTCHullWave/system/setWavesDict
@@ -0,0 +1,27 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setWavesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+U.orig
+{
+    #include "$FOAM_CASE/0/U.orig";
+}
+
+alpha   alpha.water;
+
+UMean   $U.orig.boundaryField.inlet.UMean;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/laminar/wave/0/U.orig b/tutorials/multiphase/interFoam/laminar/wave/0/U.orig
index fbe501ec7c..a5ca3c0971 100644
--- a/tutorials/multiphase/interFoam/laminar/wave/0/U.orig
+++ b/tutorials/multiphase/interFoam/laminar/wave/0/U.orig
@@ -17,7 +17,7 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0 0 0);
+internalField   uniform (2 0 0);
 
 boundaryField
 {
@@ -26,21 +26,7 @@ boundaryField
     left
     {
         type            waveVelocity;
-        origin          (0 0 0);
-        direction       (1 0 0);
-        speed           2;
-        waves
-        (
-            Airy
-            {
-                length      300;
-                amplitude   2.5;
-                phase       0;
-                angle       0;
-            }
-        );
-        scale           table ((1200 1) (1800 0));
-        crossScale      constant 1;
+        UMean           (2 0 0);
     }
     right
     {
diff --git a/tutorials/multiphase/interFoam/laminar/wave/constant/waveProperties b/tutorials/multiphase/interFoam/laminar/wave/constant/waveProperties
new file mode 100644
index 0000000000..9daba32311
--- /dev/null
+++ b/tutorials/multiphase/interFoam/laminar/wave/constant/waveProperties
@@ -0,0 +1,38 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      waveProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+origin          (0 0 0);
+
+direction       (1 0 0);
+
+waves
+(
+    Airy
+    {
+        length      300;
+        amplitude   2.5;
+        phase       0;
+        angle       0;
+    }
+);
+
+scale           table ((1200 1) (1800 0));
+
+crossScale      constant 1;
+
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/interFoam/laminar/wave/system/setWavesDict b/tutorials/multiphase/interFoam/laminar/wave/system/setWavesDict
new file mode 100644
index 0000000000..8d9d7a085b
--- /dev/null
+++ b/tutorials/multiphase/interFoam/laminar/wave/system/setWavesDict
@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  dev
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      setWavesDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+alpha   alpha.water;
+
+UMean   (2 0 0);
+
+
+// ************************************************************************* //