diff --git a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
index b5c9b0607e..257f2451af 100644
--- a/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
+++ b/applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseModel/StationaryPhaseModel/StationaryPhaseModel.C
@@ -316,11 +316,21 @@ Foam::StationaryPhaseModel<BasePhaseModel>::divq(volScalarField& he) const
 {
     const volScalarField& alpha = *this;
 
-    return -fvm::laplacian
+    const surfaceScalarField alphaKappa
     (
-        fvc::interpolate(alpha)*fvc::interpolate(this->thermo().alphahe()),
-        he
+        alpha.name() + '*' + this->thermo().kappa().name(),
+        fvc::interpolate(alpha)*fvc::interpolate(this->thermo().kappa())
     );
+
+    // Return heat flux source as an implicit energy correction
+    // to the temperature gradient flux
+    return
+       -fvc::laplacian(alphaKappa, this->thermo().T())
+       -fvm::laplacianCorrection
+        (
+            alphaKappa/fvc::interpolate(this->thermo().Cpv()),
+            he
+        );
 }
 
 
diff --git a/src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C b/src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C
index 2bfa8e5a4b..75fc64ca5e 100644
--- a/src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C
+++ b/src/ThermophysicalTransportModels/fluid/laminar/Fickian/Fickian.C
@@ -309,7 +309,7 @@ tmp<fvScalarMatrix> Fickian<BasicThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField sumJ
     (
diff --git a/src/ThermophysicalTransportModels/fluid/laminar/Fourier/Fourier.C b/src/ThermophysicalTransportModels/fluid/laminar/Fourier/Fourier.C
index 62a627ba7d..9e2c6e0a20 100644
--- a/src/ThermophysicalTransportModels/fluid/laminar/Fourier/Fourier.C
+++ b/src/ThermophysicalTransportModels/fluid/laminar/Fourier/Fourier.C
@@ -130,8 +130,12 @@ Fourier<BasicThermophysicalTransportModel>::divq(volScalarField& he) const
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-        -correction(fvm::laplacian(this->alpha()*thermo.alphahe(), he))
-        -fvc::laplacian(this->alpha()*thermo.kappa(), thermo.T());
+        -fvc::laplacian(this->alpha()*thermo.kappa(), thermo.T())
+        -fvm::laplacianCorrection
+         (
+             this->alpha()*thermo.kappa()/thermo.Cpv(),
+             he
+         );
 }
 
 
diff --git a/src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C b/src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C
index e15bca5f28..1edcfc90f4 100644
--- a/src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C
+++ b/src/ThermophysicalTransportModels/fluid/laminar/MaxwellStefan/MaxwellStefan.C
@@ -310,7 +310,7 @@ tmp<fvScalarMatrix> MaxwellStefan<BasicThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField sumJ
     (
diff --git a/src/ThermophysicalTransportModels/fluid/turbulence/eddyDiffusivity/eddyDiffusivity.C b/src/ThermophysicalTransportModels/fluid/turbulence/eddyDiffusivity/eddyDiffusivity.C
index 57824b631d..fb3e160df0 100644
--- a/src/ThermophysicalTransportModels/fluid/turbulence/eddyDiffusivity/eddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/fluid/turbulence/eddyDiffusivity/eddyDiffusivity.C
@@ -163,8 +163,8 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
-        -correction(fvm::laplacian(this->alpha()*this->alphaEff(), he))
-        -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T());
+        -fvc::laplacian(this->alpha()*this->kappaEff(), this->thermo().T())
+        -fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C b/src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
index b3e9b7a802..491d33da70 100644
--- a/src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/fluid/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
@@ -162,7 +162,7 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
     const PtrList<volScalarField>& Y = composition.Y();
 
     tmpDivq.ref() -=
-        correction(fvm::laplacian(this->alpha()*this->alphaEff(), he));
+        fvm::laplacianCorrection(this->alpha()*this->alphaEff(), he);
 
     surfaceScalarField hGradY
     (
diff --git a/src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.C b/src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.C
index f8368c9bed..9a933d7626 100644
--- a/src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.C
+++ b/src/ThermophysicalTransportModels/solid/anisotropic/anisotropic.C
@@ -401,19 +401,17 @@ Foam::solidThermophysicalTransportModels::anisotropic::divq
 {
     const solidThermo& thermo = this->thermo();
     const volSymmTensorField Kappa(this->Kappa());
+    const surfaceVectorField& Sf = thermo.mesh().Sf();
+    const surfaceScalarField& magSf = thermo.mesh().magSf();
 
     // Return heat flux source as an implicit energy correction
     // to the temperature gradient flux
     return
        -fvc::laplacian(Kappa, thermo.T())
-       -correction
+       -fvm::laplacianCorrection
         (
-            fvm::laplacian
-            (
-                Kappa/thermo.Cv(),
-                e,
-                "laplacian(alphae,e)"
-            )
+            (Sf & fvc::interpolate(Kappa/thermo.Cv()) & Sf)/sqr(magSf),
+            e
         );
 }
 
diff --git a/src/ThermophysicalTransportModels/solid/isotropic/isotropic.C b/src/ThermophysicalTransportModels/solid/isotropic/isotropic.C
index be11430a5f..91f8633023 100644
--- a/src/ThermophysicalTransportModels/solid/isotropic/isotropic.C
+++ b/src/ThermophysicalTransportModels/solid/isotropic/isotropic.C
@@ -105,15 +105,7 @@ Foam::solidThermophysicalTransportModels::isotropic::divq
     // to the temperature gradient flux
     return
        -fvc::laplacian(kappa(), thermo.T())
-       -correction
-        (
-            fvm::laplacian
-            (
-                kappa()/thermo.Cv(),
-                e,
-                "laplacian(alphae,e)"
-            )
-        );
+       -fvm::laplacianCorrection(kappa()/thermo.Cv(), e);
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
index 663b5f60c9..0c41f4835b 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.C
@@ -27,6 +27,7 @@ License
 #include "surfaceFields.H"
 #include "fvMatrix.H"
 #include "laplacianScheme.H"
+#include "gaussLaplacianScheme.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -229,9 +230,9 @@ laplacian
     const word& name
 )
 {
-    tmp<fvMatrix<Type>> Laplacian(fvm::laplacian(tgamma(), vf, name));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tgamma(), vf, name));
     tgamma.clear();
-    return Laplacian;
+    return tLaplacian;
 }
 
 
@@ -260,9 +261,9 @@ laplacian
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    tmp<fvMatrix<Type>> Laplacian(fvm::laplacian(tgamma(), vf));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tgamma(), vf));
     tgamma.clear();
-    return Laplacian;
+    return tLaplacian;
 }
 
 
@@ -325,9 +326,79 @@ laplacian
     const GeometricField<Type, fvPatchField, volMesh>& vf
 )
 {
-    tmp<fvMatrix<Type>> tfvm(fvm::laplacian(tGamma(), vf));
+    tmp<fvMatrix<Type>> tLaplacian(fvm::laplacian(tGamma(), vf));
     tGamma.clear();
-    return tfvm;
+    return tLaplacian;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const GeometricField<scalar, fvPatchField, volMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return fvm::laplacianCorrection(fvc::interpolate(gamma), vf);
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const tmp<GeometricField<scalar, fvPatchField, volMesh>>& tgamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tLaplacianCorrection
+    (
+        fvm::laplacianCorrection(tgamma(), vf)
+    );
+    tgamma.clear();
+    return tLaplacianCorrection;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const GeometricField<scalar, fvsPatchField, surfaceMesh>& gamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    return correction
+    (
+        fv::gaussLaplacianScheme<Type, scalar>::fvmLaplacianUncorrected
+        (
+            gamma*vf.mesh().magSf(),
+            vf.mesh().nonOrthDeltaCoeffs(),
+            vf
+        )
+    );
+}
+
+
+template<class Type>
+tmp<fvMatrix<Type>>
+laplacianCorrection
+(
+    const tmp<GeometricField<scalar, fvsPatchField, surfaceMesh>>& tGamma,
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+)
+{
+    tmp<fvMatrix<Type>> tLaplacianCorrection
+    (
+        fvm::laplacianCorrection(tGamma(), vf)
+    );
+    tGamma.clear();
+    return tLaplacianCorrection;
 }
 
 
diff --git a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
index e5652a7f4b..b831e65403 100644
--- a/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
+++ b/src/finiteVolume/finiteVolume/fvm/fvmLaplacian.H
@@ -176,6 +176,36 @@ namespace fvm
         const tmp<GeometricField<GType, fvsPatchField, surfaceMesh>>&,
         const GeometricField<Type, fvPatchField, volMesh>&
     );
+
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const GeometricField<scalar, fvPatchField, volMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const tmp<GeometricField<scalar, fvPatchField, volMesh>>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const GeometricField<scalar, fvsPatchField, surfaceMesh>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
+
+    template<class Type>
+    tmp<fvMatrix<Type>> laplacianCorrection
+    (
+        const tmp<GeometricField<scalar, fvsPatchField, surfaceMesh>>&,
+        const GeometricField<Type, fvPatchField, volMesh>&
+    );
 }