diff --git a/applications/solvers/faParkerFukushimaFoam/createFaFields.H b/applications/solvers/faParkerFukushimaFoam/createFaFields.H
index 15a5872e67cae1e4af767e8c2ac73621beb6096c..dc407ab3a198b70c7a9727b2f66b0d13bc22e849 100644
--- a/applications/solvers/faParkerFukushimaFoam/createFaFields.H
+++ b/applications/solvers/faParkerFukushimaFoam/createFaFields.H
@@ -111,11 +111,24 @@
         g - gn*n
     );
 
-    areaScalarField ew
+    areaScalarField geff
     (
         IOobject
         (
-            "ew",
+            "geff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        gn
+    );
+
+    areaScalarField Sd
+    (
+        IOobject
+        (
+            "Sd",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -125,11 +138,11 @@
         dimensionedScalar(dimVelocity)
     );
 
-    areaScalarField Sd
+    areaScalarField Sa
     (
         IOobject
         (
-            "Sd",
+            "Sa",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -168,18 +181,3 @@
         dimensionedScalar("one", dimless, 1)
     );
 
-    const faPatchList& patches = aMesh.boundary();
-
-    const labelList& owner = aMesh.edgeOwner();
-
-    forAll(Us.boundaryField(), patchi)
-    {
-        if (Us.boundaryField()[patchi].type() != "processor" &&
-            Us.boundaryField()[patchi].fixesValue()){
-
-            for(label faceI = patches[patchi].start(); faceI < patches[patchi].start()+patches[patchi].size(); faceI++)
-            {
-                boundaryCell[owner[faceI]] = 0;
-            }
-        }
-    }
diff --git a/applications/solvers/faParkerFukushimaFoam/createSubmodels.H b/applications/solvers/faParkerFukushimaFoam/createSubmodels.H
index ad16a59ec5e241e25cf09555c7ce5ef9a1441fe5..a19d5a0238ca6279198de7d52eced314cde9c912 100644
--- a/applications/solvers/faParkerFukushimaFoam/createSubmodels.H
+++ b/applications/solvers/faParkerFukushimaFoam/createSubmodels.H
@@ -4,7 +4,40 @@ autoPtr<suspensionEntrainmentModel> entrainment
     suspensionEntrainmentModel::New(transportProperties, Us, h, hentrain, c, tau)
 );
 
+
+autoPtr<ambientEntrainmentModel> ambientEntrainment
+(
+    ambientEntrainmentModel::New(transportProperties, Us, h, c)
+);
+
+
 autoPtr<suspensionFrictionModel> friction
 (
     suspensionFrictionModel::New(transportProperties, Us, h, c)
 );
+
+
+autoPtr<suspensionDepositionModel> sedimentation
+(
+    suspensionDepositionModel::New(transportProperties, Us, h, c, tau)
+);
+
+
+if (relaxBoundaries)
+{
+    const faPatchList& patches = aMesh.boundary();
+
+    const labelList& owner = aMesh.edgeOwner();
+
+    forAll(Us.boundaryField(), patchi)
+    {
+        if (Us.boundaryField()[patchi].type() != "processor" &&
+            Us.boundaryField()[patchi].fixesValue()){
+
+            for(label faceI = patches[patchi].start(); faceI < patches[patchi].start()+patches[patchi].size(); faceI++)
+            {
+                boundaryCell[owner[faceI]] = 0;
+            }
+        }
+    }
+}
diff --git a/applications/solvers/faParkerFukushimaFoam/faParkerFukushimaFoam.C b/applications/solvers/faParkerFukushimaFoam/faParkerFukushimaFoam.C
index efc5d6d16922aa874fa2c6373846b10a3e1c4807..ce7ff6dbac5940404e00dcaa5d8ee2155ee46675 100644
--- a/applications/solvers/faParkerFukushimaFoam/faParkerFukushimaFoam.C
+++ b/applications/solvers/faParkerFukushimaFoam/faParkerFukushimaFoam.C
@@ -44,6 +44,8 @@ Author
 #include "SolverPerformance.H"
 #include "suspensionFrictionModel.H"
 #include "suspensionEntrainmentModel.H"
+#include "suspensionDepositionModel.H"
+#include "ambientEntrainmentModel.H"
 
 #define CREATE_FIELDS createFaFields.H
 
@@ -116,6 +118,7 @@ int main(int argc, char *argv[])
 
             Info<< "Correction iteration " << iCorr << endl;
 
+            geff = max(gn - (fac::ndiv(phis, Us)&n), dimensionedScalar(dimVelocity/dimTime));
 
             const areaVectorField & tauSc = friction->tauSc();
             const areaScalarField & tauSp = friction->tauSp();
@@ -129,9 +132,16 @@ int main(int argc, char *argv[])
               + fam::Sp(tauSp, Us)  // ustar**2 impl. part
               ==
                 R*gs*c*h
-              - 1./2.*R*gn*fac::grad(c*h*h)*boundaryCell
+              - 1./2.*R*geff*fac::grad(c*h*h)*boundaryCell
             );
 
+            if (correctMomentum)
+            {
+                UsEqn += R*c*fam::ddt(h, Us)
+                       + R*h*fam::ddt(c, Us)
+                       + R*fam::div(fac::interpolate(c)*phi2s, Us);
+            }
+
             if (!final)
                 UsEqn.relax();
 
@@ -146,22 +156,29 @@ int main(int argc, char *argv[])
 
 
             //Calculate water entrainment
-            const areaScalarField Ri(R*gn*c*h/(magSqr(Us)+sqr(u0))); //Parker et al. (1986) Eq. (10)
-            ew = 0.00153/(0.0204+Ri)*mag(Us); //Parker et al. (1986) Eq. (19)
-            ew = ew*pos(hwem-h); //clip entrainment if maximum water height is reached
+            Sa = ambientEntrainment->Sm();
+            Sa = Sa*pos(hwem-h); //clip entrainment if maximum water height is reached
 
             //Additional sink term for water
-            const areaScalarField vs(R*gn*Ds*Ds/18./nu); //Terminal velocity
+            const areaScalarField vs(R*geff*Ds*Ds/18./nu); //Terminal velocity
             areaScalarField dw(pos(h-hmin*2.)*vs);
 
 
+            //Calculate sediment entrainment
+            areaScalarField Sm = entrainment->Sm();
+            //no entrainment if there is not a minimal flow thickness
+            Sm = Sm*pos(h-hentmin);
+            Sm = min(Sm, hentrain/runTime.deltaT());
+
+
             //Parker et al. (1986), Eq. (3)
             faScalarMatrix hEqn
             (
                 fam::ddt(h)
               + fam::div(phis, h)
               ==
-                ew
+                Sa
+//              + Sm/(hentrain+dimensionedScalar("hentrainmin", dimLength, 1e-5))*hentrain
             );
 
             if (waterSink)
@@ -183,29 +200,18 @@ int main(int argc, char *argv[])
                 h = max(h, hmin);
             }
 
-
-            //Calculate sediment entrainment
-            areaScalarField Sm = entrainment->Sm();
-
-            //no entrainment if there is not a minimal flow thickness
-            Sm = Sm*pos(h-hentmin);
-
-            Sm = min(Sm, hentrain/runTime.deltaT());
-
-            //Calculate sedimation rate Sd following Parker et al. 1986
-            const areaScalarField mu(sqrt(mag(tau))/vs); //Parker et al. (1986), Eq. (21)
-            const areaScalarField r0(1.+31.5*pow(mu+SMALL, -1.46)); //Parker et al. (1986), Eq. (20)
-            //See Parker et al. (1986), Eq. (4), last term
-            Sd = vs*r0*c;
+            //Calculate sedimation rate
+            Sd = sedimentation->Sd();
             Sd = min(Sd, h*c/runTime.deltaT());
 
+
             faScalarMatrix cEqn
             (
                 fam::ddt(h, c)
               + fam::div(phi2s, c)
              ==
                 Sm/(hentrain+dimensionedScalar("hentrainmin", dimLength, 1e-5))*hentrain
-              - fam::Sp(Sd/(c+dimensionedScalar("cmin", dimless, 1e-7)), c)
+              - fam::Sp(Sd/(c+dimensionedScalar("cmin", dimless, 1e-8)), c)
             );
 
             if(!final)
@@ -222,7 +228,7 @@ int main(int argc, char *argv[])
                 fam::ddt(hentrain)
              ==
               - fam::Sp(Sm/(hentrain+dimensionedScalar("heintrainmin", dimLength, 1e-5)), hentrain)
-              + Sd/(c+dimensionedScalar("cmin", dimless, 1e-7))*c
+              + Sd/(c+dimensionedScalar("cmin", dimless, 1e-8))*c
             };
 
             hentrainEqn.solve();
diff --git a/applications/solvers/faParkerFukushimaFoam/readTransportProperties.H b/applications/solvers/faParkerFukushimaFoam/readTransportProperties.H
index 96642efeaf02ea99ad1db49bf04ecd9f9ff1780c..5527aeed2240697b9ba18e0986abebdfc35d23bc 100644
--- a/applications/solvers/faParkerFukushimaFoam/readTransportProperties.H
+++ b/applications/solvers/faParkerFukushimaFoam/readTransportProperties.H
@@ -33,6 +33,9 @@ Switch waterSink(transportProperties.getOrDefault<Switch>("waterSink", false));
 
 Switch bindHeight(transportProperties.getOrDefault<Switch>("bindHeight", true));
 
+Switch correctMomentum(transportProperties.getOrDefault<Switch>("correctMomentum", true));
+
+Switch relaxBoundaries(transportProperties.getOrDefault<Switch>("relaxBoundaries", false));
 
 Info << "Running with" << endl
      << "    hwem " << hwem << endl
@@ -41,5 +44,7 @@ Info << "Running with" << endl
      << "    h0   " << h0 << endl
      << "    u0   " << u0 << endl
      << "    waterSink is " << waterSink << endl
-     << "    bindHeight is " << bindHeight << endl << endl;
+     << "    bindHeight is " << bindHeight << endl
+     << "    correctMomentum is " << correctMomentum << endl
+     << "    relaxBoundaries is " << relaxBoundaries << endl << endl;
 
diff --git a/src/avalanche/Make/files b/src/avalanche/Make/files
index e070c23b69cdaa4d135073b4beb0ac65cff79f19..52f8cdddcc3f1edd0451f77a5ba8462a9f3fc961 100644
--- a/src/avalanche/Make/files
+++ b/src/avalanche/Make/files
@@ -1,21 +1,29 @@
 deposition = deposition
 $(deposition)/depositionModel/depositionModel.C
 $(deposition)/depositionModel/depositionModelNew.C
+$(deposition)/suspensionDepositionModel/suspensionDepositionModel.C
+$(deposition)/suspensionDepositionModel/suspensionDepositionModelNew.C
 $(deposition)/depositionOff/depositionOff.C
 $(deposition)/Stoppingprofile/Stoppingprofile.C
+$(deposition)/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.C
 
 entrainment = entrainment
 $(entrainment)/entrainmentModel/entrainmentModel.C
 $(entrainment)/entrainmentModel/entrainmentModelNew.C
 $(entrainment)/suspensionEntrainmentModel/suspensionEntrainmentModel.C
 $(entrainment)/suspensionEntrainmentModel/suspensionEntrainmentModelNew.C
+$(entrainment)/ambientEntrainmentModel/ambientEntrainmentModel.C
+$(entrainment)/ambientEntrainmentModel/ambientEntrainmentModelNew.C
 $(entrainment)/entrainmentOff/entrainmentOff.C
 $(entrainment)/Erosionenergy/Erosionenergy.C
 $(entrainment)/Front/Front.C
 $(entrainment)/Medina/Medina.C
 $(entrainment)/Ramms/Ramms.C
-$(entrainment)/suspensionParkerFukushimaEntrainment/suspensionParkerFukushimaEntrainment.C \
-
+$(entrainment)/suspensionParkerFukushimaEntrainment/suspensionParkerFukushimaEntrainment.C
+$(entrainment)/suspensionEntrainmentOff/suspensionEntrainmentOff.C
+$(entrainment)/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.C
+$(entrainment)/ambientTurnerEntrainment/ambientTurnerEntrainment.C
+$(entrainment)/ambientAnceyEntrainment/ambientAnceyEntrainment.C
 
 friction = friction
 $(friction)/frictionModel/frictionModel.C
