diff --git a/dumux/discretization/cellcentered/tpfa/fickslaw.hh b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
index 4ebd1fad13fa777bf370eb44de47a22ec07334e9..202c3d90021a88e5d2e23041f533380a26beb616 100644
--- a/dumux/discretization/cellcentered/tpfa/fickslaw.hh
+++ b/dumux/discretization/cellcentered/tpfa/fickslaw.hh
@@ -79,7 +79,8 @@ public:
                        const ElementVolumeVariables& elemVolVars,
                        const SubControlVolumeFace& scvFace,
                        int phaseIdx, int compIdx,
-                       const FluxVarsCache& fluxVarsCache)
+                       const FluxVarsCache& fluxVarsCache,
+                       bool useMoles = true)
     {
         // diffusion tensors are always solution dependent
         Scalar tij = calculateTransmissibility_(problem, element, fvGeometry, elemVolVars, scvFace, phaseIdx, compIdx);
@@ -91,12 +92,24 @@ public:
         // and the outside volume variables
         const auto& outsideVolVars = elemVolVars[scvFace.outsideScvIdx()];
 
-        // compute the diffusive flux
-        const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
-        const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
-        const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
+        // compute the diffusive flux using mole fractions
+        if (useMoles)
+        {
+            const auto xInside = insideVolVars.moleFraction(phaseIdx, compIdx);
+            const auto xOutside = outsideVolVars.moleFraction(phaseIdx, compIdx);
+            const auto rho = 0.5*(insideVolVars.molarDensity(phaseIdx) + outsideVolVars.molarDensity(phaseIdx));
+
+            return rho*tij*(xInside - xOutside);
+        }
+        // compute the diffusive flux using mass fractions
+        else
+        {
+            const auto xInside = insideVolVars.massFraction(phaseIdx, compIdx);
+            const auto xOutside = outsideVolVars.massFraction(phaseIdx, compIdx);
+            const auto rho = 0.5*(insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx));
 
-        return rho*tij*(xInside - xOutside);
+            return rho*tij*(xInside - xOutside);
+        }
     }
 
     static Stencil stencil(const Problem& problem,