diff --git a/YEqn.H b/YEqn.H
index 2552408..07b0699 100644
--- a/YEqn.H
+++ b/YEqn.H
@@ -18,11 +18,18 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     // volVectorField Vc(U);
     Vc = dimensionedVector("zero", dimLength/dimTime, U[0]*0.0);
 
+    phiC = dimensionedScalar("zero", dimDensity*dimVelocity*dimArea, 0.0);
+
     forAll(Y, i)
     {
         Vc += diff.D(i) * fvc::grad(Y[i]);
     }
 
+    forAll(Y, i)
+    {
+        phiC += linearInterpolate(rho * diff.D(i)) * fvc::snGrad(Y[i]) * mesh.magSf();
+    }
+
     forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)
diff --git a/createFields.H b/createFields.H
index d6f13dc..b69d6a9 100644
--- a/createFields.H
+++ b/createFields.H
@@ -141,3 +141,17 @@ volScalarField dQ
 #include "createMRF.H"
 
 diff.correct();
+
+
+surfaceScalarField phiC
+(
+    IOobject
+    (
+        "phiC",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*U) & mesh.Sf()
+);