diff --git a/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.C b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..bb767e5c2182fe615b9b391fdb7d207ecacdda5f
--- /dev/null
+++ b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.C
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "suspensionDepositionModel.H"
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(suspensionDepositionModel, 0);
+    defineRunTimeSelectionTable(suspensionDepositionModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::suspensionDepositionModel::readDict
+(
+    const word& type,
+    const dictionary& dict
+)
+{
+    suspensionDepositionProperties_ = dict;
+    coeffDict_ = suspensionDepositionProperties_.optionalSubDict(type + "Coeffs");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::suspensionDepositionModel::suspensionDepositionModel
+(
+    const word& type,
+    const dictionary& suspensionDepositionProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c,
+    const areaVectorField& tau
+)
+:
+    suspensionDepositionProperties_(suspensionDepositionProperties),
+    coeffDict_
+    (
+        suspensionDepositionProperties_.optionalSubDict(type + "Coeffs")
+    ),
+    R_("R", dimless, suspensionDepositionProperties_),
+    Us_(Us),
+    h_(h),
+    c_(c),
+    tau_(tau),
+    Sd_
+    (
+        IOobject
+        (
+            "Sd",
+            Us_.time().timeName(),
+            Us_.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        Us_.mesh(),
+        dimensionedScalar(dimVelocity)
+    )
+{
+    Info<< "    with " << nl
+        << "    " << R_ << endl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.H b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..ecbc7640c36802d3a4d79f8efbd4492d2da663f5
--- /dev/null
+++ b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModel.H
@@ -0,0 +1,200 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::suspensionDepositionModels
+
+Description
+    A namespace for various suspensionDeposition model implementations.
+
+Class
+    Foam::suspensionDepositionModel
+
+Description
+    An abstract base class for suspension deposition models
+
+SourceFiles
+    suspensionDepositionModel.C
+    suspensionDepositionModelNew.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef suspensionDepositionModel_H
+#define suspensionDepositionModel_H
+
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+#include "areaFieldsFwd.H"
+#include "FieldFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class suspensionDepositionModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class suspensionDepositionModel
+{
+protected:
+
+    // Protected data
+
+        dictionary suspensionDepositionProperties_;
+
+        //- Model coefficients dictionary
+        dictionary coeffDict_;
+
+        //- The density ratio
+        dimensionedScalar R_;
+
+        //- Reference to the surface velocity field
+        const areaVectorField& Us_;
+
+        //- Reference to the flow height field
+        const areaScalarField& h_;
+
+        //- Reference to the sediment concentration field
+        const areaScalarField& c_;
+
+        //- Reference to the bottom stress field
+        const areaVectorField& tau_;
+
+        //- Sink term
+        mutable areaScalarField Sd_;
+
+
+    // Protected Member Functions
+
+        //- Read/update the suspensionDepositionProperties and coeffDict dictionaries
+        void readDict(const word& type, const dictionary& dict);
+
+
+        //- Disallow copy construct
+        suspensionDepositionModel(const suspensionDepositionModel&) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const suspensionDepositionModel&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("suspensionDepositionModel");
+
+
+    // Declare run-time constructor selection table
+
+#ifndef SWIG
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            suspensionDepositionModel,
+            dictionary,
+            (
+                const dictionary& suspensionDepositionProperties,
+                const areaVectorField& Us,
+                const areaScalarField& h,
+                const areaScalarField& c,
+                const areaVectorField& tau
+            ),
+            (suspensionDepositionProperties, Us, h, c, tau)
+        );
+#endif
+
+
+    // Selectors
+
+        //- Return a reference to the selected suspensionDeposition model
+        static autoPtr<suspensionDepositionModel> New
+        (
+            const dictionary& suspensionDepositionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c,
+            const areaVectorField& tau
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        suspensionDepositionModel
+        (
+            const word& type,
+            const dictionary& suspensionDepositionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c,
+            const areaVectorField& tau
+        );
+
+
+    //- Destructor
+    virtual ~suspensionDepositionModel() = default;
+
+
+    // Member Functions
+
+        //- Read suspensionDepositionProperties dictionary
+        virtual bool read(const dictionary& suspensionDepositionProperties) = 0;
+
+        //- Return the friction properties dictionary
+        const dictionary& suspensionDepositionProperties() const
+        {
+            return suspensionDepositionProperties_;
+        }
+
+        //- Const access to the model coefficients dictionary
+        virtual const dictionary& coeffDict() const
+        {
+            return coeffDict_;
+        }
+
+
+        //- Return Sink by suspensionDeposition
+        virtual const areaScalarField& Sd() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModelNew.C b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..d93c17219fdc06101c5d055bd044e2918bd38f64
--- /dev/null
+++ b/src/avalanche/deposition/suspensionDepositionModel/suspensionDepositionModelNew.C
@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "suspensionDepositionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::suspensionDepositionModel> Foam::suspensionDepositionModel::New
+(
+    const dictionary& dict,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c,
+    const areaVectorField& tau
+
+)
+{
+    const word modelName(dict.get<word>("suspensionDepositionModel"));
+
+    Info<< "Selecting suspension deposition model " << modelName << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelName);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "suspensionDepositionModel",
+            modelName,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<suspensionDepositionModel>
+    (
+        cstrIter()(dict, Us, h, c, tau)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.C b/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.C
new file mode 100644
index 0000000000000000000000000000000000000000..31b361b72815cbd5fb7ef9152531f7c6eb77c903
--- /dev/null
+++ b/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.C
@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "suspensionParkerFukushimaDeposition.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace suspensionDepositionModels
+{
+    defineTypeNameAndDebug(suspensionParkerFukushimaDeposition, 0);
+    addToRunTimeSelectionTable(suspensionDepositionModel, suspensionParkerFukushimaDeposition, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::suspensionDepositionModels::suspensionParkerFukushimaDeposition::suspensionParkerFukushimaDeposition
+(
+    const dictionary& depositionProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c,
+    const areaVectorField& tau
+)
+:
+    suspensionDepositionModel(type(), depositionProperties, Us, h, c, tau),
+    Ds_("Ds", coeffDict_),
+    nu_("nu", coeffDict_),
+    gs_(Us.db().lookupObject<areaVectorField>("gs")),
+    gn_(Us.db().lookupObject<areaScalarField>("gn")),
+    geff_(Us.db().lookupObject<areaScalarField>("geff"))
+{
+    Info<< "    " << Ds_ << nl
+        << "    " << nu_ << nl
+        << "    " << Ds_ << nl << endl;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::suspensionDepositionModels::suspensionParkerFukushimaDeposition::Sd() const
+{
+
+    const areaScalarField vs(R_*geff_*Ds_*Ds_/18./nu_); //Terminal velocity
+    const areaScalarField mu(sqrt(mag(tau_))/(vs+dimensionedScalar(dimVelocity, SMALL))); //Parker et al. (1986), Eq. (21)
+    const areaScalarField r0(1.+31.5*pow(mu+SMALL, -1.46)); //Parker et al. (1986), Eq. (20)
+    //See Parker et al. (1986), Eq. (4), last term
+    Sd_ = vs*r0*c_;
+
+    return Sd_;
+}
+
+
+bool Foam::suspensionDepositionModels::suspensionParkerFukushimaDeposition::read
+(
+    const dictionary& depositionProperties
+)
+{
+    readDict(type(), depositionProperties);
+
+    coeffDict_.readEntry("R", R_);
+    coeffDict_.readEntry("Ds", Ds_);
+    coeffDict_.readEntry("nu", nu_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.H b/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.H
new file mode 100644
index 0000000000000000000000000000000000000000..6456c7a58908b898a53f6f84fc559ae65fda6dc9
--- /dev/null
+++ b/src/avalanche/deposition/suspensionParkerFukushimaDeposition/suspensionParkerFukushimaDeposition.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2017 Matthias Rauter
+-------------------------------------------------------------------------------
+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::suspensionDepositionModels::suspensionParkerFukushimaDeposition
+
+Description
+    Calculat the sedimentation rate in a turbidity current  following
+    Parker et al., 1986
+
+SourceFiles
+    suspensionParkerFukushimaDeposition.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef suspensionParkerFukushimaDeposition_H
+#define suspensionParkerFukushimaDeposition_H
+
+#include "suspensionDepositionModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace suspensionDepositionModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class suspensionParkerFukushimaDeposition Declaration
+\*---------------------------------------------------------------------------*/
+
+class suspensionParkerFukushimaDeposition
+:
+    public suspensionDepositionModel
+{
+    // Private data
+
+        //- Form-factor for deposition curve
+        dimensionedScalar Ds_;
+
+        //- Form-factor for deposition curve
+        dimensionedScalar nu_;
+
+        //- Reference to the gravitation fields
+        const areaVectorField& gs_;
+        const areaScalarField& gn_;
+        const areaScalarField& geff_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ParkerFukushimaDeposition");
+
+
+    // Constructors
+
+        //- Construct from components
+        suspensionParkerFukushimaDeposition
+        (
+            const dictionary& frictionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c,
+            const areaVectorField& tau
+        );
+
+
+    //- Destructor
+    virtual ~suspensionParkerFukushimaDeposition() = default;
+
+
+    // Member Functions
+
+        //- Return Sink by deposition
+        virtual const areaScalarField& Sd() const;
+
+        //- Read depositionProperties dictionary
+        virtual bool read(const dictionary& depositionProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.C b/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.C
new file mode 100644
index 0000000000000000000000000000000000000000..1748e36228a5dbaaca9ab97954139eb78091e223
--- /dev/null
+++ b/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.C
@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "ambientAnceyEntrainment.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+    defineTypeNameAndDebug(ambientAnceyEntrainment, 0);
+    addToRunTimeSelectionTable(ambientEntrainmentModel, ambientAnceyEntrainment, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ambientEntrainmentModels::ambientAnceyEntrainment::ambientAnceyEntrainment
+(
+    const dictionary& entrainmentProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c
+)
+:
+    ambientEntrainmentModel(type(), entrainmentProperties, Us, h, c),
+    alpha1_("alpha1", dimless, coeffDict_),
+    alpha2_("alpha2", dimless, coeffDict_),
+    geff_(Us_.db().lookupObject<areaScalarField>("geff"))
+{
+    Info << "    " << alpha1_ << nl
+         << "    " << alpha2_ << nl << endl;
+
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::ambientEntrainmentModels::ambientAnceyEntrainment::Sm() const
+{
+    //Calculate water entrainment
+    dimensionedScalar u0(dimVelocity, 1e-5);
+    const areaScalarField Ri(R_*geff_*c_*h_/(magSqr(Us_)+sqr(u0)));
+    Sm_ = (exp(-alpha1_*Ri*Ri)*(1-pos(Ri-1))+
+           exp(-alpha1_)/(Ri+SMALL)*pos(Ri-1))
+          *mag(Us_)*alpha2_;
+
+    return Sm_;
+}
+
+
+bool Foam::ambientEntrainmentModels::ambientAnceyEntrainment::read
+(
+    const dictionary& entrainmentProperties
+)
+{
+    readDict(type(), entrainmentProperties);
+
+    coeffDict_.readEntry("R", R_);
+    coeffDict_.readEntry("alpha1", alpha1_);
+    coeffDict_.readEntry("alpha2", alpha2_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.H b/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.H
new file mode 100644
index 0000000000000000000000000000000000000000..200195ebbd81836c58cbd635ce88a17810405c1f
--- /dev/null
+++ b/src/avalanche/entrainment/ambientAnceyEntrainment/ambientAnceyEntrainment.H
@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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::ambientEntrainmentModels::ambientAnceyEntrainment
+
+Description
+    Sediment entrainment model following Ancey (2004)
+
+    Definition follows Issler et al. (2018), doi.org/10.1017/jog.2017.62
+    with an additional linear factor alpha_2
+
+
+SourceFiles
+    ambientAnceyEntrainment.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ambientAnceyEntrainment_H
+#define ambientAnceyEntrainment_H
+
+#include "ambientEntrainmentModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Erosionenergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class ambientAnceyEntrainment
+:
+    public ambientEntrainmentModel
+{
+        // Private data
+
+        //- Entrainment parameter
+        dimensionedScalar alpha1_;
+
+        //- Entrainment parameter
+        dimensionedScalar alpha2_;
+
+       //- Reference to the graviational acceleration field
+        const areaScalarField& geff_;
+
+public:
+
+    //- Runtime type information
+    TypeName("AnceyEntrainment");
+
+
+    // Constructors
+
+        //- Construct from components
+        ambientAnceyEntrainment
+        (
+            const dictionary& frictionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c
+        );
+
+
+    //- Destructor
+    virtual ~ambientAnceyEntrainment() = default;
+
+
+    // Member Functions
+
+        //- Return the Source by entrainment
+        virtual const areaScalarField& Sm() const;
+
+        //- Read entrainmentProperties dictionary
+        virtual bool read(const dictionary& entrainmentProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.C b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..8f95aea2317cda5f1baa63595924662f2aafedb3
--- /dev/null
+++ b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.C
@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "ambientEntrainmentModel.H"
+#include "fvCFD.H"
+#include "faCFD.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(ambientEntrainmentModel, 0);
+    defineRunTimeSelectionTable(ambientEntrainmentModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+void Foam::ambientEntrainmentModel::readDict
+(
+    const word& type,
+    const dictionary& dict
+)
+{
+    ambientEntrainmentProperties_ = dict;
+    coeffDict_ = ambientEntrainmentProperties_.optionalSubDict(type + "Coeffs");
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ambientEntrainmentModel::ambientEntrainmentModel
+(
+    const word& type,
+    const dictionary& ambientEntrainmentProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c
+)
+:
+    ambientEntrainmentProperties_(ambientEntrainmentProperties),
+    coeffDict_
+    (
+        ambientEntrainmentProperties_.optionalSubDict(type + "Coeffs")
+    ),
+    R_("R", dimless, ambientEntrainmentProperties_),
+    Us_(Us),
+    h_(h),
+    c_(c),
+    Sm_
+    (
+        IOobject
+        (
+            "Sm",
+            Us_.time().timeName(),
+            Us_.db(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        Us_.mesh(),
+        dimensionedScalar(dimVelocity)
+    )
+{
+    Info<< "    with " << nl
+        << "    " << R_  << nl;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.H b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..08b84cb5630c9c2e7f5be031a5b4271e45d1dc45
--- /dev/null
+++ b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModel.H
@@ -0,0 +1,193 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Namespace
+    Foam::ambientEntrainmentModels
+
+Description
+    A namespace for various entrainment model implementations for faParkerFukushimaFoam.
+
+Class
+    Foam::ambientEntrainmentModel
+
+Description
+    An abstract base class for entrainment models
+
+SourceFiles
+    ambientEntrainmentModel.C
+    ambientEntrainmentModelNew.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ambientEntrainmentModel_H
+#define ambientEntrainmentModel_H
+
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+#include "areaFieldsFwd.H"
+#include "FieldFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class entrainmentModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class ambientEntrainmentModel
+{
+protected:
+
+    // Protected data
+
+        dictionary ambientEntrainmentProperties_;
+
+        //- Model coefficients dictionary
+        dictionary coeffDict_;
+
+        //- Density ratio R = rhos/rhow
+        dimensionedScalar R_;
+
+        //- Reference to the surface velocity field
+        const areaVectorField& Us_;
+
+        //- Reference to the flow height field
+        const areaScalarField& h_;
+
+        //- Reference to the sediment concentration field
+        const areaScalarField& c_;
+
+        //- Source term
+        mutable areaScalarField Sm_;
+
+
+    // Protected Member Functions
+
+        //- Read/update the entrainmentProperties and coeffDict dictionaries
+        void readDict(const word& type, const dictionary& dict);
+
+
+        //- Disallow copy construct
+        ambientEntrainmentModel(const ambientEntrainmentModel&) = delete;
+
+        //- Disallow default bitwise assignment
+        void operator=(const ambientEntrainmentModel&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ambientEntrainmentModel");
+
+
+    // Declare run-time constructor selection table
+
+#ifndef SWIG
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            ambientEntrainmentModel,
+            dictionary,
+            (
+                const dictionary& entrainmentProperties,
+                const areaVectorField& Us,
+                const areaScalarField& h,
+                const areaScalarField& c
+            ),
+            (entrainmentProperties, Us, h, c)
+        );
+#endif
+
+
+    // Selectors
+
+        //- Return a reference to the selected entrainment model
+        static autoPtr<ambientEntrainmentModel> New
+        (
+            const dictionary& entrainmentProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        ambientEntrainmentModel
+        (
+            const word& type,
+            const dictionary& entrainmentProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c
+        );
+
+
+    //- Destructor
+    virtual ~ambientEntrainmentModel() = default;
+
+
+    // Member Functions
+
+        //- Read entrainmentProperties dictionary
+        virtual bool read(const dictionary& ambientEntrainmentProperties) = 0;
+
+        //- Return the entrainment properties dictionary
+        const dictionary& ambientEntrainmentProperties() const
+        {
+            return ambientEntrainmentProperties_;
+        }
+
+        //- Const access to the model coefficients dictionary
+        virtual const dictionary& coeffDict() const
+        {
+            return coeffDict_;
+        }
+
+        //- Return the Source by entrainment
+        virtual const areaScalarField& Sm() const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModelNew.C b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModelNew.C
new file mode 100644
index 0000000000000000000000000000000000000000..587ff1c36ae668fe6f8fa5dd41976288395fe7e3
--- /dev/null
+++ b/src/avalanche/entrainment/ambientEntrainmentModel/ambientEntrainmentModelNew.C
@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "ambientEntrainmentModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::ambientEntrainmentModel>
+Foam::ambientEntrainmentModel::New
+(
+    const dictionary& dict,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c
+)
+{
+    const word modelName(dict.get<word>("ambientEntrainmentModel"));
+
+    Info<< "Selecting ambient entrainment model " << modelName << endl;
+
+    auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelName);
+
+    if (!cstrIter.found())
+    {
+        FatalIOErrorInLookup
+        (
+            dict,
+            "ambientEntrainmentModel",
+            modelName,
+            *dictionaryConstructorTablePtr_
+        ) << exit(FatalIOError);
+    }
+
+    return autoPtr<ambientEntrainmentModel>
+    (
+        cstrIter()(dict, Us, h, c)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.C b/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.C
new file mode 100644
index 0000000000000000000000000000000000000000..7162b3c7360db08cbfe984445476331076fc6620
--- /dev/null
+++ b/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.C
@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "ambientParkerFukushimaEntrainment.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+    defineTypeNameAndDebug(ambientParkerFukushimaEntrainment, 0);
+    addToRunTimeSelectionTable(ambientEntrainmentModel, ambientParkerFukushimaEntrainment, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ambientEntrainmentModels::ambientParkerFukushimaEntrainment::ambientParkerFukushimaEntrainment
+(
+    const dictionary& entrainmentProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c
+)
+:
+    ambientEntrainmentModel(type(), entrainmentProperties, Us, h, c),
+    ewf_("ewf", dimless, coeffDict_),
+    Ri0_("Ri0", dimless, coeffDict_),
+    geff_(Us_.db().lookupObject<areaScalarField>("geff"))
+{
+    Info << "    " << ewf_ << nl
+         << "    " << Ri0_ << nl << endl;
+
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::ambientEntrainmentModels::ambientParkerFukushimaEntrainment::Sm() const
+{
+    dimensionedScalar u0(dimVelocity, 1e-5);
+    const areaScalarField Ri(R_*geff_*c_*h_/(magSqr(Us_)+sqr(u0))); //Parker et al. (1986) Eq. (10)
+    Sm_ = ewf_/(Ri0_+Ri)*mag(Us_); //Parker et al. (1986) Eq. (19)
+
+    return Sm_;
+}
+
+
+bool Foam::ambientEntrainmentModels::ambientParkerFukushimaEntrainment::read
+(
+    const dictionary& entrainmentProperties
+)
+{
+    readDict(type(), entrainmentProperties);
+
+    coeffDict_.readEntry("R", R_);
+    coeffDict_.readEntry("ewf", ewf_);
+    coeffDict_.readEntry("Ri0", Ri0_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.H b/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.H
new file mode 100644
index 0000000000000000000000000000000000000000..5dd9693f0b6ca0c41a863d6dafa9161d8bc0aacd
--- /dev/null
+++ b/src/avalanche/entrainment/ambientParkerFukushimaEntrainment/ambientParkerFukushimaEntrainment.H
@@ -0,0 +1,125 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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::ambientEntrainmentModels::ambientParkerFukushimaEntrainment
+
+Description
+    Sediment entrainment model following Parker et al. (1986)
+
+    Definition follows Parker et al. (1986), doi.org/10.1017/S0022112086001404
+
+
+SourceFiles
+    ambientParkerFukushimaEntrainment.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ambientParkerFukushimaEntrainment_H
+#define ambientParkerFukushimaEntrainment_H
+
+#include "ambientEntrainmentModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Erosionenergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class ambientParkerFukushimaEntrainment
+:
+    public ambientEntrainmentModel
+{
+    // Private data
+
+        //- Entrainment parameter
+        dimensionedScalar ewf_;
+
+        //- Entrainment parameter
+        dimensionedScalar Ri0_;
+
+       //- Reference to the graviational acceleration field
+        const areaScalarField& geff_;
+
+public:
+
+    //- Runtime type information
+    TypeName("ParkerFukushimaEntrainment");
+
+
+    // Constructors
+
+        //- Construct from components
+        ambientParkerFukushimaEntrainment
+        (
+            const dictionary& frictionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c
+        );
+
+
+    //- Destructor
+    virtual ~ambientParkerFukushimaEntrainment() = default;
+
+
+    // Member Functions
+
+        //- Return the Source by entrainment
+        virtual const areaScalarField& Sm() const;
+
+        //- Read entrainmentProperties dictionary
+        virtual bool read(const dictionary& entrainmentProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.C b/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.C
new file mode 100644
index 0000000000000000000000000000000000000000..61c9e69f2bd1a538ba619f16d7b6d7e6772d9343
--- /dev/null
+++ b/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.C
@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2020 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "ambientTurnerEntrainment.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+    defineTypeNameAndDebug(ambientTurnerEntrainment, 0);
+    addToRunTimeSelectionTable(ambientEntrainmentModel, ambientTurnerEntrainment, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::ambientEntrainmentModels::ambientTurnerEntrainment::ambientTurnerEntrainment
+(
+    const dictionary& entrainmentProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& c
+)
+:
+    ambientEntrainmentModel(type(), entrainmentProperties, Us, h, c),
+    alpha1_("alpha1", dimless, coeffDict_),
+    alpha2_("alpha2", dimless, coeffDict_),
+    Ri0_("Ri0", dimless, coeffDict_),
+    geff_(Us_.db().lookupObject<areaScalarField>("geff"))
+{
+    Info << "    " << R_  << nl
+         << "    " << alpha1_ << nl
+         << "    " << alpha2_ << nl
+         << "    " << Ri0_ << nl << endl;
+
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::ambientEntrainmentModels::ambientTurnerEntrainment::Sm() const
+{
+    //Calculate water entrainment
+    dimensionedScalar u0(dimVelocity, 1e-5);
+    const areaScalarField Ri(R_*geff_*c_*h_/(magSqr(Us_)+sqr(u0)));
+    Sm_ = (1-pos(Ri-Ri0_))*(Ri0_-Ri)/(alpha1_+alpha2_*Ri)*mag(Us_); //Turner (1986), Eq 23
+
+    return Sm_;
+}
+
+
+bool Foam::ambientEntrainmentModels::ambientTurnerEntrainment::read
+(
+    const dictionary& entrainmentProperties
+)
+{
+    readDict(type(), entrainmentProperties);
+
+    coeffDict_.readEntry("R", R_);
+    coeffDict_.readEntry("Ri0", Ri0_);
+    coeffDict_.readEntry("alpha1", alpha1_);
+    coeffDict_.readEntry("alpha2", alpha2_);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.H b/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.H
new file mode 100644
index 0000000000000000000000000000000000000000..62e49a619ad1d5a7391b1d90fca43d476fea9c93
--- /dev/null
+++ b/src/avalanche/entrainment/ambientTurnerEntrainment/ambientTurnerEntrainment.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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::ambientEntrainmentModels::ambientTurnerEntrainment
+
+Description
+    Sediment entrainment model following Turner (1986)
+
+    Definition follows Issler et al. (2018), doi.org/10.1017/jog.2017.62
+
+SourceFiles
+    ambientTurnerEntrainment.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ambientTurnerEntrainment_H
+#define ambientTurnerEntrainment_H
+
+#include "ambientEntrainmentModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace ambientEntrainmentModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Erosionenergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class ambientTurnerEntrainment
+:
+    public ambientEntrainmentModel
+{
+    // Private data
+
+
+        //- Entrainment parameter
+        dimensionedScalar alpha1_;
+
+        //- Entrainment parameter
+        dimensionedScalar alpha2_;
+
+        //- Entrainment parameter
+        dimensionedScalar Ri0_;
+
+       //- Reference to the graviational acceleration field
+        const areaScalarField& geff_;
+
+public:
+
+    //- Runtime type information
+    TypeName("TurnerEntrainment");
+
+
+    // Constructors
+
+        //- Construct from components
+        ambientTurnerEntrainment
+        (
+            const dictionary& frictionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& c
+        );
+
+
+    //- Destructor
+    virtual ~ambientTurnerEntrainment() = default;
+
+
+    // Member Functions
+
+        //- Return the Source by entrainment
+        virtual const areaScalarField& Sm() const;
+
+        //- Read entrainmentProperties dictionary
+        virtual bool read(const dictionary& entrainmentProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.C b/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.C
new file mode 100644
index 0000000000000000000000000000000000000000..426cd1a54d673ad6c0c7f8d513ae11ff9a30c0d1
--- /dev/null
+++ b/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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/>.
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "faCFD.H"
+#include "suspensionEntrainmentOff.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace suspensionEntrainmentModels
+{
+    defineTypeNameAndDebug(suspensionEntrainmentOff, 0);
+    addToRunTimeSelectionTable(suspensionEntrainmentModel, suspensionEntrainmentOff, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::suspensionEntrainmentModels::suspensionEntrainmentOff::suspensionEntrainmentOff
+(
+    const dictionary& entrainmentProperties,
+    const areaVectorField& Us,
+    const areaScalarField& h,
+    const areaScalarField& hentrain,
+    const areaScalarField& c,
+    const areaVectorField& tau
+)
+:
+    suspensionEntrainmentModel(type(), entrainmentProperties, Us, h, hentrain, c, tau)
+{
+    Info << endl;
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+const Foam::areaScalarField&
+Foam::suspensionEntrainmentModels::suspensionEntrainmentOff::Sm() const
+{
+    return Sm_;
+}
+
+
+bool Foam::suspensionEntrainmentModels::suspensionEntrainmentOff::read
+(
+    const dictionary& entrainmentProperties
+)
+{
+    readDict(type(), entrainmentProperties);
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.H b/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.H
new file mode 100644
index 0000000000000000000000000000000000000000..257c773d19641ea01bb7b3314279adc522d8eea5
--- /dev/null
+++ b/src/avalanche/entrainment/suspensionEntrainmentOff/suspensionEntrainmentOff.H
@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | avalanche module
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Matthias Rauter
+-------------------------------------------------------------------------------
+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::suspensionEntrainmentModels::suspensionEntrainmentOff
+
+Description
+    An dummy entrainment model for no entrainment.
+
+SourceFiles
+    suspensionEntrainmentOff.C
+
+Author
+    Matthias Rauter matthias@rauter.it
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef suspensionEntrainmentOff_H
+#define suspensionEntrainmentOff_H
+
+#include "suspensionEntrainmentModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+#include "IOdictionary.H"
+#include "typeInfo.H"
+#include "runTimeSelectionTables.H"
+#include "dimensionedScalar.H"
+#include "tmp.H"
+#include "autoPtr.H"
+#include "faMatrices.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace suspensionEntrainmentModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Erosionenergy Declaration
+\*---------------------------------------------------------------------------*/
+
+class suspensionEntrainmentOff
+:
+    public suspensionEntrainmentModel
+{
+    // Private data
+
+public:
+
+    //- Runtime type information
+    TypeName("entrainmentOff");
+
+
+    // Constructors
+
+        //- Construct from components
+        suspensionEntrainmentOff
+        (
+            const dictionary& frictionProperties,
+            const areaVectorField& Us,
+            const areaScalarField& h,
+            const areaScalarField& hentrain,
+            const areaScalarField& c,
+            const areaVectorField& tau
+        );
+
+
+    //- Destructor
+    virtual ~suspensionEntrainmentOff() = default;
+
+
+    // Member Functions
+
+        //- Return the Source by entrainment
+        virtual const areaScalarField& Sm() const;
+
+        //- Read entrainmentProperties dictionary
+        virtual bool read(const dictionary& entrainmentProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace frictionModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //