diff --git a/flowexperimental/comp/flowexp_comp.hpp b/flowexperimental/comp/flowexp_comp.hpp
index 213d63c07..44591c560 100644
--- a/flowexperimental/comp/flowexp_comp.hpp
+++ b/flowexperimental/comp/flowexp_comp.hpp
@@ -20,7 +20,7 @@
 #define FLOWEXP_COMP_HPP
 
 #include <opm/material/constraintsolvers/PTFlash.hpp>
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/models/discretization/common/baseauxiliarymodule.hh>
 #include <opm/models/ptflash/flashmodel.hh>
@@ -169,7 +169,7 @@ private:
     static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
 
 public:
-    using type = Opm::GenericOilGasFluidSystem<Scalar, num_comp>;
+    using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp>;
 };
 template<class TypeTag, int NumComp>
 struct EnableMech<TypeTag, TTag::FlowExpCompProblem<NumComp>> {
diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh
index 578b62d76..9027e4334 100644
--- a/opm/models/ptflash/flashindices.hh
+++ b/opm/models/ptflash/flashindices.hh
@@ -51,13 +51,16 @@ class FlashIndices
     using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
 
 public:
-    static constexpr bool waterEnabled = false;
+    //! All phases active (note: immiscible/"dummy" water phase)
+    static constexpr bool waterEnabled = true;
     static constexpr bool gasEnabled = true;
     static constexpr bool oilEnabled = true;
-    static constexpr int waterPhaseIdx = -1;
-    static constexpr int numPhases = 2;
+
+    //! number of active phases
+    static constexpr int numPhases = 3;
+
     //! number of equations/primary variables
-    static const int numEq = numComponents + EnergyIndices::numEq_;
+    static const int numEq = numComponents + EnergyIndices::numEq_ + 1;
 
     // Primary variable indices
 
@@ -66,6 +69,9 @@ public:
 
     //! Index of the molefraction of the first component
     static constexpr int z0Idx = pressure0Idx + 1;
+
+    //! Index of water saturation
+    static constexpr int water0Idx = z0Idx + numComponents - 1;
     
     // equation indices
 
diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh
index 7e71ef5a3..a7038898d 100644
--- a/opm/models/ptflash/flashintensivequantities.hh
+++ b/opm/models/ptflash/flashintensivequantities.hh
@@ -76,6 +76,7 @@ class FlashIntensiveQuantities
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { dimWorld = GridView::dimensionworld };
     enum { pressure0Idx = Indices::pressure0Idx };
+    enum { water0Idx = Indices::water0Idx};
 
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
@@ -215,20 +216,22 @@ public:
         (R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
 
 
-        // Update saturation
-        // \Note: the current implementation assume oil-gas system.
+        // Update oil/gas saturation; water saturation is a primary variable
+        Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx);
         Evaluation L = fluidState_.L();
-        Evaluation So = Opm::max((L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
-        Evaluation Sg = Opm::max(1 - So, 0.0);
-        Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg);
+        Evaluation So = Opm::max((1 - Sw) * (L * Z_L / ( L * Z_L + (1 - L) * Z_V)), 0.0);
+        Evaluation Sg = Opm::max(1 - So - Sw, 0.0);
+        Scalar sumS = Opm::getValue(So) + Opm::getValue(Sg) + Opm::getValue(Sw);
         So /= sumS;
         Sg /= sumS;
+        Sw /= sumS;
 
-        fluidState_.setSaturation(0, So);
-        fluidState_.setSaturation(1, Sg);
+        fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So);
+        fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg);
+        fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
 
-        fluidState_.setCompressFactor(0, Z_L);
-        fluidState_.setCompressFactor(1, Z_V);
+        fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L);
+        fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);
 
         // Print saturation
         if (flashVerbosity >= 5) {
diff --git a/opm/models/ptflash/flashlocalresidual.hh b/opm/models/ptflash/flashlocalresidual.hh
index c528189e7..8d6172fc1 100644
--- a/opm/models/ptflash/flashlocalresidual.hh
+++ b/opm/models/ptflash/flashlocalresidual.hh
@@ -49,12 +49,16 @@ class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResid
     using Indices = GetPropType<TypeTag, Properties::Indices>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
 
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+    enum { water0Idx = Indices::water0Idx };
     enum { conti0EqIdx = Indices::conti0EqIdx };
 
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     using DiffusionModule = Opm::DiffusionModule<TypeTag, enableDiffusion>;
 
@@ -77,15 +81,25 @@ public:
         const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
         const auto& fs = intQuants.fluidState();
 
-        // compute storage term of all components within all phases
-        for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-            unsigned eqIdx = conti0EqIdx + compIdx;
-            storage[eqIdx] +=
-            Toolbox::template decay<LhsEval>(fs.massFraction(phaseIdx, compIdx))
-                * Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
+        // compute water storage term
+        if (phaseIdx == waterPhaseIdx) {
+            unsigned eqIdx = conti0EqIdx + numComponents;
+            storage[eqIdx] = 
+                Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
                 * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
                 * Toolbox::template decay<LhsEval>(intQuants.porosity());
         }
+        else {
+            // compute storage term of all components within oil/gas phases
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                unsigned eqIdx = conti0EqIdx + compIdx;
+                storage[eqIdx] +=
+                Toolbox::template decay<LhsEval>(fs.massFraction(phaseIdx, compIdx))
+                    * Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
+                    * Toolbox::template decay<LhsEval>(fs.saturation(phaseIdx))
+                    * Toolbox::template decay<LhsEval>(intQuants.porosity());
+            }
+        }
 
         EnergyModule::addPhaseStorage(storage, elemCtx.intensiveQuantities(dofIdx, timeIdx), phaseIdx);
     }
@@ -146,19 +160,31 @@ public:
                     up.fluidState().density(phaseIdx)
                     * extQuants.volumeFlux(phaseIdx);
 
-                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    flux[conti0EqIdx + compIdx] +=
-                        tmp*up.fluidState().massFraction(phaseIdx, compIdx);
+                if (phaseIdx == waterPhaseIdx) {
+                    unsigned eqIdx = conti0EqIdx + numComponents;
+                    flux[eqIdx] = tmp;
+                }
+                else {
+                    for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                        flux[conti0EqIdx + compIdx] +=
+                            tmp*up.fluidState().massFraction(phaseIdx, compIdx);
+                    }
                 }
             }
             else {
                 Evaluation tmp =
                     Toolbox::value(up.fluidState().density(phaseIdx))
                      * extQuants.volumeFlux(phaseIdx);
-
-                for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                    flux[conti0EqIdx + compIdx] +=
-                        tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
+                
+                if (phaseIdx == waterPhaseIdx) {
+                    unsigned eqIdx = conti0EqIdx + numComponents;
+                    flux[eqIdx] = tmp;
+                }
+                else {
+                    for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                        flux[conti0EqIdx + compIdx] +=
+                            tmp*Toolbox::value(up.fluidState().massFraction(phaseIdx, compIdx));
+                    }
                 }
             }
         }
diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh
index 32500f38b..093f4a6f7 100644
--- a/opm/models/ptflash/flashnewtonmethod.hh
+++ b/opm/models/ptflash/flashnewtonmethod.hh
@@ -125,6 +125,12 @@ protected:
         for (unsigned compIdx = 0; compIdx < numComponents - 1; ++compIdx) {
             clampValue_(nextValue[z0Idx + compIdx], tol, 1-tol);
         }
+
+        // limit change in water saturation to 0.2
+        constexpr Scalar dSwMax = 0.2;
+        if (update[Indices::water0Idx] > dSwMax) {
+            nextValue[Indices::water0Idx] = currentValue[Indices::water0Idx] - dSwMax;
+        }
     }
 private:
     void clampValue_(Scalar& val, Scalar minVal, Scalar maxVal) const
diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh
index 2df1c67c3..7f9914327 100644
--- a/opm/models/ptflash/flashprimaryvariables.hh
+++ b/opm/models/ptflash/flashprimaryvariables.hh
@@ -61,6 +61,12 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     // primary variable indices
     enum { z0Idx = Indices::z0Idx };
     enum { pressure0Idx = Indices::pressure0Idx };
+    enum { water0Idx = Indices::water0Idx };
+
+    // phase indices
+    enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
+    enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
+    enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
 
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -108,10 +114,15 @@ public:
         // the energy module
         EnergyModule::setPriVarTemperatures(*this, fluidState);
 
+        // assign components total fraction
         for (int i = 0; i < numComponents - 1; ++i)
             (*this)[z0Idx + i] = getValue(fluidState.moleFraction(i));
 
-        (*this)[pressure0Idx] = getValue(fluidState.pressure(0));
+        // assign pressure
+        (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx));
+
+        // assign water saturation
+        (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
     }
 
     /*!
@@ -121,12 +132,13 @@ public:
      */
     void print(std::ostream& os = std::cout) const
     {
-        os << "(p_" << FluidSystem::phaseName(0) << " = "
+        os << "(p_" << FluidSystem::phaseName(FluidSystem::oilPhaseIdx) << " = "
            << this->operator[](pressure0Idx);
         for (unsigned compIdx = 0; compIdx < numComponents - 2; ++compIdx) {
             os << ", z_" << FluidSystem::componentName(compIdx) << " = "
                << this->operator[](z0Idx + compIdx);
         }
+        os << ", S_w = " << this->operator[](water0Idx);
         os << ")" << std::flush;
     }
 };
diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp
index be6a2616b..575423d32 100644
--- a/opm/simulators/flow/FlowProblemComp.hpp
+++ b/opm/simulators/flow/FlowProblemComp.hpp
@@ -362,6 +362,9 @@ public:
             Dune::FieldVector<Scalar, numComponents> z(0.0);
             Scalar sumMoles = 0.0;
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
+                if (phaseIdx == waterPhaseIdx){ 
+                    continue;
+                }
                 const auto saturation = fs.saturation(phaseIdx);
                 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
                     Scalar tmp = fs.molarity(phaseIdx, compIdx) * saturation;
@@ -527,6 +530,10 @@ protected:
                                             - waterSaturationData[dofIdx]
                                             - gasSaturationData[dofIdx]);
             }
+            if (water_active) {
+                dofFluidState.setSaturation(FluidSystem::waterPhaseIdx,
+                                            waterSaturationData[dofIdx]);
+            }
 
             //////
             // set phase pressures
diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp
index abd8f3cc3..2d6847fad 100644
--- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp
+++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp
@@ -30,7 +30,7 @@
 #include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
 #include <opm/material/fluidsystems/BlackOilDefaultIndexTraits.hpp>
 
-#include <opm/material/fluidsystems/GenericOilGasFluidSystem.hpp>
+#include <opm/material/fluidsystems/GenericOilGasWaterFluidSystem.hpp>
 
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
 #include <opm/input/eclipse/EclipseState/Runspec.hpp>
@@ -1677,7 +1677,7 @@ INSTANTIATE_TYPE(float)
 #endif
 
 #define INSTANTIATE_COMP(NUM) \
-    template<class T> using FS##NUM = GenericOilGasFluidSystem<T, NUM>; \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM>; \
     template class GenericOutputBlackoilModule<FS##NUM<double>>;
 
 INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput