diff --git a/applications/solvers/combustion/fireFoam/YEEqn.H b/applications/solvers/combustion/fireFoam/YEEqn.H
index 479ae4c2eb..6726a1964f 100644
--- a/applications/solvers/combustion/fireFoam/YEEqn.H
+++ b/applications/solvers/combustion/fireFoam/YEEqn.H
@@ -23,7 +23,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
               ==
                 parcels.SYi(i, Yi)
               + surfaceFilm.Srho(i)
@@ -63,7 +63,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             )
           : -dpdt
         )
-      - fvm::laplacian(thermophysicalTransport->alphaEff(), he)
+      + thermophysicalTransport->divq(he)
      ==
         combustion->Qdot()
       + radiation->Sh(thermo, he)
diff --git a/applications/solvers/combustion/reactingFoam/YEqn.H b/applications/solvers/combustion/reactingFoam/YEqn.H
index 9346dd3f97..8149acd6a0 100644
--- a/applications/solvers/combustion/reactingFoam/YEqn.H
+++ b/applications/solvers/combustion/reactingFoam/YEqn.H
@@ -23,7 +23,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
              ==
                 reaction->R(Yi)
               + fvOptions(rho, Yi)
diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
index 7d886e7fa5..0dfa63e3bb 100644
--- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
+++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/YEqn.H
@@ -42,7 +42,7 @@ if (Y.size())
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport.alphaEff(), Yi)
+              + thermophysicalTransport.divj(Yi)
              ==
                 reaction.R(Yi)
               + fvOptions(rho, Yi)
diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
index 3824190673..04ff49768b 100644
--- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
+++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H
@@ -24,7 +24,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
               ==
                 coalParcels.SYi(i, Yi)
               + combustion->R(Yi)
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
index c2aa62d37e..1f97c7236d 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H
@@ -24,7 +24,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
               ==
                 parcels.SYi(i, Yi)
               + fvOptions(rho, Yi)
diff --git a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
index c8e6cabc00..35b68bcc07 100644
--- a/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
+++ b/applications/solvers/lagrangian/reactingParcelFoam/simpleReactingParcelFoam/YEqn.H
@@ -22,7 +22,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             fvScalarMatrix YEqn
             (
                 mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
              ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
diff --git a/applications/solvers/lagrangian/sprayFoam/YEqn.H b/applications/solvers/lagrangian/sprayFoam/YEqn.H
index 1ee3dd74b9..b365b0bc4a 100644
--- a/applications/solvers/lagrangian/sprayFoam/YEqn.H
+++ b/applications/solvers/lagrangian/sprayFoam/YEqn.H
@@ -23,7 +23,7 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             (
                 fvm::ddt(rho, Yi)
               + mvConvection->fvmDiv(phi, Yi)
-              - fvm::laplacian(thermophysicalTransport->alphaEff(), Yi)
+              + thermophysicalTransport->divj(Yi)
              ==
                 parcels.SYi(i, Yi)
               + combustion->R(Yi)
diff --git a/src/ThermophysicalTransportModels/LES/LESThermophysicalTransportModel/LESThermophysicalTransportModel.H b/src/ThermophysicalTransportModels/LES/LESThermophysicalTransportModel/LESThermophysicalTransportModel.H
index 2ad4fcf982..69854c0fed 100644
--- a/src/ThermophysicalTransportModels/LES/LESThermophysicalTransportModel/LESThermophysicalTransportModel.H
+++ b/src/ThermophysicalTransportModels/LES/LESThermophysicalTransportModel/LESThermophysicalTransportModel.H
@@ -200,12 +200,6 @@ public:
             return this->thermo().alphahe(patchi);
         }
 
-        //- Return the heat flux
-        virtual tmp<volVectorField> q() const = 0;
-
-        //- Return the source term for the energy equation
-        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const = 0;
-
         //- Correct the LES transport
         virtual void correct();
 
diff --git a/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.C b/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.C
index 5d69333fed..65155275b0 100644
--- a/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.C
@@ -122,7 +122,7 @@ tmp<volVectorField>eddyDiffusivity<BasicThermophysicalTransportModel>::q() const
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->alphaEff()
+       -this->alpha()*this->alphaEff()
        *fvc::grad(this->thermo().he())
     );
 }
@@ -135,7 +135,36 @@ eddyDiffusivity<BasicThermophysicalTransportModel>::divq
     volScalarField& he
 ) const
 {
-    return -fvm::laplacian(this->alphaEff(), he);
+    return -fvm::laplacian(this->alpha()*this->alphaEff(), he);
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<volVectorField>eddyDiffusivity<BasicThermophysicalTransportModel>::j
+(
+    const volScalarField& Yi
+) const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "j(" + Yi.name() + ')',
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+        -this->alpha()*this->alphaEff()*fvc::grad(Yi)
+    );
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+eddyDiffusivity<BasicThermophysicalTransportModel>::divj
+(
+    volScalarField& Yi
+) const
+{
+    return -fvm::laplacian(this->alpha()*this->alphaEff(), Yi);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.H b/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.H
index 035cd66a09..d36c8ba891 100644
--- a/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.H
+++ b/src/ThermophysicalTransportModels/LES/LESeddyDiffusivity/LESeddyDiffusivity.H
@@ -157,6 +157,12 @@ public:
         //- Return the source term for the energy equation
         virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
         //- Correct the eddyDiffusivity viscosity
         virtual void correct();
 };
