diff --git a/src/waves/waveModels/Airy/Airy.C b/src/waves/waveModels/Airy/Airy.C
index bea8cb9856..deca8efe81 100644
--- a/src/waves/waveModels/Airy/Airy.C
+++ b/src/waves/waveModels/Airy/Airy.C
@@ -153,6 +153,20 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::Airy::velocity
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::waveModels::Airy::pressure
+(
+    const scalar t,
+    const scalar u,
+    const vector2DField& xz
+) const
+{
+    // It is a fluke of the formulation that the time derivative of the velocity
+    // potential equals the x-derivative multipled by the celerity. This allows
+    // for this shortcut in evaluating the unsteady pressure.
+    return celerity()*velocity(t, u, xz)->component(0);
+}
+
+
 void Foam::waveModels::Airy::write(Ostream& os) const
 {
     waveModel::write(os);
diff --git a/src/waves/waveModels/Airy/Airy.H b/src/waves/waveModels/Airy/Airy.H
index 1ee70783d5..79307c2850 100644
--- a/src/waves/waveModels/Airy/Airy.H
+++ b/src/waves/waveModels/Airy/Airy.H
@@ -169,6 +169,16 @@ public:
             const vector2DField& xz
         ) const;
 
+        //- Get the wave pressure at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, and z with
+        //  negative gravity.
+        virtual tmp<scalarField> pressure
+        (
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const;
+
         //- Write
         virtual void write(Ostream& os) const;
 };
diff --git a/src/waves/waveModels/solitary/solitary.C b/src/waves/waveModels/solitary/solitary.C
index 5143c4034d..b2a22d8a28 100644
--- a/src/waves/waveModels/solitary/solitary.C
+++ b/src/waves/waveModels/solitary/solitary.C
@@ -153,6 +153,18 @@ Foam::tmp<Foam::vector2DField> Foam::waveModels::solitary::velocity
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::waveModels::solitary::pressure
+(
+    const scalar t,
+    const scalar u,
+    const vector2DField& xz
+) const
+{
+    NotImplemented;
+    return tmp<scalarField>(nullptr);
+}
+
+
 void Foam::waveModels::solitary::write(Ostream& os) const
 {
     waveModel::write(os);
diff --git a/src/waves/waveModels/solitary/solitary.H b/src/waves/waveModels/solitary/solitary.H
index 5b7880f0b2..11cf45020b 100644
--- a/src/waves/waveModels/solitary/solitary.H
+++ b/src/waves/waveModels/solitary/solitary.H
@@ -157,6 +157,16 @@ public:
             const vector2DField& xz
         ) const;
 
+        //- Get the wave pressure at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, and z with
+        //  negative gravity.
+        virtual tmp<scalarField> pressure
+        (
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const;
+
         //- Write
         virtual void write(Ostream& os) const;
 };
diff --git a/src/waves/waveModels/waveModel/waveModel.H b/src/waves/waveModels/waveModel/waveModel.H
index cccafe0d6a..bff683753a 100644
--- a/src/waves/waveModels/waveModel/waveModel.H
+++ b/src/waves/waveModels/waveModel/waveModel.H
@@ -147,6 +147,16 @@ public:
             const vector2DField& xz
         ) const = 0;
 
+        //- Get the wave pressure at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, and z with
+        //  negative gravity.
+        virtual tmp<scalarField> pressure
+        (
+            const scalar t,
+            const scalar u,
+            const vector2DField& xz
+        ) const = 0;
+
         //- Write
         virtual void write(Ostream& os) const;
 };
diff --git a/src/waves/waveSuperposition/waveSuperposition.C b/src/waves/waveSuperposition/waveSuperposition.C
index 67c068acdc..90aeb15024 100644
--- a/src/waves/waveSuperposition/waveSuperposition.C
+++ b/src/waves/waveSuperposition/waveSuperposition.C
@@ -108,6 +108,38 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::velocity
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pressure
+(
+    const scalar t,
+    const vectorField& xyz
+) const
+{
+    scalarField result(xyz.size(), 0);
+
+    forAll(waveModels_, wavei)
+    {
+        const vector2D d(cos(waveAngles_[wavei]), sin(waveAngles_[wavei]));
+        const vector2DField xz
+        (
+            zip
+            (
+                d & zip(xyz.component(0), xyz.component(1)),
+                tmp<scalarField>(xyz.component(2))
+            )
+        );
+        const vector2DField uw
+        (
+            waveModels_[wavei].velocity(t, d.x()*speed_, xz)
+        );
+        result += waveModels_[wavei].pressure(t, d.x()*speed_, xz);
+    }
+
+    tmp<scalarField> s = scale(zip(xyz.component(0), xyz.component(1)));
+
+    return s*result;
+}
+
+
 Foam::tmp<Foam::scalarField> Foam::waveSuperposition::scale
 (
     const vector2DField& xy
@@ -266,6 +298,31 @@ Foam::tmp<Foam::vectorField> Foam::waveSuperposition::UGas
 }
 
 
+Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pLiquid
+(
+    const scalar t,
+    const vectorField& p
+) const
+{
+    tensor axes;
+    scalar u;
+    vectorField xyz(p.size());
+    transformation(p, axes, u, xyz);
+
+    return pressure(t, xyz);
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::waveSuperposition::pGas
+(
+    const scalar t,
+    const vectorField& p
+) const
+{
+    return - pLiquid(t, p);
+}
+
+
 void Foam::waveSuperposition::write(Ostream& os) const
 {
     os.writeKeyword("origin") << origin_ << token::END_STATEMENT << nl;
diff --git a/src/waves/waveSuperposition/waveSuperposition.H b/src/waves/waveSuperposition/waveSuperposition.H
index b2fefe7887..e5f802141d 100644
--- a/src/waves/waveSuperposition/waveSuperposition.H
+++ b/src/waves/waveSuperposition/waveSuperposition.H
@@ -99,11 +99,12 @@ class waveSuperposition
         //- Get the wave velocity at a given time, mean velocity and local
         //  coordinates. Local x is aligned with the mean velocity, z with
         //  negative gravity, and y is perpendicular to both.
-        tmp<vectorField> velocity
-        (
-            const scalar t,
-            const vectorField& xyz
-        ) const;
+        tmp<vectorField> velocity(const scalar t, const vectorField& xyz) const;
+
+        //- Get the wave pressure at a given time, mean velocity and local
+        //  coordinates. Local x is aligned with the mean velocity, z with
+        //  negative gravity, and y is perpendicular to both.
+        tmp<scalarField> pressure(const scalar t, const vectorField& xyz) const;
 
         //- Get the scaling factor, calculated from the optional scaling
         //  functions. X and y are the same as for the elevation method.
@@ -139,6 +140,12 @@ public:
         //- Get the gas velocity at a given time and global positions
         tmp<vectorField> UGas(const scalar t, const vectorField& p) const;
 
+        //- Get the liquid pressure at a given time and global positions
+        tmp<scalarField> pLiquid(const scalar t, const vectorField& p) const;
+
+        //- Get the gas pressure at a given time and global positions
+        tmp<scalarField> pGas(const scalar t, const vectorField& p) const;
+
         //- Get the mean flow velocity
         inline vector UMean() const
         {