diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
index 88dce166..f7d3d217 100644
--- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
+++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C
@@ -92,7 +92,7 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        if(hasEvolved)
+        if(hasEvolved && smoothenForces)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
         }
diff --git a/applications/solvers/cfdemSolverRhoPimple/createFields.H b/applications/solvers/cfdemSolverRhoPimple/createFields.H
index 35364cbd..23ddbab4 100644
--- a/applications/solvers/cfdemSolverRhoPimple/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/createFields.H
@@ -177,6 +177,17 @@ Info<< "Reading thermophysical properties\n" << endl;
         )
     );
 
+    bool smoothenForces
+    (
+        pimple.dict().lookupOrDefault<bool>
+        (
+            "smoothenForces",
+            false
+        )
+    );
+    if (smoothenForces) Info << "Smoothening implicit particle forces.\n" << endl;
+    else Info << "Not smoothening implicit particle forces.\n" << endl;
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
index 8aa2dea3..45af8a5c 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.C
@@ -68,6 +68,18 @@ twoWayMPI::twoWayMPI
 :
     dataExchangeModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    DEMVariableNames_(propsDict_.lookupOrDefault<wordList>("DEMVariables",wordList(0))),
+    DEMVariables_
+    (
+        IOobject
+        (
+            "DEMVariables",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+         )
+    ),
     lmp(NULL)
 {
     Info << "Starting up LIGGGHTS for first time execution" << endl;
@@ -100,7 +112,30 @@ twoWayMPI::~twoWayMPI()
     delete lmp;
 }
 
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+void twoWayMPI::updateDEMVariables()
+{
+    scalar variablevalue = 0.0;
+    word variablename("");
+    DEMVariables_.clear();
+    for (label i=0; i<DEMVariableNames_.size(); i++)
+    {
+        variablename = DEMVariableNames_[i];
+        variablevalue = DEMVariableValue(variablename);
+        DEMVariables_.append(variablevalue);
+    }
+}
+
+double twoWayMPI::DEMVariableValue(word variablename)
+{
+    return liggghts_get_variable(lmp,variablename.c_str());
+}
+
+
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
+
 void twoWayMPI::getData
 (
     word name,
@@ -355,6 +390,9 @@ bool twoWayMPI::couple(int i)
         particleCloud_.clockM().start(4,"LIGGGHTS_reallocArrays");
         particleCloud_.reAllocArrays();
         particleCloud_.clockM().stop("LIGGGHTS_reallocArrays");
+
+        // retrieve DEM variables if present
+        updateDEMVariables();
     }
 
     return coupleNow;
diff --git a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
index 9bd66fd5..9801ec52 100644
--- a/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
+++ b/src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMPI/twoWayMPI.H
@@ -44,6 +44,7 @@ SourceFiles
 #include <mpi.h>
 
 #include "dataExchangeModel.H"
+#include "scalarIOList.H"
 
 //=================================//
 //LAMMPS/LIGGGHTS
@@ -71,7 +72,14 @@ private:
     dictionary propsDict_;
     MPI_Comm comm_liggghts;
 
+    wordList DEMVariableNames_;
+
+    scalarIOList DEMVariables_; 
+
   // private member functions
+    double DEMVariableValue(word);
+
+    void updateDEMVariables();
 
 protected:
     LAMMPS_NS::LAMMPS *lmp;
@@ -146,7 +154,6 @@ public:
         int getNumberOfTypes() const;
         double* getTypeVol() const;
 
-
         scalar getCG() const { return lmp->force->cg(); }
 };
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C
index 3e5e2ac3..e4eef410 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.C
@@ -218,8 +218,6 @@ void BeetstraDrag::setForce() const
                     voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
