diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
index ad5d475ab..021297ebe 100644
--- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
+++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
@@ -181,10 +181,7 @@ public:
         typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
 
         const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(waterPhaseIdx));
-        Valgrind::CheckDefined(Sw);
-        const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
-        Valgrind::CheckDefined(result);
-        return result;
+        return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
     }
 
     /*!
@@ -327,18 +324,35 @@ public:
     {
         typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
 
-        Scalar Swco = params.Swl();
-
         Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+        Scalar So = FsToolbox::value(fluidState.saturation(oilPhaseIdx));
         Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
-        Sw = std::min(1.0, std::max(Swco, Sw));
-        Sg = std::min(1.0, std::max(0.0, Sg));
 
-        Scalar Sw_ow = Sg + Sw;
-        Scalar So_go = 1 - Sw_ow + Swco;
+        if (params.inconsistentHysteresisUpdate()) {
+            Sg = std::min(1.0, std::max(0.0, Sg));
+            // NOTE: the saturations which are passed to update the hysteresis curves are
+            // inconsistent with the ones used to calculate the relative permabilities. We do
+            // it like this anyway because (a) the saturation functions of opm-core do it
+            // this way (b) the simulations seem to converge better (which is not too much
+            // surprising actually, because the time step does not start on a kink in the
+            // solution) and (c) the Eclipse 100 simulator may do the same.
+            //
+            // Though be aware that from a physical perspective this is definitively
+            // incorrect!
+            params.oilWaterParams().update(/*pcSw=*/1 - So, /*krwSw=*/1 - So, /*krn_Sw=*/1 - So);
+            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
+        }
+        else {
+            Scalar Swco = params.Swl();
+            Sw = std::min(1.0, std::max(0.0, Sw));
+            Sg = std::min(1.0, std::max(0.0, Sg));
 
-        params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
-        params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
+            Scalar Sw_ow = Sg + std::max(Swco, Sw);
+            Scalar So_go = 1 + Swco - Sw_ow;
+
+            params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
+            params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
+        }
     }
 };
 } // namespace Opm
diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp
index 8e331ce7f..dc43f9e71 100644
--- a/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp
+++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp
@@ -124,6 +124,18 @@ public:
     Scalar Swl() const
     { assertFinalized_(); return Swl_; }
 
+    /*!
+     * \brief Specify whether inconsistent saturations should be used to update the
+     *        hysteresis parameters.
+     *
+     * Returning 'true' is wrong from a physical point of view because the saturations
+     * which are used to update the hysteresis parameters are calculated differently than
+     * the ones used to calculate the relperms and capillary pressures. Since Eclipse
+     * E100 probably uses inconsistent saturations, we return true here anyway.
+     */
+    bool inconsistentHysteresisUpdate() const
+    { return true; }
+
 private:
 #ifndef NDEBUG
     void assertFinalized_() const