diff --git a/examples/problems/co2ptflashproblem.hh b/examples/problems/co2ptflashproblem.hh
index a75036acb..e7df55c89 100644
--- a/examples/problems/co2ptflashproblem.hh
+++ b/examples/problems/co2ptflashproblem.hh
@@ -77,6 +77,14 @@ struct NumComp<TypeTag, TTag::CO2PTBaseProblem> {
     static constexpr int value = 3;
 };
 
+template <class TypeTag, class MyTypeTag>
+struct EnableDummyWater { using type = UndefinedProperty; };
+
+template <class TypeTag>
+struct EnableDummyWater<TypeTag, TTag::CO2PTBaseProblem> {
+    static constexpr bool value = true;
+};
+
 // Set the grid type: --->2D
 template <class TypeTag>
 struct Grid<TypeTag, TTag::CO2PTBaseProblem> { using type = Dune::YaspGrid</*dim=*/2>; };
@@ -105,9 +113,10 @@ struct FluidSystem<TypeTag, TTag::CO2PTBaseProblem>
 private:
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     static constexpr int num_comp = getPropValue<TypeTag, Properties::NumComp>();
+    static constexpr bool enable_water = getPropValue<TypeTag, Properties::EnableDummyWater>();
 
 public:
-    using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp>;
+    using type = Opm::GenericOilGasWaterFluidSystem<Scalar, num_comp, enable_water>;
 };
 
 // Set the material Law
diff --git a/opm/models/common/multiphasebaseproperties.hh b/opm/models/common/multiphasebaseproperties.hh
index f067612f9..968a750ba 100644
--- a/opm/models/common/multiphasebaseproperties.hh
+++ b/opm/models/common/multiphasebaseproperties.hh
@@ -83,6 +83,8 @@ struct EnableDispersion { using type = UndefinedProperty; };
 //! Enable convective mixing?
 template<class TypeTag, class MyTypeTag>
 struct EnableConvectiveMixing { using type = UndefinedProperty; };
+template <class TypeTag, class MyTypeTag>
+struct EnableWater { using type = UndefinedProperty; };
 
 } // namespace Opm::Properties
 
diff --git a/opm/models/ptflash/flashindices.hh b/opm/models/ptflash/flashindices.hh
index 9027e4334..ac4c0ebf4 100644
--- a/opm/models/ptflash/flashindices.hh
+++ b/opm/models/ptflash/flashindices.hh
@@ -48,19 +48,20 @@ class FlashIndices
 {
     static constexpr int numComponents = getPropValue<TypeTag, Properties::NumComponents>();
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableWater = getPropValue<TypeTag, Properties::EnableWater>() };
     using EnergyIndices = Opm::EnergyIndices<PVOffset + numComponents, enableEnergy>;
 
 public:
     //! All phases active (note: immiscible/"dummy" water phase)
-    static constexpr bool waterEnabled = true;
+    static constexpr bool waterEnabled = enableWater;
     static constexpr bool gasEnabled = true;
     static constexpr bool oilEnabled = true;
 
     //! number of active phases
-    static constexpr int numPhases = 3;
+    static constexpr int numPhases = enableWater ? 3 : 2;
 
     //! number of equations/primary variables
-    static const int numEq = numComponents + EnergyIndices::numEq_ + 1;
+    static const int numEq = numComponents + EnergyIndices::numEq_ + (enableWater ? 1 : 0);
 
     // Primary variable indices
 
@@ -71,7 +72,7 @@ public:
     static constexpr int z0Idx = pressure0Idx + 1;
 
     //! Index of water saturation
-    static constexpr int water0Idx = z0Idx + numComponents - 1;
+    static constexpr int water0Idx = enableWater ? z0Idx + numComponents - 1 : -1000;
     
     // equation indices
 
diff --git a/opm/models/ptflash/flashintensivequantities.hh b/opm/models/ptflash/flashintensivequantities.hh
index a7038898d..4fad4b706 100644
--- a/opm/models/ptflash/flashintensivequantities.hh
+++ b/opm/models/ptflash/flashintensivequantities.hh
@@ -78,6 +78,8 @@ class FlashIntensiveQuantities
     enum { pressure0Idx = Indices::pressure0Idx };
     enum { water0Idx = Indices::water0Idx};
 
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
@@ -216,8 +218,12 @@ public:
         (R * fluidState_.temperature(FluidSystem::gasPhaseIdx));
 
 
-        // Update oil/gas saturation; water saturation is a primary variable
-        Evaluation Sw = priVars.makeEvaluation(water0Idx, timeIdx);
+        // Update saturation
+        Evaluation Sw = 0.0;
+        if constexpr (waterEnabled) {
+            Sw = priVars.makeEvaluation(water0Idx, timeIdx);
+            fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
+        }
         Evaluation L = fluidState_.L();
         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);
@@ -228,7 +234,6 @@ public:
 
         fluidState_.setSaturation(FluidSystem::oilPhaseIdx, So);
         fluidState_.setSaturation(FluidSystem::gasPhaseIdx, Sg);
-        fluidState_.setSaturation(FluidSystem::waterPhaseIdx, Sw);
 
         fluidState_.setCompressFactor(FluidSystem::oilPhaseIdx, Z_L);
         fluidState_.setCompressFactor(FluidSystem::gasPhaseIdx, Z_V);
@@ -253,7 +258,10 @@ public:
 
         // set the phase viscosity and density
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-            paramCache.updatePhase(fluidState_, phaseIdx);
+            if (phaseIdx == static_cast<unsigned int>(FluidSystem::oilPhaseIdx) 
+             || phaseIdx == static_cast<unsigned int>(FluidSystem::gasPhaseIdx)) {
+                paramCache.updatePhase(fluidState_, phaseIdx);
+            }
 
             const Evaluation& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
 
diff --git a/opm/models/ptflash/flashlocalresidual.hh b/opm/models/ptflash/flashlocalresidual.hh
index 8d6172fc1..29ac021db 100644
--- a/opm/models/ptflash/flashlocalresidual.hh
+++ b/opm/models/ptflash/flashlocalresidual.hh
@@ -65,6 +65,8 @@ class FlashLocalResidual: public GetPropType<TypeTag, Properties::DiscLocalResid
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     using EnergyModule = Opm::EnergyModule<TypeTag, enableEnergy>;
 
+    static const bool waterEnabled = Indices::waterEnabled;
+
     using Toolbox = Opm::MathToolbox<Evaluation>;
 
 public:
@@ -82,7 +84,7 @@ public:
         const auto& fs = intQuants.fluidState();
 
         // compute water storage term
