From c0da2a5effadca4b4409a888ecf4b38ab798ded9 Mon Sep 17 00:00:00 2001
From: Will Bainbridge <http://cfd.direct>
Date: Fri, 13 Oct 2023 12:37:01 +0100
Subject: [PATCH] waveModels: New irregular wave model

This model builds up phase fraction and velocity fields from multiple
first-order waves, sampled from a selectable wave spectrum.

Usage:

    Property  | Description                   | Required? | Default
    ----------+-------------------------------+-----------+-------------
    depth     | The water depth [m]           | no        | great
    spectrum  | The wave spectrum             | yes       |
    n         | The number of times to sample | yes       |
              | the spectrum                  |           |
    span      | The fractional range across   | no        | (0.01 0.99)
              | which to sample the spectrum  |           |
    setFormat | The format with which to plot | no        | none
              | the spectrum                  |           |

Example specification in constant/waveProperties:

    waves
    (
        irregular
        {
            spectrum    PiersonMoskowitz; // or JONSWAP

            PiersonMoskowitzCoeffs
            {
                U19_5       15;
            }

            JONSWAPCoeffs
            {
                U10         10;
                F           200e3;
            }

            n           16;
            span        (0.01 0.99);
        }
    );
---
 src/waves/Make/files                          |   5 +
 src/waves/Make/options                        |   6 +-
 .../waveAlpha/waveAlphaFvPatchScalarField.H   |   6 +-
 .../waveInletOutletFvPatchField.H             |   8 +-
 .../waveVelocityFvPatchVectorField.H          |   4 +-
 src/waves/waveModels/Airy/Airy.C              | 115 ++-----
 src/waves/waveModels/Airy/Airy.H              | 120 +++----
 src/waves/waveModels/Airy/AiryCoeffs.C        | 182 +++++++++++
 src/waves/waveModels/Airy/AiryCoeffs.H        | 172 ++++++++++
 src/waves/waveModels/Airy/AiryI.H             |  77 +++++
 src/waves/waveModels/Stokes2/Stokes2.C        |  54 ++--
 src/waves/waveModels/Stokes2/Stokes2.H        |  24 +-
 src/waves/waveModels/Stokes5/Stokes5.C        |  65 ++--
 src/waves/waveModels/Stokes5/Stokes5.H        |  22 +-
 src/waves/waveModels/irregular/irregular.C    | 301 ++++++++++++++++++
 src/waves/waveModels/irregular/irregular.H    | 199 ++++++++++++
 .../irregular/waveSpectra/JONSWAP/JONSWAP.C   | 109 +++++++
 .../irregular/waveSpectra/JONSWAP/JONSWAP.H   | 152 +++++++++
 .../PiersonMoskowitz/PiersonMoskowitz.C       | 122 +++++++
 .../PiersonMoskowitz/PiersonMoskowitz.H       | 166 ++++++++++
 .../waveSpectra/waveSpectrum/waveSpectrum.C   | 211 ++++++++++++
 .../waveSpectra/waveSpectrum/waveSpectrum.H   | 172 ++++++++++
 src/waves/waveModels/solitary/solitary.H      |  21 ++
 src/waves/waveModels/waveModel/waveModel.C    |   8 +-
 .../waveAtmBoundaryLayerSuperposition.H       |  18 +-
 .../waveSuperposition/waveSuperposition.C     |   2 +-
 .../waveSuperposition/waveSuperposition.H     |  30 +-
 27 files changed, 2094 insertions(+), 277 deletions(-)
 create mode 100644 src/waves/waveModels/Airy/AiryCoeffs.C
 create mode 100644 src/waves/waveModels/Airy/AiryCoeffs.H
 create mode 100644 src/waves/waveModels/Airy/AiryI.H
 create mode 100644 src/waves/waveModels/irregular/irregular.C
 create mode 100644 src/waves/waveModels/irregular/irregular.H
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.C
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.H
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.C
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.H
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.C
 create mode 100644 src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.H

diff --git a/src/waves/Make/files b/src/waves/Make/files
index b999601919..9fd9b25367 100644
--- a/src/waves/Make/files
+++ b/src/waves/Make/files
@@ -1,9 +1,14 @@
 waveModels/waveModel/waveModel.C
 waveModels/waveModel/waveModelNew.C
+waveModels/Airy/AiryCoeffs.C
 waveModels/Airy/Airy.C
 waveModels/Stokes2/Stokes2.C
 waveModels/Stokes5/Stokes5.C
 waveModels/solitary/solitary.C
+waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.C
+waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.C
+waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.C
+waveModels/irregular/irregular.C
 
 waveSuperpositions/waveSuperposition/waveSuperposition.C
 waveSuperpositions/waveSuperposition/waveSuperpositionNew.C
diff --git a/src/waves/Make/options b/src/waves/Make/options
index 4709d60eae..c39852a3bd 100644
--- a/src/waves/Make/options
+++ b/src/waves/Make/options
@@ -2,10 +2,12 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/atmosphericModels/lnInclude
+    -I$(LIB_SRC)/atmosphericModels/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -ldynamicMesh \
-    -latmosphericModels
+    -latmosphericModels \
+    -lsampling
diff --git a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
index b2daff1949..a650cef434 100644
--- a/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
+++ b/src/waves/derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.H
@@ -32,9 +32,9 @@ Description
 
 Usage
     \table
-        Property  | Description                            | Req'd? | Default
-        phi       | Name of the flux field                 | no     | phi
-        liquid    | Is the alpha field that of the liquid? | no     | true
+        Property  | Description                            | Required? | Default
+        phi       | Name of the flux field                 | no        | phi
+        liquid    | Is the alpha field that of the liquid? | no        | true
     \endtable
 
     Example of the boundary condition specification:
diff --git a/src/waves/derivedFvPatchFields/waveInletOutlet/waveInletOutletFvPatchField.H b/src/waves/derivedFvPatchFields/waveInletOutlet/waveInletOutletFvPatchField.H
index c80421f9a9..8880da2f1b 100644
--- a/src/waves/derivedFvPatchFields/waveInletOutlet/waveInletOutletFvPatchField.H
+++ b/src/waves/derivedFvPatchFields/waveInletOutlet/waveInletOutletFvPatchField.H
@@ -32,10 +32,10 @@ Description
 
 Usage
     \table
-        Property        | Description                | Req'd? | Default
-        phi             | Name of the flux field     | no     | phi
-        inletValueAbove | inlet value above the wave | no     | None
-        inletValueBelow | inlet value below the wave | no     | None
+        Property        | Description                | Required? | Default
+        phi             | Name of the flux field     | no        | phi
+        inletValueAbove | inlet value above the wave | no        | None
+        inletValueBelow | inlet value below the wave | no        | None
     \endtable
 
     Example of the boundary condition specification:
diff --git a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
index 500d073221..60d4b1ebf1 100644
--- a/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
+++ b/src/waves/derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.H
@@ -32,8 +32,8 @@ Description
 
 Usage
     \table
-        Property  | Description                            | Req'd? | Default
-        phi       | Name of the flux field                 | no     | phi
+        Property  | Description                            | Required? | Default
+        phi       | Name of the flux field                 | no        | phi
     \endtable
 
     Example of the boundary condition specification:
diff --git a/src/waves/waveModels/Airy/Airy.C b/src/waves/waveModels/Airy/Airy.C
index c6f58ecef7..830c420d5b 100644
--- a/src/waves/waveModels/Airy/Airy.C
+++ b/src/waves/waveModels/Airy/Airy.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,16 +38,15 @@ namespace waveModels
 }
 }
 
-
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::scalar Foam::waveModels::Airy::length
+Foam::scalar Foam::waveModels::Airy::readLength
 (
     const dictionary& dict,
     const scalar depth,
     const scalar amplitude,
     const scalar g,
-    scalar (*modelCelerity)(scalar, scalar, scalar, scalar)
+    scalar (*celerityPtr)(const AiryCoeffs&)
 )
 {
     const bool haveLength = dict.found("length");
@@ -66,90 +65,24 @@ Foam::scalar Foam::waveModels::Airy::length
     }
     else
     {
-        const scalar period = dict.lookup<scalar>("period");
-
-        scalar length0 = 0;
-        scalar length1 = 1.5*g*sqr(period)/(2*constant::mathematical::pi);
-
-        // Bisect down to round off error to find the wave length
-        for (label i = 0; i < ceil(std::log2(1/small)); ++ i)
-        {
-            const scalar length = (length0 + length1)/2;
-
-            const scalar t = length/modelCelerity(depth, amplitude, length, g);
-
-            (t < period ? length0 : length1) = length;
-        }
-
-        return (length0 + length1)/2;
+        return
+            AiryCoeffs
+            (
+                depth,
+                amplitude,
+                dict.lookup<scalar>("period"),
+                g,
+                celerityPtr
+            ).length;
     }
 }
 
 
