diff --git a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
index 12fe71cc51..f82e512e06 100644
--- a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
+++ b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.C
@@ -105,18 +105,18 @@ Fourier<TurbulenceThermophysicalTransportModel>::DEff
 
 
 template<class BasicThermophysicalTransportModel>
-tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::q() const
+tmp<surfaceScalarField> Fourier<BasicThermophysicalTransportModel>::q() const
 {
     const thermoModel& thermo = this->thermo();
 
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -(this->alpha()*thermo.kappa())*fvc::grad(thermo.T())
+        -fvc::interpolate(this->alpha()*thermo.kappa())*fvc::snGrad(thermo.T())
     );
 }
 
@@ -136,7 +136,7 @@ Fourier<BasicThermophysicalTransportModel>::divq(volScalarField& he) const
 
 
 template<class BasicThermophysicalTransportModel>
-tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::j
+tmp<surfaceScalarField> Fourier<BasicThermophysicalTransportModel>::j
 (
     const volScalarField& Yi
 ) const
@@ -147,7 +147,7 @@ tmp<volVectorField>Fourier<BasicThermophysicalTransportModel>::j
            " unityLewisFourier"
         << exit(FatalError);
 
-    return tmp<volVectorField>(nullptr);
+    return tmp<surfaceScalarField>(nullptr);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
index 0fc68cc55b..7d2ef66db0 100644
--- a/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
+++ b/src/ThermophysicalTransportModels/laminar/Fourier/Fourier.H
@@ -113,13 +113,13 @@ public:
         ) const;
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const;
+        virtual tmp<surfaceScalarField> q() const;
 
         //- 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;
+        virtual tmp<surfaceScalarField> j(const volScalarField& Yi) const;
 
         //- Return the source term for the given specie mass-fraction equation
         virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
diff --git a/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.C b/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.C
index 7f68a751e9..f2d613135d 100644
--- a/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.C
+++ b/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.C
@@ -25,6 +25,8 @@ License
 
 #include "unityLewisFourier.H"
 #include "fvmLaplacian.H"
+#include "fvcSnGrad.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,17 +71,18 @@ bool unityLewisFourier<BasicThermophysicalTransportModel>::read()
 
 
 template<class BasicThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 unityLewisFourier<BasicThermophysicalTransportModel>::q() const
 {
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->thermo().alpha()*this->alpha()*fvc::grad(this->thermo().he())
+       -fvc::interpolate(this->thermo().alpha()*this->alpha())
+       *fvc::snGrad(this->thermo().he())
     );
 }
 
@@ -94,19 +97,20 @@ divq(volScalarField& he) const
 
 
 template<class BasicThermophysicalTransportModel>
-tmp<volVectorField>unityLewisFourier<BasicThermophysicalTransportModel>::j
+tmp<surfaceScalarField>unityLewisFourier<BasicThermophysicalTransportModel>::j
 (
     const volScalarField& Yi
 ) const
 {
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "j(" + Yi.name() + ')',
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->thermo().alpha()*this->alpha()*fvc::grad(Yi)
+       -fvc::interpolate(this->thermo().alpha()*this->alpha())
+       *fvc::snGrad(Yi)
     );
 }
 
diff --git a/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.H b/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.H
index a80a039cd5..30ee9a7cba 100644
--- a/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.H
+++ b/src/ThermophysicalTransportModels/laminar/unityLewisFourier/unityLewisFourier.H
@@ -119,13 +119,13 @@ public:
         }
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const;
+        virtual tmp<surfaceScalarField> q() const;
 
         //- 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;
+        virtual tmp<surfaceScalarField> j(const volScalarField& Yi) const;
 
         //- Return the source term for the given specie mass-fraction equation
         virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
diff --git a/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H b/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
index e7fd06c587..3961e35c25 100644
--- a/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
+++ b/src/ThermophysicalTransportModels/thermophysicalTransportModel/thermophysicalTransportModel.H
@@ -128,13 +128,13 @@ public:
         ) const = 0;
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const = 0;
+        virtual tmp<surfaceScalarField> q() const = 0;
 
         //- 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;
+        virtual tmp<surfaceScalarField> 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;
diff --git a/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.C b/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.C
index 1fd4f177ed..96d54dccda 100644
--- a/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.C
@@ -137,17 +137,18 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::DEff
 
 
 template<class TurbulenceThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 eddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
 {
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -(this->alpha()*this->kappaEff()*fvc::grad(this->thermo().T()))
+       -fvc::interpolate(this->alpha()*this->kappaEff())
+       *fvc::snGrad(this->thermo().T())
     );
 }
 
@@ -168,7 +169,7 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
 
 
 template<class TurbulenceThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 eddyDiffusivity<TurbulenceThermophysicalTransportModel>::j
 (
     const volScalarField& Yi
@@ -180,7 +181,7 @@ eddyDiffusivity<TurbulenceThermophysicalTransportModel>::j
            " nonUnityLewisEddyDiffusivity or unityLewisEddyDiffusivity"
         << exit(FatalError);
 
-    return tmp<volVectorField>(nullptr);
+    return tmp<surfaceScalarField>(nullptr);
 }
 
 
diff --git a/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.H b/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.H
index 28bdc4136b..b4c9c615d0 100644
--- a/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.H
+++ b/src/ThermophysicalTransportModels/turbulence/eddyDiffusivity/eddyDiffusivity.H
@@ -182,13 +182,13 @@ public:
         ) const;
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const;
+        virtual tmp<surfaceScalarField> q() const;
 
         //- 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;
+        virtual tmp<surfaceScalarField> j(const volScalarField& Yi) const;
 
         //- Return the source term for the given specie mass-fraction equation
         virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
diff --git a/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C b/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
index 018261cf53..3ba5212454 100644
--- a/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.C
@@ -91,19 +91,20 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
 
 
 template<class TurbulenceThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
 {
-    tmp<volVectorField> tmpq
+    tmp<surfaceScalarField> tmpq
     (
-        volVectorField::New
+        surfaceScalarField::New
         (
             IOobject::groupName
             (
                 "q",
                 this->momentumTransport().alphaRhoPhi().group()
             ),
-           -(this->alpha()*this->kappaEff()*fvc::grad(this->thermo().T()))
+           -fvc::interpolate(this->alpha()*this->kappaEff())
+           *fvc::snGrad(this->thermo().T())
         )
     );
 
@@ -115,9 +116,12 @@ nonUnityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
         forAll(Y, i)
         {
             tmpq.ref() -=
-                this->alpha()*DEff(Y[i])
-               *composition.HE(i, this->thermo().p(), this->thermo().T())
-               *fvc::grad(Y[i]);
+                fvc::interpolate
+                (
+                    this->alpha()*DEff(Y[i])
+                   *composition.HE(i, this->thermo().p(), this->thermo().T())
+                )
+               *fvc::snGrad(Y[i]);
         }
     }
 
diff --git a/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.H b/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.H
index 1b36ffbcc2..30e75e4e90 100644
--- a/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.H
+++ b/src/ThermophysicalTransportModels/turbulence/nonUnityLewisEddyDiffusivity/nonUnityLewisEddyDiffusivity.H
@@ -144,7 +144,7 @@ public:
         }
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const;
+        virtual tmp<surfaceScalarField> q() const;
 
         //- Return the source term for the energy equation
         virtual tmp<fvScalarMatrix> divq(volScalarField& he) const;
diff --git a/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.C b/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.C
index 364d3d958a..c29c7a012e 100644
--- a/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.C
+++ b/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.C
@@ -138,17 +138,18 @@ bool unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::read()
 
 
 template<class TurbulenceThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::q() const
 {
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "q",
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->alphaEff()*this->alpha()*fvc::grad(this->thermo().he())
+       -fvc::interpolate(this->alphaEff()*this->alpha())
+       *fvc::snGrad(this->thermo().he())
     );
 }
 
@@ -165,20 +166,20 @@ unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::divq
 
 
 template<class TurbulenceThermophysicalTransportModel>
-tmp<volVectorField>
+tmp<surfaceScalarField>
 unityLewisEddyDiffusivity<TurbulenceThermophysicalTransportModel>::j
 (
     const volScalarField& Yi
 ) const
 {
-    return volVectorField::New
+    return surfaceScalarField::New
     (
         IOobject::groupName
         (
             "j(" + Yi.name() + ')',
             this->momentumTransport().alphaRhoPhi().group()
         ),
-       -this->DEff(Yi)*this->alpha()*fvc::grad(Yi)
+       -fvc::interpolate(this->DEff(Yi)*this->alpha())*fvc::snGrad(Yi)
     );
 }
 
diff --git a/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.H b/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.H
index 9f9c864f93..727f2f38a1 100644
--- a/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.H
+++ b/src/ThermophysicalTransportModels/turbulence/unityLewisEddyDiffusivity/unityLewisEddyDiffusivity.H
@@ -200,13 +200,13 @@ public:
         }
 
         //- Return the heat flux
-        virtual tmp<volVectorField> q() const;
+        virtual tmp<surfaceScalarField> q() const;
 
         //- 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;
+        virtual tmp<surfaceScalarField> j(const volScalarField& Yi) const;
 
         //- Return the source term for the given specie mass-fraction equation
         virtual tmp<fvScalarMatrix> divj(volScalarField& Yi) const;
diff --git a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
index 89b0c0b492..0aa981bada 100644
--- a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
+++ b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.C
@@ -58,52 +58,6 @@ void Foam::functionObjects::wallHeatFlux::writeFileHeader(const label i)
 }
 
 
-Foam::tmp<Foam::volScalarField>
-Foam::functionObjects::wallHeatFlux::calcWallHeatFlux(const volVectorField& q)
-{
-    tmp<volScalarField> twallHeatFlux
-    (
-        volScalarField::New
-        (
-            type(),
-            mesh_,
-            dimensionedScalar(dimMass/pow3(dimTime), 0)
-        )
-    );
-
-    volScalarField::Boundary& wallHeatFluxBf =
-        twallHeatFlux.ref().boundaryFieldRef();
-
-    const volVectorField::Boundary& qBf = q.boundaryField();
-
-    forAllConstIter(labelHashSet, patchSet_, iter)
-    {
-        const label patchi = iter.key();
-
-        const vectorField& Sfp = mesh_.Sf().boundaryField()[patchi];
-        const scalarField& magSfp = mesh_.magSf().boundaryField()[patchi];
-
-        wallHeatFluxBf[patchi] = (-Sfp/magSfp) & qBf[patchi];
-    }
-
-    if (foundObject<volScalarField>("qr"))
-    {
-        const volScalarField& qr = lookupObject<volScalarField>("qr");
-
-        const volScalarField::Boundary& radHeatFluxBf = qr.boundaryField();
-
-        forAllConstIter(labelHashSet, patchSet_, iter)
-        {
-            const label patchi = iter.key();
-
-            wallHeatFluxBf[patchi] -= radHeatFluxBf[patchi];
-        }
-    }
-
-    return twallHeatFlux;
-}
-
-
 Foam::tmp<Foam::volScalarField>
 Foam::functionObjects::wallHeatFlux::calcWallHeatFlux
 (
diff --git a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
index aaaaca21d4..33799d2a1f 100644
--- a/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
+++ b/src/functionObjects/field/wallHeatFlux/wallHeatFlux.H
@@ -108,9 +108,6 @@ protected:
         //- File header information
         virtual void writeFileHeader(const label i);
 
-        //- Calculate the heat-flux
-        tmp<volScalarField> calcWallHeatFlux(const volVectorField& q);
-
         //- Calculate the heat-flux
         tmp<volScalarField> calcWallHeatFlux(const surfaceScalarField& q);
 
diff --git a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
index 72dfe70bba..06003555d1 100644
--- a/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
+++ b/src/thermophysicalModels/solidThermo/solidThermo/heSolidThermo.C
@@ -325,6 +325,9 @@ template<class BasicSolidThermo, class MixtureType>
 Foam::tmp<Foam::surfaceScalarField>
 Foam::heSolidThermo<BasicSolidThermo, MixtureType>::q() const
 {
+    const fvMesh& mesh = this->T_.mesh();
+    mesh.setFluxRequired(this->T_.name());
+
     return
       - (
             isotropic()