diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
index 43c0f9b2a9df4708c2e2386ec8a640acf5413a56..e754fa07b8be5ed0e770e4edf1204faa792a168b 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/Make/files
@@ -12,4 +12,8 @@ heatTransferModels/heatTransferModel/heatTransferModel.C
 heatTransferModels/heatTransferModel/newHeatTransferModel.C
 heatTransferModels/RanzMarshall/RanzMarshall.C
 
+liftModels/liftModel/liftModel.C
+liftModels/liftModel/newLiftModel.C
+liftModels/constantCoefficient/constantCoefficient.C
+
 LIB = $(FOAM_LIBBIN)/libcompressibleEulerianInterfacialModels
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.C
new file mode 100644
index 0000000000000000000000000000000000000000..069a9838849412833e6121f2da3a6fc783feeb4c
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.C
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "constantCoefficient.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace liftModels
+{
+    defineTypeNameAndDebug(constantCoefficient, 0);
+
+    addToRunTimeSelectionTable
+    (
+        liftModel,
+        constantCoefficient,
+        dictionary
+    );
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::liftModels::constantCoefficient::constantCoefficient
+(
+    const dictionary& dict,
+    const volScalarField& alpha1,
+    const phaseModel& phase1,
+    const phaseModel& phase2
+)
+:
+    liftModel(dict, alpha1, phase1, phase2),
+    Cl_("Cl", dimless, dict.lookup("Cl"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::liftModels::constantCoefficient::~constantCoefficient()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volVectorField> Foam::liftModels::constantCoefficient::F
+(
+    const volVectorField& U
+) const
+{
+    return
+        Cl_
+       *(phase1_*phase1_.rho() + phase2_*phase2_.rho())
+       *(
+            (phase1_.U() - phase2_.U())
+          ^ fvc::curl(U)
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.H
new file mode 100644
index 0000000000000000000000000000000000000000..054dc858bea59b662f4276ed750595ea5b859b96
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/constantCoefficient/constantCoefficient.H
@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::liftModels::constantCoefficient
+
+Description
+
+SourceFiles
+    constantCoefficient.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef constantCoefficient_H
+#define constantCoefficient_H
+
+#include "liftModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace liftModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class constantCoefficient Declaration
+\*---------------------------------------------------------------------------*/
+
+class constantCoefficient
+:
+    public liftModel
+{
+    // Private data
+
+        //- Constant lift coefficient
+        dimensionedScalar Cl_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("constantCoefficient");
+
+
+    // Constructors
+
+        //- Construct from components
+        constantCoefficient
+        (
+            const dictionary& dict,
+            const volScalarField& alpha1,
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        );
+
+
+    //- Destructor
+    virtual ~constantCoefficient();
+
+
+    // Member Functions
+
+        //- Lift force
+        tmp<volVectorField> F(const volVectorField& U) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace liftModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..36dbde08ee1edd74883e36d4431a21b7240aff92
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C
@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "liftModel.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(liftModel, 0);
+    defineRunTimeSelectionTable(liftModel, dictionary);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::liftModel::liftModel
+(
+    const dictionary& dict,
+    const volScalarField& alpha1,
+    const phaseModel& phase1,
+    const phaseModel& phase2
+)
+:
+    dict_(dict),
+    alpha1_(alpha1),
+    phase1_(phase1),
+    phase2_(phase2)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::liftModel::~liftModel()
+{}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H
new file mode 100644
index 0000000000000000000000000000000000000000..ef16d1660582884a10e1de18b4a8ac923a2e4e2b
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H
@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+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::liftModel
+
+Description
+
+SourceFiles
+    liftModel.C
+    newLiftModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef liftModel_H
+#define liftModel_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "dictionary.H"
+#include "phaseModel.H"
+#include "runTimeSelectionTables.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class liftModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class liftModel
+{
+protected:
+
+    // Protected data
+
+        const dictionary& dict_;
+        const volScalarField& alpha1_;
+        const phaseModel& phase1_;
+        const phaseModel& phase2_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("liftModel");
+
+
+    // Declare runtime construction
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            liftModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                const volScalarField& alpha1,
+                const phaseModel& phase1,
+                const phaseModel& phase2
+            ),
+            (dict, alpha1, phase1, phase2)
+        );
+
+
+    // Constructors
+
+        liftModel
+        (
+            const dictionary& dict,
+            const volScalarField& alpha1,
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        );
+
+
+    //- Destructor
+    virtual ~liftModel();
+
+
+    // Selectors
+
+        static autoPtr<liftModel> New
+        (
+            const dictionary& dict,
+            const volScalarField& alpha1,
+            const phaseModel& phase1,
+            const phaseModel& phase2
+        );
+
+
+    // Member Functions
+
+        //- Lift force
+        virtual tmp<volVectorField> F(const volVectorField& U) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/newLiftModel.C b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/newLiftModel.C
new file mode 100644
index 0000000000000000000000000000000000000000..21aad9034ef6a03acb819c1f6540e4a10a02e3d6
--- /dev/null
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/newLiftModel.C
@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "liftModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::autoPtr<Foam::liftModel> Foam::liftModel::New
+(
+    const dictionary& dict,
+    const volScalarField& alpha1,
+    const phaseModel& phase1,
+    const phaseModel& phase2
+)
+{
+    word liftModelType
+    (
+        dict.subDict(phase1.name()).lookup("type")
+    );
+
+    Info << "Selecting liftModel for phase "
+        << phase1.name()
+        << ": "
+        << liftModelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(liftModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn("liftModel::New")
+            << "Unknown liftModelType type "
+            << liftModelType << endl << endl
+            << "Valid liftModel types are : " << endl
+            << dictionaryConstructorTablePtr_->sortedToc()
+            << exit(FatalError);
+    }
+
+    return
+        cstrIter()
+        (
+            dict.subDict(phase1.name()).subDict(liftModelType + "Coeffs"),
+            alpha1,
+            phase1,
+            phase2
+        );
+}
+
+
+// ************************************************************************* //
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
index 8ba3382b86967b3d0d74522d8b0bcff64b8c215e..e1f9b221a35d38f187f33942878164f2201dd465 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -111,13 +111,6 @@ Foam::twoPhaseSystem::twoPhaseSystem
         lookup("Cvm")
     ),
 
-    Cl_
-    (
-        "Cl",
-        dimless,
-        lookup("Cl")
-    ),
-
     drag1_
     (
         dragModel::New
@@ -162,6 +155,28 @@ Foam::twoPhaseSystem::twoPhaseSystem
         )
     ),
 
+    lift1_
+    (
+        liftModel::New
+        (
+            subDict("lift"),
+            phase1_,
+            phase1_,
+            phase2_
+        )
+    ),
+
+    lift2_
+    (
+        liftModel::New
+        (
+            subDict("lift"),
+            phase2_,
+            phase2_,
+            phase1_
+        )
+    ),
+
     dispersedPhase_(lookup("dispersedPhase")),
 
     residualPhaseFraction_
@@ -311,11 +326,28 @@ Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::liftForce
     );
     volVectorField& liftForce = tliftForce();
 
-    volVectorField Ur(phase1_.U() - phase2_.U());
-
-    liftForce =
-        Cl_*(phase1_*phase1_.rho() + phase2_*phase2_.rho())
-       *(Ur ^ fvc::curl(U));
+    if (dispersedPhase_ == phase1_.name())
+    {
+        liftForce = lift1().F(U);
+    }
+    else if (dispersedPhase_ == phase2_.name())
+    {
+        liftForce = lift2().F(U);
+    }
+    else if (dispersedPhase_ == "both")
+    {
+        liftForce =
+        (
+            phase2_*lift1().F(U)
+          + phase1_*lift2().F(U)
+        );
+    }
+    else
+    {
+        FatalErrorIn("twoPhaseSystem::liftForce()")
+            << "dispersedPhase: " << dispersedPhase_ << " is incorrect"
+            << exit(FatalError);
+    }
 
     // Remove lift at fixed-flux boundaries
     forAll(phase1_.phi().boundaryField(), patchi)
@@ -631,7 +663,6 @@ bool Foam::twoPhaseSystem::read()
 
         lookup("sigma") >> sigma_;
         lookup("Cvm") >> Cvm_;
-        lookup("Cl") >> Cl_;
 
         // drag1_->read(*this);
         // drag2_->read(*this);
diff --git a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
index 0b501152c15f410169d8617da1c1cff494dce192..967ffc222641a77f4e79f33167b715fde48262f5 100644
--- a/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
+++ b/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ SourceFiles
 #include "IOdictionary.H"
 #include "phaseModel.H"
 #include "dragModel.H"
+#include "liftModel.H"
 #include "heatTransferModel.H"
 #include "volFields.H"
 #include "surfaceFields.H"
@@ -57,7 +58,7 @@ namespace Foam
 // Forward declarations
 class dragModel;
 class heatTransferModel;
-
+class liftModel;
 
 /*---------------------------------------------------------------------------*\
                       Class twoPhaseSystem Declaration
@@ -69,31 +70,57 @@ class twoPhaseSystem
 {
     // Private data
 
+        //- Reference to the mesh
         const fvMesh& mesh_;
 
+        //- Phase model 1
         phaseModel phase1_;
+
+        //- Phase model 2
         phaseModel phase2_;
 
+        //- Total volumetric flux
         surfaceScalarField phi_;
 
+        //- 
         volScalarField dgdt_;
 
+        //- Surface tension coefficient
         dimensionedScalar sigma_;
 
+        //- Virtual mass coefficient
         dimensionedScalar Cvm_;
-        dimensionedScalar Cl_;
 
+        //- Drag model for phase 1
         autoPtr<dragModel> drag1_;
+
+        //- Drag model for phase 2
         autoPtr<dragModel> drag2_;
 
+        //- Heat transfer model for phase 1 
         autoPtr<heatTransferModel> heatTransfer1_;
+
+        //- Heat transfer model for phase 2
         autoPtr<heatTransferModel> heatTransfer2_;
 
+        //- Lift model for phase 1
+        autoPtr<liftModel> lift1_;
+
+        //- Lift model for phase 2
+        autoPtr<liftModel> lift2_;
+
+        //- Name of the dispersed phase, or "both"
         word dispersedPhase_;
+
+        //- Residual phase fraction
         scalar residualPhaseFraction_;
+
+        //- Redisual slip
         dimensionedScalar residualSlip_;
 
 
+    // Private member functions
+
         //- Return the mixture flux
         tmp<surfaceScalarField> calcPhi() const;
 
@@ -113,147 +140,174 @@ public:
 
     // Member Functions
 
-        const fvMesh& mesh() const
-        {
-            return mesh_;
-        }
+        //- Return the drag coefficient
+        tmp<volScalarField> dragCoeff() const;
+
+        //- Return the heat transfer coefficient
+        tmp<volScalarField> heatTransferCoeff() const;
+
+        //- Return the lift force
+        tmp<volVectorField> liftForce(const volVectorField& U) const;
+
+        //- Return the mixture density
+        tmp<volScalarField> rho() const;
+
+        //- Return the mixture velocity
+        tmp<volVectorField> U() const;
 
-        const phaseModel& phase1() const
-        {
-            return phase1_;
-        }
+        //- Solve for the two-phase-fractions
+        void solve();
 
-        const phaseModel& phase2() const
-        {
-            return phase2_;
-        }
+        //- Correct two-phase properties other than turbulence
+        void correct();
 
-        const phaseModel& otherPhase(const phaseModel& phase) const
-        {
-            if (&phase == &phase1_)
+        //- Correct two-phase turbulence
+        void correctTurbulence();
+
+        //- Read base phaseProperties dictionary
+        bool read();
+
+        // Access
+
+            //- Return the mesh
+            const fvMesh& mesh() const
             {
-                return phase2_;
+                return mesh_;
             }
-            else
+
+            //- Return phase model 1
+            const phaseModel& phase1() const
             {
                 return phase1_;
             }
-        }
 
-        phaseModel& phase1()
-        {
-            return phase1_;
-        }
-
-        phaseModel& phase2()
-        {
-            return phase2_;
-        }
+            //- Return non-const access to phase model 1
+            phaseModel& phase1()
+            {
+                return phase1_;
+            }
 
-        //- Return the mixture flux
-        const surfaceScalarField& phi() const
-        {
-            return phi_;
-        }
+            //- Return phase model 2
+            const phaseModel& phase2() const
+            {
+                return phase2_;
+            }
 
-        //- Return the mixture flux
-        surfaceScalarField& phi()
-        {
-            return phi_;
-        }
-
-        const volScalarField& dgdt() const
-        {
-            return dgdt_;
-        }
-
-        volScalarField& dgdt()
-        {
-            return dgdt_;
-        }
-
-        const dragModel& drag1() const
-        {
-            return drag1_();
-        }
-
-        const dragModel& drag2() const
-        {
-            return drag2_();
-        }
-
-        const dragModel& drag(const phaseModel& phase) const
-        {
-            if (&phase == &phase1_)
+            //- Return non-const access to phase model 2
+            phaseModel& phase2()
             {
-                return drag1_();
+                return phase2_;
             }
-            else
+
+            //- Return the phase not given as an argument
+            const phaseModel& otherPhase(const phaseModel& phase) const
             {
-                return drag2_();
+                if (&phase == &phase1_)
+                {
+                    return phase2_;
+                }
+                else
+                {
+                    return phase1_;
+                }
             }
-        }
 
-        scalar residualPhaseFraction() const
-        {
-            return residualPhaseFraction_;
-        }
+            //- Return the mixture flux
+            const surfaceScalarField& phi() const
+            {
+                return phi_;
+            }
 
-        const dimensionedScalar& residualSlip() const
-        {
-            return residualSlip_;
-        }
+            //- Return non-const access to the the mixture flux
+            surfaceScalarField& phi()
+            {
+                return phi_;
+            }
 
-        tmp<volScalarField> dragCoeff() const;
-        tmp<volVectorField> liftForce(const volVectorField& U) const;
+            //- Return
+            const volScalarField& dgdt() const
+            {
+                return dgdt_;
+            }
 
-        const heatTransferModel& heatTransfer1() const
-        {
-            return heatTransfer1_();
-        }
+            //- Return non-const access to the
+            volScalarField& dgdt()
+            {
+                return dgdt_;
+            }
 
-        const heatTransferModel& heatTransfer2() const
-        {
-            return heatTransfer2_();
-        }
+            //- Return the drag model for phase 1
+            const dragModel& drag1() const
+            {
+                return drag1_();
+            }
 
-        tmp<volScalarField> heatTransferCoeff() const;
+            //- Return the drag model for phase 2
+            const dragModel& drag2() const
+            {
+                return drag2_();
+            }
 
-        //- Return the mixture density
-        tmp<volScalarField> rho() const;
+            //- Return the drag model for the supplied phase
+            const dragModel& drag(const phaseModel& phase) const
+            {
+                if (&phase == &phase1_)
+                {
+                    return drag1_();
+                }
+                else
+                {
+                    return drag2_();
+                }
+            }
 
-        //- Return the mixture velocity
-        tmp<volVectorField> U() const;
+            //- Return non-const access to the residual phase fraction
+            scalar residualPhaseFraction() const
+            {
+                return residualPhaseFraction_;
+            }
 
-        //- Return the surface tension coefficient
-        dimensionedScalar sigma() const
-        {
-            return sigma_;
-        }
+            //- Return the residual slip
+            const dimensionedScalar& residualSlip() const
+            {
+                return residualSlip_;
+            }
 
-        //- Return the virtual-mass coefficient
-        dimensionedScalar Cvm() const
-        {
-            return Cvm_;
-        }
+            //- Return the heat transfer model for phase 1
+            const heatTransferModel& heatTransfer1() const
+            {
+                return heatTransfer1_();
+            }
 
-        //- Return the lift coefficient
-        dimensionedScalar Cl() const
-        {
-            return Cl_;
-        }
+            //- Return the heat transfer model for phase 2
+            const heatTransferModel& heatTransfer2() const
+            {
+                return heatTransfer2_();
+            }
 
-        //- Solve for the two-phase-fractions
-        void solve();
+            //- Return the lift model for phase 1
+            const liftModel& lift1() const
+            {
+                return lift1_();
+            }
 
-        //- Correct two-phase properties other than turbulence
-        void correct();
+            //- Return the lift model for phase 2
+            const liftModel& lift2() const
+            {
+                return lift2_();
+            }
 
-        //- Correct two-phase turbulence
-        void correctTurbulence();
+            //- Return the surface tension coefficient
+            dimensionedScalar sigma() const
+            {
+                return sigma_;
+            }
 
-        //- Read base phaseProperties dictionary
-        bool read();
+            //- Return the virtual-mass coefficient
+            dimensionedScalar Cvm() const
+            {
+                return Cvm_;
+            }
 };
 
 
diff --git a/applications/utilities/postProcessing/sampling/sample/sample.C b/applications/utilities/postProcessing/sampling/sample/sample.C
index e357f2f7e535e1fc8f028399c80ae22e24d52cd2..59abc2e2f18ec888ad66ac078f9d29a90f250a0a 100644
--- a/applications/utilities/postProcessing/sampling/sample/sample.C
+++ b/applications/utilities/postProcessing/sampling/sample/sample.C
@@ -132,6 +132,11 @@ int main(int argc, char *argv[])
             )
         );
 
+        // Note: both IOsampledSets and IOsampledSurfaces read the
+        //       same dictionary. Unregister one to make sure no duplicates
+        //       trying to register
+        sSetsPtr().checkOut();
+
         sSurfsPtr.reset
         (
             new IOsampledSurfaces
@@ -160,6 +165,8 @@ int main(int argc, char *argv[])
             )
         );
 
+        sSetsPtr().checkOut();
+
         sSurfsPtr.reset
         (
             new IOsampledSurfaces
diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
index 86e2a29e6109b9857d0f1fdebd108f1fe1432dc6..6d5f4d2d47fc743f53127a19a3ebf5a5571747d3 100644
--- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
+++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "fvSchemes.H"
 #include "Time.H"
+#include "steadyStateDdtScheme.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -108,6 +109,11 @@ void Foam::fvSchemes::read(const dictionary& dict)
     )
     {
         defaultDdtScheme_ = ddtSchemes_.lookup("default");
+        steady_ =
+        (
+            word(defaultDdtScheme_)
+         == fv::steadyStateDdtScheme<scalar>::typeName
+        );
     }
 
 
@@ -364,7 +370,8 @@ Foam::fvSchemes::fvSchemes(const objectRegistry& obr)
             tokenList()
         )()
     ),
-    defaultFluxRequired_(false)
+    defaultFluxRequired_(false),
+    steady_(false)
 {
     // persistent settings across reads is incorrect
     clear();
diff --git a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H
index cebbc8142ab6720a04decb209fa73669aeebeff9..d16dd7a8168eb8061828310a734d584250148156 100644
--- a/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H
+++ b/src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,6 +78,10 @@ class fvSchemes
         dictionary fluxRequired_;
         bool defaultFluxRequired_;
 
+        //- Steady-state run indicator
+        //  Set true if the default ddtScheme is steadyState
+        bool steady_;
+
 
     // Private Member Functions
 
@@ -128,6 +132,18 @@ public:
 
             bool fluxRequired(const word& name) const;
 
+            //- Return true if the default ddtScheme is steadyState
+            bool steady() const
+            {
+                return steady_;
+            }
+
+            //- Return true if the default ddtScheme is not steadyState
+            bool transient() const
+            {
+                return !steady_;
+            }
+
 
         // Read
 
diff --git a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
index 3ca9ddd62c525faf6985ed6366698f702e6f5b7e..8283924496cb739fe189c9a3539581abd376241d 100644
--- a/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
+++ b/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -2348,12 +2348,19 @@ Foam::operator&
     GeometricField<Type, fvPatchField, volMesh>& Mphi = tMphi();
 
     // Loop over field components
-    for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+    if (M.hasDiag())
+    {
+        for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
+        {
+            scalarField psiCmpt(psi.field().component(cmpt));
+            scalarField boundaryDiagCmpt(M.diag());
+            M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
+            Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
+        }
+    }
+    else
     {
-        scalarField psiCmpt(psi.field().component(cmpt));
-        scalarField boundaryDiagCmpt(M.diag());
-        M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
-        Mphi.internalField().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
+        Mphi.internalField() = pTraits<Type>::zero;
     }
 
     Mphi.internalField() += M.lduMatrix::H(psi.field()) + M.source();
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index d536dc047fdb9a3ebd33fbe0687d898c7ddc1d8e..4851de3e2e64e9d0bc6979e3958fec69e9b598f6 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -198,16 +198,6 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
 
         alphaEqn.solve();
 
-        //// updated stress
-        //tauPrime.internalField() =
-        //    this->particleStressModel_->tauPrime
-        //    (
-        //        alpha_.internalField(),
-        //        rho.internalField(),
-        //        uSqrAverage.internalField()
-        //    )();
-        //tauPrime.correctBoundaryConditions();
-
 
         // Generate correction fields
         // ~~~~~~~~~~~~~~~~~
@@ -229,8 +219,6 @@ void Foam::PackingModels::Implicit<CloudType>::cacheFields(const bool store)
             (
                 cloudName + ":uCorrect",
                 fvc::reconstruct(phiCorrect_())
-            //  - deltaT*fvc::grad(tauPrime)/(rho*alpha_)
-            //  + (applyGravity_ ? deltaT*g*(1.0 - rhoc/rho) : 0.0)
             )
 
         );
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
index a3bec3a5778874f03b2b236e3064f8e8dd1156e9..47cc489688743c2d1c05de60d36df7efd4835c31 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/LES/bubbleColumn/constant/phaseProperties
@@ -51,6 +51,26 @@ heatTransfer
     water   RanzMarshall;
 }
 
+lift
+{
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    water
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
@@ -59,9 +79,6 @@ residualSlip            1e-2;
 // Virtual-mass coefficient
 Cvm             0.5;
 
-// Lift coefficient
-Cl              0;
-
 // Dispersed-phase turbulence coefficient
 Ct              1;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
index a3bec3a5778874f03b2b236e3064f8e8dd1156e9..47cc489688743c2d1c05de60d36df7efd4835c31 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/bubbleColumn/constant/phaseProperties
@@ -51,6 +51,26 @@ heatTransfer
     water   RanzMarshall;
 }
 
+lift
+{
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    water
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
@@ -59,9 +79,6 @@ residualSlip            1e-2;
 // Virtual-mass coefficient
 Cvm             0.5;
 
-// Lift coefficient
-Cl              0;
-
 // Dispersed-phase turbulence coefficient
 Ct              1;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
index b5eb1e24eee49e4b66cb673a6065fd15d6ee3ae0..b44f6c59349130c320b3a5c76caf5525f6571366 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/RAS/fluidisedBed/constant/phaseProperties
@@ -50,6 +50,26 @@ heatTransfer
     air         RanzMarshall;
 }
 
+lift
+{
+    particles
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase     particles;
 
 residualPhaseFraction   1e-3;
@@ -58,9 +78,6 @@ residualSlip            1e-2;
 // Virtual-mass coefficient
 Cvm             0;
 
-// Lift coefficient
-Cl              0;
-
 // Minimum allowable pressure
 pMin            10000;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
index b7325116e31f9bcf61fd7177b45779058dd024ce..99146f4bc03087eee69649b8e3165707e57ab937 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumn/constant/phaseProperties
@@ -51,6 +51,26 @@ heatTransfer
     water   RanzMarshall;
 }
 
+lift
+{
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    water
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
@@ -59,9 +79,6 @@ residualSlip            1e-2;
 // Virtual-mass ceofficient
 Cvm             0.5;
 
-// Lift coefficient
-Cl              0;
-
 // Minimum allowable pressure
 pMin            10000;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties
index d9304132f75d573a0a262e64d6cba880d083b95c..d91a9d1da71893010a3a41d1212643722478a855 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/bubbleColumnIATE/constant/phaseProperties
@@ -71,6 +71,26 @@ heatTransfer
     water   RanzMarshall;
 }
 
+lift
+{
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    water
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
@@ -79,9 +99,6 @@ residualSlip            1e-2;
 // Virtual-mass ceofficient
 Cvm             0.5;
 
-// Lift coefficient
-Cl              0;
-
 // Minimum allowable pressure
 pMin            10000;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
index cb95638812c8df4060622783826d56924dcc110f..91671832abd38462e5eeed6a418397eaf9ba3ab3 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed/constant/phaseProperties
@@ -50,6 +50,26 @@ heatTransfer
     air         RanzMarshall;
 }
 
+lift
+{
+    particles
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase     particles;
 
 residualPhaseFraction   1e-3;
@@ -58,9 +78,6 @@ residualSlip            1e-2;
 // Virtual-mass ceofficient
 Cvm             0;
 
-// Lift coefficient
-Cl              0;
-
 // Minimum allowable pressure
 pMin            10000;
 
diff --git a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
index c31927d84d98018694a365bfe22eac8e5c4e75ca..427e6ef28f85bd6cca3c30ec191205119e86ec2d 100644
--- a/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
+++ b/tutorials/multiphase/twoPhaseEulerFoam/laminar/mixerVessel2D/constant/phaseProperties
@@ -51,6 +51,26 @@ heatTransfer
     water   RanzMarshall;
 }
 
+lift
+{
+    air
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+    water
+    {
+        type constantCoefficient;
+        constantCoefficientCoeffs
+        {
+            Cl          0;
+        }
+    }
+}
+
 dispersedPhase          both;
 
 residualPhaseFraction   1e-3;
@@ -59,9 +79,6 @@ residualSlip            1e-2;
 // Virtual-mass coefficient
 Cvm             0.5;
 
-// Lift coefficient
-Cl              0;
-
 // Minimum allowable pressure
 pMin            10000;