-// * * * * * * * * * * Static Protected Member Functions  * * * * * * ** * * //
-
-Foam::scalar Foam::waveModels::Airy::k(const scalar length)
-{
-    return 2*Foam::constant::mathematical::pi/length;
-}
-
-
-bool Foam::waveModels::Airy::deep(const scalar depth, const scalar length)
-{
-    return depth*k(length) > log(great);
-}
-
-
-Foam::scalar Foam::waveModels::Airy::celerity
-(
-    const scalar depth,
-    const scalar amplitude,
-    const scalar length,
-    const scalar g
-)
-{
-    return sqrt(g/k(length)*tanh(k(length)*depth));
-}
-
-
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::angle
-(
-    const scalar t,
-    const scalarField& x
-) const
+Foam::scalar Foam::waveModels::Airy::celerity(const AiryCoeffs& coeffs)
 {
-    return phase_ + k()*(x - celerity()*t);
-}
-
-
-Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::vi
-(
-    const label i,
-    const scalar t,
-    const vector2DField& xz
-) const
-{
-    const scalarField x(xz.component(0));
-    const scalarField z(xz.component(1));
-
-    const scalarField phi(angle(t, x));
-
-    const scalar kzGreat = log(i*great);
-    const scalarField kz(min(max(k()*z, - kzGreat), kzGreat));
-
-    if (deep())
-    {
-        return i*exp(kz)*zip(cos(i*phi), sin(i*phi));
-    }
-    else
-    {
-        const scalar kd = k()*depth();
-        const scalarField kdz(max(scalar(0), kd + kz));
-
-        return i*zip(cosh(i*kdz)*cos(i*phi), sinh(i*kdz)*sin(i*phi))/sinh(kd);
-    }
+    return coeffs.celerity();
 }
 
 
@@ -170,16 +103,16 @@ Foam::waveModels::Airy::Airy
     const dictionary& dict,
     const scalar g,
     const word& modelName,
-    scalar (*modelCelerity)(scalar, scalar, scalar, scalar)
+    scalar (*celerityPtr)(const AiryCoeffs&)
 )
 :
     waveModel(dict, g),
     depth_(dict.lookupOrDefault<scalar>("depth", great)),
     amplitude_(Function1<scalar>::New("amplitude", dict)),
-    length_(length(dict, depth_, amplitude(), g, modelCelerity)),
+    length_(readLength(dict, depth_, amplitude(), g, celerityPtr)),
     phase_(dict.lookup<scalar>("phase"))
 {
-    const scalar c = modelCelerity(depth(), amplitude(), length(), g);
+    const scalar c = celerityPtr(coeffs());
 
     Info<< waveModel::typeName << ": " << modelName
         << ": period = " << length_/c
@@ -195,13 +128,19 @@ Foam::waveModels::Airy::~Airy()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::waveModels::Airy::celerity() const
+{
+    return celerity(coeffs());
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::elevation
 (
     const scalar t,
     const scalarField& x
 ) const
 {
-    return amplitude(t)*cos(angle(t, x));
+    return amplitude(t)*cos(coeffs().angle(phase_, t, x));
 }
 
 
@@ -211,9 +150,9 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::velocity
     const vector2DField& xz
 ) const
 {
-    const scalar ka = k()*amplitude(t);
+    const scalar ka = coeffs().k()*amplitude(t);
 
-    return Airy::celerity()*ka*vi(1, t, xz);
+    return Airy::celerity()*ka*coeffs().vi(1, phase_, t, xz);
 }
 
 
@@ -221,7 +160,7 @@ void Foam::waveModels::Airy::write(Ostream& os) const
 {
     waveModel::write(os);
 
-    if (!deep())
+    if (!coeffs().deep())
     {
         writeEntry(os, "depth", depth_);
     }
diff --git a/src/waves/waveModels/Airy/Airy.H b/src/waves/waveModels/Airy/Airy.H
index 322c15e239..0d3249003a 100644
--- a/src/waves/waveModels/Airy/Airy.H
+++ b/src/waves/waveModels/Airy/Airy.H
@@ -36,6 +36,29 @@ Description
 
     See the leading terms of equations 18 and 19.
 
+Usage
+    \table
+        Property  | Description            | Required?         | Default
+        depth     | The water depth [m]    | no                | great
+        amplitude | The amplitude [m]      | yes               |
+        length    | The wave length [m]    | if period not set |
+        period    | The wave period [s]    | if length not set |
+        phase     | The phase offset [rad] | yes               |
+    \endtable
+
+    Example specification in constant/waveProperties:
+    \verbatim
+    waves
+    (
+        Airy
+        {
+            length      40;
+            amplitude   0.5;
+            phase       0;
+        }
+    );
+    \endverbatim
+
 SourceFiles
     Airy.C
 
@@ -45,6 +68,7 @@ SourceFiles
 #define Airy_H
 
 #include "waveModel.H"
+#include "AiryCoeffs.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,81 +90,43 @@ class Airy
         //- Depth [m]
         const scalar depth_;
 
-        //- Peak-to-peak mean amplitude [m]
+        //- Amplitude [m]
         autoPtr<Function1<scalar>> amplitude_;
 
-        //- Peak-to-peak length [m]
+        //- Wavelength [m]
         const scalar length_;
 
         //- Phase offset [rad]
         const scalar phase_;
 
 
-private:
-
     // Private Member Functions
 
         //- Read and return the wave length from the dictionary. Either reads
         //  the length directly, or reads the period and depth and calculates
         //  the length.
-        static scalar length
+        static scalar readLength
         (
             const dictionary& dict,
             const scalar depth,
             const scalar amplitude,
             const scalar g,
-            scalar (*celerityPtr)(scalar, scalar, scalar, scalar)
+            scalar (*celerityPtr)(const AiryCoeffs&)
         );
 
 
 protected:
 
-    // Static Protected Member Functions
-
-        //- The angular wavenumber [rad/m]
-        static scalar k(const scalar length);
-
-        //- Return whether shallow and intermediate effects are to be omitted
-        static bool deep(const scalar length, const scalar depth);
-
-        //- The wave celerity [m/s]
-        static scalar celerity
-        (
-            const scalar depth,
-            const scalar amplitude,
-            const scalar length,
-            const scalar g
-        );
-
-
     // Protected Member Functions
 
-        //- The angular wavenumber [rad/m]
-        scalar k() const
-        {
-            return k(length_);
-        }
+        //- Return the wave coefficients
+        AiryCoeffs coeffs(const scalar t) const;
 
-        //- Return whether shallow and intermediate effects are to be omitted
-        bool deep() const
-        {
-            return deep(depth(), length());
-        }
+        //- Return the wave coefficients at steady state
+        AiryCoeffs coeffs() const;
 
-        //- Angle of the oscillation [rad]
-        tmp<scalarField> angle
-        (
-            const scalar t,
-            const scalarField& x
-        ) const;
-
-        //- Return the non-dimensionalised i-th harmonic of the velocity
-        tmp<vector2DField> vi
-        (
-            const label i,
-            const scalar t,
-            const vector2DField& xz
-        ) const;
+        //- The wave celerity [m/s]
+        static scalar celerity(const AiryCoeffs&);
 
 
 public:
@@ -160,8 +146,7 @@ public:
             const dictionary& dict,
             const scalar g,
             const word& modelName = Airy::typeName,
-            scalar (*modelCelerity)(scalar, scalar, scalar, scalar) =
-                &Airy::celerity
+            scalar (*celerityPtr)(const AiryCoeffs&) = &AiryCoeffs::celerity
         );
 
         //- Construct a clone
@@ -180,40 +165,23 @@ public:
         // Access
 
             //- Get the depth
-            scalar depth() const
-            {
-                return depth_;
-            }
+            scalar depth() const;
 
             //- Get the amplitude
-            scalar amplitude(const scalar t) const
-            {
-                return amplitude_->value(t);
-            }
+            scalar amplitude(const scalar t) const;
 
             //- Get the amplitude at steady state
-            scalar amplitude() const
-            {
-                return amplitude_->value(great);
-            }
+            scalar amplitude() const;
 
             //- Get the length
-            scalar length() const
-            {
-                return length_;
-            }
+            scalar length() const;
 
             //- Get the phase
-            scalar phase() const
-            {
-                return phase_;
-            }
+            scalar phase() const;
+
 
         //- The wave celerity [m/s]
-        virtual scalar celerity() const
-        {
-            return celerity(depth(), amplitude(), length(), g());
-        }
+        virtual scalar celerity() const;
 
         //- Get the wave elevation at a given time and local coordinates. Local
         //  x is aligned with the direction of propagation.
@@ -244,6 +212,10 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "AiryI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 #endif
 
 // ************************************************************************* //
diff --git a/src/waves/waveModels/Airy/AiryCoeffs.C b/src/waves/waveModels/Airy/AiryCoeffs.C
new file mode 100644
index 0000000000..7630990739
--- /dev/null
+++ b/src/waves/waveModels/Airy/AiryCoeffs.C
@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2023 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 "AiryCoeffs.H"
+#include "mathematicalConstants.H"
+
+// * * * * * * * * * * * Private Static Member Functions * * * * * * * * * * //
+
+Foam::scalar Foam::waveModels::AiryCoeffs::calcLength
+(
+    const scalar depth,
+    const scalar amplitude,
+    const scalar period,
+    const scalar g,
+    scalar (*celerityPtr)(const AiryCoeffs&)
+)
+{
+    scalar length0 = 0;
+    scalar length1 = 1.5*g*sqr(period)/(2*constant::mathematical::pi);
+
+    // Bisect down to round off error to find the wave length
+    for (label i = 0; i < ceil(std::log2(1/small)); ++ i)
+    {
+        const scalar length = (length0 + length1)/2;
+
+        const scalar t =
+            length/celerityPtr(AiryCoeffs(depth, amplitude, length, g));
+
+        (t < period ? length0 : length1) = length;
+    }
+
+    return (length0 + length1)/2;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveModels::AiryCoeffs::AiryCoeffs
+(
+    const scalar depth,
+    const scalar amplitude,
+    const scalar length,
+    const scalar g
+)
+:
+    depth(depth),
+    amplitude(amplitude),
+    length(length),
+    g(g)
+{}
+
+
+Foam::waveModels::AiryCoeffs::AiryCoeffs
+(
+    const scalar depth,
+    const scalar amplitude,
+    const scalar period,
+    const scalar g,
+    scalar (*celerityPtr)(const AiryCoeffs&)
+)
+:
+    depth(depth),
+    amplitude(amplitude),
+    length(calcLength(depth, amplitude, period, g, celerityPtr)),
+    g(g)
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::waveModels::AiryCoeffs::k() const
+{
+    return 2*Foam::constant::mathematical::pi/length;
+}
+
+
+bool Foam::waveModels::AiryCoeffs::deep() const
+{
+    return depth*k() > log(great);
+}
+
+
+Foam::scalar Foam::waveModels::AiryCoeffs::celerity(const AiryCoeffs& coeffs)
+{
+    return sqrt(coeffs.g/coeffs.k()*tanh(coeffs.k()*coeffs.depth));
+}
+
+
+Foam::scalar Foam::waveModels::AiryCoeffs::celerity() const
+{
+    return celerity(*this);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::AiryCoeffs::angle
+(
+    const scalar phase,
+    const scalar t,
+    const scalarField& x
+) const
+{
+    return phase + k()*(x - celerity()*t);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::AiryCoeffs::elevation
+(
+    const scalar phase,
+    const scalar t,
+    const scalarField& x
+) const
+{
+    return amplitude*cos(angle(phase, t, x));
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::AiryCoeffs::vi
+(
+    const label i,
+    const scalar phase,
+    const scalar t,
+    const vector2DField& xz
+) const
+{
+    const scalarField x(xz.component(0));
+    const scalarField z(xz.component(1));
+
+    const scalarField phi(angle(phase, t, x));
+
+    const scalar kzGreat = log(i*great);
+    const scalarField kz(min(max(k()*z, - kzGreat), kzGreat));
+
+    if (deep())
+    {
+        return i*exp(kz)*zip(cos(i*phi), sin(i*phi));
+    }
+    else
+    {
+        const scalar kd = k()*depth;
+        const scalarField kdz(max(scalar(0), kd + kz));
+
+        return i*zip(cosh(i*kdz)*cos(i*phi), sinh(i*kdz)*sin(i*phi))/sinh(kd);
+    }
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::AiryCoeffs::velocity
+(
+    const scalar phase,
+    const scalar t,
+    const vector2DField& xz
+) const
+{
+    const scalar ka = k()*amplitude;
+
+    return celerity()*ka*vi(1, phase, t, xz);
+}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/Airy/AiryCoeffs.H b/src/waves/waveModels/Airy/AiryCoeffs.H
new file mode 100644
index 0000000000..dafdef6e5c
--- /dev/null
+++ b/src/waves/waveModels/Airy/AiryCoeffs.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2023 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/>.
+
+Class
+    Foam::waveModels::AiryCoeffs
+
+Description
+    Calculation engine for the Airy wave model and other models that are a
+    correction on top of the Airy model or a superposition of Airy models
+
+SourceFiles
+    AiryCoeffs.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AiryCoeffs_H
+#define AiryCoeffs_H
+
+#include "waveModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class AiryCoeffs Declaration
+\*---------------------------------------------------------------------------*/
+
+class AiryCoeffs
+{
+private:
+
+    // Private Static Member Functions
+
+        //- Calculate the length from the depth, amplitude, period and g.
+        //  Requires a function which returns celerity given the depth,
+        //  amplitude, length and g. Iterates to solve this function for the
+        //  specified period.
+        static scalar calcLength
+        (
+            const scalar depth,
+            const scalar amplitude,
+            const scalar period,
+            const scalar g,
+            scalar (*celerityPtr)(const AiryCoeffs&)
+        );
+
+
+public:
+
+    // Public Data
+
+        //- Depth [m]
+        const scalar depth;
+
+        //- Amplitude [m]
+        const scalar amplitude;
+
+        //- Wavelength [m]
+        const scalar length;
+
+        //- Gravitational acceleration [m/s^2]
+        const scalar g;
+
+
+    // Constructors
+
+        //- Construct from components
+        AiryCoeffs
+        (
+            const scalar depth,
+            const scalar amplitude,
+            const scalar length,
+            const scalar g
+        );
+
+        //- Construct from components but with period instead of length
+        AiryCoeffs
+        (
+            const scalar depth,
+            const scalar amplitude,
+            const scalar period,
+            const scalar g,
+            scalar (*celerityPtr)(const AiryCoeffs&)
+        );
+
+
+    // Member Functions
+
+        //- The angular wavenumber [rad/m]
+        scalar k() const;
+
+        //- Return whether shallow and intermediate effects are to be omitted
+        bool deep() const;
+
+        //- The wave celerity [m/s]
+        static scalar celerity(const AiryCoeffs& coeffs);
+
+        //- The wave celerity [m/s]
+        scalar celerity() const;
+
+        //- Angle of the oscillation [rad]
+        tmp<scalarField> angle
+        (
+            const scalar phase,
+            const scalar t,
+            const scalarField& x
+        ) const;
+
+        //- Get the wave elevation at a given time and local coordinates. Local
+        //  x is aligned with the direction of propagation.
+        tmp<scalarField> elevation
+        (
+            const scalar phase,
+            const scalar t,
+            const scalarField& x
+        ) const;
+
+        //- Return the non-dimensionalised i-th harmonic of the velocity
+        tmp<vector2DField> vi
+        (
+            const label i,
+            const scalar phase,
+            const scalar t,
+            const vector2DField& xz
+        ) const;
+
+        //- 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.
+        tmp<vector2DField> velocity
+        (
+            const scalar phase,
+            const scalar t,
+            const vector2DField& xz
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace waveModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/Airy/AiryI.H b/src/waves/waveModels/Airy/AiryI.H
new file mode 100644
index 0000000000..aaef8d7831
--- /dev/null
+++ b/src/waves/waveModels/Airy/AiryI.H
@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2017-2023 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 "Airy.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+inline Foam::waveModels::AiryCoeffs Foam::waveModels::Airy::coeffs
+(
+    const scalar t
+) const
+{
+    return AiryCoeffs(depth_, amplitude_->value(t), length_, g());
+}
+
+
+inline Foam::waveModels::AiryCoeffs Foam::waveModels::Airy::coeffs() const
+{
+    return AiryCoeffs(depth_, amplitude(), length_, g());
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::waveModels::Airy::depth() const
+{
+    return depth_;
+}
+
+
+inline Foam::scalar Foam::waveModels::Airy::amplitude(const scalar t) const
+{
+    return amplitude_->value(t);
+}
+
+
+inline Foam::scalar Foam::waveModels::Airy::amplitude() const
+{
+    return amplitude_->value(great);
+}
+
+
+inline Foam::scalar Foam::waveModels::Airy::length() const
+{
+    return length_;
+}
+
+
+inline Foam::scalar Foam::waveModels::Airy::phase() const
+{
+    return phase_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/Stokes2/Stokes2.C b/src/waves/waveModels/Stokes2/Stokes2.C
index 76721bfa1f..94407b74cf 100644
--- a/src/waves/waveModels/Stokes2/Stokes2.C
+++ b/src/waves/waveModels/Stokes2/Stokes2.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,21 +40,15 @@ namespace waveModels
 
 // * * * * * * * * * * Static Protected Member Functions  * * * * * * ** * * //
 
-Foam::scalar Foam::waveModels::Stokes2::celerity
-(
-    const scalar depth,
-    const scalar amplitude,
-    const scalar length,
-    const scalar g
-)
+Foam::scalar Foam::waveModels::Stokes2::celerity(const AiryCoeffs& coeffs)
 {
     static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k(length)*depth, - kdGreat), kdGreat);
-    const scalar ka = k(length)*amplitude;
+    const scalar kd = min(max(coeffs.k()*coeffs.depth, - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*coeffs.amplitude;
 
-    const scalar S = deep(depth, length) ? 0 : 1/cosh(2*kd);
+    const scalar S = coeffs.deep() ? 0 : 1/cosh(2*kd);
 
-    const scalar C0 = Airy::celerity(depth, amplitude, length, g);
+    const scalar C0 = coeffs.celerity();
     const scalar C2ByC0 = (2 + 7*sqr(S))/4/sqr(1 - S);
 
     if (debug)
@@ -62,7 +56,7 @@ Foam::scalar Foam::waveModels::Stokes2::celerity
         Info<< "C2 = " << C2ByC0*C0 << endl;
     }
 
-    return Airy::celerity(depth, amplitude, length, g) + sqr(ka)*C2ByC0*C0;
+    return Airy::celerity(coeffs) + sqr(ka)*C2ByC0*C0;
 }
 
 
@@ -73,10 +67,10 @@ Foam::waveModels::Stokes2::Stokes2
     const dictionary& dict,
     const scalar g,
     const word& modelName,
-    scalar (*modelCelerity)(scalar, scalar, scalar, scalar)
+    scalar (*celerityPtr)(const AiryCoeffs&)
 )
 :
-    Airy(dict, g, modelName, modelCelerity)
+    Airy(dict, g, modelName, celerityPtr)
 {}
 
 
@@ -88,17 +82,25 @@ Foam::waveModels::Stokes2::~Stokes2()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::waveModels::Stokes2::celerity() const
+{
+    return celerity(coeffs());
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes2::elevation
 (
     const scalar t,
     const scalarField& x
 ) const
 {
-    static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k()*depth(), - kdGreat), kdGreat);
-    const scalar ka = k()*amplitude(t);
+    const AiryCoeffs coeffs = this->coeffs();
 
-    const scalar T = deep() ? 1 : tanh(kd);
+    static const scalar kdGreat = log(great);
+    const scalar kd = min(max(coeffs.k()*depth(), - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*amplitude(t);
+
+    const scalar T = coeffs.deep() ? 1 : tanh(kd);
 
     const scalar B22 = (3/sqr(T) - 1)/T/4;
 
@@ -109,7 +111,7 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes2::elevation
 
     return
         Airy::elevation(t, x)
-      + (1/k())*sqr(ka)*B22*cos(2*angle(t, x));
+      + (1/coeffs.k())*sqr(ka)*B22*cos(2*coeffs.angle(phase(), t, x));
 }
 
 
@@ -119,11 +121,13 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes2::velocity
     const vector2DField& xz
 ) const
 {
-    static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k()*depth(), - kdGreat), kdGreat);
-    const scalar ka = k()*amplitude(t);
+    const AiryCoeffs coeffs = this->coeffs();
 
-    const scalar A22ByA11 = deep() ? 0 : 0.375/pow3(sinh(kd));
+    static const scalar kdGreat = log(great);
+    const scalar kd = min(max(coeffs.k()*depth(), - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*amplitude(t);
+
+    const scalar A22ByA11 = coeffs.deep() ? 0 : 0.375/pow3(sinh(kd));
 
     if (debug)
     {
@@ -133,7 +137,7 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes2::velocity
 
     return
         Airy::velocity(t, xz)
-      + Airy::celerity()*sqr(ka)*A22ByA11*vi(2, t, xz);
+      + Airy::celerity()*sqr(ka)*A22ByA11*coeffs.vi(2, phase(), t, xz);
 }
 
 
diff --git a/src/waves/waveModels/Stokes2/Stokes2.H b/src/waves/waveModels/Stokes2/Stokes2.H
index 5dc31f7e74..58b5e87517 100644
--- a/src/waves/waveModels/Stokes2/Stokes2.H
+++ b/src/waves/waveModels/Stokes2/Stokes2.H
@@ -36,6 +36,12 @@ Description
 
     See equations 18 and 19.
 
+Usage
+    The parameters of this model are identical to those used by the Airy model
+
+See also
+    Foam::waveModels::Airy
+
 SourceFiles
     Stokes2.C
 
@@ -66,18 +72,12 @@ protected:
     // Protected Member Functions
 
         //- The wave celerity [m/s]
-        static scalar celerity
-        (
-            const scalar depth,
-            const scalar amplitude,
-            const scalar length,
-            const scalar g
-        );
+        static scalar celerity(const AiryCoeffs&);
 
 
 public:
 
-   //- Runtime type information
+    //- Runtime type information
     TypeName("Stokes2");
 
 
@@ -89,8 +89,7 @@ public:
             const dictionary& dict,
             const scalar g,
             const word& modelName = Stokes2::typeName,
-            scalar (*modelCelerity)(scalar, scalar, scalar, scalar) =
-                &Stokes2::celerity
+            scalar (*celerityPtr)(const AiryCoeffs&) = &Stokes2::celerity
         );
 
         //- Construct a clone
@@ -107,10 +106,7 @@ public:
     // Member Functions
 
         //- The wave celerity [m/s]
-        virtual scalar celerity() const
-        {
-            return celerity(depth(), amplitude(), length(), g());
-        }
+        virtual scalar celerity() const;
 
         //- Get the wave elevation at a given time and local coordinates. Local
         //  x is aligned with the direction of propagation.
diff --git a/src/waves/waveModels/Stokes5/Stokes5.C b/src/waves/waveModels/Stokes5/Stokes5.C
index 6471bb085b..853b865a8e 100644
--- a/src/waves/waveModels/Stokes5/Stokes5.C
+++ b/src/waves/waveModels/Stokes5/Stokes5.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,21 +40,15 @@ namespace waveModels
 
 // * * * * * * * * * * Static Protected Member Functions  * * * * * * ** * * //
 
-Foam::scalar Foam::waveModels::Stokes5::celerity
-(
-    const scalar depth,
-    const scalar amplitude,
-    const scalar length,
-    const scalar g
-)
+Foam::scalar Foam::waveModels::Stokes5::celerity(const AiryCoeffs& coeffs)
 {
     static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k(length)*depth, - kdGreat), kdGreat);
-    const scalar ka = k(length)*amplitude;
+    const scalar kd = min(max(coeffs.k()*coeffs.depth, - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*coeffs.amplitude;
 
-    const scalar S = deep(depth, length) ? 0 : 1/cosh(2*kd);
+    const scalar S = coeffs.deep() ? 0 : 1/cosh(2*kd);
 
-    const scalar C0 = Airy::celerity(depth, amplitude, length, g);
+    const scalar C0 = coeffs.celerity();
     const scalar C4ByC0 =
         1.0/32/pow5(1 - S)
        *(4 + 32*S - 116*sqr(S) - 400*pow3(S) - 71*pow4(S) + 146*pow5(S));
@@ -64,7 +58,7 @@ Foam::scalar Foam::waveModels::Stokes5::celerity
         Info<< "C4 = " << C4ByC0*C0 << endl;
     }
 
-    return Stokes2::celerity(depth, amplitude, length, g) + pow4(ka)*C4ByC0*C0;
+    return Stokes2::celerity(coeffs) + pow4(ka)*C4ByC0*C0;
 }
 
 
@@ -75,10 +69,10 @@ Foam::waveModels::Stokes5::Stokes5
     const dictionary& dict,
     const scalar g,
     const word& modelName,
-    scalar (*modelCelerity)(scalar, scalar, scalar, scalar)
+    scalar (*celerityPtr)(const AiryCoeffs&)
 )
 :
-    Stokes2(dict, g, modelName, modelCelerity)
+    Stokes2(dict, g, modelName, celerityPtr)
 {}
 
 
@@ -90,17 +84,26 @@ Foam::waveModels::Stokes5::~Stokes5()
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
+Foam::scalar Foam::waveModels::Stokes5::celerity() const
+{
+    return celerity(coeffs());
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
 (
     const scalar t,
     const scalarField& x
 ) const
 {
-    static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k()*depth(), - kdGreat), kdGreat);
-    const scalar ka = k()*amplitude(t);
+    const AiryCoeffs coeffs = this->coeffs();
 
-    const scalar S = deep() ? 0 : 1/cosh(2*kd), T = deep() ? 1 : tanh(kd);
+    static const scalar kdGreat = log(great);
+    const scalar kd = min(max(coeffs.k()*depth(), - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*amplitude(t);
+
+    const scalar S = coeffs.deep() ? 0 : 1/cosh(2*kd);
+    const scalar T = coeffs.deep() ? 1 : tanh(kd);
 
     const scalar B31 =
       - 3.0/8/pow3(1 - S)
@@ -143,11 +146,11 @@ Foam::tmp<Foam::scalarField> Foam::waveModels::Stokes5::elevation
             << "B55 = " << B55 << endl;
     }
 
-    const scalarField phi(angle(t, x));
+    const scalarField phi(coeffs.angle(phase(), t, x));
 
     return
         Stokes2::elevation(t, x)
-      + (1/k())
+      + (1/coeffs.k())
        *(
            pow3(ka)*B31*(cos(phi) - cos(3*phi))
          + pow4(ka)*(B42*cos(2*phi) + B44*cos(4*phi))
@@ -162,12 +165,14 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes5::velocity
     const vector2DField& xz
 ) const
 {
-    static const scalar kdGreat = log(great);
-    const scalar kd = min(max(k()*depth(), - kdGreat), kdGreat);
-    const scalar ka = k()*amplitude(t);
+    const AiryCoeffs coeffs = this->coeffs();
 
-    const scalar S = deep() ? 0 : 1/cosh(2*kd);
-    const scalar SByA11 = deep() ? 0 : S*sinh(kd);
+    static const scalar kdGreat = log(great);
+    const scalar kd = min(max(coeffs.k()*depth(), - kdGreat), kdGreat);
+    const scalar ka = coeffs.k()*amplitude(t);
+
+    const scalar S = coeffs.deep() ? 0 : 1/cosh(2*kd);
+    const scalar SByA11 = coeffs.deep() ? 0 : S*sinh(kd);
 
     const scalar A31ByA11 =
         1.0/8/pow3(1 - S)
@@ -226,15 +231,17 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Stokes5::velocity
             << "A55 = " << A55ByA11*A11 << endl;
     }
 
-    const vector2DField v1(vi(1, t, xz)), v3(vi(3, t, xz));
+    auto v = [&](const label i) { return coeffs.vi(i, phase(), t, xz); };
+
+    const vector2DField v1(v(1)), v3(v(3));
 
     return
         Stokes2::velocity(t, xz)
       + Airy::celerity()
        *(
             pow3(ka)*(A31ByA11*v1 + A33ByA11*v3)
-          + pow4(ka)*(A42ByA11*vi(2, t, xz) + A44ByA11*vi(4, t, xz))
-          + pow5(ka)*(A51ByA11*v1 + A53ByA11*v3 + A55ByA11*vi(5, t, xz))
+          + pow4(ka)*(A42ByA11*v(2) + A44ByA11*v(4))
+          + pow5(ka)*(A51ByA11*v1 + A53ByA11*v3 + A55ByA11*v(5))
         );
 }
 
diff --git a/src/waves/waveModels/Stokes5/Stokes5.H b/src/waves/waveModels/Stokes5/Stokes5.H
index d173b40881..125621d5d3 100644
--- a/src/waves/waveModels/Stokes5/Stokes5.H
+++ b/src/waves/waveModels/Stokes5/Stokes5.H
@@ -35,6 +35,12 @@ Description
         216-234.
     \endverbatim
 
+Usage
+    The parameters of this model are identical to those used by the Airy model
+
+See also:
+    Foam::waveModels::Airy
+
 SourceFiles
     Stokes5.C
 
@@ -65,13 +71,7 @@ protected:
     // Protected Member Functions
 
         //- The wave celerity [m/s]
-        static scalar celerity
-        (
-            const scalar depth,
-            const scalar amplitude,
-            const scalar length,
-            const scalar g
-        );
+        static scalar celerity(const AiryCoeffs&);
 
 
 public:
@@ -88,8 +88,7 @@ public:
             const dictionary& dict,
             const scalar g,
             const word& modelName = Stokes5::typeName,
-            scalar (*modelCelerity)(scalar, scalar, scalar, scalar) =
-                &Stokes5::celerity
+            scalar (*celerityPtr)(const AiryCoeffs&) = &Stokes5::celerity
         );
 
         //- Construct a clone
@@ -106,10 +105,7 @@ public:
     // Member Functions
 
         //- The wave celerity [m/s]
-        virtual scalar celerity() const
-        {
-            return celerity(depth(), amplitude(), length(), g());
-        }
+        virtual scalar celerity() const;
 
         //- Get the wave elevation at a given time and local coordinates. Local
         //  x is aligned with the direction of propagation.
diff --git a/src/waves/waveModels/irregular/irregular.C b/src/waves/waveModels/irregular/irregular.C
new file mode 100644
index 0000000000..ee0a9fe6dc
--- /dev/null
+++ b/src/waves/waveModels/irregular/irregular.C
@@ -0,0 +1,301 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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 "irregular.H"
+#include "mathematicalConstants.H"
+#include "OSspecific.H"
+#include "Random.H"
+#include "rawSetWriter.H"
+#include "writeFile.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+    defineTypeNameAndDebug(irregular, 0);
+    addToRunTimeSelectionTable(waveModel, irregular, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::waveModels::AiryCoeffs
+Foam::waveModels::irregular::coeffs(const label i) const
+{
+    return AiryCoeffs(depth_, amplitudes_[i], lengths_[i], g());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveModels::irregular::irregular(const irregular& wave)
+:
+    waveModel(wave),
+    depth_(wave.depth_),
+    spectrum_(wave.spectrum_, false),
+    n_(wave.n_),
+    span_(wave.span_),
+    formatter_(wave.formatter_, false),
+    amplitudes_(wave.amplitudes_),
+    lengths_(wave.lengths_),
+    phases_(wave.phases_)
+{}
+
+
+Foam::waveModels::irregular::irregular
+(
+    const dictionary& dict,
+    const scalar g,
+    const word& modelName
+)
+:
+    waveModel(dict, g),
+    depth_(dict.lookupOrDefault<scalar>("depth", great)),
+    spectrum_(waveSpectrum::New(dict, g)),
+    n_(dict.lookup<label>("n")),
+    span_(dict.lookupOrDefault("span", Pair<scalar>(0.01, 0.99))),
+    formatter_
+    (
+        dict.found("setFormat")
+      ? setWriter::New(dict.lookup("setFormat"), dict).ptr()
+      : nullptr
+    ),
+    amplitudes_(n_),
+    lengths_(n_),
+    phases_(n_)
+{
+    // Get the bounds of the frequency space
+    const Pair<scalar> f01
+    (
+        spectrum_->fFraction(span_.first()),
+        spectrum_->fFraction(span_.second())
+    );
+
+    // Create a discretisation in frequency space
+    const scalarField x(scalarField(scalarList(identityMap(n_ + 1)))/n_);
+    const scalarField f((1 - x)*f01.first() + x*f01.second());
+
+    // Evaluate the integrals of the distribution
+    const scalarField integralS(spectrum_->integralS(f));
+    const scalarField integralFS(spectrum_->integralFS(f));
+
+    // Set the wave parameters using means constructed from the integrals over
+    // the frequency intervals
+    scalarList periods(n_);
+    forAll(amplitudes_, i)
+    {
+        const label j0 = n_ - 1 - i, j1 = n_ - i;
+
+        // Integrating the of power spectral density gives the square of the
+        // root-mean-squared value. So, a factor of two is needed to convert to
+        // peak amplitude.
+        amplitudes_[i] = sqrt(2*(integralS[j1] - integralS[j0]));
+
+        // Frequency is the integrated average, Integral(f*S)/Integral(S), for
+        // this interval. Period is the reciprocal of this value.
+        periods[i] =
+            (integralS[j1] - integralS[j0])/(integralFS[j1] - integralFS[j0]);
+
+        // Let the Airy calculation engine convert the above to wavelength
+        lengths_[i] =
+            AiryCoeffs
+            (
+                depth_,
+                amplitudes_[i],
+                periods[i],
+                g,
+                AiryCoeffs::celerity
+            ).length;
+    }
+
+    // Set random phases
+    Random rndGen(0);
+    forAll(phases_, i)
+    {
+        phases_[i] = rndGen.scalarAB(0, constant::mathematical::twoPi);
+    }
+
+    // Report
+    const scalar avgAmplitude =
+        sqrt(2*(integralS.last() - integralS.first())/(f.last() - f.first()));
+    const scalar avgPeroid =
+        (integralS.last() - integralS.first())
+       /(integralFS.last() - integralFS.first());
+    const scalar avgLength =
+        AiryCoeffs(depth_, avgAmplitude, avgPeroid, g, AiryCoeffs::celerity)
+       .length;
+    const string::size_type pad =
+        waveModel::typeName.length() + modelName.length() + 2;
+    Info<< waveModel::typeName << ": " << modelName
+        << ": min/average/max period = "
+        << periods.first() << '/' << avgPeroid << '/' << periods.last()
+        << endl << string(pad, ' ').c_str()
+        << "  min/average/max length = "
+        << lengths_.first() << '/' << avgLength << '/' << lengths_.last()
+        << endl;
+
+    // Write out the continuous distribution
+    if (Pstream::master() && formatter_.valid())
+    {
+        static const scalar t = 0.1;
+
+        const Pair<scalar> g01
+        (
+            max(rootVSmall, (1 + t)*(f01.first() - t*f01.second())),
+            - t*f01.first() + (1 + t)*f01.second()
+        );
+
+        const label nn = 1024;
+        const scalarField xx(scalarField(scalarList(identityMap(nn + 1)))/nn);
+        const scalarField ff((1 - xx)*g01.first() + xx*g01.second());
+
+        formatter_->write
+        (
+            cwd()
+           /functionObjects::writeFile::outputPrefix
+           /waveModel::typeName + "s"
+           /modelName + spectrum_->type().capitalise(),
+            "continuous",
+            coordSet(true, "f", ff),
+            "S",
+            spectrum_->S(ff)()
+        );
+    }
+
+    // Write out the sampled distribution
+    if (Pstream::master() && formatter_.valid())
+    {
+        scalarField ff(3*n_ + 1), SS(3*n_ + 1);
+
+        forAll(f, i)
+        {
+            if (i != 0)
+            {
+                ff[3*i - 1] = f[i];
+                SS[3*i - 1] = sqr(amplitudes_[n_ - i])/2/(f[i] - f[i - 1]);
+            }
+
+            ff[3*i] = f[i];
+            SS[3*i] = 0;
+
+            if (i != n_)
+            {
+                ff[3*i + 1] = f[i];
+                SS[3*i + 1] = sqr(amplitudes_[n_ - 1 - i])/2/(f[i + 1] - f[i]);
+            }
+        }
+
+        formatter_->write
+        (
+            cwd()
+           /functionObjects::writeFile::outputPrefix
+           /waveModel::typeName + "s"
+           /modelName + spectrum_->type().capitalise(),
+            "sampled",
+            coordSet(true, "f", ff),
+            "S",
+            SS
+        );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveModels::irregular::~irregular()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::scalar Foam::waveModels::irregular::celerity() const
+{
+    return coeffs(amplitudes_.size() - 1).celerity();
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveModels::irregular::elevation
+(
+    const scalar t,
+    const scalarField& x
+) const
+{
+    tmp<scalarField> tResult(new scalarField(x.size(), scalar(0)));
+    scalarField& result = tResult.ref();
+
+    forAll(amplitudes_, i)
+    {
+        result += coeffs(i).elevation(phases_[i], t, x);
+    }
+
+    return tResult;
+}
+
+
+Foam::tmp<Foam::vector2DField> Foam::waveModels::irregular::velocity
+(
+    const scalar t,
+    const vector2DField& xz
+) const
+{
+    tmp<vector2DField> tResult(new vector2DField(xz.size(), vector2D::zero));
+    vector2DField& result = tResult.ref();
+
+    forAll(amplitudes_, i)
+    {
+        result += coeffs(i).velocity(phases_[i], t, xz);
+    }
+
+    return tResult;
+}
+
+
+void Foam::waveModels::irregular::write(Ostream& os) const
+{
+    waveModel::write(os);
+
+    if (!coeffs(amplitudes_.size() - 1).deep())
+    {
+        writeEntry(os, "depth", depth_);
+    }
+
+    writeEntry(os, "spectrum", spectrum_->type());
+
+    os  << indent << word(spectrum_->type() + "Coeffs") << nl
+        << indent << token::BEGIN_BLOCK << nl << incrIndent;
+    spectrum_->write(os);
+    os  << decrIndent
+        << indent << token::END_BLOCK << nl;
+
+    writeEntry(os, "n", n_);
+    writeEntryIfDifferent(os, "span", Pair<scalar>(0.01, 0.99), span_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/irregular.H b/src/waves/waveModels/irregular/irregular.H
new file mode 100644
index 0000000000..090eb0458c
--- /dev/null
+++ b/src/waves/waveModels/irregular/irregular.H
@@ -0,0 +1,199 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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/>.
+
+Class
+    Foam::waveModels::irregular
+
+Description
+    Irregular wave model. Phase fraction and velocity field are built up from
+    multiple first-order waves, sampled from a selectable wave spectrum.
+
+Usage
+    \table
+        Property  | Description                     | Required? | Default
+        depth     | The water depth [m]             | no        | great
+        spectrum  | The wave spectrum               | yes       |
+        n         | The number of times to sample \\
+                    the spectrum                    | yes       |
+        span      | The fractional range across \\
+                    which to sample the spectrum    | no        | (0.01 0.99)
+        setFormat | The format with which to plot \\
+                    the spectrum                    | no        | none
+    \endtable
+
+    Example specification in constant/waveProperties:
+    \verbatim
+    waves
+    (
+        irregular
+        {
+            spectrum    PiersonMoskowitz;
+
+            PiersonMoskowitzCoeffs
+            {
+                U19_5       15;
+            }
+
+            n           16;
+            span        (0.01 0.99);
+        }
+    );
+    \endverbatim
+
+SourceFiles
+    irregular.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef irregular_H
+#define irregular_H
+
+#include "waveModel.H"
+#include "waveSpectrum.H"
+#include "AiryCoeffs.H"
+#include "setWriter.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveModels
+{
+
+/*---------------------------------------------------------------------------*\
+                            Class irregular Declaration
+\*---------------------------------------------------------------------------*/
+
+class irregular
+:
+    public waveModel
+{
+    // Private Data
+
+        //- Depth [m]
+        const scalar depth_;
+
+        //- Spectrum
+        autoPtr<waveSpectrum> spectrum_;
+
+        //- The number of times to sample the spectrum
+        const label n_;
+
+        //- The fractional range across which to sample the spectrum
+        const Pair<scalar> span_;
+
+        //- Formatter with which to plot the spectrum
+        autoPtr<setWriter> formatter_;
+
+        //- Amplitudes [m]
+        scalarList amplitudes_;
+
+        //- Wavelengths [m]
+        scalarList lengths_;
+
+        //- Phase offsets [rad]
+        scalarList phases_;
+
+
+    // Private Member Function
+
+        //- Return the i-th wave coefficients
+        AiryCoeffs coeffs(const label i) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("irregular");
+
+
+    // Constructors
+
+        //- Construct a copy
+        irregular(const irregular& wave);
+
+        //- Construct from a dictionary and gravity
+        irregular
+        (
+            const dictionary& dict,
+            const scalar g,
+            const word& modelName = irregular::typeName
+        );
+
+        //- Construct a clone
+        virtual autoPtr<waveModel> clone() const
+        {
+            return autoPtr<waveModel>(new irregular(*this));
+        }
+
+
+    //- Destructor
+    virtual ~irregular();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Get the depth
+            scalar depth() const
+            {
+                return depth_;
+            }
+
+
+        //- The wave celerity [m/s]
+        virtual scalar celerity() const;
+
+        //- 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 scalarField& x
+        ) const;
+
+        //- 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 vector2DField& xz
+        ) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace waveModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.C b/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.C
new file mode 100644
index 0000000000..114b5f2a8b
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.C
@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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 "JONSWAP.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveSpectra
+{
+    defineTypeNameAndDebug(JONSWAP, 0);
+    addToRunTimeSelectionTable(waveSpectrum, JONSWAP, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveSpectra::JONSWAP::JONSWAP
+(
+    const JONSWAP& spectrum
+)
+:
+    waveSpectrum(spectrum),
+    U10_(spectrum.U10_),
+    F_(spectrum.F_)
+{}
+
+
+Foam::waveSpectra::JONSWAP::JONSWAP
+(
+    const dictionary& dict,
+    const scalar g
+)
+:
+    waveSpectrum(dict, g),
+    U10_(dict.lookup<scalar>("U10")),
+    F_(dict.lookup<scalar>("F"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveSpectra::JONSWAP::~JONSWAP()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectra::JONSWAP::S
+(
+    const scalarField& f
+) const
+{
+    const scalar alpha = 0.076*pow(sqr(U10_)/(F_*g()), 0.22);
+    const scalarField w(twoPi*f);
+    const scalar wp = 22*pow(sqr(g())/(U10_*F_), 1.0/3.0);
+    const scalar gamma = 3.3;
+    const scalarField sigma(0.07 + 0.02*pos(w - wp));
+    const scalarField r(exp(- sqr(w - wp)/(2*sqr(sigma)*sqr(wp))));
+    return twoPi*alpha*sqr(g())/pow5(w)*exp(- 1.2*pow4(wp/w))*pow(gamma, r);
+}
+
+
+Foam::scalar Foam::waveSpectra::JONSWAP::fFraction(const scalar fraction) const
+{
+    const scalar wp = 22*pow(sqr(g())/(U10_*F_), 1.0/3.0);
+    const scalar f1 = wp/pow025(rootSmall/1.2)/twoPi;
+    return waveSpectrum::fFraction(fraction, f1);
+}
+
+
+void Foam::waveSpectra::JONSWAP::write(Ostream& os) const
+{
+    waveSpectrum::write(os);
+
+    writeEntry(os, "U10", U10_);
+    writeEntry(os, "F", F_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.H b/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.H
new file mode 100644
index 0000000000..b522a562e4
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/JONSWAP/JONSWAP.H
@@ -0,0 +1,152 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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/>.
+
+Class
+    Foam::JONSWAP
+
+Description
+    JONSWAP wave spectrum. This is similar to the Pierson-Moskowitz spectrum,
+    but with an additional empirical correction to account for the fetch
+    (distance to the lee shore).
+
+    References:
+    \verbatim
+        Hasselmann, K., Barnett, T. P., Bouws, E., Carlson, H., Cartwright, \\
+        D. E., Enke, K., ... & Walden, H. (1973).
+        Measurements of wind-wave growth and swell decay during the Joint \\
+        North Sea Wave Project (JONSWAP).
+        Ergaenzungsheft zur Deutschen Hydrographischen Zeitschrift, Reihe A.
+    \endverbatim
+
+    \verbatim
+        Stewart, R. H. (2008).
+        Introduction to physical oceanography.
+        Robert H. Stewart.
+    \endverbatim
+
+Usage
+    \table
+        Property | Description                             | Required? | Default
+        U10      | The air speed 10 metres above the \\
+                   surface [m/s]                           | yes       |
+        F        | The fetch (distance from the lee \\
+                   shore) [m]                              | yes       |
+    \endtable
+
+    Example specification:
+    \verbatim
+    spectrum    JONSWAP;
+
+    JONSWAPCoeffs
+    {
+        U10         10;
+        F           200e3;
+    }
+    \endverbatim
+
+See also
+    Foam::waveSpectra::PiersonMoskowitz
+
+SourceFiles
+    JONSWAP.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef JONSWAP_H
+#define JONSWAP_H
+
+#include "waveSpectrum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveSpectra
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class JONSWAP Declaration
+\*---------------------------------------------------------------------------*/
+
+class JONSWAP
+:
+    public waveSpectrum
+{
+    // Private Data
+
+        //- Reference velocity magnitude [m/s]
+        const scalar U10_;
+
+        //- Fetch (distance from the lee shore) [m]
+        const scalar F_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("JONSWAP");
+
+
+    // Constructors
+
+        //- Construct a copy
+        JONSWAP(const JONSWAP& spectrum);
+
+        //- Construct from a dictionary and gravity
+        JONSWAP(const dictionary& dict, const scalar g);
+
+        //- Construct a clone
+        virtual autoPtr<waveSpectrum> clone() const
+        {
+            return autoPtr<waveSpectrum>(new JONSWAP(*this));
+        }
+
+
+    //- Destructor
+    virtual ~JONSWAP();
+
+
+    // Member Functions
+
+        //- Evaluate the wave spectral density at the given frequencies [m^2/Hz]
+        virtual tmp<scalarField> S(const scalarField& f) const;
+
+        //- Return the frequency below which a given fraction of the spectrum's
+        //  total energy falls []
+        virtual scalar fFraction(const scalar fraction) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace waveSpectra
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.C b/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.C
new file mode 100644
index 0000000000..0458dd469c
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.C
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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 "PiersonMoskowitz.H"
+#include "mathematicalConstants.H"
+#include "addToRunTimeSelectionTable.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveSpectra
+{
+    defineTypeNameAndDebug(PiersonMoskowitz, 0);
+    addToRunTimeSelectionTable(waveSpectrum, PiersonMoskowitz, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveSpectra::PiersonMoskowitz::PiersonMoskowitz
+(
+    const PiersonMoskowitz& spectrum
+)
+:
+    waveSpectrum(spectrum),
+    U19_5_(spectrum.U19_5_),
+    alpha_(spectrum.alpha_),
+    beta_(spectrum.beta_)
+{}
+
+
+Foam::waveSpectra::PiersonMoskowitz::PiersonMoskowitz
+(
+    const dictionary& dict,
+    const scalar g
+)
+:
+    waveSpectrum(dict, g),
+    U19_5_(dict.lookup<scalar>("U19_5")),
+    alpha_(dict.lookupOrDefault<scalar>("alpha", 8.1e-3)),
+    beta_(dict.lookupOrDefault<scalar>("beta", 0.74))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveSpectra::PiersonMoskowitz::~PiersonMoskowitz()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectra::PiersonMoskowitz::S
+(
+    const scalarField& f
+) const
+{
+    const scalarField w(twoPi*f);
+    const scalar w0 = g()/U19_5_;
+    return twoPi*alpha_*sqr(g())/pow5(w)*exp(- beta_*pow4(w0/w));
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectra::PiersonMoskowitz::integralS
+(
+    const scalarField& f
+) const
+{
+    const scalarField w(twoPi*f);
+    const scalar w0 = g()/U19_5_;
+    const scalar integralSInf = twoPi*alpha_*sqr(g())/(4*beta_*pow4(w0))/twoPi;
+    return integralSInf*exp(- beta_*pow4(w0/w));
+}
+
+
+Foam::scalar Foam::waveSpectra::PiersonMoskowitz::fFraction
+(
+    const scalar fraction
+) const
+{
+    const scalar w0 = g()/U19_5_;
+    return w0/pow025(- log(fraction)/beta_)/twoPi;
+}
+
+
+void Foam::waveSpectra::PiersonMoskowitz::write(Ostream& os) const
+{
+    waveSpectrum::write(os);
+
+    writeEntry(os, "U19_5", U19_5_);
+    writeEntryIfDifferent<scalar>(os, "alpha", 8.1e-3, alpha_);
+    writeEntryIfDifferent<scalar>(os, "beta", 0.74, beta_);
+}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.H b/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.H
new file mode 100644
index 0000000000..20238a3cfb
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/PiersonMoskowitz/PiersonMoskowitz.H
@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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/>.
+
+Class
+    Foam::PiersonMoskowitz
+
+Description
+    Pierson-Moskowitz wave spectrum. This spectrum has the following form:
+
+        \f[
+            S(\omega) = \frac{\alpha g^2}{\omega^5} \exp \left(- \beta \\
+            \left( \frac{\omega_0}{\omega} \right)^4 \right)
+        \f]
+
+    \vartable
+        \omega    | angular frequency [rad/s], equal to \f$2 \pi f\f$
+        f         | frequency [Hz]
+        S(\omega) | spectral density [m^2/Hz]
+        \alpha    | coefficient, equal to 8.1e3
+        \beta     | coefficient, equal to 0.74
+        \omega_0  | reference angular frequency, equal to \f$g/U_{19.5}\f$
+        U_{19.5}  | reference velocity magnitude, 19.5 metres above the sea \\
+                    surface
+    \endvartable
+
+    References:
+    \verbatim
+        Pierson Jr, W. J., & Moskowitz, L. (1964).
+        A proposed spectral form for fully developed wind seas based on the \\
+        similarity theory of SA Kitaigorodskii.
+        Journal of geophysical research, 69(24), 5181-5190.
+    \endverbatim
+
+    \verbatim
+        Stewart, R. H. (2008).
+        Introduction to physical oceanography.
+        Robert H. Stewart.
+    \endverbatim
+
+Usage
+    \table
+        Property | Description                             | Required? | Default
+        U19_5    | The air speed 19.5 metres above the \\
+                   surface [m/s]                           | yes       |
+    \endtable
+
+    Example specification:
+    \verbatim
+    spectrum    PiersonMoskowitz;
+
+    PiersonMoskowitzCoeffs
+    {
+        U19_5       15;
+    }
+    \endverbatim
+
+SourceFiles
+    PiersonMoskowitz.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef PiersonMoskowitz_H
+#define PiersonMoskowitz_H
+
+#include "waveSpectrum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace waveSpectra
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class PiersonMoskowitz Declaration
+\*---------------------------------------------------------------------------*/
+
+class PiersonMoskowitz
+:
+    public waveSpectrum
+{
+    // Private Data
+
+        //- Reference velocity magnitude
+        const scalar U19_5_;
+
+        //- Alpha coefficient
+        const scalar alpha_;
+
+        //- Beta coefficient
+        const scalar beta_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("PiersonMoskowitz");
+
+
+    // Constructors
+
+        //- Construct a copy
+        PiersonMoskowitz(const PiersonMoskowitz& spectrum);
+
+        //- Construct from a dictionary and gravity
+        PiersonMoskowitz(const dictionary& dict, const scalar g);
+
+        //- Construct a clone
+        virtual autoPtr<waveSpectrum> clone() const
+        {
+            return autoPtr<waveSpectrum>(new PiersonMoskowitz(*this));
+        }
+
+
+    //- Destructor
+    virtual ~PiersonMoskowitz();
+
+
+    // Member Functions
+
+        //- Evaluate the wave spectral density at the given frequencies [m^2/Hz]
+        virtual tmp<scalarField> S(const scalarField& f) const;
+
+        //- Evaluate the integral of the wave spectral density at the given
+        //  frequencies [m^2]
+        virtual tmp<scalarField> integralS(const scalarField& f) const;
+
+        //- Return the frequency below which a given fraction of the spectrum's
+        //  total power falls []
+        virtual scalar fFraction(const scalar fraction) const;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace waveSpectra
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.C b/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.C
new file mode 100644
index 0000000000..ab7a6e64e3
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.C
@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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 "waveSpectrum.H"
+#include "unintegrable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(waveSpectrum, 0);
+    defineRunTimeSelectionTable(waveSpectrum, dictionary);
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectrum::refined
+(
+    const label n,
+    const scalarField& x
+)
+{
+    if (n == x.size())
+    {
+        return x;
+    }
+
+    tmp<scalarField> txx(new scalarField(n, -vGreat));
+    scalarField& xx = txx.ref();
+
+    if (n > x.size())
+    {
+        for (label i = 0; i < x.size() - 1; ++ i)
+        {
+            const label j0 = i*(n - 1)/(x.size() - 1);
+            const label j1 = (i + 1)*(n - 1)/(x.size() - 1);
+
+            xx[j0] = x[i];
+
+            for (label j = j0 + 1; j < j1; ++ j)
+            {
+                const scalar f = scalar(j - j0)/(j1 - j0);
+
+                xx[j] = (1 - f)*x[i] + f*x[i + 1];
+            }
+        }
+
+        xx.last() = x.last();
+    }
+    else
+    {
+        forAll(xx, j)
+        {
+            const label i = j*(x.size() - 1)/(n - 1);
+
+            xx[j] = x[i];
+        }
+    }
+
+    return txx;
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectrum::refined
+(
+    const label n,
+    const tmp<scalarField>& x
+)
+{
+    if (n == x().size())
+    {
+        return x;
+    }
+
+    return refined(n, x());
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+Foam::scalar Foam::waveSpectrum::fFraction
+(
+    const scalar fraction,
+    const scalar f1
+) const
+{
+    // Integrate from zero to the maximum frequency
+    const scalarField ff(refined(n_, scalarField(scalarList({small*f1, f1}))));
+    const scalarField integralSS(integralS(ff));
+
+    // Sample the integrated values for the frequency at the given fraction
+    return
+        distributions::unintegrable::sample
+        (
+            ff,
+            integralSS/integralSS.last(),
+            fraction
+        );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::waveSpectrum::waveSpectrum(const waveSpectrum& spectrum)
+:
+    g_(spectrum.g_),
+    n_(spectrum.n_)
+{}
+
+
+Foam::waveSpectrum::waveSpectrum(const dictionary& dict, const scalar g)
+:
+    g_(g),
+    n_((1<<dict.lookupOrDefault<label>("level", 16)) + 1)
+{}
+
+
+// * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::waveSpectrum>
+Foam::waveSpectrum::New(const dictionary& dict, const scalar g)
+{
+    const word type = dict.lookup<word>("spectrum");
+
+    if (debug)
+    {
+        Info<< "Selecting " << waveSpectrum::typeName << " " << type << endl;
+    }
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(type);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorInFunction
+            << "Unknown " << waveSpectrum::typeName << " " << type << nl << nl
+            << "Valid spectrum types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return cstrIter()(dict.optionalSubDict(type + "Coeffs"), g);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::waveSpectrum::~waveSpectrum()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectrum::integralS
+(
+    const scalarField& f
+) const
+{
+    const scalarField ff(refined(n_, f));
+    return
+        refined
+        (
+            f.size(),
+            distributions::unintegrable::integrate(ff, S(ff))
+        );
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveSpectrum::integralFS
+(
+    const scalarField& f
+) const
+{
+    const scalarField ff(refined(n_, f));
+    return
+        refined
+        (
+            f.size(),
+            distributions::unintegrable::integrateX(ff, S(ff))
+        );
+}
+
+
+void Foam::waveSpectrum::write(Ostream& os) const
+{}
+
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.H b/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.H
new file mode 100644
index 0000000000..2d7795c770
--- /dev/null
+++ b/src/waves/waveModels/irregular/waveSpectra/waveSpectrum/waveSpectrum.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2023 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/>.
+
+Class
+    Foam::waveSpectrum
+
+Description
+    Base class for wave spectra
+
+SourceFiles
+    waveSpectrum.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef waveSpectrum_H
+#define waveSpectrum_H
+
+#include "dictionary.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class waveSpectrum Declaration
+\*---------------------------------------------------------------------------*/
+
+class waveSpectrum
+{
+    // Private Data
+
+        //- Gravitational acceleration [m/s^2]
+        const scalar g_;
+
+        //- Number of intervals to use for integration
+        const label n_;
+
+
+    // Private Member Functions
+
+        //- Refine or coarsen the given field to the given size, using linear
+        //  interpolation to construct intermediate values. Guaranteed
+        //  symmetric, so refined(x.size(), refined(n, x)) will return x (if n
+        //  is larger than x.size()).
+        static tmp<scalarField> refined(const label n, const scalarField& x);
+
+        //- Tmp variant of the above
+        static tmp<scalarField> refined
+        (
+            const label n,
+            const tmp<scalarField>& x
+        );
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Return the frequency below which a given fraction of the spectrum's
+        //  total energy falls, given also the maximum frequency with
+        //  appreciable energy
+        scalar fFraction(const scalar fraction, const scalar f1) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("waveSpectrum");
+
+
+    // Declare runtime construction
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        waveSpectrum,
+        dictionary,
+        (const dictionary& dict, const scalar g),
+        (dict, g)
+    );
+
+
+    // Constructors
+
+        //- Construct a copy
+        waveSpectrum(const waveSpectrum& spectrum);
+
+        //- Construct from a dictionary and gravity
+        waveSpectrum(const dictionary& dict, const scalar g);
+
+        //- Construct a clone
+        virtual autoPtr<waveSpectrum> clone() const = 0;
+
+
+    // Selectors
+
+        //- Select given a dictionary and gravity
+        static autoPtr<waveSpectrum> New
+        (
+            const dictionary& dict,
+            const scalar g
+        );
+
+
+    //- Destructor
+    virtual ~waveSpectrum();
+
+
+    // Member Functions
+
+        //- Access the gravitation acceleration [m/s^2]
+        inline scalar g() const
+        {
+            return g_;
+        }
+
+        //- Evaluate the wave spectral density at the given frequencies [m^2/Hz]
+        virtual tmp<scalarField> S(const scalarField& f) const = 0;
+
+        //- Evaluate the integral of the wave spectral density at the given
+        //  frequencies [m^2]
+        virtual tmp<scalarField> integralS(const scalarField& f) const;
+
+        //- Evaluate the integral of the wave spectral density multiplied by
+        //  the frequency at the given frequencies [Hz m^2]
+        virtual tmp<scalarField> integralFS(const scalarField& f) const;
+
+        //- Return the frequency below which a given fraction of the spectrum's
+        //  total energy falls
+        virtual scalar fFraction(const scalar fraction) const = 0;
+
+        //- Write
+        virtual void write(Ostream& os) const;
+
+
+    // Member Operators
+
+        //- Disallow default bitwise assignment
+        void operator=(const waveSpectrum&) = delete;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waves/waveModels/solitary/solitary.H b/src/waves/waveModels/solitary/solitary.H
index 7776d199e6..e3faa84090 100644
--- a/src/waves/waveModels/solitary/solitary.H
+++ b/src/waves/waveModels/solitary/solitary.H
@@ -36,6 +36,27 @@ Description
 
     See pages 314-317.
 
+Usage
+    \table
+        Property  | Description               | Required? | Default
+        depth     | The water depth [m]       | no        | great
+        amplitude | Peak amplitude [m]        | yes       |
+        offset    | The positional offset [m] | yes       |
+    \endtable
+
+    Example specification in constant/waveProperties:
+    \verbatim
+    waves
+    (
+        solitary
+        {
+            length      40;
+            amplitude   0.5;
+            offset      0;
+        }
+    );
+    \endverbatim
+
 SourceFiles
     solitary.C
 
diff --git a/src/waves/waveModels/waveModel/waveModel.C b/src/waves/waveModels/waveModel/waveModel.C
index 6d3ed6b8ca..78e4b4d906 100644
--- a/src/waves/waveModels/waveModel/waveModel.C
+++ b/src/waves/waveModels/waveModel/waveModel.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,11 +42,7 @@ Foam::waveModel::waveModel(const waveModel& wave)
 {}
 
 
-Foam::waveModel::waveModel
-(
-    const dictionary& dict,
-    const scalar g
-)
+Foam::waveModel::waveModel(const dictionary& dict, const scalar g)
 :
     g_(g)
 {}
diff --git a/src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H b/src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H
index ef2da0116c..3027207904 100644
--- a/src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H
+++ b/src/waves/waveSuperpositions/waveAtmBoundaryLayerSuperposition/waveAtmBoundaryLayerSuperposition.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,21 +37,27 @@ Description
 
 Usage
     \table
-        Property | Description                              | Req'd? | Default
-        UGasRef  | The gas velocity at the reference height | yes    |
+        Property | Description                             | Required? | Default
+        UGasRef  | The gas velocity at the reference \\
+                   height                                  | yes       |
         hRef     | The reference height relative to the \\
-                   origin of the wave coordinate system     | yes    |
-        hWaveMin | The minimum wave elevation               | yes    |
-        hWaveMax | The maximum wave elevation               | yes    |
+                   origin of the wave coordinate system    | yes       |
+        hWaveMin | The minimum wave elevation              | yes       |
+        hWaveMax | The maximum wave elevation              | yes       |
     \endtable
 
     Example specification:
     \verbatim
     type        waveAtmBoundaryLayer;
+
     // parameters for waveSuperposition ...
+
     UGasRef     (10 0 0);
+
     hRef        10;
+
     hWaveMin    -2;
+
     hWaveMax    3;
     \endverbatim
 
diff --git a/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C b/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C
index 9b144e78e9..5d2c4c6476 100644
--- a/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C
+++ b/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.C
@@ -256,7 +256,7 @@ Foam::scalar Foam::waveSuperposition::maxWaveSpeed(const scalar t) const
         maxCelerity = max(waveModels_[wavei].celerity(), maxCelerity);
     }
 
-    return mag(maxCelerity + (direction_& UMean_->value(t)));
+    return mag(maxCelerity + (direction_ & UMean_->value(t)));
 }
 
 
diff --git a/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H b/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H
index 447b9400a0..358e7606ae 100644
--- a/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H
+++ b/src/waves/waveSuperpositions/waveSuperposition/waveSuperposition.H
@@ -28,25 +28,31 @@ Description
     A wrapper around a list of wave models. Superimposes the modelled values of
     elevation and velocity. The New method looks up or constructs an
     instance of this class on demand and returns a reference. Properties are
-    read from a dictionary in constant.
+    read from a waveProperties 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    |
-        UMean      | velocity of the mean flow               | 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
+        Property        | Description                      | Required? | 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       |
+        UMean           | velocity of the mean flow        | 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
@@ -64,9 +70,13 @@ Usage
             angle       0;
         }
     );
+
     UMean       (2 0 0);
+
     scale       table ((100 1) (200 0));
+
     crossScale  constant 1;
+
     heightAboveWave no;
     \endverbatim