diff --git a/src/ThermophysicalTransportModels/RAS/RASThermophysicalTransportModel/RASThermophysicalTransportModel.H b/src/ThermophysicalTransportModels/RAS/RASThermophysicalTransportModel/RASThermophysicalTransportModel.H
index 27a6b557dc..8507631800 100644
--- a/src/ThermophysicalTransportModels/RAS/RASThermophysicalTransportModel/RASThermophysicalTransportModel.H
+++ b/src/ThermophysicalTransportModels/RAS/RASThermophysicalTransportModel/RASThermophysicalTransportModel.H
@@ -200,12 +200,6 @@ public:
             return this->thermo().alphahe(patchi);
         }
 
-        //- Return the heat flux
-        virtual tmp<volVectorField> q() const = 0;
-
-        //- Return the source term for the energy equation
-        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const = 0;
-
         //- Correct the RAS transport
         virtual void correct();
 
diff --git a/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.C b/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.C
index ae850e91b3..f6a0e7fc1d 100644
--- a/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.C
@@ -122,8 +122,7 @@ tmp<volVectorField>eddyDiffusivity<BasicThermophysicalTransportModel>::q() const
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->alphaEff()
-       *fvc::grad(this->thermo().he())
+       -this->alpha()*this->alphaEff()*fvc::grad(this->thermo().he())
     );
 }
 
@@ -135,7 +134,36 @@ eddyDiffusivity<BasicThermophysicalTransportModel>::divq
     volScalarField& he
 ) const
 {
-    return -fvm::laplacian(this->alphaEff(), he);
+    return -fvm::laplacian(this->alpha()*this->alphaEff(), he);
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<volVectorField>eddyDiffusivity<BasicThermophysicalTransportModel>::j
+(
+    const volScalarField& Yi
+) const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "j(" + Yi.name() + ')',
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+        -this->alpha()*this->alphaEff()*fvc::grad(Yi)
+    );
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+eddyDiffusivity<BasicThermophysicalTransportModel>::divj
+(
+    volScalarField& Yi
+) const
+{
+    return -fvm::laplacian(this->alpha()*this->alphaEff(), Yi);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.H b/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.H
index f4995f46da..a247807615 100644
--- a/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.H
+++ b/src/ThermophysicalTransportModels/RAS/RASeddyDiffusivity/RASeddyDiffusivity.H
@@ -157,6 +157,12 @@ public:
         //- Return the source term for the energy equation
         virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
         //- Correct the eddyDiffusivity viscosity
         virtual void correct();
 };
diff --git a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
index db6a59dafd..d4b7a4d856 100644
--- a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
+++ b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
@@ -76,7 +76,7 @@ tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::q() const
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-        -this->alpha()*fvc::grad(this->thermo().he())
+        -this->alpha()*this->thermo().alpha()*fvc::grad(this->thermo().he())
     );
 }
 
@@ -85,7 +85,33 @@ template<class BasicThermophysicalTransportModel>
 tmp<fvScalarMatrix>
 Fourier<BasicThermophysicalTransportModel>::divq(volScalarField& he) const
 {
-    return -fvm::laplacian(this->alpha(), he);
+    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), he);
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::j
+(
+    const volScalarField& Yi
+) const
+{
+    return volVectorField::New
+    (
+        IOobject::groupName
+        (
+            "j(" + Yi.name() + ')',
+            this->momentumTransport().alphaRhoPhi().group()
+        ),
+        -this->alpha()*this->thermo().alpha()*fvc::grad(Yi)
+    );
+}
+
+
+template<class BasicThermophysicalTransportModel>
+tmp<fvScalarMatrix>
+Fourier<BasicThermophysicalTransportModel>::divj(volScalarField& Yi) const
+{
+    return -fvm::laplacian(this->alpha()*this->thermo().alpha(), Yi);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
index 6dfbd9697b..dd3d4f8ed2 100644
--- a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
+++ b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
@@ -98,6 +98,12 @@ public:
         //- Return the source term for the energy equation
         virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
 
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
+
         //- Correct the Fourier viscosity
         virtual void correct();
 };
diff --git a/src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.H b/src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.H
index 6b907bce48..cedd38d512 100644
--- a/src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.H
+++ b/src/ThermophysicalTransportModels/laminar/laminarThermophysicalTransportModel/laminarThermophysicalTransportModel.H
@@ -200,12 +200,6 @@ public:
             return this->thermo().alphahe(patchi);
         }
 
-        //- Return the heat flux
-        virtual tmp<volVectorField> q() const = 0;
-
-        //- Return the source term for the energy equation
-        virtual tmp<fvScalarMatrix> divq(volScalarField& he) const = 0;
-
         //- Correct the laminar transport
         virtual void correct();
 
diff --git a/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H b/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
index 7576e0395f..266aa3bd3a 100644
--- a/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
+++ b/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
@@ -122,6 +122,12 @@ public:
         //- Return the source term for the energy equation
         virtual tmp<fvScalarMatrix> divq(volScalarField& he) const = 0;
 
+        //- Return the specie flux for the given specie mass-fraction
+        virtual tmp<volVectorField> j(const volScalarField& Yi) const = 0;
+
+        //- Return the source term for the given specie mass-fraction equation
+        virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const = 0;
+
         //- Solve the thermophysical transport model equations
         //  and correct the transport coefficients
         virtual void correct() = 0;