diff --git a/src/waveModels/Make/files b/src/waveModels/Make/files
index 76946b3c77a41a83b156749f4ae8e9d13b3ff1fe..b6f5f017a06bc348c220fbab732967254b4ff02d 100644
--- a/src/waveModels/Make/files
+++ b/src/waveModels/Make/files
@@ -15,11 +15,13 @@ waveGenerationModels/derived/StokesI/StokesIWaveModel.C
 waveGenerationModels/derived/StokesV/StokesVWaveModel.C
 waveGenerationModels/derived/irregularMultiDirectional/irregularMultiDirectionalWaveModel.C
 
-
 waveAbsorptionModels/base/waveAbsorptionModel/waveAbsorptionModel.C
 waveAbsorptionModels/derived/shallowWaterAbsorption/shallowWaterAbsorption.C
 
 derivedFvPatchFields/waveVelocity/waveVelocityFvPatchVectorField.C
 derivedFvPatchFields/waveAlpha/waveAlphaFvPatchScalarField.C
 
+fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
+fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
+
 LIB = $(FOAM_LIBBIN)/libwaveModels
diff --git a/src/waveModels/Make/options b/src/waveModels/Make/options
index 71b7873964d544eddf96d22aa40f4c3372c23c9c..c51882a5360477a5ecb79b12ab3f5799c70f038d 100644
--- a/src/waveModels/Make/options
+++ b/src/waveModels/Make/options
@@ -1,5 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
 
 LIB_LIBS = \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
new file mode 100644
index 0000000000000000000000000000000000000000..ba5818bb9616a407c9d3084fab9bd08b4ca1f1e4
--- /dev/null
+++ b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.C
@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
+-------------------------------------------------------------------------------
+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 "multiphaseMangrovesSource.H"
+#include "mathematicalConstants.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "fvmDdt.H"
+#include "fvmSup.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(multiphaseMangrovesSource, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        multiphaseMangrovesSource,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::multiphaseMangrovesSource::multiphaseMangrovesSource
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(name, modelType, dict, mesh),
+    aZone_(),
+    NZone_(),
+    CmZone_(),
+    CdZone_()
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::multiphaseMangrovesSource::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<vector>& eqn,
+    const label fieldi
+)
+{
+    const volVectorField& U = eqn.psi();
+
+    volScalarField DragForceMangroves
+    (
+        IOobject
+        (
+            typeName + ":DragForceMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("DragForceMangroves", dimDensity/dimTime, 0)
+    );
+
+    volScalarField InertiaForceMangroves
+    (
+        IOobject
+        (
+            typeName + ":InertiaForceMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("InertiaForceMangroves", dimless, 0)
+    );
+
+    const scalar pi = constant::mathematical::pi;
+
+    forAll(zoneIDs_, i)
+    {
+        const labelList& zones = zoneIDs_[i];
+
+        forAll(zones, j)
+        {
+            const label zonei = zones[j];
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            forAll(cz, k)
+            {
+                const label celli = cz[k];
+                const scalar a = aZone_[i];
+                const scalar N = NZone_[i];
+                const scalar Cm = CmZone_[i];
+                const scalar Cd = CdZone_[i];
+
+                const vector& Uc = U[celli];
+
+                DragForceMangroves[celli] = 0.5*Cd*a*N*mag(Uc);
+                InertiaForceMangroves[celli] = 0.25*(Cm+1)*pi*a*a;
+            }
+        }
+    }
+
+    DragForceMangroves.correctBoundaryConditions();
+    InertiaForceMangroves.correctBoundaryConditions();
+
+    fvMatrix<vector> MangrovesEqn
+    (
+      - fvm::Sp(rho*DragForceMangroves, U)
+      - rho*InertiaForceMangroves*fvm::ddt(rho, U)
+    );
+
+    // Contributions are added to RHS of momentum equation
+    eqn += MangrovesEqn;
+}
+
+
+bool Foam::fv::multiphaseMangrovesSource::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        if (coeffs_.found("UNames"))
+        {
+            coeffs_.lookup("UNames") >> fieldNames_;
+        }
+        else if (coeffs_.found("U"))
+        {
+            word UName(coeffs_.lookup("U"));
+            fieldNames_ = wordList(1, UName);
+        }
+        else
+        {
+            fieldNames_ = wordList(1, "U");
+        }
+
+        applied_.setSize(fieldNames_.size(), false);
+
+        // Create the Mangroves models - 1 per region
+        const dictionary& regionsDict(coeffs_.subDict("regions"));
+        const wordList regionNames(regionsDict.toc());
+        aZone_.setSize(regionNames.size(), 1);
+        NZone_.setSize(regionNames.size(), 1);
+        CmZone_.setSize(regionNames.size(), 1);
+        CdZone_.setSize(regionNames.size(), 1);
+        zoneIDs_.setSize(regionNames.size());
+
+        forAll(zoneIDs_, i)
+        {
+            const word& regionName = regionNames[i];
+            const dictionary& modelDict = regionsDict.subDict(regionName);
+
+            const word& zoneName = modelDict.lookup("cellZone");
+            zoneIDs_[i] = mesh_.cellZones().findIndices(zoneName);
+            if (zoneIDs_[i].empty())
+            {
+                FatalErrorInFunction
+                    << "Unable to find cellZone " << zoneName << nl
+                    << "Valid cellZones are:" << mesh_.cellZones().names()
+                    << exit(FatalError);
+            }
+
+            modelDict.lookup("a") >> aZone_[i];
+            modelDict.lookup("N") >> NZone_[i];
+            modelDict.lookup("Cm") >> CmZone_[i];
+            modelDict.lookup("Cd") >> CdZone_[i];
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H
new file mode 100644
index 0000000000000000000000000000000000000000..0e8ec5b03a176e56bd83ab6ad8ae8ce2456a6ae6
--- /dev/null
+++ b/src/waveModels/fvOptions/multiphaseMangrovesSource/multiphaseMangrovesSource.H
@@ -0,0 +1,154 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
+-------------------------------------------------------------------------------
+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::multiphaseMangrovesSource
+
+Group
+    grpFvOptionsSources
+
+Description
+
+Usage
+
+SourceFiles
+    multiphaseMangrovesSource.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMangrovesSource_H
+#define multiphaseMangrovesSource_H
+
+#include "fvOption.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class MangrovesModel;
+
+namespace fv
+{
+
+
+/*---------------------------------------------------------------------------*\
+                   Class multiphaseMangrovesSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMangrovesSource
+:
+    public option
+{
+
+protected:
+
+    // Protected data
+
+        // Coefficients per cell zone
+
+            //- Width of the vegetation element
+            scalarList aZone_;
+
+            //- Number of plants per unit of area
+            scalarList NZone_;
+
+            //- Inertia coefficient
+            scalarList CmZone_;
+
+            //- Drag coefficient
+            scalarList CdZone_;
+
+            //- Porosity coefficient
+            //scalarList MangrovesZone_;
+
+            //- Zone indices
+            labelListList zoneIDs_;
+
+        // Field properties
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multiphaseMangrovesSource(const multiphaseMangrovesSource&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const multiphaseMangrovesSource&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseMangrovesSource");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMangrovesSource
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMangrovesSource()
+    {}
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Add implicit contribution to phase momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<vector>& eqn,
+                const label fieldi
+            );
+
+
+        // IO
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..ddbfce77fe83bc08cf79c83ad2c1275373a3024c
--- /dev/null
+++ b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.C
@@ -0,0 +1,335 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
+-------------------------------------------------------------------------------
+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 "multiphaseMangrovesTurbulenceModel.H"
+#include "mathematicalConstants.H"
+#include "fvMesh.H"
+#include "fvMatrices.H"
+#include "fvmDdt.H"
+#include "fvcGrad.H"
+#include "fvmDiv.H"
+#include "fvmLaplacian.H"
+#include "fvmSup.H"
+#include "surfaceInterpolate.H"
+#include "surfaceFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace fv
+{
+    defineTypeNameAndDebug(multiphaseMangrovesTurbulenceModel, 0);
+    addToRunTimeSelectionTable
+    (
+        option,
+        multiphaseMangrovesTurbulenceModel,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::fv::multiphaseMangrovesTurbulenceModel::multiphaseMangrovesTurbulenceModel
+(
+    const word& name,
+    const word& modelType,
+    const dictionary& dict,
+    const fvMesh& mesh
+)
+:
+    option(name, modelType, dict, mesh),
+    aZone_(),
+    NZone_(),
+    CkpZone_(),
+    CepZone_(),
+    CdZone_(),
+    UName_("U"),
+    kName_("k"),
+    epsilonName_("epsilon")
+{
+    read(dict);
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
+(
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+
+    // Set alpha as zero-gradient
+    const volScalarField& epsilon =
+        mesh_.lookupObject<volScalarField>(epsilonName_);
+    const volScalarField& k = mesh_.lookupObject<volScalarField>(kName_);
+
+    volScalarField epsilonMangroves
+    (
+        IOobject
+        (
+            typeName + ":epsilonMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("epsilonMangroves", dimMass/dimMass/dimTime, 0)
+    );
+
+    volScalarField kMangroves
+    (
+        IOobject
+        (
+            typeName + ":kMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("kMangroves", dimTime/dimTime/dimTime, 0)
+    );
+
+    forAll(zoneIDs_, i)
+    {
+        const labelList& zones = zoneIDs_[i];
+
+        forAll(zones, j)
+        {
+            const label zonei = zones[j];
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            forAll(cz, kk)
+            {
+                const label celli = cz[kk];
+                const scalar a = aZone_[i];
+                const scalar N = NZone_[i];
+                const scalar Ckp = CkpZone_[i];
+                const scalar Cep = CepZone_[i];
+                const scalar Cd = CdZone_[i];
+
+                const vector& Uc = U[celli];
+
+                const scalar f = Cd*a*N*mag(Uc);
+                epsilonMangroves[celli] = Cep*f;
+                kMangroves[celli] = Ckp*f;
+            }
+        }
+    }
+
+    if (eqn.psi().name() == epsilonName_)
+    {
+    	epsilonMangroves.correctBoundaryConditions();
+	    fvMatrix<scalar> epsilonEqn
+	    (
+      	    - fvm::Sp(epsilonMangroves, epsilon)
+    	);
+	    eqn += epsilonEqn;
+    }
+    else if (eqn.psi().name() == kName_)
+    {
+    	kMangroves.correctBoundaryConditions();
+	    fvMatrix<scalar> kEqn
+	    (
+      	    - fvm::Sp(kMangroves, k)
+    	);
+	    eqn += kEqn;
+    }
+}
+
+
+bool Foam::fv::multiphaseMangrovesTurbulenceModel::read(const dictionary& dict)
+{
+    if (option::read(dict))
+    {
+        if (coeffs_.found("epsilonNames"))
+        {
+            coeffs_.lookup("epsilonNames") >> fieldNames_;
+        }
+        else if (coeffs_.found("epsilon"))
+        {
+            word UName(coeffs_.lookup("epsilon"));
+            fieldNames_ = wordList(1, UName);
+        }
+        else
+        {
+	        fieldNames_.setSize(2);
+	        fieldNames_[0] = "epsilon";
+    	    fieldNames_[1] = "k";
+        }
+
+        applied_.setSize(fieldNames_.size(), false);
+
+        // Create the Mangroves models - 1 per region
+        const dictionary& regionsDict(coeffs_.subDict("regions"));
+        const wordList regionNames(regionsDict.toc());
+        aZone_.setSize(regionNames.size(), 1);
+        NZone_.setSize(regionNames.size(), 1);
+        CkpZone_.setSize(regionNames.size(), 1);
+        CepZone_.setSize(regionNames.size(), 1);
+        CdZone_.setSize(regionNames.size(), 1);
+        zoneIDs_.setSize(regionNames.size());
+
+        forAll(zoneIDs_, i)
+        {
+            const word& regionName = regionNames[i];
+            const dictionary& modelDict = regionsDict.subDict(regionName);
+
+            const word& zoneName = modelDict.lookup("cellZone");
+            zoneIDs_[i] = mesh_.cellZones().findIndices(zoneName);
+            if (zoneIDs_[i].empty())
+            {
+                FatalErrorInFunction
+                    << "Unable to find cellZone " << zoneName << nl
+                    << "Valid cellZones are:" << mesh_.cellZones().names()
+                    << exit(FatalError);
+            }
+
+            modelDict.lookup("a") >> aZone_[i];
+            modelDict.lookup("N") >> NZone_[i];
+            modelDict.lookup("Ckp") >> CkpZone_[i];
+            modelDict.lookup("Cep") >> CepZone_[i];
+            modelDict.lookup("Cd") >> CdZone_[i];
+        }
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
+(
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    const volVectorField& U = mesh_.lookupObject<volVectorField>(UName_);
+
+    // Set alpha as zero-gradient
+    const volScalarField& epsilon =
+        mesh_.lookupObject<volScalarField>(epsilonName_);
+    const volScalarField& k = mesh_.lookupObject<volScalarField>(kName_);
+
+    volScalarField epsilonMangroves
+    (
+        IOobject
+        (
+            typeName + ":epsilonMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("epsilonMangroves", dimMass/dimMass/dimTime, 0)
+    );
+
+    volScalarField kMangroves
+    (
+        IOobject
+        (
+            typeName + ":kMangroves",
+            mesh_.time().timeName(),
+            mesh_.time(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("kMangroves", dimTime/dimTime/dimTime, 0)
+    );
+
+    forAll(zoneIDs_, i)
+    {
+        const labelList& zones = zoneIDs_[i];
+
+        forAll(zones, j)
+        {
+            const label zonei = zones[j];
+            const cellZone& cz = mesh_.cellZones()[zonei];
+
+            forAll(cz, kk)
+            {
+                const label celli = cz[kk];
+                const scalar a = aZone_[i];
+                const scalar N = NZone_[i];
+                const scalar Ckp = CkpZone_[i];
+                const scalar Cep = CepZone_[i];
+                const scalar Cd = CdZone_[i];
+
+                const vector& Uc = U[celli];
+
+    epsilonMangroves[celli] = Cep * Cd * a * N * sqrt ( Uc.x()*Uc.x() + Uc.y()*Uc.y() + Uc.z()*Uc.z());
+    kMangroves[celli] = Ckp * Cd * a * N * sqrt ( Uc.x()*Uc.x() + Uc.y()*Uc.y() + Uc.z()*Uc.z());
+            }
+        }
+    }
+
+    if (eqn.psi().name() == "epsilon")
+    {
+	Info<< " applying source to epsilon stuff.... " << endl;
+    	epsilonMangroves.correctBoundaryConditions();
+	fvMatrix<scalar> epsilonEqn
+	(
+      	    - fvm::Sp(rho*epsilonMangroves, epsilon)
+    	);
+	eqn += epsilonEqn;
+    }
+    else if (eqn.psi().name() == "k")
+    {
+	Info<< " applying source to k stuff.... " << endl;
+    	kMangroves.correctBoundaryConditions();
+	fvMatrix<scalar> kEqn
+	(
+      	    - fvm::Sp(rho*kMangroves, k)
+    	);
+	eqn += kEqn;
+    }
+
+}
+
+void Foam::fv::multiphaseMangrovesTurbulenceModel::addSup
+(
+    const volScalarField& alpha,
+    const volScalarField& rho,
+    fvMatrix<scalar>& eqn,
+    const label fieldi
+)
+{
+    NotImplemented;
+}
+
+// ************************************************************************* //
diff --git a/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..713106218316f0fae4a9a5e43a55b259b3b0e257
--- /dev/null
+++ b/src/waveModels/fvOptions/multiphaseMangrovesTurbulenceModel/multiphaseMangrovesTurbulenceModel.H
@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2017 OpenCFD Ltd.
+     \\/     M anipulation  | Copyright (C) 2017 IH-Cantabria
+-------------------------------------------------------------------------------
+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::multiphaseMangrovesTurbulenceModel
+
+Group
+    grpFvOptionsSources
+
+Description
+
+Usage
+
+SourceFiles
+    multiphaseMangrovesTurbulenceModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef multiphaseMangrovesTurbulenceModel_H
+#define multiphaseMangrovesTurbulenceModel_H
+
+#include "fvOption.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+namespace fv
+{
+
+
+/*---------------------------------------------------------------------------*\
+              Class multiphaseMangrovesTurbulenceModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class multiphaseMangrovesTurbulenceModel
+:
+    public option
+{
+
+protected:
+
+    // Protected data
+
+        // Coefficients per cell zone
+
+            //- Width of the vegetation element
+            scalarList aZone_;
+
+            //- Number of plants per unit of area
+            scalarList NZone_;
+
+            //- Ckp
+            scalarList CkpZone_;
+
+            //- Cep
+            scalarList CepZone_;
+
+            //- Drag coefficient
+            scalarList CdZone_;
+
+            //- Zone indices
+            labelListList zoneIDs_;
+
+
+        // Field properties
+
+            //- Name of U; default = U
+            word UName_;
+
+            //- Name of k; default = k
+            word kName_;
+
+            //- Name of epsilon; default = epsilon
+            word epsilonName_;
+
+
+private:
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        multiphaseMangrovesTurbulenceModel
+        (
+            const multiphaseMangrovesTurbulenceModel&
+        );
+
+        //- Disallow default bitwise assignment
+        void operator=(const multiphaseMangrovesTurbulenceModel&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("multiphaseMangrovesTurbulenceModel");
+
+
+    // Constructors
+
+        //- Construct from components
+        multiphaseMangrovesTurbulenceModel
+        (
+            const word& name,
+            const word& modelType,
+            const dictionary& dict,
+            const fvMesh& mesh
+        );
+
+
+    //- Destructor
+    virtual ~multiphaseMangrovesTurbulenceModel()
+    {}
+
+
+    // Member Functions
+
+        // Add explicit and implicit contributions
+
+            //- Add implicit contribution to momentum equation
+            virtual void addSup
+            (
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            );
+
+            //- Add implicit contribution to compressible momentum equation
+            virtual void addSup
+            (
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            );
+
+            //- Add implicit contribution to phase momentum equation
+            virtual void addSup
+            (
+                const volScalarField& alpha,
+                const volScalarField& rho,
+                fvMatrix<scalar>& eqn,
+                const label fieldi
+            );
+
+        // IO
+
+            //- Read dictionary
+            virtual bool read(const dictionary& dict);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fv
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //