diff --git a/src/fvOptions/Make/files b/src/fvOptions/Make/files
index 2e91b0709551ab989f9bb81c86e355945545336f..085a1af9450da604231126a6fbbf36a5f7e5150d 100644
--- a/src/fvOptions/Make/files
+++ b/src/fvOptions/Make/files
@@ -18,6 +18,7 @@ $(derivedSources)/explicitPorositySource/explicitPorositySource.C
$(derivedSources)/jouleHeatingSource/jouleHeatingSource.C
$(derivedSources)/meanVelocityForce/meanVelocityForce.C
$(derivedSources)/meanVelocityForce/patchMeanVelocityForce/patchMeanVelocityForce.C
+$(derivedSources)/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C
$(derivedSources)/phaseLimitStabilization/phaseLimitStabilization.C
$(derivedSources)/radialActuationDiskSource/radialActuationDiskSource.C
$(derivedSources)/rotorDiskSource/rotorDiskSource.C
@@ -57,4 +58,5 @@ $(derivedConstraints)/velocityDampingConstraint/velocityDampingConstraint.C
corrections/limitTemperature/limitTemperature.C
corrections/limitVelocity/limitVelocity.C
+
LIB = $(FOAM_LIBBIN)/libfvOptions
diff --git a/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C b/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C
new file mode 100644
index 0000000000000000000000000000000000000000..ae76698548280db155087bd51c024e83df0f0de4
--- /dev/null
+++ b/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.C
@@ -0,0 +1,241 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "multiphaseStabilizedTurbulence.H"
+#include "fvMatrices.H"
+#include "turbulentTransportModel.H"
+#include "gravityMeshObject.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+ defineTypeNameAndDebug(multiphaseStabilizedTurbulence, 0);
+
+ addToRunTimeSelectionTable
+ (
+ option,
+ multiphaseStabilizedTurbulence,
+ dictionary
+ );
+}
+}
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::fv::multiphaseStabilizedTurbulence::multiphaseStabilizedTurbulence
+(
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+)
+:
+ option(sourceName, modelType, dict, mesh),
+ rhoName_(coeffs_.getOrDefault<word>("rho", "rho")),
+ Cmu_
+ (
+ dimensionedScalar::lookupOrAddToDict
+ (
+ "Cmu",
+ coeffs_,
+ 0.09
+ )
+ ),
+ C_
+ (
+ dimensionedScalar::lookupOrAddToDict
+ (
+ "C",
+ coeffs_,
+ 1.51
+ )
+ ),
+ lambda2_
+ (
+ dimensionedScalar::lookupOrAddToDict
+ (
+ "lambda2",
+ coeffs_,
+ 0.1
+ )
+ ),
+ alpha_
+ (
+ dimensionedScalar::lookupOrAddToDict
+ (
+ "alpha",
+ coeffs_,
+ 1.36
+ )
+ )
+{
+ fieldNames_.setSize(2, "undefined");
+
+ // Note: incompressible only
+ const auto* turbPtr =
+ mesh_.findObject<incompressible::turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ if (turbPtr)
+ {
+ const tmp<volScalarField>& tk = turbPtr->k();
+ fieldNames_[0] = tk().name();
+
+ const tmp<volScalarField>& tnut = turbPtr->nut();
+ fieldNames_[1] = tnut().name();
+
+ Log << " Applying model to " << fieldNames_[0]
+ << " and " << fieldNames_[1] << endl;
+ }
+ else
+ {
+ FatalErrorInFunction
+ << "Unable to find incompressible turbulence model"
+ << exit(FatalError);
+ }
+
+ applied_.setSize(fieldNames_.size(), false);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::fv::multiphaseStabilizedTurbulence::addSup
+(
+ const volScalarField& rho,
+ fvMatrix<scalar>& eqn,
+ const label fieldi
+)
+{
+ // Not applicable to compressible cases
+ NotImplemented;
+}
+
+
+void Foam::fv::multiphaseStabilizedTurbulence::addSup
+(
+ fvMatrix<scalar>& eqn,
+ const label fieldi
+)
+{
+ if (fieldi != 0)
+ {
+ return;
+ }
+
+ Log << this->name() << ": applying bouyancy production term to "
+ << eqn.psi().name() << endl;
+
+ // Buoyancy production in k eqn
+
+ const auto* turbPtr =
+ mesh_.findObject<incompressible::turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ if (!turbPtr)
+ {
+ FatalErrorInFunction
+ << "Unable to find incompressible turbulence model"
+ << exit(FatalError);
+ }
+
+
+ tmp<volScalarField> tepsilon = turbPtr->epsilon();
+ const volScalarField& epsilon = tepsilon();
+ const volScalarField& k = eqn.psi();
+
+ // Note: using solver density field for incompressible multiphase cases
+ const auto& rho = mesh_.lookupObject<volScalarField>(rhoName_);
+
+ const auto& g = meshObjects::gravity::New(mesh_.time());
+
+ const dimensionedScalar eps0("eps0", epsilon.dimensions(), SMALL);
+
+ // Note: differing from reference by replacing nut/k by Cmu*k/epsilon
+ const volScalarField GbyK
+ (
+ "GbyK",
+ Cmu_*k/(epsilon + eps0)*alpha_*(g & fvc::grad(rho))/rho
+ );
+
+ return eqn -= fvm::SuSp(GbyK, k);
+}
+
+
+void Foam::fv::multiphaseStabilizedTurbulence::correct(volScalarField& field)
+{
+ if (field.name() != fieldNames_[1])
+ {
+ return;
+ }
+
+ Log << this->name() << ": correcting " << field.name() << endl;
+
+ const auto* turbPtr =
+ mesh_.findObject<incompressible::turbulenceModel>
+ (
+ turbulenceModel::propertiesName
+ );
+
+ // nut correction
+ const auto& U = turbPtr->U();
+ tmp<volScalarField> tepsilon = turbPtr->epsilon();
+ const auto& epsilon = tepsilon();
+ tmp<volScalarField> tk = turbPtr->k();
+ const auto& k = tk();
+
+ tmp<volTensorField> tgradU = fvc::grad(U);
+ const auto& gradU = tgradU();
+ const dimensionedScalar pSmall("pSmall", dimless/sqr(dimTime), SMALL);
+ const volScalarField pRat
+ (
+ magSqr(symm(gradU))/(magSqr(skew(gradU)) + pSmall)
+ );
+
+ const volScalarField epsilonTilde
+ (
+ max
+ (
+ epsilon,
+ lambda2_*C_*pRat*epsilon
+ )
+ );
+
+ field = Cmu_*sqr(k)/epsilonTilde;
+ field.correctBoundaryConditions();
+}
+
+
+// ************************************************************************* //
diff --git a/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.H b/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.H
new file mode 100644
index 0000000000000000000000000000000000000000..5543cc28e3f0358e1e953ac7f7b63d3e6657b38c
--- /dev/null
+++ b/src/fvOptions/sources/derived/multiphaseStabilizedTurbulence/multiphaseStabilizedTurbulence.H
@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2019 OpenCFD Ltd
+-------------------------------------------------------------------------------
+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::fv::multiphaseStabilizedTurbulence
+
+Group
+ grpFvOptionsSources
+
+Description
+ Applies corrections to the turbulence kinetic energy equation and turbulence
+ viscosity field for incompressible multiphase flow cases.
+
+ Turbulence kinetic energy is over-predicted in VOF solvers at the phase
+ interface and throughout the water column in nearly-potential flow regions
+ beneath surface waves.
+
+ This fvOption applies corrections based on the references:
+ \verbatim
+ Buoyancy source term in turbulence kinetic energy equation:
+ Devolder, B., Rauwoens, P., and Troch, P. (2017).
+ Application of a buoyancy-modified k-w SST turbulence model to
+ simulate wave run-up around a monopile subjected to regular waves
+ using OpenFOAM.
+ Coastal Engineering, 125, 81-94.
+
+ Correction to turbulence viscosity:
+ Larsen, B.E. and Fuhrman, D.R. (2018).
+ On the over-production of turbulence beneath surface waves in
+ Reynolds-averaged Navier-Stokes models
+ J. Fluid Mech, 853, 419-460
+ \endverbatim
+
+Usage
+ Example usage:
+
+ \verbatim
+ multiphaseStabilizedTurbulence1
+ {
+ type multiphaseStabilizedTurbulence;
+ active yes;
+
+ multiphaseStabilizedTurbulenceCoeffs
+ {
+ // Optional coefficients
+ lambda2 0.1; // A value of 0 sets the nut correction to 0
+ Cmu 0.09; // from k-epsilon model
+ C 1.51; // from k-omega model
+ alpha 1.36; // 1/Prt
+ }
+ }
+ \endverbatim
+
+ The model C coefficient for the k-epsilon model equates to C2/C1 = 1.33;
+ the (default) value of 1.51 comes from the k-omega model and is more
+ conservative.
+
+SourceFiles
+ multiphaseStabilizedTurbulence.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fv_multiphaseStabilizedTurbulence_H
+#define fv_multiphaseStabilizedTurbulence_H
+
+#include "fvOption.H"
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+
+/*---------------------------------------------------------------------------*\
+ Class multiphaseStabilizedTurbulence Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseStabilizedTurbulence
+:
+ public option
+{
+ // Private data
+
+ //- Name of density field
+ const word rhoName_;
+
+
+ // Model coefficients
+
+ //- k-epsilon model Cmu coefficient
+ dimensionedScalar Cmu_;
+
+ //- Model coefficient
+ // For k-epsilon models this equates to C2/C1 = 1.33 and for
+ // k-omega = 1.51. Here the default is the more conservative 1.51
+ dimensionedScalar C_;
+
+ //- lambda2 coefficient; default = 0.1
+ dimensionedScalar lambda2_;
+
+ //- Buoyancy coefficient
+ dimensionedScalar alpha_;
+
+
+ // Private Member Functions
+
+ //- No copy construct
+ multiphaseStabilizedTurbulence
+ (
+ const multiphaseStabilizedTurbulence&
+ ) = delete;
+
+ //- No copy assignment
+ void operator=(const multiphaseStabilizedTurbulence&) = delete;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("multiphaseStabilizedTurbulence");
+
+
+ // Constructors
+
+ //- Construct from explicit source name and mesh
+ multiphaseStabilizedTurbulence
+ (
+ const word& sourceName,
+ const word& modelType,
+ const dictionary& dict,
+ const fvMesh& mesh
+ );
+
+
+ // Member Functions
+
+ //- Add explicit contribution to compressible k equation
+ virtual void addSup
+ (
+ const volScalarField& rho,
+ fvMatrix<scalar>& eqn,
+ const label fieldi
+ );
+
+ //- Add explicit contribution to incompressible k equation
+ virtual void addSup
+ (
+ fvMatrix<scalar>& eqn,
+ const label fieldi
+ );
+
+ //- Correct the turbulence viscosity
+ virtual void correct(volScalarField& field);
+
+ //- Read source dictionary
+ virtual bool read(const dictionary& dict)
+ {
+ return true;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //