diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.cpp b/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
index db7ab622..48b1feb3 100644
--- a/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
+++ b/opm/core/props/satfunc/SaturationPropsFromDeck.cpp
@@ -136,9 +136,11 @@ namespace Opm
                                            double* pc,
                                            double* dpcds) const
     {
-        assert(cells != 0);
+        assert(cells != 0);        
+        assert(phaseUsage_.phase_used[BlackoilPhases::Liquid]);
 
         const int np = numPhases();
+
         if (dpcds) {
             ExplicitArraysSatDerivativesFluidState fluidState(phaseUsage_);
             typedef ExplicitArraysSatDerivativesFluidState::Evaluation Evaluation;
@@ -151,12 +153,26 @@ namespace Opm
                 MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
 
                 // copy the values calculated using opm-material to the target arrays
-                for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
-                    double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
-                    pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].value;
+                for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
+                    // skip unused phases
+                    if ( ! phaseUsage_.phase_used[canonicalPhaseIdx]) {
+                        continue;
+                    }
+                    const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
 
-                    for (int satPhaseIdx = 0; satPhaseIdx < np; ++satPhaseIdx)
-                        dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx].derivatives[satPhaseIdx];
+                    const double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
+                    // in opm-material the wetting phase is the reference phase
+                    // for two-phase problems i.e water for oil-water system,
+                    // but for flow it is always oil. Add oil (liquid) capillary pressure value
+                    // to shift the reference phase to oil
+                    pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].value + sign * capillaryPressures[canonicalPhaseIdx].value;
+                    for (int canonicalSatPhaseIdx = 0; canonicalSatPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalSatPhaseIdx) {
+                        if ( ! phaseUsage_.phase_used[canonicalSatPhaseIdx])
+                            continue;
+
+                        const int satPhaseIdx = phaseUsage_.phase_pos[canonicalSatPhaseIdx];
+                        dpcds[np*np*i + satPhaseIdx*np + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid].derivatives[canonicalSatPhaseIdx] + sign * capillaryPressures[canonicalPhaseIdx].derivatives[canonicalSatPhaseIdx];
+                    }
                 }
             }
         } else {
@@ -170,9 +186,18 @@ namespace Opm
                 MaterialLaw::capillaryPressures(capillaryPressures, params, fluidState);
 
                 // copy the values calculated using opm-material to the target arrays
-                for (int pcPhaseIdx = 0; pcPhaseIdx < np; ++pcPhaseIdx) {
-                    double sign = (pcPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
-                    pc[np*i + pcPhaseIdx] = sign*capillaryPressures[pcPhaseIdx];
+                for (int canonicalPhaseIdx = 0; canonicalPhaseIdx < BlackoilPhases::MaxNumPhases; ++canonicalPhaseIdx) {
+                    // skip unused phases
+                    if ( ! phaseUsage_.phase_used[canonicalPhaseIdx])
+                        continue;
+
+                    const int pcPhaseIdx = phaseUsage_.phase_pos[canonicalPhaseIdx];
+                    double sign = (canonicalPhaseIdx == BlackoilPhases::Aqua)? -1.0 : 1.0;
+                    // in opm-material the wetting phase is the reference phase
+                    // for two-phase problems i.e water for oil-water system,
+                    // but for flow it is always oil. Add oil (liquid) capillary pressure value
+                    // to shift the reference phase to oil
+                    pc[np*i + pcPhaseIdx] = capillaryPressures[BlackoilPhases::Liquid] + sign * capillaryPressures[canonicalPhaseIdx];
                 }
             }
         }
diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp
index b3d750c3..20f56c5a 100644
--- a/opm/core/simulator/initStateEquil_impl.hpp
+++ b/opm/core/simulator/initStateEquil_impl.hpp
@@ -767,28 +767,34 @@ namespace Opm
                 double sat[BlackoilPhases::MaxNumPhases];
                 double threshold_sat = 1.0e-6;
 
-                sat[waterpos] = smax[waterpos];
-                sat[gaspos] = smax[gaspos];
-                sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
-                if (sw > smax[waterpos]-threshold_sat ) {
-                    sat[waterpos] = smax[waterpos];
-                    props.capPress(1, sat, &cell, pc, 0);                   
-                    phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
-                } else if (sg > smax[gaspos]-threshold_sat) {
-                    sat[gaspos] = smax[gaspos];
-                    props.capPress(1, sat, &cell, pc, 0);                   
-                    phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
-                }
-                if (sg < smin[gaspos]+threshold_sat) {
-                    sat[gaspos] = smin[gaspos];
-                    props.capPress(1, sat, &cell, pc, 0);
-                    phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
-                }
-                if (sw < smin[waterpos]+threshold_sat) {
-                    sat[waterpos] = smin[waterpos];
-                    props.capPress(1, sat, &cell, pc, 0);
-                    phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
-                }
+                sat[oilpos] = 1.0;
+                  if (water) {
+                      sat[waterpos] = smax[waterpos];
+                      sat[oilpos] -= sat[waterpos];
+                  }
+                  if (gas) {
+                      sat[gaspos] = smax[gaspos];
+                      sat[oilpos] -= sat[gaspos];
+                  }
+                  if (water && sw > smax[waterpos]-threshold_sat ) {
+                      sat[waterpos] = smax[waterpos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
+                  } else if (gas && sg > smax[gaspos]-threshold_sat) {
+                      sat[gaspos] = smax[gaspos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
+                  }
+                  if (gas && sg < smin[gaspos]+threshold_sat) {
+                      sat[gaspos] = smin[gaspos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
+                  }
+                  if (water && sw < smin[waterpos]+threshold_sat) {
+                      sat[waterpos] = smin[waterpos];
+                      props.capPress(1, sat, &cell, pc, 0);
+                      phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
+                  }
             }
             return phase_saturations;
         }