-                // in case a fines phase is present, void fraction needs to be adapted
-                adaptVoidfraction(voidfraction, cellI);
 
                 if (multiTypes_)
                 {
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H
index 21245471..d66b81f3 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDrag/BeetstraDrag.H
@@ -90,8 +90,6 @@ protected:
 
     bool usePC_;
 
-    virtual void adaptVoidfraction(double&, label) const {}
-
     virtual scalar effDiameter(double d, label cellI, label index) const {return d;}
 
     virtual scalar meanSauterDiameter(double d, label cellI) const {return d;}
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
index d0580349..ae6e05e7 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.C
@@ -59,19 +59,15 @@ BeetstraDragPoly::BeetstraDragPoly
     {
         dFine_ = readScalar(propsDict_.lookup("dFine"));
         volScalarField& alphaP(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaP")));
-     //   alphaP_.set(&alphaP);
         alphaP_ = &alphaP;
         volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
-     //   alphaSt_.set(&alphaSt);
         alphaSt_ = &alphaSt;
         volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
-//        dSauter_.set(&dSauter);
         dSauter_ = &dSauter;
     }
     else
     {
         volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
-  //      dSauter_.set(&dSauter);
         dSauter_ = &dSauter;
     }
 }
@@ -84,13 +80,6 @@ BeetstraDragPoly::~BeetstraDragPoly()
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) const
-{
-    if (fines_) voidfraction -= (*alphaSt_)[cellI];
-    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
-}
-
 scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
 {
     scalar dS = (*dSauter_)[cellI];
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
index 67096b8b..8c91b7cf 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/BeetstraDragPoly/BeetstraDragPoly.H
@@ -52,8 +52,6 @@ protected:
 
     volScalarField* dSauter_;
 
-    void adaptVoidfraction(double&, label) const;
-
     scalar effDiameter(double, label, label) const;
     
     scalar meanSauterDiameter(double, label) const;
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.C
index 5c4e6efa..ce779474 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/Fines.C
@@ -50,6 +50,11 @@ Fines::Fines
     forceModel(dict,sm),
     finesFields_(dict,sm)
 {
+    // safety check: fines should be first force model in list because they correct the voidfraction field
+    if (particleCloud_.forceModels()[0] != "dSauter" || particleCloud_.forceModels()[1] != "Fines")
+    {
+        FatalError <<"Force models dSauter and Fines need to be first in list.\n" << abort(FatalError);
+    }
 }
 
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
index 1147bbe6..eb8c536f 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.C
@@ -45,11 +45,15 @@ FinesFields::FinesFields
     propsDict_(dict.subDict("FinesFieldsProps")),
     smoothing_(propsDict_.lookupOrDefault<bool>("smoothing",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    clogKin_(propsDict_.lookupOrDefault<bool>("kineticClogging",false)),
+    clogStick_(propsDict_.lookupOrDefault<bool>("stickyClogging",false)),
+    movingBed_(propsDict_.lookupOrDefault<bool>("movingBed",true)),
+    fluxFieldName_(propsDict_.lookupOrDefault<word>("fluxFieldName","phi")),
+    phi_(sm.mesh().lookupObject<surfaceScalarField> (fluxFieldName_)),
     velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
-    // this is probably really bad
-    voidfraction_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_))),
+    voidfraction_(sm.mesh().lookupObjectRef<volScalarField> (voidfractionFieldName_)),
     UsFieldName_(propsDict_.lookupOrDefault<word>("granVelFieldName","Us")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     pFieldName_(propsDict_.lookupOrDefault<word>("pFieldName","p")),
@@ -200,7 +204,6 @@ FinesFields::FinesFields
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(1,0,-1,0,0), 0)
-        //dimensionedVector("zero", dimensionSet(1,-2,-1,0,0), vector::zero)
     ),
     uDyn_
     (   IOobject
@@ -219,54 +222,114 @@ FinesFields::FinesFields
     rhoFine_("rhoFine",dimensionSet(1,-3,0,0,0),0.0),
     g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81)),
     alphaDynMax_(0.1),
-    alphaMax_(readScalar(propsDict_.lookup ("alphaMax"))),
-    critVoidfraction_(readScalar(propsDict_.lookup ("critVoidfraction"))),
-    depRate_(readScalar(propsDict_.lookup ("depRate"))),
+    alphaMax_(propsDict_.lookupOrDefault<scalar>("alphaMax",0.95)),
+    alphaMinClog_(propsDict_.lookupOrDefault<scalar>("alphaMinClog",0.3)),
+    critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.05)),
+    deltaT_(voidfraction_.mesh().time().deltaTValue()),
+    depositionLength_(readScalar(propsDict_.lookup ("depositionLength"))),
     exponent_(-1.33),
-    nCrit_(readScalar(propsDict_.lookup ("nCrit"))),
-    poresizeWidth_(readScalar(propsDict_.lookup ("poresizeWidth"))),
+    nCrit_(0.0),
+    poresizeWidth_(0.0),
     prefactor_(10.5),
-    ratioHydraulicPore_(1.5)
+    ratioHydraulicPore_(1.5),
+    tauRelease_(readScalar(propsDict_.lookup ("tauRelease"))),
+    uBindHigh_(propsDict_.lookupOrDefault<scalar>("uBindHigh",0.0)),
+    uBindLow_(propsDict_.lookupOrDefault<scalar>("uBindLow",0.0)),
+    uMin_(0.001)
 {
     Sds_.write();
 
     if (propsDict_.found("prefactor"))
-        prefactor_=readScalar(propsDict_.lookup ("prefactor"));
-    if (propsDict_.found("exponent"))
-        exponent_=readScalar(propsDict_.lookup ("exponent"));
-    if (propsDict_.found("dFine"))
-        dFine_.value()=readScalar(propsDict_.lookup ("dFine"));
-    else
-        FatalError <<"Please specify dFine.\n" << abort(FatalError);
-    if (propsDict_.found("diffCoeff"))
-        diffCoeff_.value()=readScalar(propsDict_.lookup ("diffCoeff"));
-    if (propsDict_.found("rhoFine"))
-        rhoFine_.value()=readScalar(propsDict_.lookup ("rhoFine"));
-    else
-        FatalError <<"Please specify rhoFine.\n" << abort(FatalError);
-    if (propsDict_.found("nuAve"))
-        nuAve_.value()=readScalar(propsDict_.lookup ("nuAve"));
-    if (propsDict_.found("alphaDynMax"))
-        alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax"));
-
-    if(verbose_)
     {
-      alphaG_.writeOpt() = IOobject::AUTO_WRITE;
-      alphaG_.write();
-      alphaP_.writeOpt() = IOobject::AUTO_WRITE;
-      alphaP_.write();
-      deltaAlpha_.writeOpt() = IOobject::AUTO_WRITE;
-      deltaAlpha_.write();
-      dHydMix_.writeOpt() = IOobject::AUTO_WRITE;
-      dHydMix_.write();
-      DragCoeff_.writeOpt() = IOobject::AUTO_WRITE;
-      DragCoeff_.write();
-      dSauterMix_.writeOpt() = IOobject::AUTO_WRITE;
-      dSauterMix_.write();
-      FanningCoeff_.writeOpt() = IOobject::AUTO_WRITE;
-      FanningCoeff_.write();
-      Froude_.writeOpt() = IOobject::AUTO_WRITE;
-      Froude_.write();
+        prefactor_=readScalar(propsDict_.lookup ("prefactor"));
+    }
+    if (propsDict_.found("exponent"))
+    {
+        exponent_=readScalar(propsDict_.lookup ("exponent"));
+    }
+    if (propsDict_.found("dFine"))
+    {
+        dFine_.value()=readScalar(propsDict_.lookup ("dFine"));
+    }
+    else
+    {
+        FatalError <<"Please specify dFine.\n" << abort(FatalError);
+    }
+    if (propsDict_.found("diffCoeff"))
+    {
+        diffCoeff_.value()=readScalar(propsDict_.lookup ("diffCoeff"));
+    }
+    if (propsDict_.found("rhoFine"))
+    {
+        rhoFine_.value()=readScalar(propsDict_.lookup ("rhoFine"));
+    }
+    else
+    {
+        FatalError <<"Please specify rhoFine.\n" << abort(FatalError);
+    }
+    if (propsDict_.found("nuAve"))
+    {
+        nuAve_.value()=readScalar(propsDict_.lookup ("nuAve"));
+    }
+    if (propsDict_.found("alphaDynMax"))
+    {
+        alphaDynMax_=readScalar(propsDict_.lookup ("alphaDynMax"));
+    }
+
+    if (clogKin_)
+    {
+        nCrit_ = readScalar(propsDict_.lookup ("nCrit"));
+        poresizeWidth_ = readScalar(propsDict_.lookup ("poresizeWidth"));
+    }
+
+    if (clogStick_)
+    {
+        if (uBindHigh_ - uBindLow_ < SMALL)
+        {
+            FatalError <<"No reasonable values for critical binding velocities.\n" << abort(FatalError);
+        }
+
+        uReconstructed_.set
+        (
+            new volVectorField
+            (
+                IOobject
+                (
+                    "uReconstructed",
+                    sm.mesh().time().timeName(),
+                    sm.mesh(),
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                U_//sm.mesh()
+            )
+        );
+
+    }
+
+    if (verbose_)
+    {
+        alphaG_.writeOpt() = IOobject::AUTO_WRITE;
+        alphaG_.write();
+        alphaP_.writeOpt() = IOobject::AUTO_WRITE;
+        alphaP_.write();
+        deltaAlpha_.writeOpt() = IOobject::AUTO_WRITE;
+        deltaAlpha_.write();
+        dHydMix_.writeOpt() = IOobject::AUTO_WRITE;
+        dHydMix_.write();
+        DragCoeff_.writeOpt() = IOobject::AUTO_WRITE;
+        DragCoeff_.write();
+        dSauterMix_.writeOpt() = IOobject::AUTO_WRITE;
+        dSauterMix_.write();
+        FanningCoeff_.writeOpt() = IOobject::AUTO_WRITE;
+        FanningCoeff_.write();
+        Froude_.writeOpt() = IOobject::AUTO_WRITE;
+        Froude_.write();
+        if (clogStick_)
+        {
+            uReconstructed_().writeOpt() = IOobject::AUTO_WRITE;
+            uReconstructed_().write();
+        }
     }
 }
 
@@ -281,76 +344,93 @@ FinesFields::~FinesFields()
 
 void FinesFields::update()
 {
-    if(verbose_)  Info << "FinesFields: Updating alphaP.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating alphaP.\n" << endl;
     updateAlphaP();
-    if(verbose_)  Info << "FinesFields: Updating alphaG.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating alphaG.\n" << endl;
     updateAlphaG();
-    if(verbose_)  Info << "FinesFields: Calculating source terms.\n" << endl;
+    if (clogStick_)
+    {
+        if (verbose_)  Info << "FinesFields: Updating uReconstructed.\n" << endl;
+        updateUReconstructed();
+    }
+    if (verbose_)  Info << "FinesFields: Calculating source terms.\n" << endl;
     calcSource();
-    if(verbose_)  Info << "FinesFields: Updating dSauter.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating dSauter.\n" << endl;
     updateDSauter();
-    if(verbose_)  Info << "FinesFields: Updating dHydMix.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating dHydMix.\n" << endl;
     updateDHydMix();
-    if(verbose_)  Info << "FinesFields: Updating Froude.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating Froude.\n" << endl;
     updateFroude();
-    if(verbose_)  Info << "FinesFields: Updating FanningCoeff.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating FanningCoeff.\n" << endl;
     updateFanningCoeff();
-    if(verbose_)  Info << "FinesFields: Updating DragCoeff.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating DragCoeff.\n" << endl;
     updateDragCoeff();
-    if(verbose_)  Info << "FinesFields: Updating uDyn.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Updating uDyn.\n" << endl;
     updateUDyn();
-    if(verbose_)  Info << "FinesFields: Integrating alphas.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Integrating alphas.\n" << endl;
     integrateFields();
-    if(verbose_)  Info << "FinesFields: Update finished.\n" << endl;
+    if (verbose_)  Info << "FinesFields: Update finished.\n" << endl;
 }
 
 
 void FinesFields::calcSource()
 {
+    if (!clogKin_ & !clogStick_) return;
     Sds_.primitiveFieldRef() = 0;
     deltaAlpha_.primitiveFieldRef() = 0.0;
-    scalar f(0.0);
-    scalar critpore(0.0);
-    scalar dmean(0.0);
-    scalar d1(0.0);
-    scalar d2(0.0);
+    scalar fKin = 0.0;
+    scalar fStick = 0.0;
+    scalar critpore = 0.0;
+    scalar dmean = 0.0;
+    scalar d1 = 0.0;
+    scalar d2 = 0.0;
+    scalar magU = 0.0;
+    scalar tauDeposition = 0.0;
 
     forAll(Sds_,cellI)
     {
-        // calculate everything in units auf dSauter
-        critpore = nCrit_*dFine_.value()/dSauter_[cellI];
-        // pore size from hydraulic radius
-        dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
-        // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
-        dmean /= ratioHydraulicPore_;
-        d1 = dmean * (1 - poresizeWidth_);
-        d2 = dmean * (1 + poresizeWidth_);
+        fKin = 0.0;
+        fStick = 0.0;
+        if (clogKin_ && alphaP_[cellI] > alphaMinClog_) // no kinetic cloggig in dilute regions
+        {
+            // calculate everything in units auf dSauter
+            critpore = nCrit_*dFine_.value()/dSauter_[cellI];
+            // pore size from hydraulic radius
+            dmean = 2 * (1 - alphaP_[cellI]) / ( (1 + poresizeWidth_*poresizeWidth_/3) * 3 * alphaP_[cellI] );
+            // Sweeney and Martin, Acta Materialia 51 (2003): ratio of hydraulic to pore throat radius
+            dmean /= ratioHydraulicPore_;
+            d1 = dmean * (1 - poresizeWidth_);
+            d2 = dmean * (1 + poresizeWidth_);
 
-        f = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
-        if (f < 0)
-        {
-            f = 0.0;
+            fKin = (critpore*critpore*critpore - d1 * d1 * d1) / (d2 * d2 * d2 - d1 * d1 * d1);
+            if (fKin < 0) fKin = 0.0;
+            else if (fKin > 1.0) fKin = 1.0;
         }
-        else if (f > 1.0)
+
+        if (clogStick_)
         {
-            f = 1.0;
+            magU = mag(uReconstructed_()[cellI]); // use U reconstructed from phi to suppress oscillations at interfaces
+           // fStick = 1.0 / ( 1.0 + magU/uBind_) * alphaP_[cellI] / 0.65;
+            if (magU < uBindLow_) fStick = 1.0;
+            else if (magU > uBindHigh_) fStick = 0.0;
+            else fStick = 1.0 - (magU - uBindLow_) / (uBindHigh_ - uBindLow_);
+      //      if (fStick > 1.0) fStick = 1.0;
+            fStick *= alphaP_[cellI] / 0.65;
         }
 
         // at this point, voidfraction is still calculated from the true particle sizes
-        deltaAlpha_[cellI] = f * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
-        // too much volume occupied: release it (50% per time step)
+        deltaAlpha_[cellI] = max(fKin,fStick) * (alphaMax_ - alphaP_[cellI]) - alphaSt_[cellI];
+        // too much volume occupied: release it
         if (deltaAlpha_[cellI] < 0.0)
         {
-            Sds_[cellI] = 0.5*deltaAlpha_[cellI];
-        }
-        // volume too occupy available: deposit at most 80% of dyn hold up
-        else if (depRate_ * deltaAlpha_[cellI] > 0.8 * alphaDyn_[cellI])
-        {
-            Sds_[cellI] = 0.8 * alphaDyn_[cellI];
+            Sds_[cellI] = deltaAlpha_[cellI] / tauRelease_;
         }
+        // volume to occupy available
         else
         {
-            Sds_[cellI] = depRate_ * deltaAlpha_[cellI];
+            tauDeposition = depositionLength_ / max(mag(U_[cellI]),uMin_);
+            if (tauDeposition < 4.0*deltaT_) tauDeposition = 4.0*deltaT_; // do not deposit all dyn hold-up within less than 3 time steps
+            Sds_[cellI] = alphaDyn_[cellI] / tauDeposition;
         }
     }
 }
@@ -364,11 +444,12 @@ void FinesFields::integrateFields()
 
     fvScalarMatrix alphaStEqn
     (
-          fvm::ddt(alphaSt_)
-        + fvm::div(phiSt,alphaSt_)
+        fvm::ddt(alphaSt_)
         ==
         Sds_
     );
+    if (movingBed_) alphaStEqn += fvm::div(phiSt,alphaSt_);
+
     fvScalarMatrix alphaDynEqn
     (
         fvm::ddt(alphaDyn_)
@@ -377,10 +458,11 @@ void FinesFields::integrateFields()
         ==
         -Sds_
     );
+
     alphaStEqn.solve();
     alphaDynEqn.solve();
 
-    if(smoothing_)
+    if (smoothing_)
         particleCloud_.smoothingM().smoothen(alphaDyn_);
 
     // limit hold-ups, should be done more elegantly
@@ -422,13 +504,14 @@ void FinesFields::integrateFields()
 }
 
 
-void FinesFields::updateAlphaG()
+void FinesFields::updateAlphaG() // called after updateAlphaP() - correct voidfraction by removing space occupied by fines
 {
     alphaG_ = max(voidfraction_ - alphaSt_ - alphaDyn_, critVoidfraction_);
+    voidfraction_ = alphaG_;
 }
 
 
-void FinesFields::updateAlphaP()
+void FinesFields::updateAlphaP() // called first in the update cycle - voidfraction_ is current with DEM data
 {
     alphaP_ = 1.0 - voidfraction_ + SMALL;
 }
@@ -439,7 +522,7 @@ void FinesFields::updateDHydMix()
     forAll(dHydMix_,cellI)
     {
         scalar aPSt =  alphaP_[cellI] + alphaSt_[cellI];
-        if(aPSt < SMALL || aPSt > 1 - SMALL)
+        if (aPSt < SMALL || aPSt > 1 - SMALL)
             dHydMix_[cellI] = SMALL;
         else
             dHydMix_[cellI] = 2*(1 - aPSt) / (3*aPSt ) * dSauterMix_[cellI];
@@ -459,11 +542,11 @@ void FinesFields::updateDragCoeff()
     forAll(DragCoeff_,cellI)
     {
         Ref1 = Ref[cellI];
-        if(Ref1 <= SMALL)
+        if (Ref1 <= SMALL)
             Cd = 24.0 / SMALL;
-        else if(Ref1 <= 1.0)
+        else if (Ref1 <= 1.0)
             Cd = 24.0 / Ref1;
-        else if(Ref1 <= 1000)
+        else if (Ref1 <= 1000)
             Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
         else
             Cd = 0.44;
@@ -475,11 +558,11 @@ void FinesFields::updateDragCoeff()
         forAll(DragCoeff_.boundaryField()[patchI], faceI)
         {
             Ref1 = Ref.boundaryField()[patchI][faceI];
-            if(Ref1 <= SMALL)
+            if (Ref1 <= SMALL)
                 Cd = 24.0 / SMALL;
-            else if(Ref1 <= 1.0)
+            else if (Ref1 <= 1.0)
                 Cd = 24.0 / Ref1;
-            else if(Ref1 <= 1000)
+            else if (Ref1 <= 1000)
                 Cd = 24 * (1.0 + 0.15 * Foam::pow(Ref1,0.687) ) / Ref1;
             else
                 Cd = 0.44;
@@ -496,9 +579,9 @@ void FinesFields::updateDSauter()
     {
         scalar aP = alphaP_[cellI];
         scalar aSt = alphaSt_[cellI];
-        if(aSt < SMALL)
+        if (aSt < SMALL)
             dSauterMix_[cellI] = dSauter_[cellI];
-        else if(aP < SMALL)
+        else if (aP < SMALL)
             dSauterMix_[cellI] = dFine_.value();
         else
             dSauterMix_[cellI] = (aP + aSt) / (aP / dSauter_[cellI] + aSt / dFine_.value() );
@@ -541,7 +624,7 @@ void FinesFields::updateUDyn()
     {
         scalar mU(mag(U_[cellI]));
         scalar muDyn(mag(uDyn_[cellI]));
-        if(muDyn > mU && muDyn > SMALL)
+        if (muDyn > mU && muDyn > SMALL)
         {
             uDyn_[cellI] *= mU / muDyn;
         }
@@ -552,13 +635,24 @@ void FinesFields::updateUDyn()
         {
             scalar mU(mag(U_.boundaryField()[patchI][faceI]));
             scalar muDyn(mag(uDyn_.boundaryField()[patchI][faceI]));
-            if(muDyn > mU && muDyn > SMALL)
+            if (muDyn > mU && muDyn > SMALL)
             {
                 uDyn_.boundaryFieldRef()[patchI][faceI] *= mU / muDyn;
             }
         }
 }
 
+void FinesFields::updateUReconstructed()
+{
+    if (phi_.dimensions() == dimensionSet(1, 0, -1, 0, 0)) // compressible
+    {
+        uReconstructed_() = fvc::reconstruct(phi_) / (rhoG_ * voidfraction_);
+    }
+    else
+    {
+        uReconstructed_() = fvc::reconstruct(phi_) / voidfraction_;
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H
index 3098ff2f..b35fc5f3 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FinesFields.H
@@ -46,6 +46,16 @@ private:
 
     bool verbose_;
 
+    bool clogKin_;
+
+    bool clogStick_;
+
+    bool movingBed_;
+
+    word fluxFieldName_;
+
+    const surfaceScalarField& phi_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -90,11 +100,12 @@ private:
 
     volScalarField Sds_;
 
-    //volVectorField massFluxDyn_;
     surfaceScalarField massFluxDyn_;
 
     volVectorField uDyn_;
 
+    autoPtr<volVectorField> uReconstructed_; // velocity field can show oscillations at porosity steps
+
     dimensionedScalar dFine_;
 
     dimensionedScalar diffCoeff_;
@@ -109,9 +120,13 @@ private:
 
     scalar alphaMax_;
 
+    scalar alphaMinClog_;
+
     scalar critVoidfraction_;
 
-    scalar depRate_;
+    scalar deltaT_;
+
+    scalar depositionLength_;
 
     scalar exponent_;
 
@@ -123,6 +138,14 @@ private:
 
     scalar ratioHydraulicPore_;
 
+    scalar tauRelease_;
+
+    scalar uBindHigh_;
+
+    scalar uBindLow_;
+
+    scalar uMin_;
+
     void calcSource();
 
     void integrateFields();
@@ -143,6 +166,8 @@ private:
 
     void updateUDyn();
 
+    void updateUReconstructed();
+
 public:
 
     //- Runtime type information
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C
index 3bca651c..8c3154ae 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/forceModel/forceModel.C
@@ -69,7 +69,8 @@ forceModel::forceModel
             IOobject::AUTO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
+        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero),
+        "zeroGradient"
     ),
     expParticleForces_
     (   IOobject
@@ -81,7 +82,8 @@ forceModel::forceModel
             IOobject::AUTO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
+        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero),
+       "zeroGradient"
     ),
     coupleForce_(true),
     modelType_(sm.modelType()),
diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
index 3f1f907e..c97fc9b0 100644
--- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
+++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.C
@@ -65,14 +65,70 @@ constDiffSmoothing::constDiffSmoothing
     propsDict_(dict.subDict(typeName + "Props")),
     lowerLimit_(readScalar(propsDict_.lookup("lowerLimit"))),
     upperLimit_(readScalar(propsDict_.lookup("upperLimit"))),
-    smoothingLength_(dimensionedScalar("smoothingLength", dimLength, readScalar(propsDict_.lookup("smoothingLength")))),
-    smoothingLengthReferenceField_(dimensionedScalar("smoothingLengthReferenceField", dimLength, readScalar(propsDict_.lookup("smoothingLength")))),
-    DT_("DT", dimensionSet(0,2,-1,0,0), 0.),
+    smoothingLength_(propsDict_.lookupOrDefault<scalar>("smoothingLength", -1.0)),
+    smoothingLengthReference_(propsDict_.lookupOrDefault<scalar>("smoothingLengthReference",smoothingLength_)),
+    smoothingLengthFieldName_(propsDict_.lookupOrDefault<word>("smoothingLengthFieldName","smoothingLengthField")),
+    smoothingLengthField_
+    (   IOobject
+        (
+            smoothingLengthFieldName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("smoothingLength", dimensionSet(0,1,0,0,0,0,0), smoothingLength_),
+        "zeroGradient"
+    ),
+    smoothingLengthReferenceFieldName_(propsDict_.lookupOrDefault<word>("smoothingLengthReferenceFieldName","smoothingLengthReferenceField")),
+    smoothingLengthReferenceField_
+    (   IOobject
+        (
+            smoothingLengthReferenceFieldName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("smoothingLengthReference", dimensionSet(0,1,0,0,0,0,0), smoothingLengthReference_),
+        "zeroGradient"
+    ),
+    DT_
+    (   IOobject
+        (
+            "DT",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("DT", dimensionSet(0,2,-1,0,0,0,0), 0.0),
+        "zeroGradient"
+    ),
     verbose_(propsDict_.found("verbose"))
 {
-    if(propsDict_.found("smoothingLengthReferenceField"))
+    // either use scalar or field parameters for smoothing
+    if (smoothingLength_ > 0.0 || smoothingLengthReference_ > 0.0)
     {
-       smoothingLengthReferenceField_.value() = double(readScalar(propsDict_.lookup("smoothingLengthReferenceField")));
+        if (smoothingLengthField_.headerOk() || smoothingLengthReferenceField_.headerOk())
+        {
+            FatalError <<"constDiffSmoothing: Either use scalar or field parameter for smoothing.\n" << abort(FatalError);
+        }
+    }
+
+    if (smoothingLength_ < 0.0 && !smoothingLengthField_.headerOk())
+    {
+        FatalError <<"constDiffSmoothing: Provide scalar or field parameter for smoothing.\n" << abort(FatalError);
+    }
+  
+    // if no scalar length for smoothing wrt reference field is provided and no
+    // such field, use smoothingLengthField
+    if (smoothingLengthReference_ < 0.0 && !smoothingLengthReferenceField_.headerOk())
+    {
+        smoothingLengthReferenceField_ = smoothingLengthField_;
     }
 
     checkFields(sSmoothField_);
@@ -104,8 +160,8 @@ void constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
     sSmoothField.oldTime()=fieldSrc;
     sSmoothField.oldTime().correctBoundaryConditions();
 
-    double deltaT = sSmoothField.mesh().time().deltaTValue();
-    DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
+    dimensionedScalar deltaT = sSmoothField.mesh().time().deltaT();
+    DT_ = smoothingLengthField_ * smoothingLengthField_ / deltaT;
 
     // do smoothing
     solve
@@ -145,8 +201,8 @@ void constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
     vSmoothField.oldTime()=fieldSrc;
     vSmoothField.oldTime().correctBoundaryConditions();
 
-    double deltaT = vSmoothField_.mesh().time().deltaTValue();
-    DT_.value() = smoothingLength_.value() * smoothingLength_.value() / deltaT;
+    dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
+    DT_ = smoothingLengthField_ * smoothingLengthField_ / deltaT;
 
     // do smoothing
     solve
@@ -183,9 +239,8 @@ void constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc) const
     double sourceStrength = 1e5; //large number to keep reference values constant
 
     dimensionedScalar deltaT = vSmoothField.mesh().time().deltaT();
-    DT_.value() = smoothingLengthReferenceField_.value()
-                         * smoothingLengthReferenceField_.value() / deltaT.value();
-
+    DT_ = smoothingLengthReferenceField_ * smoothingLengthReferenceField_ / deltaT;
+        
     tmp<volScalarField> NLarge
     (
         new volScalarField
diff --git a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.H b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.H
index dc6ec5aa..bb7e5ea2 100644
--- a/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.H
+++ b/src/lagrangian/cfdemParticle/subModels/smoothingModel/constDiffSmoothing/constDiffSmoothing.H
@@ -60,11 +60,25 @@ class constDiffSmoothing
 private:
 
     dictionary propsDict_;
+
     const scalar lowerLimit_;
+
     const scalar upperLimit_;
-    const dimensionedScalar smoothingLength_;
-    dimensionedScalar smoothingLengthReferenceField_;
-    mutable dimensionedScalar DT_;
+
+    const scalar smoothingLength_;
+
+    const scalar smoothingLengthReference_;
+
+    word smoothingLengthFieldName_;
+
+    mutable volScalarField smoothingLengthField_;
+
+    word smoothingLengthReferenceFieldName_;
+
+    mutable volScalarField smoothingLengthReferenceField_;
+
+    mutable volScalarField DT_;
+
     const bool verbose_;
 
 public:
diff --git a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/constant/couplingProperties b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/constant/couplingProperties
index 6569a646..df0e2572 100755
--- a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/constant/couplingProperties
+++ b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/constant/couplingProperties
@@ -38,7 +38,7 @@ meshMotionModel noMeshMotion;
 
 regionModel allRegion;
 
-IOModel "basicIO";
+IOModel "off";
 
 dataExchangeModel twoWayMPI;//twoWayM2M;//twoWayFiles;//oneWayVTK;//
 
@@ -55,7 +55,7 @@ forceModels
     dSauter
     Fines
     FanningDynFines
-    ErgunStatFines
+    BeetstraDragPoly
     gradPForce
     viscForce
 );
@@ -140,13 +140,14 @@ viscForceProps
     interpolation;
 }
 
-ErgunStatFinesProps
+BeetstraDragProps
 {
+    fines true;
+    dFine 0.000388;
     velFieldName "U";
     granVelFieldName "Us";
     densityFieldName "rho";
     voidfractionFieldName "voidfraction";
-    phi 1;
 }
 
 FanningDynFinesProps
@@ -207,7 +208,7 @@ dividedProps
 {
     alphaMin 0.25;
     scaleUpVol 1.0;
-    weight 1.0; //1.33;
+    weight 1.0;
     verbose;
 }
 
diff --git a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/blockMeshDict b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/blockMeshDict
index 6f51c7d8..f429a14f 100644
--- a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/blockMeshDict
+++ b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/blockMeshDict
@@ -13,7 +13,7 @@ FoamFile
     object      blockMeshDict;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#include "../../geometry"
+#include "../geometry"
 convertToMeters 1;
 
 vertices #codeStream
diff --git a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/controlDict b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/controlDict
index 40c67cd6..8b4d7592 100644
--- a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/controlDict
+++ b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/controlDict
@@ -59,7 +59,7 @@ functions
 
         volInt
 	{
-        	type            volRegion;
+        	type            volFieldValue;
 	 	functionObjectLibs ("libfieldFunctionObjects.so");
 		writeControl   timeStep;
                 log             true;
@@ -76,7 +76,7 @@ functions
 
         inflow
         {
-            type            surfaceRegion;
+            type            surfaceFieldValue;
             libs ("libfieldFunctionObjects.so");
             writeControl   timeStep;
             log             true;
diff --git a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/fvOptions b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/fvOptions
index 9f7f44e2..fd8f8b41 100644
--- a/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/fvOptions
+++ b/tutorials/cfdemSolverRhoPimple/FinesColumn/CFD/system/fvOptions
@@ -15,20 +15,6 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-/*
-source1
-{
-    type            temperatureLimitsConstraint;
-    selectionMode   all;
-    active          true;
-
-        temperatureLimitsConstraintCoeffs
-        {
-            Tmin     288;
-            Tmax     298;
-        }
-}
-*/
 source1
 {
     type        limitTemperature;
@@ -37,7 +23,9 @@ source1
     {
         active          yes;
         selectionMode   all;
-        Tmin     288;
-        Tmax     298;
+        Tmin     288; // OF4
+        Tmax     298; // OF4
+        min      $Tmin; // OF5,OF6
+        max      $Tmax; // OF5,OF6
     }
 }