-        if (phaseIdx == waterPhaseIdx) {
+        if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
             unsigned eqIdx = conti0EqIdx + numComponents;
             storage[eqIdx] = 
                 Toolbox::template decay<LhsEval>(fs.density(phaseIdx))
@@ -160,7 +162,7 @@ public:
                     up.fluidState().density(phaseIdx)
                     * extQuants.volumeFlux(phaseIdx);
 
-                if (phaseIdx == waterPhaseIdx) {
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
                     unsigned eqIdx = conti0EqIdx + numComponents;
                     flux[eqIdx] = tmp;
                 }
@@ -176,7 +178,7 @@ public:
                     Toolbox::value(up.fluidState().density(phaseIdx))
                      * extQuants.volumeFlux(phaseIdx);
                 
-                if (phaseIdx == waterPhaseIdx) {
+                if (waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)) {
                     unsigned eqIdx = conti0EqIdx + numComponents;
                     flux[eqIdx] = tmp;
                 }
diff --git a/opm/models/ptflash/flashmodel.hh b/opm/models/ptflash/flashmodel.hh
index f348523af..8a192f50e 100644
--- a/opm/models/ptflash/flashmodel.hh
+++ b/opm/models/ptflash/flashmodel.hh
@@ -135,6 +135,10 @@ template<class TypeTag>
 struct EnableEnergy<TypeTag, TTag::FlashModel>
 { static constexpr bool value = false; };
 
+template<class TypeTag>
+struct EnableWater<TypeTag, TTag::MultiPhaseBaseModel> 
+{ static constexpr int value = GetPropType<TypeTag, Properties::FluidSystem>::waterEnabled; };
+
 } // namespace Opm::Properties
 
 namespace Opm {
diff --git a/opm/models/ptflash/flashnewtonmethod.hh b/opm/models/ptflash/flashnewtonmethod.hh
index 093f4a6f7..5b6fa1ddb 100644
--- a/opm/models/ptflash/flashnewtonmethod.hh
+++ b/opm/models/ptflash/flashnewtonmethod.hh
@@ -63,6 +63,8 @@ class FlashNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonMeth
     enum { z0Idx = Indices::z0Idx };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
 
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+
 public:
     /*!
      * \copydoc FvBaseNewtonMethod::FvBaseNewtonMethod(Problem& )
@@ -126,10 +128,12 @@ protected:
             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;
+        if constexpr (waterEnabled) {
+            // 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:
diff --git a/opm/models/ptflash/flashprimaryvariables.hh b/opm/models/ptflash/flashprimaryvariables.hh
index 7f9914327..897c322ff 100644
--- a/opm/models/ptflash/flashprimaryvariables.hh
+++ b/opm/models/ptflash/flashprimaryvariables.hh
@@ -63,6 +63,8 @@ class FlashPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { pressure0Idx = Indices::pressure0Idx };
     enum { water0Idx = Indices::water0Idx };
 
+    static constexpr bool waterEnabled = Indices::waterEnabled;
+
     // phase indices
     enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
     enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
@@ -122,7 +124,9 @@ public:
         (*this)[pressure0Idx] = getValue(fluidState.pressure(oilPhaseIdx));
 
         // assign water saturation
-        (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
+        if constexpr (waterEnabled) {
+            (*this)[water0Idx] = getValue(fluidState.saturation(waterPhaseIdx));
+        }
     }
 
     /*!
@@ -138,7 +142,9 @@ public:
             os << ", z_" << FluidSystem::componentName(compIdx) << " = "
                << this->operator[](z0Idx + compIdx);
         }
-        os << ", S_w = " << this->operator[](water0Idx);
+        if constexpr (waterEnabled) {
+            os << ", S_w = " << this->operator[](water0Idx);
+        }
         os << ")" << std::flush;
     }
 };
diff --git a/opm/simulators/flow/FlowProblemComp.hpp b/opm/simulators/flow/FlowProblemComp.hpp
index 575423d32..d6e6c817f 100644
--- a/opm/simulators/flow/FlowProblemComp.hpp
+++ b/opm/simulators/flow/FlowProblemComp.hpp
@@ -362,7 +362,7 @@ public:
             Dune::FieldVector<Scalar, numComponents> z(0.0);
             Scalar sumMoles = 0.0;
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
-                if (phaseIdx == waterPhaseIdx){ 
+                if (Indices::waterEnabled && phaseIdx == static_cast<unsigned int>(waterPhaseIdx)){ 
                     continue;
                 }
                 const auto saturation = fs.saturation(phaseIdx);
@@ -494,7 +494,7 @@ protected:
         const bool gas_active = FluidSystem::phaseIsActive(gasPhaseIdx);
         const bool oil_active = FluidSystem::phaseIsActive(oilPhaseIdx);
 
-        if (water_active && Indices::numPhases > 1)
+        if (water_active && Indices::numPhases > 2)
             waterSaturationData = fp.get_double("SWAT");
         else
             waterSaturationData.resize(numDof);
diff --git a/opm/simulators/flow/GenericOutputBlackoilModule.cpp b/opm/simulators/flow/GenericOutputBlackoilModule.cpp
index 2d6847fad..19e11d4cd 100644
--- a/opm/simulators/flow/GenericOutputBlackoilModule.cpp
+++ b/opm/simulators/flow/GenericOutputBlackoilModule.cpp
@@ -1676,17 +1676,27 @@ INSTANTIATE_TYPE(double)
 INSTANTIATE_TYPE(float)
 #endif
 
-#define INSTANTIATE_COMP(NUM) \
-    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM>; \
+#define INSTANTIATE_COMP_THREEPHASE(NUM) \
+    template<class T> using FS##NUM = GenericOilGasWaterFluidSystem<T, NUM, true>; \
     template class GenericOutputBlackoilModule<FS##NUM<double>>;
 
-INSTANTIATE_COMP(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput
+INSTANTIATE_COMP_THREEPHASE(0) // \Note: to register the parameter ForceDisableFluidInPlaceOutput
+INSTANTIATE_COMP_THREEPHASE(2)
+INSTANTIATE_COMP_THREEPHASE(3)
+INSTANTIATE_COMP_THREEPHASE(4)
+INSTANTIATE_COMP_THREEPHASE(5)
+INSTANTIATE_COMP_THREEPHASE(6)
+INSTANTIATE_COMP_THREEPHASE(7)
 
-INSTANTIATE_COMP(2)
-INSTANTIATE_COMP(3)
-INSTANTIATE_COMP(4)
-INSTANTIATE_COMP(5)
-INSTANTIATE_COMP(6)
-INSTANTIATE_COMP(7)
+#define INSTANTIATE_COMP_TWOPHASE(NUM) \
+    template<class T> using GFS##NUM = GenericOilGasWaterFluidSystem<T, NUM, false>; \
+    template class GenericOutputBlackoilModule<GFS##NUM<double>>;
+
+INSTANTIATE_COMP_TWOPHASE(2)
+INSTANTIATE_COMP_TWOPHASE(3)
+INSTANTIATE_COMP_TWOPHASE(4)
+INSTANTIATE_COMP_TWOPHASE(5)
+INSTANTIATE_COMP_TWOPHASE(6)
+INSTANTIATE_COMP_TWOPHASE(7)
 
 } // namespace Opm