diff --git a/dumux/material/constraintsolvers/compositionfromfugacities2pncmin.hh b/dumux/material/constraintsolvers/compositionfromfugacities2pncmin.hh
index 1596b9ef18fffa288b10f2f90a95992c6d9aa626..dcb295fcdcde07acd2d7f8740cee9adc2a3ca02d 100644
--- a/dumux/material/constraintsolvers/compositionfromfugacities2pncmin.hh
+++ b/dumux/material/constraintsolvers/compositionfromfugacities2pncmin.hh
@@ -144,8 +144,8 @@ protected:
 };
 
 template <class Scalar, class FluidSystem>
-DUNE_DEPRECATED_MSG("compositionFromFugacities2pncmin is deprecated. Use CompositionFromFugacities2pncmin (capital C) instead.")
-class compositionFromFugacities2pncmin
+class DUNE_DEPRECATED_MSG("compositionFromFugacities2pncmin is deprecated. Use CompositionFromFugacities2pncmin (capital C) instead.")
+  compositionFromFugacities2pncmin
   : public CompositionFromFugacities2pncmin<Scalar, FluidSystem>
 { };
 } // end namespace Dumux
diff --git a/dumux/porousmediumflow/2pnc/implicit/model.hh b/dumux/porousmediumflow/2pnc/implicit/model.hh
index 8b7c2fef46bd20e6aaf29db1c37eccfdb529ec8d..cfea62d30d5b08f5e40a03fff1fbb5c3f8dd8682 100644
--- a/dumux/porousmediumflow/2pnc/implicit/model.hh
+++ b/dumux/porousmediumflow/2pnc/implicit/model.hh
@@ -275,15 +275,15 @@ public:
         // get the number of degrees of freedom
         auto numDofs = this->numDofs();
 
-        ScalarField *Sg            = writer.allocateManagedBuffer (numDofs);
-        ScalarField *Sl            = writer.allocateManagedBuffer (numDofs);
-        ScalarField *pg            = writer.allocateManagedBuffer (numDofs);
-        ScalarField *pl            = writer.allocateManagedBuffer (numDofs);
+        ScalarField *sn            = writer.allocateManagedBuffer (numDofs);
+        ScalarField *sw            = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pn            = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pw            = writer.allocateManagedBuffer (numDofs);
         ScalarField *pc            = writer.allocateManagedBuffer (numDofs);
-        ScalarField *rhoL          = writer.allocateManagedBuffer (numDofs);
-        ScalarField *rhoG          = writer.allocateManagedBuffer (numDofs);
-        ScalarField *mobL          = writer.allocateManagedBuffer (numDofs);
-        ScalarField *mobG          = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoW          = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoN          = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobW          = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobN          = writer.allocateManagedBuffer (numDofs);
         ScalarField *temperature   = writer.allocateManagedBuffer (numDofs);
         ScalarField *poro          = writer.allocateManagedBuffer (numDofs);
         VectorField *velocityN = writer.template allocateManagedBuffer<double, dim>(numDofs);
@@ -334,15 +334,15 @@ public:
                 auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
                 GlobalPosition globalPos = fvGeometry.subContVol[scvIdx].global;
-                (*Sg)[dofIdxGlobal]             = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                (*Sl)[dofIdxGlobal]             = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                (*pg)[dofIdxGlobal]             = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                (*pl)[dofIdxGlobal]             = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                (*sn)[dofIdxGlobal]             = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                (*sw)[dofIdxGlobal]             = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                (*pn)[dofIdxGlobal]             = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                (*pw)[dofIdxGlobal]             = elemVolVars[scvIdx].pressure(wPhaseIdx);
                 (*pc)[dofIdxGlobal]             = elemVolVars[scvIdx].capillaryPressure();
-                (*rhoL)[dofIdxGlobal]           = elemVolVars[scvIdx].density(wPhaseIdx);
-                (*rhoG)[dofIdxGlobal]           = elemVolVars[scvIdx].density(nPhaseIdx);
-                (*mobL)[dofIdxGlobal]           = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                (*mobG)[dofIdxGlobal]           = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                (*rhoW)[dofIdxGlobal]           = elemVolVars[scvIdx].density(wPhaseIdx);
+                (*rhoN)[dofIdxGlobal]           = elemVolVars[scvIdx].density(nPhaseIdx);
+                (*mobW)[dofIdxGlobal]           = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                (*mobN)[dofIdxGlobal]           = elemVolVars[scvIdx].mobility(nPhaseIdx);
                 (*poro)[dofIdxGlobal]           = elemVolVars[scvIdx].porosity();
                 (*temperature)[dofIdxGlobal]    = elemVolVars[scvIdx].temperature();
 
@@ -367,15 +367,15 @@ public:
 
         } // loop over element
 
-        writer.attachDofData(*Sg, "Sg", isBox);
-        writer.attachDofData(*Sl, "Sl", isBox);
-        writer.attachDofData(*pg, "pg", isBox);
-        writer.attachDofData(*pl, "pl", isBox);
+        writer.attachDofData(*sn, "Sn", isBox);
+        writer.attachDofData(*sw, "Sw", isBox);
+        writer.attachDofData(*pn, "pn", isBox);
+        writer.attachDofData(*pw, "pw", isBox);
         writer.attachDofData(*pc, "pc", isBox);
-        writer.attachDofData(*rhoL, "rhoL", isBox);
-        writer.attachDofData(*rhoG, "rhoG", isBox);
-        writer.attachDofData(*mobL, "mobL", isBox);
-        writer.attachDofData(*mobG, "mobG", isBox);
+        writer.attachDofData(*rhoW, "rhoW", isBox);
+        writer.attachDofData(*rhoN, "rhoN", isBox);
+        writer.attachDofData(*mobW, "mobW", isBox);
+        writer.attachDofData(*mobN, "mobN", isBox);
         writer.attachDofData(*poro, "porosity", isBox);
         writer.attachDofData(*temperature, "temperature", isBox);
         writer.attachDofData(*perm[0], "Kxx", isBox);
@@ -583,93 +583,92 @@ protected:
                 if (staticDat_[dofIdxGlobal].wasSwitched)
                     Smin = -0.01;
 
-                //if saturation of liquid phase is smaller 0 switch
+                //if saturation of wetting phase is smaller 0 switch
                 if (volVars.saturation(wPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
-                    //liquid phase has to disappear
-                    std::cout << "Liquid Phase disappears at vertex " << dofIdxGlobal
-                                << ", coordinated: " << globalPos << ", Sl: "
+                    //wetting phase has to disappear
+                    std::cout << "Wetting Phase disappears at vertex " << dofIdxGlobal
+                                << ", coordinated: " << globalPos << ", Sw: "
                                 << volVars.saturation(wPhaseIdx) << std::endl;
                     newPhasePresence = nPhaseOnly;
 
-                    //switch not depending on formulation
-                    //switch "Sl" to "xgH20"
+                    //switch saturation to xnH20 (not depending on formulation)
                     globalSol[dofIdxGlobal][switchIdx]
                             = volVars.moleFraction(nPhaseIdx, wCompIdx /*H2O*/);
 
-                    //switch all secondary components to mole fraction in gas phase
+                    //switch all secondary components to mole fraction in nonwetting phase
                     for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
                         globalSol[dofIdxGlobal][compIdx] = volVars.moleFraction(nPhaseIdx,compIdx);
                 }
-                //if saturation of gas phase is smaller than 0 switch
+                //if saturation of nonwetting phase is smaller than 0 switch
                 else if (volVars.saturation(nPhaseIdx) <= Smin)
                 {
                     wouldSwitch = true;
-                    //gas phase has to disappear
-                    std::cout << "Gas Phase disappears at vertex " << dofIdxGlobal
-                                << ", coordinated: " << globalPos << ", Sg: "
+                    //nonwetting phase has to disappear
+                    std::cout << "Nonwetting Phase disappears at vertex " << dofIdxGlobal
+                                << ", coordinated: " << globalPos << ", Sn: "
                                 << volVars.saturation(nPhaseIdx) << std::endl;
                     newPhasePresence = wPhaseOnly;
 
-                    //switch "Sl" to "xlN2"
+                    //switch saturation to xwN2, not depending on formulation
                     globalSol[dofIdxGlobal][switchIdx]
                             = volVars.moleFraction(wPhaseIdx, nCompIdx /*N2*/);
                 }
             }
             else if (phasePresence == nPhaseOnly)
             {
-                Scalar xlmax = 1;
-                Scalar sumxl = 0;
-                //Calculate sum of mole fractions in the hypothetical liquid phase
+                Scalar xwmax = 1;
+                Scalar sumxw = 0;
+                //Calculate sum of mole fractions in the hypothetical wetting phase
                 for (int compIdx = 0; compIdx < numComponents; compIdx++)
                 {
-                    sumxl += volVars.moleFraction(wPhaseIdx, compIdx);
+                    sumxw += volVars.moleFraction(wPhaseIdx, compIdx);
                 }
-                if (sumxl > xlmax)
+                if (sumxw > xwmax)
                     wouldSwitch = true;
                 if (staticDat_[dofIdxGlobal].wasSwitched)
-                    xlmax *=1.02;
-                //liquid phase appears if sum is larger than one
-                if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
+                    xwmax *=1.02;
+                //wetting phase appears if sum is larger than one
+                if (sumxw/*sum of mole fractions*/ > xwmax/*1*/)
                 {
-                    std::cout << "Liquid Phase appears at vertex " << dofIdxGlobal
-                            << ", coordinated: " << globalPos << ", sumxl: "
-                            << sumxl << std::endl;
+                    std::cout << "Wetting Phase appears at vertex " << dofIdxGlobal
+                            << ", coordinated: " << globalPos << ", sumxw: "
+                            << sumxw << std::endl;
                     newPhasePresence = bothPhases;
 
-                    //saturation of the liquid phase set to 0.0001 (if formulation pnsw and vice versa)
+                    //saturation of the wetting phase set to 0.0001
                     if (formulation == pnsw)
                         globalSol[dofIdxGlobal][switchIdx] = 0.0001;
                     else if (formulation == pwsn)
                         globalSol[dofIdxGlobal][switchIdx] = 0.9999;
 
-                    //switch all secondary components back to liquid mole fraction
+                    //switch all secondary components back to wetting mole fraction
                     for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
                         globalSol[dofIdxGlobal][compIdx] = volVars.moleFraction(wPhaseIdx,compIdx);
                 }
             }
             else if (phasePresence == wPhaseOnly)
             {
-                Scalar xgmax = 1;
-                Scalar sumxg = 0;
-                //Calculate sum of mole fractions in the hypothetical liquid phase
+                Scalar xnmax = 1;
+                Scalar sumxn = 0;
+                //Calculate sum of mole fractions in the hypothetical nonwetting phase
                 for (int compIdx = 0; compIdx < numComponents; compIdx++)
                 {
-                    sumxg += volVars.moleFraction(nPhaseIdx, compIdx);
+                    sumxn += volVars.moleFraction(nPhaseIdx, compIdx);
                 }
-                if (sumxg > xgmax)
+                if (sumxn > xnmax)
                     wouldSwitch = true;
                 if (staticDat_[dofIdxGlobal].wasSwitched)
-                    xgmax *=1.02;
-                //liquid phase appears if sum is larger than one
-                if (sumxg > xgmax)
+                    xnmax *=1.02;
+                //nonwetting phase appears if sum is larger than one
+                if (sumxn > xnmax)
                 {
-                    std::cout << "Gas Phase appears at vertex " << dofIdxGlobal
-                            << ", coordinated: " << globalPos << ", sumxg: "
-                            << sumxg << std::endl;
+                    std::cout << "Nonwetting Phase appears at vertex " << dofIdxGlobal
+                            << ", coordinated: " << globalPos << ", sumxn: "
+                            << sumxn << std::endl;
                     newPhasePresence = bothPhases;
-                    //saturation of the liquid phase set to 0.9999 (if formulation pnsw and vice versa)
+                    //saturation of the wetting phase set to 0.9999
                     if (formulation == pnsw)
                         globalSol[dofIdxGlobal][switchIdx] = 0.9999;
                     else if (formulation == pwsn)
diff --git a/dumux/porousmediumflow/2pnc/implicit/propertydefaults.hh b/dumux/porousmediumflow/2pnc/implicit/propertydefaults.hh
index 36c5f83bc91f409a772af3aaea9c5c687d447c5d..3cf9509b4ede48f8f2638e6f50631e48796dc469 100644
--- a/dumux/porousmediumflow/2pnc/implicit/propertydefaults.hh
+++ b/dumux/porousmediumflow/2pnc/implicit/propertydefaults.hh
@@ -72,7 +72,7 @@ private:
 public:
     static const int value = FluidSystem::numComponents;
 };
-//! The major components belonging to the existing phases are mentioned here e.g., 2 for water and air being the major component in the liquid and gas phases in a 2 phase system
+//! The major components belonging to the existing phases, e.g. 2 for water and air being the major components in a liquid-gas-phase system
 SET_PROP(TwoPNC, NumMajorComponents)
 {
 private:
@@ -81,7 +81,7 @@ private:
 public:
     static const int value = FluidSystem::numPhases;
     static_assert(value == 2,
-                  "The model is restricted to two-phases, thus number of major components must also be two.");
+                  "The model is restricted to two phases, thus number of major components must also be two.");
 };
 
 /*!
diff --git a/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh b/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh
index 8867840502e757b648461aafc8681bd515b234c0..23dc3348f763ab039bd070e6eec9cef86c9e450f 100644
--- a/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2pnc/implicit/volumevariables.hh
@@ -142,7 +142,7 @@ public:
             Scalar kr;
             if (phaseIdx == wPhaseIdx)
                 kr = MaterialLaw::krw(materialParams, saturation(wPhaseIdx));
-            else // ATTENTION: krn requires the liquid saturation // as parameter!
+            else // ATTENTION: krn requires the wetting-phase saturation as parameter!
                 kr = MaterialLaw::krn(materialParams, saturation(wPhaseIdx));
 
             mobility_[phaseIdx] = kr / fluidState_.viscosity(phaseIdx);
@@ -254,7 +254,7 @@ public:
         // now comes the tricky part: calculate phase composition
         if (phasePresence == bothPhases) {
             // both phases are present, phase composition results from
-            // the gas <-> liquid equilibrium. This is
+            // the nonwetting <-> wetting equilibrium. This is
             // the job of the "MiscibleMultiPhaseComposition"
             // constraint solver
 
@@ -282,12 +282,12 @@ public:
                     moleFrac[compIdx] = priVars[compIdx];
 
 
-            Scalar sumMoleFracNotGas = 0;
+            Scalar sumMoleFracOtherComponents = 0;
             for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
-                    sumMoleFracNotGas+=moleFrac[compIdx];
+                    sumMoleFracOtherComponents+=moleFrac[compIdx];
 
-            sumMoleFracNotGas += moleFrac[wCompIdx];
-            moleFrac[nCompIdx] = 1 - sumMoleFracNotGas;
+            sumMoleFracOtherComponents += moleFrac[wCompIdx];
+            moleFrac[nCompIdx] = 1 - sumMoleFracOtherComponents;
 
 
             // Set fluid state mole fractions
@@ -306,9 +306,9 @@ public:
 
         }
         else if (phasePresence == wPhaseOnly){
-        // only the liquid phase is present, i.e. liquid phase
+        // only the wetting phase is present, i.e. wetting phase
         // composition is stored explicitly.
-        // extract _mass_ fractions in the gas phase
+        // extract _mass_ fractions in the nonwetting phase
             Dune::FieldVector<Scalar, numComponents> moleFrac;
 
             for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
diff --git a/dumux/porousmediumflow/2pncmin/implicit/model.hh b/dumux/porousmediumflow/2pncmin/implicit/model.hh
index e7d3be552b9616b7da84941ab9e92ca8e9c31670..59884fc075e526bd44d9c6c5be3b26b1763e1bf6 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/model.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/model.hh
@@ -186,15 +186,15 @@ public:
         auto numDofs = this->numDofs();
 
         // create the required scalar fields
-        ScalarField *Sg           = writer.allocateManagedBuffer(numDofs);
-        ScalarField *Sl           = writer.allocateManagedBuffer(numDofs);
-        ScalarField *pg           = writer.allocateManagedBuffer (numDofs);
-        ScalarField *pl           = writer.allocateManagedBuffer (numDofs);
+        ScalarField *sn           = writer.allocateManagedBuffer(numDofs);
+        ScalarField *sw           = writer.allocateManagedBuffer(numDofs);
+        ScalarField *pn           = writer.allocateManagedBuffer (numDofs);
+        ScalarField *pw           = writer.allocateManagedBuffer (numDofs);
         ScalarField *pc           = writer.allocateManagedBuffer (numDofs);
-        ScalarField *rhoL         = writer.allocateManagedBuffer (numDofs);
-        ScalarField *rhoG         = writer.allocateManagedBuffer (numDofs);
-        ScalarField *mobL         = writer.allocateManagedBuffer (numDofs);
-        ScalarField *mobG         = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoW         = writer.allocateManagedBuffer (numDofs);
+        ScalarField *rhoN         = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobW         = writer.allocateManagedBuffer (numDofs);
+        ScalarField *mobN         = writer.allocateManagedBuffer (numDofs);
         ScalarField *phasePresence = writer.allocateManagedBuffer (numDofs);
         ScalarField *temperature  = writer.allocateManagedBuffer (numDofs);
         ScalarField *poro         = writer.allocateManagedBuffer (numDofs);
@@ -251,15 +251,15 @@ public:
             {
                 auto dofIdxGlobal = this->dofMapper().subIndex(element, scvIdx, dofCodim);
 
-                (*Sg)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
-                (*Sl)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
-                (*pg)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
-                (*pl)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
+                (*sn)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(nPhaseIdx);
+                (*sw)[dofIdxGlobal] = elemVolVars[scvIdx].saturation(wPhaseIdx);
+                (*pn)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(nPhaseIdx);
+                (*pw)[dofIdxGlobal] = elemVolVars[scvIdx].pressure(wPhaseIdx);
                 (*pc)[dofIdxGlobal] = elemVolVars[scvIdx].capillaryPressure();
-                (*rhoL)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
-                (*rhoG)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
-                (*mobL)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
-                (*mobG)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
+                (*rhoW)[dofIdxGlobal] = elemVolVars[scvIdx].density(wPhaseIdx);
+                (*rhoN)[dofIdxGlobal] = elemVolVars[scvIdx].density(nPhaseIdx);
+                (*mobW)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(wPhaseIdx);
+                (*mobN)[dofIdxGlobal] = elemVolVars[scvIdx].mobility(nPhaseIdx);
                 (*poro)[dofIdxGlobal] = elemVolVars[scvIdx].porosity();
 
                 for (int sPhaseIdx = 0; sPhaseIdx < numSPhases; ++sPhaseIdx)
@@ -289,15 +289,15 @@ public:
             }
         }  // loop over element
 
-        writer.attachDofData(*Sg, "Sg", isBox);
-        writer.attachDofData(*Sl, "Sl", isBox);
-        writer.attachDofData(*pg, "pg", isBox);
-        writer.attachDofData(*pl, "pl", isBox);
+        writer.attachDofData(*sn, "Sn", isBox);
+        writer.attachDofData(*sw, "Sw", isBox);
+        writer.attachDofData(*pn, "pn", isBox);
+        writer.attachDofData(*pw, "pw", isBox);
         writer.attachDofData(*pc, "pc", isBox);
-        writer.attachDofData(*rhoL, "rhoL", isBox);
-        writer.attachDofData(*rhoG, "rhoG", isBox);
-        writer.attachDofData(*mobL, "mobL", isBox);
-        writer.attachDofData(*mobG, "mobG", isBox);
+        writer.attachDofData(*rhoW, "rhoW", isBox);
+        writer.attachDofData(*rhoN, "rhoN", isBox);
+        writer.attachDofData(*mobW, "mobW", isBox);
+        writer.attachDofData(*mobN, "mobN", isBox);
         writer.attachDofData(*poro, "porosity", isBox);
         writer.attachDofData(*permeabilityFactor, "permeabilityFactor", isBox);
         writer.attachDofData(*temperature, "temperature", isBox);
@@ -412,87 +412,86 @@ protected:
             if (this->staticDat_[dofIdxGlobal].wasSwitched)
                 Smin = -0.01;
 
-            //if saturation of liquid phase is smaller 0 switch
+            //if saturation of wetting phase is smaller 0 switch
             if (volVars.saturation(wPhaseIdx) <= Smin)
             {
                 wouldSwitch = true;
-                //liquid phase has to disappear
-                std::cout << "Liquid Phase disappears at vertex " << dofIdxGlobal
-                            << ", coordinated: " << globalPos << ", Sl: "
+                //wetting phase has to disappear
+                std::cout << "Wetting Phase disappears at vertex " << dofIdxGlobal
+                            << ", coordinated: " << globalPos << ", Sw: "
                             << volVars.saturation(wPhaseIdx) << std::endl;
                 newPhasePresence = nPhaseOnly;
 
-                //switch not depending on formulation
-                //switch "Sl" to "xgH20"
+                //switch saturation to xnH20 (not depending on formulation)
                 globalSol[dofIdxGlobal][switchIdx]
                         = volVars.moleFraction(nPhaseIdx, wCompIdx /*H2O*/);
-                //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+                //Here unlike 2pnc model we do not switch all components to to mole fraction in nonwetting phase
             }
-            //if saturation of gas phase is smaller than 0 switch
+            //if saturation of nonwetting phase is smaller than 0 switch
             else if (volVars.saturation(nPhaseIdx) <= Smin)
             {
                 wouldSwitch = true;
-                //gas phase has to disappear
-                std::cout << "Gas Phase disappears at vertex " << dofIdxGlobal
-                            << ", coordinated: " << globalPos << ", Sg: "
+                //nonwetting phase has to disappear
+                std::cout << "Nonwetting Phase disappears at vertex " << dofIdxGlobal
+                            << ", coordinated: " << globalPos << ", Sn: "
                             << volVars.saturation(nPhaseIdx) << std::endl;
                 newPhasePresence = wPhaseOnly;
 
-                //switch "Sl" to "xlN2"
+                //switch saturation to xwN2 (not depending on formulation)
                 globalSol[dofIdxGlobal][switchIdx] = volVars.moleFraction(wPhaseIdx, nCompIdx /*N2*/);
             }
         }
         else if (phasePresence == nPhaseOnly)
         {
-            Scalar sumxl = 0;
-            //Calculate sum of mole fractions (water and air) in the hypothetical liquid phase
+            Scalar sumxw = 0;
+            //Calculate sum of mole fractions (water and air) in the hypothetical wetting phase
             for (int compIdx = 0; compIdx < numComponents; compIdx++)
             {
-                sumxl += volVars.moleFraction(wPhaseIdx, compIdx);
+                sumxw += volVars.moleFraction(wPhaseIdx, compIdx);
             }
-            Scalar xlmax = 1.0;
-            if (sumxl > xlmax)
+            Scalar xwmax = 1.0;
+            if (sumxw > xwmax)
                 wouldSwitch = true;
             if (this->staticDat_[dofIdxGlobal].wasSwitched)
-                xlmax *=1.02;
+                xwmax *=1.02;
 
             //if the sum of the mole fractions would be larger than
             //1, wetting phase appears
-            if (sumxl/*sum of mole fractions*/ > xlmax/*1*/)
+            if (sumxw/*sum of mole fractions*/ > xwmax/*1*/)
             {
-                // liquid phase appears
-                std::cout << "Liquid Phase appears at vertex " << dofIdxGlobal
-                          << ", coordinated: " << globalPos << ", sumxl: "
-                          << sumxl << std::endl;
+                // wetting phase appears
+                std::cout << "Wetting Phase appears at vertex " << dofIdxGlobal
+                          << ", coordinated: " << globalPos << ", sumxw: "
+                          << sumxw << std::endl;
                 newPhasePresence = bothPhases;
                 if (formulation == pnsw)
                     globalSol[dofIdxGlobal][switchIdx] = 0.0;
                 else if (formulation == pwsn)
                     globalSol[dofIdxGlobal][switchIdx] = 1.0;
-                //Here unlike 2pnc model we do not switch all components to to mole fraction in gas phase
+                //Here unlike 2pnc model we do not switch all components to to mole fraction in nonwetting phase
             }
         }
         else if (phasePresence == wPhaseOnly)
         {
-            Scalar xgmax = 1;
-            Scalar sumxg = 0;
-            //Calculate sum of mole fractions in the hypothetical gas phase
+            Scalar xnmax = 1;
+            Scalar sumxn = 0;
+            //Calculate sum of mole fractions in the hypothetical nonwetting phase
             for (int compIdx = 0; compIdx < numComponents; compIdx++)
             {
-                sumxg += volVars.moleFraction(nPhaseIdx, compIdx);
+                sumxn += volVars.moleFraction(nPhaseIdx, compIdx);
             }
-            if (sumxg > xgmax)
+            if (sumxn > xnmax)
                 wouldSwitch = true;
             if (this->staticDat_[dofIdxGlobal].wasSwitched)
-                xgmax *=1.02;
-            //liquid phase appears if sum is larger than one
-            if (sumxg > xgmax)
+                xnmax *=1.02;
+            //nonwetting phase appears if sum is larger than one
+            if (sumxn > xnmax)
             {
-                std::cout << "Gas Phase appears at vertex " << dofIdxGlobal
-                          << ", coordinated: " << globalPos << ", sumxg: "
-                          << sumxg << std::endl;
+                std::cout << "Nonwetting Phase appears at vertex " << dofIdxGlobal
+                          << ", coordinated: " << globalPos << ", sumxn: "
+                          << sumxn << std::endl;
                 newPhasePresence = bothPhases;
-                //saturation of the liquid phase set to 0.9999 (if formulation pnsw and vice versa)
+                //saturation of the wetting phase set to 0.9999 (if formulation pnsw and vice versa)
                 if (formulation == pnsw)
                     globalSol[dofIdxGlobal][switchIdx] = 0.999;
                 else if (formulation == pwsn)
diff --git a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
index 720c012deaed50aa54968f26fb8205128fc2bdbe..c96f161e13d0ad532d761f781b4b8de067b8e66e 100644
--- a/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
+++ b/dumux/porousmediumflow/2pncmin/implicit/volumevariables.hh
@@ -275,7 +275,7 @@ public:
         if (phasePresence == bothPhases)
         {
             // both phases are present, phase composition results from
-            // the gas <-> liquid equilibrium. This is
+            // the nonwetting <-> wetting equilibrium. This is
             // the job of the "MiscibleMultiPhaseComposition"
             // constraint solver
 
@@ -314,13 +314,13 @@ public:
                 /(fugCoeffG[compIdx]*fluidState.pressure(nPhaseIdx));
 
             moleFrac[wCompIdx] =  priVars[switchIdx];
-            Scalar sumMoleFracNotGas = 0;
+            Scalar sumMoleFracOtherComponents = 0;
             for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
             {
-                    sumMoleFracNotGas+=moleFrac[compIdx];
+                    sumMoleFracOtherComponents+=moleFrac[compIdx];
             }
-            sumMoleFracNotGas += moleFrac[wCompIdx];
-            moleFrac[nCompIdx] = 1 - sumMoleFracNotGas;
+            sumMoleFracOtherComponents += moleFrac[wCompIdx];
+            moleFrac[nCompIdx] = 1 - sumMoleFracOtherComponents;
 
             // Set fluid state mole fractions
             for (int compIdx=0; compIdx<numComponents; ++compIdx)
@@ -342,9 +342,9 @@ public:
         else if (phasePresence == wPhaseOnly)
         {
 
-            // only the liquid phase is present, i.e. liquid phase
+            // only the wetting phase is present, i.e. wetting phase
             // composition is stored explicitly.
-            // extract _mass_ fractions in the gas phase
+            // extract _mass_ fractions in the nonwetting phase
             Dune::FieldVector<Scalar, numComponents> moleFrac;
 
             for (int compIdx=numMajorComponents; compIdx<numComponents; ++compIdx)
@@ -414,7 +414,7 @@ public:
     { return permeabilityFactor_; }
 
     /*!
-     * \brief Returns the mole fraction of the salinity in the liquid phase
+     * \brief Returns the mole fraction of the salinity in the wetting phase
      */
     Scalar moleFracSalinity() const
     {
@@ -422,7 +422,7 @@ public:
     }
 
     /*!
-     * \brief Returns the salinity (mass fraction) in the liquid phase
+     * \brief Returns the salinity (mass fraction) in the wetting phase
      */
     Scalar salinity() const
     {
diff --git a/test/references/fuelcell2pncbox-reference.vtu b/test/references/fuelcell2pncbox-reference.vtu
index 80abd5597f0a53e04d8b3ba3d256a80b15c63be5..51349dfa14bf560ec3dcdf083b9860edafb5c54e 100644
--- a/test/references/fuelcell2pncbox-reference.vtu
+++ b/test/references/fuelcell2pncbox-reference.vtu
@@ -2,8 +2,8 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="126" NumberOfPoints="154">
-      <PointData Scalars="Sg">
-        <DataArray type="Float32" Name="Sg" NumberOfComponents="1" format="ascii">
+      <PointData Scalars="Sn">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0.699989 0.699989 0.699991 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991
           0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991
           0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991 0.699989 0.699991
@@ -18,7 +18,7 @@
           0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
           0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7 0.7
         </DataArray>
-        <DataArray type="Float32" Name="Sl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           0.300011 0.300011 0.300009 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009
           0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009
           0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009 0.300011 0.300009
@@ -33,7 +33,7 @@
           0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
           0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3 0.3
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
           99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8
           99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8
           99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8
@@ -48,7 +48,7 @@
           100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
           100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
         </DataArray>
-        <DataArray type="Float32" Name="pl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
@@ -78,7 +78,7 @@
           -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1
           -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1 -11440.1
         </DataArray>
-        <DataArray type="Float32" Name="rhoL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
@@ -93,7 +93,7 @@
           993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403
           993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403
         </DataArray>
-        <DataArray type="Float32" Name="rhoG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
           1.0976 1.0976 1.10151 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151
           1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151
           1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151 1.0976 1.10151
@@ -108,7 +108,7 @@
           1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907
           1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907 1.11907
         </DataArray>
-        <DataArray type="Float32" Name="mobL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
           129.817 129.817 129.816 129.816 129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816
           129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816
           129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816 129.817 129.816
@@ -123,7 +123,7 @@
           129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808
           129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808 129.808
         </DataArray>
-        <DataArray type="Float32" Name="mobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           18761 18761 18685.6 18685.6 18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6
           18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6
           18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6 18761 18685.6
@@ -168,51 +168,6 @@
           310 310 310 310 310 310 310 310 310 310 310 310
           310 310 310 310 310 310 310 310 310 310
         </DataArray>
-        <DataArray type="Float32" Name="reactionSourceH2O [mol/(sm^2)]" NumberOfComponents="1" format="ascii">
-          747.822 747.822 0 0 747.822 0 747.822 0 747.822 0 747.822 0
-          747.822 0 747.822 0 747.822 0 747.822 0 747.822 0 747.822 0
-          747.822 0 747.822 0 747.822 0 747.822 0 747.822 0 747.822 0
-          747.822 0 747.822 0 747.822 0 747.822 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0
-        </DataArray>
-        <DataArray type="Float32" Name="reactionSourceO2 [mol/(sm^2)]" NumberOfComponents="1" format="ascii">
-          -255.16 -255.16 0 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
-          -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
-          -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
-          -255.16 0 -255.16 0 -255.16 0 -255.16 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0
-        </DataArray>
-        <DataArray type="Float32" Name="currentDensity [A/cm^2]" NumberOfComponents="1" format="ascii">
-          0.451352 0.451352 0 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
-          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
-          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
-          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0 0 0
-          0 0 0 0 0 0 0 0 0 0
-        </DataArray>
         <DataArray type="Float32" Name="Kxx" NumberOfComponents="1" format="ascii">
           5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11
           5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11 5e-11
@@ -378,6 +333,51 @@
           0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375
           0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375 0.383375
         </DataArray>
+        <DataArray type="Float32" Name="reactionSourceH2O [mol/(sm^2)]" NumberOfComponents="1" format="ascii">
+          747.822 747.822 0 0 747.822 0 747.822 0 747.822 0 747.822 0
+          747.822 0 747.822 0 747.822 0 747.822 0 747.822 0 747.822 0
+          747.822 0 747.822 0 747.822 0 747.822 0 747.822 0 747.822 0
+          747.822 0 747.822 0 747.822 0 747.822 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="reactionSourceO2 [mol/(sm^2)]" NumberOfComponents="1" format="ascii">
+          -255.16 -255.16 0 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
+          -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
+          -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0 -255.16 0
+          -255.16 0 -255.16 0 -255.16 0 -255.16 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0
+        </DataArray>
+        <DataArray type="Float32" Name="currentDensity [A/cm^2]" NumberOfComponents="1" format="ascii">
+          0.451352 0.451352 0 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
+          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
+          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0 0.451352 0
+          0.451352 0 0.451352 0 0.451352 0 0.451352 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0 0 0
+          0 0 0 0 0 0 0 0 0 0
+        </DataArray>
       </PointData>
       <CellData Scalars="process rank">
         <DataArray type="Float32" Name="process rank" NumberOfComponents="1" format="ascii">
diff --git a/test/references/fuelcell2pnccc-reference.vtu b/test/references/fuelcell2pnccc-reference.vtu
index 2071e1e1a87c785ac31d40349c80b134daf9c9f8..fa2279344b1ef04ac32e8b17556a4d0a6554d675 100644
--- a/test/references/fuelcell2pnccc-reference.vtu
+++ b/test/references/fuelcell2pnccc-reference.vtu
@@ -2,8 +2,8 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="126" NumberOfPoints="154">
-      <CellData Scalars="Sg">
-        <DataArray type="Float32" Name="Sg" NumberOfComponents="1" format="ascii">
+      <CellData Scalars="Sn">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999
           0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.69999 0.699992 0.699992 0.699992
           0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992 0.699992
@@ -16,7 +16,7 @@
           0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999 0.699999
           0.699999 0.699999 0.699999 0.699999 0.699999 0.699999
         </DataArray>
-        <DataArray type="Float32" Name="Sl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001
           0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.30001 0.300008 0.300008 0.300008
           0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008 0.300008
@@ -29,7 +29,7 @@
           0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001 0.300001
           0.300001 0.300001 0.300001 0.300001 0.300001 0.300001
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
           99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8
           99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.8 99999.9 99999.9 99999.9
           99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9 99999.9
@@ -42,7 +42,7 @@
           100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000 100000
           100000 100000 100000 100000 100000 100000
         </DataArray>
-        <DataArray type="Float32" Name="pl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
           111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440 111440
@@ -68,7 +68,7 @@
           -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2
           -11440.2 -11440.2 -11440.2 -11440.2 -11440.2 -11440.2
         </DataArray>
-        <DataArray type="Float32" Name="rhoL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
           993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402 993.402
@@ -81,7 +81,7 @@
           993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403 993.403
           993.403 993.403 993.403 993.403 993.403 993.403
         </DataArray>
-        <DataArray type="Float32" Name="rhoG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
           1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924
           1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.09924 1.10315 1.10315 1.10315
           1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315 1.10315
@@ -94,7 +94,7 @@
           1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742 1.11742
           1.11742 1.11742 1.11742 1.11742 1.11742 1.11742
         </DataArray>
-        <DataArray type="Float32" Name="mobL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
           129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817
           129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.817 129.815 129.815 129.815
           129.815 129.815 129.815 129.815 129.815 129.815 129.815 129.815 129.815 129.815 129.815 129.815
@@ -107,7 +107,7 @@
           129.809 129.809 129.809 129.809 129.809 129.809 129.809 129.809 129.809 129.809 129.809 129.809
           129.809 129.809 129.809 129.809 129.809 129.809
         </DataArray>
-        <DataArray type="Float32" Name="mobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3
           18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18729.3 18654.1 18654.1 18654.1
           18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1 18654.1
diff --git a/test/references/saltflushbox2pncmin-reference.vtu b/test/references/saltflushbox2pncmin-reference.vtu
index 8d1ce3d9e2261db49764c0afdc773d378c614301..e2b9e4b3f36f3822987654e497b5af9a59e55bc3 100644
--- a/test/references/saltflushbox2pncmin-reference.vtu
+++ b/test/references/saltflushbox2pncmin-reference.vtu
@@ -2,8 +2,8 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="200" NumberOfPoints="231">
-      <PointData Scalars="Sg" Vectors="velocityW">
-        <DataArray type="Float32" Name="Sg" NumberOfComponents="1" format="ascii">
+      <PointData Scalars="Sn" Vectors="velocityW">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0.05 0.0302391 0.05 0.0450697 0.02192 0.0374154 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 0.0258503 0.0571775 0.061978 0.0899974 0.0899913 0.137245 0.19658 0.379178
@@ -25,7 +25,7 @@
           0 0 0 0 0 0 0 0.042327 0.177972 0.445703 0.643397 0.648713
           0.648832 0.650173 0.8
         </DataArray>
-        <DataArray type="Float32" Name="Sl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           0.95 0.969761 0.95 0.95493 0.97808 0.962585 1 1 1 1 1 1
           1 1 1 1 1 1 1 1 1 1 1 1
           1 1 1 1 0.97415 0.942823 0.938022 0.910003 0.910009 0.862755 0.80342 0.620822
@@ -47,7 +47,7 @@
           1 1 1 1 1 1 1 0.957673 0.822028 0.554297 0.356603 0.351287
           0.351168 0.349827 0.2
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
           1.3e+07 1.29986e+07 1.3e+07 1.29901e+07 1.29823e+07 1.29734e+07 1.29646e+07 1.29539e+07 1.28541e+07 1.2846e+07 1.2633e+07 1.26254e+07
           1.23741e+07 1.23646e+07 1.21013e+07 1.2091e+07 1.18222e+07 1.18115e+07 1.154e+07 1.15291e+07 1.13776e+07 1.13667e+07 1.13358e+07 1.13248e+07
           1.12936e+07 1.12827e+07 1.12509e+07 1.12401e+07 1.12061e+07 1.11965e+07 1.11532e+07 1.11418e+07 1.1092e+07 1.10809e+07 1.10305e+07 1.10188e+07
@@ -69,7 +69,7 @@
           1.23139e+07 1.20168e+07 1.17246e+07 1.14347e+07 1.12691e+07 1.12272e+07 1.11855e+07 1.11438e+07 1.11032e+07 1.10422e+07 1.10218e+07 1.10128e+07
           1.10072e+07 1.10025e+07 1.1e+07
         </DataArray>
-        <DataArray type="Float32" Name="pl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           1.29995e+07 1.29981e+07 1.29995e+07 1.29896e+07 1.29818e+07 1.29729e+07 1.29641e+07 1.29534e+07 1.28536e+07 1.28455e+07 1.26325e+07 1.26249e+07
           1.23736e+07 1.23641e+07 1.21008e+07 1.20905e+07 1.18217e+07 1.1811e+07 1.15395e+07 1.15286e+07 1.13771e+07 1.13662e+07 1.13353e+07 1.13243e+07
           1.12931e+07 1.12822e+07 1.12504e+07 1.12396e+07 1.12056e+07 1.1196e+07 1.11526e+07 1.11413e+07 1.10915e+07 1.10804e+07 1.10299e+07 1.10181e+07
@@ -113,7 +113,7 @@
           499.687 499.687 499.687 499.687 499.687 499.687 499.687 513.455 566.681 750.862 1129.39 1149.06
           1149.51 1154.64 7500
         </DataArray>
-        <DataArray type="Float32" Name="rhoL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           928.515 933.143 928.515 929.964 940.667 935.54 1068.71 1062.76 1125.7 1126 1125.23 1125.54
           1124.66 1124.95 1123.96 1124.22 1123.15 1123.36 1122.22 1122.38 1120.54 1120.62 1119.06 1119.05
           1117.88 1117.79 1117.04 1116.89 1116.62 1116.62 1116.6 1116.6 1116.58 1116.57 1116.55 1116.55
@@ -135,7 +135,7 @@
           1125.98 1125.39 1124.61 1123.63 1121.61 1119.59 1117.79 1116.6 1116.58 1116.56 1116.55 1116.54
           1116.54 1116.54 1218.97
         </DataArray>
-        <DataArray type="Float32" Name="rhoG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
           107.068 107.058 107.068 106.986 106.926 106.849 106.784 106.692 105.867 105.796 104.03 103.963
           101.878 101.796 99.6119 99.5229 97.2945 97.203 94.9533 94.8607 93.6188 93.5266 93.2872 93.1962
           92.9496 92.8598 92.604 92.5152 92.2354 92.1554 91.7944 91.7001 91.2848 91.1926 90.7728 90.675
@@ -157,7 +157,7 @@
           101.362 98.8917 96.4651 94.0596 92.7018 92.3759 92.0493 91.7166 91.378 90.8701 90.7003 90.625
           90.5781 90.5397 90.6021
         </DataArray>
-        <DataArray type="Float32" Name="mobL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
           4371.14 4533.53 4371.14 4373.87 4370.03 4247.87 3020.03 3040.21 2955.05 2955.3 2954.76 2955.01
           2954.42 2954.63 2954 2954.19 2953.53 2953.67 2953.02 2953.12 2952.13 2952.17 2951.47 2951.47
           2951.05 2951.02 2950.8 2950.76 2600.38 2204.34 2147.9 1839.83 1839.89 1396.85 959.881 227.056
@@ -179,7 +179,7 @@
           2955.46 2955.09 2954.58 2953.96 2952.73 2951.72 2951.03 2385.97 1083.87 114.079 4.35451 3.79262
           3.78069 3.64836 0
         </DataArray>
-        <DataArray type="Float32" Name="mobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           19.0304 4.09502 19.0304 13.8887 1.50781 7.87457 0 0 0 0 0 0
           0 0 0 0 0 0 0 0 0 0 0 0
           0 0 0 0 2.52098 28.5449 36.3922 111.078 111.056 386.054 1095.61 6889.56
diff --git a/test/references/saltflushcc2pncmin-reference.vtu b/test/references/saltflushcc2pncmin-reference.vtu
index ebf73404aedb65bcec08c218a5b5d3dd6f717f80..68c61c057b35af984a667f3cfef244811c58af95 100644
--- a/test/references/saltflushcc2pncmin-reference.vtu
+++ b/test/references/saltflushcc2pncmin-reference.vtu
@@ -2,8 +2,8 @@
 <VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">
   <UnstructuredGrid>
     <Piece NumberOfCells="200" NumberOfPoints="231">
-      <CellData Scalars="Sg" Vectors="velocityW">
-        <DataArray type="Float32" Name="Sg" NumberOfComponents="1" format="ascii">
+      <CellData Scalars="Sn" Vectors="velocityW">
+        <DataArray type="Float32" Name="Sn" NumberOfComponents="1" format="ascii">
           0.0367877 0.0287734 0 0 0 0 0 0 0 0 0 0
           0 0.0647832 0.0896264 0.148761 0.486796 0.620356 0.621248 0.626118 0.0491666 0.0401371 0 0
           0 0 0 0 0 0 0 0 0 0.0798296 0.113544 0.215626
@@ -22,7 +22,7 @@
           0.0604441 0.0673876 0 0 0 0 0 0 0 0 0 0
           0.0446795 0.18798 0.499829 0.639674 0.641078 0.641178 0.641345 0.644413
         </DataArray>
-        <DataArray type="Float32" Name="Sl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="Sw" NumberOfComponents="1" format="ascii">
           0.963212 0.971227 1 1 1 1 1 1 1 1 1 1
           1 0.935217 0.910374 0.851239 0.513204 0.379644 0.378752 0.373882 0.950833 0.959863 1 1
           1 1 1 1 1 1 1 1 1 0.92017 0.886456 0.784374
@@ -41,7 +41,7 @@
           0.939556 0.932612 1 1 1 1 1 1 1 1 1 1
           0.95532 0.81202 0.500171 0.360326 0.358922 0.358822 0.358655 0.355587
         </DataArray>
-        <DataArray type="Float32" Name="pg" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pn" NumberOfComponents="1" format="ascii">
           1.30022e+07 1.29913e+07 1.29731e+07 1.28615e+07 1.26502e+07 1.24178e+07 1.21743e+07 1.19248e+07 1.1672e+07 1.14174e+07 1.1271e+07 1.12331e+07
           1.11948e+07 1.11596e+07 1.11099e+07 1.10568e+07 1.10164e+07 1.10116e+07 1.10077e+07 1.1003e+07 1.29968e+07 1.29831e+07 1.29641e+07 1.28549e+07
           1.26444e+07 1.24096e+07 1.21647e+07 1.19145e+07 1.16614e+07 1.14066e+07 1.12601e+07 1.12223e+07 1.11838e+07 1.11494e+07 1.10991e+07 1.10443e+07
@@ -60,7 +60,7 @@
           1.29726e+07 1.29315e+07 1.28945e+07 1.28402e+07 1.26506e+07 1.23723e+07 1.2104e+07 1.18412e+07 1.15812e+07 1.13227e+07 1.1175e+07 1.11379e+07
           1.11018e+07 1.10696e+07 1.10266e+07 1.10152e+07 1.10105e+07 1.10066e+07 1.10032e+07 1.10005e+07
         </DataArray>
-        <DataArray type="Float32" Name="pl" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="pw" NumberOfComponents="1" format="ascii">
           1.30016e+07 1.29908e+07 1.29726e+07 1.2861e+07 1.26497e+07 1.24173e+07 1.21738e+07 1.19243e+07 1.16715e+07 1.14169e+07 1.12705e+07 1.12326e+07
           1.11943e+07 1.11591e+07 1.11094e+07 1.10563e+07 1.10156e+07 1.10105e+07 1.10066e+07 1.10019e+07 1.29963e+07 1.29826e+07 1.29636e+07 1.28544e+07
           1.26439e+07 1.24091e+07 1.21642e+07 1.1914e+07 1.16609e+07 1.14061e+07 1.12596e+07 1.12218e+07 1.11833e+07 1.11489e+07 1.10986e+07 1.10437e+07
@@ -98,7 +98,7 @@
           519.706 522.163 499.687 499.687 499.687 499.687 499.687 499.687 499.687 499.687 499.687 499.687
           514.254 571.295 815.754 1116.2 1121.12 1121.47 1122.06 1133.07
         </DataArray>
-        <DataArray type="Float32" Name="rhoL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoW" NumberOfComponents="1" format="ascii">
           935.631 948.289 1125.66 1125.23 1124.68 1123.99 1123.18 1122.26 1121.25 1120.2 1118.6 1117.44
           1116.66 1116.6 1116.58 1116.56 1116.55 1116.54 1116.54 1116.54 929.842 936.225 1096.56 1125.9
           1125.36 1124.66 1123.8 1122.81 1121.72 1120.58 1118.82 1117.51 1116.62 1116.6 1116.58 1116.56
@@ -117,7 +117,7 @@
           928.504 928.678 951.637 1126.76 1126.49 1126.02 1125.32 1124.37 1123.22 1121.92 1119.69 1117.8
           1116.58 1116.57 1116.55 1116.54 1116.54 1116.54 1116.54 1116.54
         </DataArray>
-        <DataArray type="Float32" Name="rhoG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="rhoN" NumberOfComponents="1" format="ascii">
           107.089 107.005 106.859 105.933 104.179 102.25 100.229 98.1594 96.064 93.9543 92.7525 92.4507
           92.1406 91.8482 91.434 90.9918 90.6547 90.615 90.5827 90.5431 107.042 106.931 106.778 105.871
           104.123 102.173 100.142 98.0679 95.9704 93.86 92.6595 92.3595 92.0493 91.7632 91.344 90.8878
@@ -136,7 +136,7 @@
           106.84 106.497 106.191 105.739 104.163 101.848 99.6195 97.4395 95.2854 93.1452 91.9396 91.6529
           91.3665 91.0981 90.7399 90.6453 90.6063 90.5734 90.545 90.5227
         </DataArray>
-        <DataArray type="Float32" Name="mobL" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobW" NumberOfComponents="1" format="ascii">
           4258.3 3947.06 2954.98 2954.69 2954.34 2953.93 2953.45 2952.95 2952.44 2951.96 2951.31 2950.92
           2950.71 2115.43 1843.68 1302.26 69.6703 7.54034 7.39172 6.6185 4287.8 4156.76 2960.07 2955.21
           2954.86 2954.4 2953.87 2953.29 2952.7 2952.15 2951.39 2950.94 2950.7 1947.5 1607.64 844.308
@@ -155,7 +155,7 @@
           4132.31 3961.93 4428.75 2955.93 2955.77 2955.47 2955 2954.37 2953.64 2952.9 2951.77 2951.04
           2356.47 1015.78 58.7773 4.78358 4.6182 4.6066 4.58728 4.24261
         </DataArray>
-        <DataArray type="Float32" Name="mobG" NumberOfComponents="1" format="ascii">
+        <DataArray type="Float32" Name="mobN" NumberOfComponents="1" format="ascii">
           7.47737 3.51284 0 0 0 0 0 0 0 0 0 0
           0 41.5731 109.722 488.569 13315.4 24280.2 24363 24817.2 18.0855 9.75845 0 0
           0 0 0 0 0 0 0 0 0 77.711 221.121 1428.22