diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
index 832cefa73d09298b893f0e0c700a88dc5259ba53..6a88d1d9ceb045d97bfc4edd6cf6d800536e5228 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/options
@@ -4,8 +4,8 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/transportModel \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
index 9d502f7dd48d882c9c771de2719fb139c6bddc7c..bbf37904ca83743a749d6936ce5d54ad4229295b 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
@@ -41,7 +41,13 @@ Foam::RASModels::kineticTheoryModel::kineticTheoryModel
     const word& type
 )
 :
-    eddyViscosity<RASModel<PhaseCompressibleTurbulenceModel<phaseModel> > >
+    eddyViscosity
+    <
+        RASModel<EddyDiffusivity<ThermalDiffusivity
+        <
+            PhaseCompressibleTurbulenceModel<phaseModel>
+        > > >
+    >
     (
         type,
         alpha,
@@ -183,7 +189,10 @@ bool Foam::RASModels::kineticTheoryModel::read()
     (
         eddyViscosity
         <
-            RASModel<PhaseCompressibleTurbulenceModel<phaseModel> >
+            RASModel<EddyDiffusivity<ThermalDiffusivity
+            <
+                PhaseCompressibleTurbulenceModel<phaseModel>
+            > > >
         >::read()
     )
     {
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
index bfa1f872cb877ecdfd0e15e5e2d0faf8cbc5ef25..9aadb83d12738a7f3fb91b54cff381d5772ab15a 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H
@@ -48,6 +48,8 @@ SourceFiles
 #include "RASModel.H"
 #include "eddyViscosity.H"
 #include "PhaseCompressibleTurbulenceModel.H"
+#include "ThermalDiffusivity.H"
+#include "EddyDiffusivity.H"
 #include "phaseModel.H"
 #include "dragModel.H"
 #include "viscosityModel.H"
@@ -72,7 +74,10 @@ class kineticTheoryModel
 :
     public eddyViscosity
     <
-        RASModel<PhaseCompressibleTurbulenceModel<phaseModel> >
+        RASModel<EddyDiffusivity<ThermalDiffusivity
+        <
+            PhaseCompressibleTurbulenceModel<phaseModel>
+        > > >
     >
 {
     // Private data
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C
index 65d857d9c842f58f147769165d86a13685ddcfc2..eb0a07004439a5b028680abf4fd024e99cb1ebe0 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C
@@ -29,6 +29,9 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
+#include "ThermalDiffusivity.H"
+#include "EddyDiffusivity.H"
+
 #include "laminar.H"
 #include "RASModel.H"
 #include "LESModel.H"
@@ -39,6 +42,7 @@ makeBaseTurbulenceModel
     volScalarField,
     compressibleTurbulenceModel,
     PhaseCompressibleTurbulenceModel,
+    ThermalDiffusivity,
     phaseModel
 );
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
index 8cca6f493f098cf93b841a128ff96286df2dbb46..adb23d40d2c63d264adb8d5ac2aa2cb501d04474 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C
@@ -40,7 +40,13 @@ Foam::RASModels::phasePressureModel::phasePressureModel
     const word& type
 )
 :
-    eddyViscosity<RASModel<PhaseCompressibleTurbulenceModel<phaseModel> > >
+    eddyViscosity
+    <
+        RASModel<EddyDiffusivity<ThermalDiffusivity
+        <
+            PhaseCompressibleTurbulenceModel<phaseModel>
+        > > >
+    >
     (
         type,
         alpha,
@@ -87,7 +93,10 @@ bool Foam::RASModels::phasePressureModel::read()
     (
         eddyViscosity
         <
-            RASModel<PhaseCompressibleTurbulenceModel<phaseModel> >
+            RASModel<EddyDiffusivity<ThermalDiffusivity
+            <
+                PhaseCompressibleTurbulenceModel<phaseModel>
+            > > >
         >::read()
     )
     {
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
index a2b17dc1c7660edd7cd938faa83167d8d8b79fa2..68734e5b4f792e2f80ff503a6fc2b2d819e05a6e 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H
@@ -54,6 +54,8 @@ SourceFiles
 #include "RASModel.H"
 #include "eddyViscosity.H"
 #include "PhaseCompressibleTurbulenceModel.H"
+#include "ThermalDiffusivity.H"
+#include "EddyDiffusivity.H"
 #include "phaseModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -71,7 +73,10 @@ class phasePressureModel
 :
     public eddyViscosity
     <
-        RASModel<PhaseCompressibleTurbulenceModel<phaseModel> >
+        RASModel<EddyDiffusivity<ThermalDiffusivity
+        <
+            PhaseCompressibleTurbulenceModel<phaseModel>
+        > > >
     >
 {
     // Private data
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C
index 2befe1585afb66487fcd56f38a52c895cbed32ce..e79a9469fd568bf737d25afb0bbc0e23f2738e32 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.C
@@ -110,10 +110,10 @@ template<class BasePhaseModel, class ThermoType>
 Foam::tmp<Foam::scalarField>
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::mu
 (
-    const label patchI
+    const label patchi
 ) const
 {
-    return thermo_->mu(patchI);
+    return thermo_->mu(patchi);
 }
 
 
@@ -129,10 +129,10 @@ template<class BasePhaseModel, class ThermoType>
 Foam::tmp<Foam::scalarField>
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::nu
 (
-    const label patchI
+    const label patchi
 ) const
 {
-    return thermo_->nu(patchI);
+    return thermo_->nu(patchi);
 }
 
 
@@ -148,10 +148,75 @@ template<class BasePhaseModel, class ThermoType>
 Foam::tmp<Foam::scalarField>
 Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::kappa
 (
-    const label patchI
+    const label patchi
 ) const
 {
-    return thermo_->kappa(patchI);
+    return thermo_->kappa(patchi);
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::kappaEff
+(
+    const volScalarField& alphat
+) const
+{
+    return thermo_->kappaEff(alphat);
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::scalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::kappaEff
+(
+    const scalarField& alphat,
+    const label patchi
+) const
+{
+    return thermo_->kappaEff(alphat, patchi);
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::alpha() const
+{
+    return thermo_->alpha();
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::scalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::alpha
+(
+    const label patchi
+) const
+{
+    return thermo_->alpha(patchi);
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::volScalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::alphaEff
+(
+    const volScalarField& alphat
+) const
+{
+    return thermo_->alphaEff(alphat);
+}
+
+
+template<class BasePhaseModel, class ThermoType>
+Foam::tmp<Foam::scalarField>
+Foam::ThermoPhaseModel<BasePhaseModel, ThermoType>::alphaEff
+(
+    const scalarField& alphat,
+    const label patchi
+) const
+{
+    return thermo_->alphaEff(alphat, patchi);
 }
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H
index 0c2b243a83c2d7e2613ceded67c82cf479194e1c..923187931ee2093b6f5ea3b5496ba4cdf0f2f4bd 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/ThermoPhaseModel/ThermoPhaseModel.H
@@ -105,19 +105,51 @@ public:
             virtual tmp<volScalarField> mu() const;
 
             //- Access the laminar dynamic viscosity
-            virtual tmp<scalarField> mu(const label patchI) const;
+            virtual tmp<scalarField> mu(const label patchi) const;
 
             //- Return the laminar kinematic viscosity
             virtual tmp<volScalarField> nu() const;
 
             //- Access the laminar kinematic viscosity
-            virtual tmp<scalarField> nu(const label patchI) const;
+            virtual tmp<scalarField> nu(const label patchi) const;
 
             //- Return the laminar thermal conductivity
             virtual tmp<volScalarField> kappa() const;
 
             //- Access the laminar thermal conductivity
-            virtual tmp<scalarField> kappa(const label patchI) const;
+            virtual tmp<scalarField> kappa(const label patchi) const;
+
+            //- Return the laminar thermal conductivity
+            virtual tmp<volScalarField> kappaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Access the laminar thermal conductivity
+            virtual tmp<scalarField> kappaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
+
+            //- Return the thermal diffusivity for enthalpy
+            virtual tmp<volScalarField> alpha() const;
+
+            //- Return the thermal diffusivity for enthalpy on a patch
+            virtual tmp<scalarField> alpha(const label patchi) const;
+
+            //- Return the thermal diffusivity for enthalpy
+            virtual tmp<volScalarField> alphaEff
+            (
+                const volScalarField& alphat
+            ) const;
+
+            //- Return the thermal diffusivity for enthalpy on a patch
+            virtual tmp<scalarField> alphaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const;
 };
 
 
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
index 082bb6c019169e69038b8fe7870ef5519601d609..34e371a04073f15a3b2c48b1b0f0e22793c6af93 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/phaseModel/phaseModel/phaseModel.H
@@ -246,19 +246,51 @@ public:
             virtual tmp<volScalarField> mu() const = 0;
 
             //- Return the laminar dynamic viscosity on a patch
-            virtual tmp<scalarField> mu(const label patchI) const = 0;
+            virtual tmp<scalarField> mu(const label patchi) const = 0;
 
             //- Return the laminar kinematic viscosity
             virtual tmp<volScalarField> nu() const = 0;
 
             //- Return the laminar kinematic viscosity on a patch
-            virtual tmp<scalarField> nu(const label patchI) const = 0;
+            virtual tmp<scalarField> nu(const label patchi) const = 0;
 
             //- Return the laminar thermal conductivity
             virtual tmp<volScalarField> kappa() const = 0;
 
             //- Return the laminar thermal conductivity on a patch
-            virtual tmp<scalarField> kappa(const label patchI) const = 0;
+            virtual tmp<scalarField> kappa(const label patchi) const = 0;
+
+            //- Return the laminar thermal conductivity
+            virtual tmp<volScalarField> kappaEff
+            (
+                const volScalarField& alphat
+            ) const = 0;
+
+            //- Access the laminar thermal conductivity
+            virtual tmp<scalarField> kappaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const = 0;
+
+            //- Return the thermal diffusivity for enthalpy
+            virtual tmp<volScalarField> alpha() const = 0;
+
+            //- Return the thermal diffusivity for enthalpy on a patch
+            virtual tmp<scalarField> alpha(const label patchi) const = 0;
+
+            //- Return the thermal diffusivity for enthalpy
+            virtual tmp<volScalarField> alphaEff
+            (
+                const volScalarField& alphat
+            ) const = 0;
+
+            //- Return the thermal diffusivity for enthalpy on a patch
+            virtual tmp<scalarField> alphaEff
+            (
+                const scalarField& alphat,
+                const label patchi
+            ) const = 0;
 
 
         // Turbulence
diff --git a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
index 25bc49fd489146aad2a4a4125f59714bc54b0e76..47ae733b14e5cea2b9b7d2b8b02e67e277f60ffd 100644
--- a/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
+++ b/applications/solvers/multiphase/reactingTwoPhaseEulerFoam/phaseSystems/reactionThermo/hRefConstThermos.C
@@ -33,30 +33,18 @@ License
 #include "perfectGas.H"
 #include "perfectFluid.H"
 #include "rhoConst.H"
-#include "incompressiblePerfectGas.H"
-#include "hConstThermo.H"
-#include "janafThermo.H"
+
 #include "sensibleEnthalpy.H"
-#include "absoluteEnthalpy.H"
-#include "absoluteInternalEnergy.H"
-#include "thermo.H"
+
+#include "hRefConstThermo.H"
 
 #include "constTransport.H"
-#include "sutherlandTransport.H"
 
 #include "pureMixture.H"
-#include "homogeneousMixture.H"
-#include "inhomogeneousMixture.H"
-#include "veryInhomogeneousMixture.H"
 #include "multiComponentMixture.H"
-#include "reactingMixture.H"
-#include "singleStepReactingMixture.H"
 
 #include "thermoPhysicsTypes.H"
 
-#include "hRefConstThermo.H"
-
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -119,6 +107,19 @@ constTransport
     >
 > constRefFluidEThermoPhysics;
 
+typedef
+constTransport
+<
+    species::thermo
+    <
+        hRefConstThermo
+        <
+            rhoConst<specie>
+        >,
+        sensibleInternalEnergy
+    >
+> constRefRhoConstEThermoPhysics;
+
 
 // pureMixture, sensibleEnthalpy:
 
@@ -218,6 +219,15 @@ makeReactionMixtureThermo
     constRefFluidEThermoPhysics
 );
 
+makeReactionMixtureThermo
+(
+    rhoThermo,
+    rhoReactionThermo,
+    heRhoThermo,
+    multiComponentMixture,
+    constRefRhoConstEThermoPhysics
+);
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
diff --git a/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.C b/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.C
index 8e66570e1617d267390dedfa893eaeb1db31f9ac..eec1be8c8ee81c98fcfb2c909cfd83cd1d80e514 100644
--- a/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.C
+++ b/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.C
@@ -30,6 +30,14 @@ License
 template<class BasicTurbulenceModel>
 void Foam::EddyDiffusivity<BasicTurbulenceModel>::correctNut()
 {
+    // Read Prt if provided
+    Prt_ = dimensioned<scalar>::lookupOrDefault
+    (
+        "Prt",
+        this->coeffDict(),
+        1.0
+    );
+
     alphat_ = this->rho_*this->nut()/Prt_;
     alphat_.correctBoundaryConditions();
 }
@@ -41,7 +49,7 @@ template<class BasicTurbulenceModel>
 Foam::EddyDiffusivity<BasicTurbulenceModel>::EddyDiffusivity
 (
     const word& type,
-    const geometricOneField& alpha,
+    const alphaField& alpha,
     const volScalarField& rho,
     const volVectorField& U,
     const surfaceScalarField& alphaRhoPhi,
@@ -62,23 +70,14 @@ Foam::EddyDiffusivity<BasicTurbulenceModel>::EddyDiffusivity
         propertiesName
     ),
 
-    // Prt_
-    // (
-    //     dimensioned<scalar>::lookupOrAddToDict
-    //     (
-    //         "Prt",
-    //         this->coeffDict_,
-    //         1.0
-    //     )
-    // ),
-
+    // Cannot read Prt yet
     Prt_("Prt", dimless, 1.0),
 
     alphat_
     (
         IOobject
         (
-            "alphat",
+            IOobject::groupName("alphat", U.group()),
             this->runTime_.timeName(),
             this->mesh_,
             IOobject::MUST_READ,
@@ -89,34 +88,6 @@ Foam::EddyDiffusivity<BasicTurbulenceModel>::EddyDiffusivity
 {}
 
 
-// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-Foam::autoPtr<Foam::EddyDiffusivity<BasicTurbulenceModel> >
-Foam::EddyDiffusivity<BasicTurbulenceModel>::New
-(
-    const volScalarField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& phi,
-    const transportModel& transport,
-    const word& propertiesName
-)
-{
-    return autoPtr<EddyDiffusivity>
-    (
-        static_cast<EddyDiffusivity*>(
-        BasicTurbulenceModel::New
-        (
-            rho,
-            U,
-            phi,
-            transport,
-            propertiesName
-        ).ptr())
-    );
-}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
diff --git a/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.H b/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.H
index a01107f6301d261be71f17642008663e46059931..a6693fac10bab2bf50950cf11ca82066f4574457 100644
--- a/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.H
+++ b/src/TurbulenceModels/compressible/EddyDiffusivity/EddyDiffusivity.H
@@ -82,7 +82,7 @@ public:
         EddyDiffusivity
         (
             const word& type,
-            const geometricOneField& alpha,
+            const alphaField& alpha,
             const volScalarField& rho,
             const volVectorField& U,
             const surfaceScalarField& alphaRhoPhi,
@@ -92,19 +92,6 @@ public:
         );
 
 
-    // Selectors
-
-        //- Return a reference to the selected turbulence model
-        static autoPtr<EddyDiffusivity> New
-        (
-            const volScalarField& rho,
-            const volVectorField& U,
-            const surfaceScalarField& phi,
-            const transportModel& trasportModel,
-            const word& propertiesName = turbulenceModel::propertiesName
-        );
-
-
     //- Destructor
     virtual ~EddyDiffusivity()
     {}
diff --git a/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.C b/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.C
index ad6de84714511139b8bb4bc9d00f0ce5318feda1..4a85a59d074422c300e3389db6a4211aebc82ae3 100644
--- a/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.C
+++ b/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.C
@@ -31,7 +31,7 @@ template<class BasicTurbulenceModel>
 Foam::ThermalDiffusivity<BasicTurbulenceModel>::ThermalDiffusivity
 (
     const word& type,
-    const geometricOneField& alpha,
+    const alphaField& alpha,
     const volScalarField& rho,
     const volVectorField& U,
     const surfaceScalarField& alphaRhoPhi,
@@ -56,6 +56,34 @@ Foam::ThermalDiffusivity<BasicTurbulenceModel>::ThermalDiffusivity
 
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+Foam::autoPtr<Foam::ThermalDiffusivity<BasicTurbulenceModel> >
+Foam::ThermalDiffusivity<BasicTurbulenceModel>::New
+(
+    const alphaField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+{
+    return autoPtr<ThermalDiffusivity>
+    (
+        static_cast<ThermalDiffusivity*>(
+        BasicTurbulenceModel::New
+        (
+            alpha,
+            rho,
+            U,
+            phi,
+            transport,
+            propertiesName
+        ).ptr())
+    );
+}
+
+
 template<class BasicTurbulenceModel>
 Foam::autoPtr<Foam::ThermalDiffusivity<BasicTurbulenceModel> >
 Foam::ThermalDiffusivity<BasicTurbulenceModel>::New
diff --git a/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.H b/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.H
index 5996798d871a8af4cd917fd44e12b8de736c13e3..2066e44a1daf176f4e74a088359358f4636bee30 100644
--- a/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.H
+++ b/src/TurbulenceModels/compressible/ThermalDiffusivity/ThermalDiffusivity.H
@@ -53,7 +53,7 @@ class ThermalDiffusivity
 
 public:
 
-    typedef geometricOneField alphaField;
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
     typedef volScalarField rhoField;
     typedef typename BasicTurbulenceModel::transportModel transportModel;
 
@@ -64,7 +64,7 @@ public:
         ThermalDiffusivity
         (
             const word& type,
-            const geometricOneField& alpha,
+            const alphaField& alpha,
             const volScalarField& rho,
             const volVectorField& U,
             const surfaceScalarField& alphaRhoPhi,
@@ -76,6 +76,18 @@ public:
 
     // Selectors
 
+        //- Return a reference to the selected turbulence model
+        static autoPtr<ThermalDiffusivity> New
+        (
+            const alphaField& alpha,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& phi,
+            const transportModel& trasportModel,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
         //- Return a reference to the selected turbulence model
         static autoPtr<ThermalDiffusivity> New
         (
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
index 02e87c7a776e08424112330dc5db2a793334a28c..8ace068f487521f17f25fb5ba8d32b8fd30bc4e7 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/convectiveHeatTransfer/convectiveHeatTransferFvPatchScalarField.C
@@ -110,14 +110,15 @@ void convectiveHeatTransferFvPatchScalarField::updateCoeffs()
 
     const label patchi = patch().index();
 
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const compressible::turbulenceModel& turbModel =
+        db().lookupObject<compressible::turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            dimensionedInternalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                compressible::turbulenceModel::propertiesName,
+                dimensionedInternalField().group()
+            )
+        );
 
     const scalarField alphaEffw(turbModel.alphaEff(patchi));
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
index 63d700d3a2e52638d63168a7657448e11f69ccf1..d3d5116348a1b46d51ab227379794b508dae23dc 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/totalFlowRateAdvectiveDiffusive/totalFlowRateAdvectiveDiffusiveFvPatchScalarField.C
@@ -32,6 +32,7 @@ License
 #include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
 Foam::totalFlowRateAdvectiveDiffusiveFvPatchScalarField::
 totalFlowRateAdvectiveDiffusiveFvPatchScalarField
 (
@@ -154,10 +155,11 @@ void Foam::totalFlowRateAdvectiveDiffusiveFvPatchScalarField::updateCoeffs()
     const label patchI = patch().index();
 
     const LESModel<EddyDiffusivity<compressible::turbulenceModel> >& turbModel =
-    db().lookupObject
-    <
-        LESModel<EddyDiffusivity<compressible::turbulenceModel> >
-    >   (
+        db().lookupObject
+        <
+            LESModel<EddyDiffusivity<compressible::turbulenceModel> >
+        >
+        (
             IOobject::groupName
             (
                 turbulenceModel::propertiesName,
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
index 98a5a05f5e75c94032da7d005f00f69c67dba5ae..7a9983fc20f987f6cff4b507d75bfcbe5f4b142e 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
@@ -200,15 +200,15 @@ void alphatJayatillekeWallFunctionFvPatchScalarField::updateCoeffs()
     const label patchi = patch().index();
 
     // Retrieve turbulence properties from model
-
-    const turbulenceModel& turbModel = db().lookupObject<turbulenceModel>
-    (
-        IOobject::groupName
+    const compressible::turbulenceModel& turbModel =
+        db().lookupObject<compressible::turbulenceModel>
         (
-            turbulenceModel::propertiesName,
-            dimensionedInternalField().group()
-        )
-    );
+            IOobject::groupName
+            (
+                compressible::turbulenceModel::propertiesName,
+                dimensionedInternalField().group()
+            )
+        );
 
     const scalar Cmu25 = pow025(Cmu_);
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
index 90d2236005691bcda5e309cc602d7b304826703c..ece3ac9e36dd897524af892cd6dc08360ccd3a78 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "alphatWallFunctionFvPatchScalarField.H"
+#include "compressibleTurbulenceModel.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "addToRunTimeSelectionTable.H"
@@ -44,8 +45,6 @@ alphatWallFunctionFvPatchScalarField::alphatWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF),
-    rhoName_("rho"),
-    nutName_("nut"),
     Prt_(0.85)
 {}
 
@@ -59,8 +58,6 @@ alphatWallFunctionFvPatchScalarField::alphatWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(ptf, p, iF, mapper),
-    rhoName_(ptf.rhoName_),
-    nutName_(ptf.nutName_),
     Prt_(ptf.Prt_)
 {}
 
@@ -73,8 +70,6 @@ alphatWallFunctionFvPatchScalarField::alphatWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(p, iF, dict),
-    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
-    nutName_(dict.lookupOrDefault<word>("nut", "nut")),
     Prt_(dict.lookupOrDefault<scalar>("Prt", 0.85))
 {}
 
@@ -85,8 +80,6 @@ alphatWallFunctionFvPatchScalarField::alphatWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(awfpsf),
-    rhoName_(awfpsf.rhoName_),
-    nutName_(awfpsf.nutName_),
     Prt_(awfpsf.Prt_)
 {}
 
@@ -98,8 +91,6 @@ alphatWallFunctionFvPatchScalarField::alphatWallFunctionFvPatchScalarField
 )
 :
     fixedValueFvPatchScalarField(awfpsf, iF),
-    rhoName_(awfpsf.rhoName_),
-    nutName_(awfpsf.nutName_),
     Prt_(awfpsf.Prt_)
 {}
 
@@ -113,13 +104,23 @@ void alphatWallFunctionFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const scalarField& rhow =
-        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
+    const label patchi = patch().index();
 
-    const scalarField& nutw =
-        patch().lookupPatchField<volScalarField, scalar>(nutName_);
+    // Retrieve turbulence properties from model
+    const compressibleTurbulenceModel& turbModel =
+        db().lookupObject<compressibleTurbulenceModel>
+        (
+            IOobject::groupName
+            (
+                compressibleTurbulenceModel::propertiesName,
+                dimensionedInternalField().group()
+            )
+        );
 
-    operator==(rhow*nutw/Prt_);
+    const scalarField& rhow = turbModel.rho().boundaryField()[patchi];
+    const tmp<scalarField> tnutw = turbModel.nut(patchi);
+
+    operator==(rhow*tnutw/Prt_);
 
     fixedValueFvPatchScalarField::updateCoeffs();
 }
@@ -128,8 +129,6 @@ void alphatWallFunctionFvPatchScalarField::updateCoeffs()
 void alphatWallFunctionFvPatchScalarField::write(Ostream& os) const
 {
     fvPatchField<scalar>::write(os);
-    writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
-    writeEntryIfDifferent<word>(os, "nut", "nut", nutName_);
     os.writeKeyword("Prt") << Prt_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H
index f44ee6b8cf13c048843de3e88f813883a16334a9..308128ecce8a13f67e86cf8eedc2132824c73200 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H
@@ -95,12 +95,6 @@ class alphatWallFunctionFvPatchScalarField
 {
     // Private data
 
-        //- Name of density field (default = rho)
-        word rhoName_;
-
-        //- Name of turbulent viscosity field (default = nut)
-        word nutName_;
-
         //- Turbulent Prandtl number (default = 0.85)
         scalar Prt_;
 
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/makeTurbulenceModel.H b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/makeTurbulenceModel.H
index 91c3a2e71b896bb3b5f51d2da480ae7cf97678de..2bf0664c1631acd35c2f7fb4bdeca6b4805dccea 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/makeTurbulenceModel.H
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/makeTurbulenceModel.H
@@ -23,7 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#define makeBaseTurbulenceModel(Alpha, Rho, baseModel, BaseModel, Transport)   \
+#define makeBaseTurbulenceModel(                                               \
+    Alpha, Rho, baseModel, BaseModel, TDModel, Transport)                      \
                                                                                \
     namespace Foam                                                             \
     {                                                                          \
@@ -41,7 +42,7 @@ License
             dictionary                                                         \
         );                                                                     \
                                                                                \
-        typedef BaseModel<CompressibleTurbulenceModel<Transport> >             \
+        typedef TDModel<BaseModel<Transport> >                                 \
             Transport##BaseModel;                                              \
                                                                                \
                                                                                \
diff --git a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
index d5342a15d7787736a99579eee81fc3b9cbb2e46c..db2d4140700320b205c91b2c77a538a17314b15d 100644
--- a/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
+++ b/src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C
@@ -23,12 +23,12 @@ License
 
 \*---------------------------------------------------------------------------*/
 
+#include "CompressibleTurbulenceModel.H"
 #include "compressibleTransportModel.H"
 #include "fluidThermo.H"
 #include "addToRunTimeSelectionTable.H"
 #include "makeTurbulenceModel.H"
 
-#include "CompressibleTurbulenceModel.H"
 #include "ThermalDiffusivity.H"
 #include "EddyDiffusivity.H"
 
@@ -43,17 +43,18 @@ makeBaseTurbulenceModel
     geometricOneField,
     volScalarField,
     compressibleTurbulenceModel,
+    CompressibleTurbulenceModel,
     ThermalDiffusivity,
     fluidThermo
 );
 
 #define makeRASModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
-    (fluidThermoThermalDiffusivity, RAS, Type)
+    (fluidThermoCompressibleTurbulenceModel, RAS, Type)
 
 #define makeLESModel(Type)                                                     \
     makeTemplatedTurbulenceModel                                               \
-    (fluidThermoThermalDiffusivity, LES, Type)
+    (fluidThermoCompressibleTurbulenceModel, LES, Type)
 
 
 // -------------------------------------------------------------------------- //
diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
index e3a3fcd3327468270110e23200f0208e7ca1f595..1bbbedc91df55a2e51dc9f80214520f217564d2d 100644
--- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
+++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H
@@ -141,8 +141,6 @@ protected:
             dimensionedScalar alphaOmega1_;
             dimensionedScalar alphaOmega2_;
 
-            dimensionedScalar Prt_;
-
             dimensionedScalar gamma1_;
             dimensionedScalar gamma2_;
 
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.air b/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.air
new file mode 100644
index 0000000000000000000000000000000000000000..60979da94f4cd22aaa91059d363efb400f581977
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.air
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.water b/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.water
new file mode 100644
index 0000000000000000000000000000000000000000..3ae0a65592e12a3e4de29a4f74fbdbe0d359f0fd
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/LES/bubbleColumn/0/alphat.water
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.air b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.air
new file mode 100644
index 0000000000000000000000000000000000000000..60979da94f4cd22aaa91059d363efb400f581977
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.air
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.water b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.water
new file mode 100644
index 0000000000000000000000000000000000000000..3ae0a65592e12a3e4de29a4f74fbdbe0d359f0fd
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumn/0/alphat.water
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.water;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/0/alphat.gas b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/0/alphat.gas
new file mode 100644
index 0000000000000000000000000000000000000000..de0de301914bd2350f7be7483c62bfd33b74e418
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/bubbleColumnEvaporatingReacting/0/alphat.gas
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.gas;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.air b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.air
new file mode 100644
index 0000000000000000000000000000000000000000..60979da94f4cd22aaa91059d363efb400f581977
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.air
@@ -0,0 +1,48 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.air;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            compressible::alphatWallFunction;
+        Prt             0.85;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.particles b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.particles
new file mode 100644
index 0000000000000000000000000000000000000000..b047d8d119ee3e73504c1aae0b99a3a3258975e3
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/RAS/fluidisedBed/0/alphat.particles
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.particles;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //
diff --git a/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/0/alphat.particles b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/0/alphat.particles
new file mode 100644
index 0000000000000000000000000000000000000000..b047d8d119ee3e73504c1aae0b99a3a3258975e3
--- /dev/null
+++ b/tutorials/multiphase/reactingTwoPhaseEulerFoam/laminar/fluidisedBed/0/alphat.particles
@@ -0,0 +1,47 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      alphat.particles;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -1 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    inlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    outlet
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    walls
+    {
+        type            calculated;
+        value           $internalField;
+    }
+
+    defaultFaces
+    {
+        type            empty;
+    }
+}
+
+// ************************************************************************* //