diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
index 848bb6b10..0850da7e2 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
@@ -257,7 +257,7 @@ namespace Opm {
 
         /// Computing the water velocity without shear-thinning for the cell faces.
         /// The water velocity will be used for shear-thinning calculation.
-        void computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
+        void computeWaterShearVelocityFaces(const V& transi, const std::vector<typename Base::ReservoirResidualQuant>& rq,
                                             const std::vector<ADB>& phasePressure, const SolutionState& state,
                                             std::vector<double>& water_vel, std::vector<double>& visc_mult);
 
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
index f202ec05a..159542ebf 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel_impl.hpp
@@ -291,14 +291,19 @@ namespace Opm {
         // Set up the common parts of the mass balance equations
         // for each active phase.
         const V transi = subset(geo_.transmissibility(), ops_.internal_faces);
-        const std::vector<ADB> kr = computeRelPerm(state);
+        {
+            const std::vector<ADB> kr = computeRelPerm(state);
+            for (int phaseIdx=0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
+                sd_.rq[phaseIdx].kr = kr[canph_[phaseIdx]];
+            }
+        }
 
 
         if (has_plyshlog_) {
             std::vector<double> water_vel;
             std::vector<double> visc_mult;
 
-            computeWaterShearVelocityFaces(transi, kr, state.canonical_phase_pressures, state, water_vel, visc_mult);
+            computeWaterShearVelocityFaces(transi, sd_.rq, state.canonical_phase_pressures, state, water_vel, visc_mult);
             if ( !polymer_props_ad_.computeShearMultLog(water_vel, visc_mult, shear_mult_faces_) ) {
                 // std::cerr << " failed in calculating the shear-multiplier " << std::endl;
                 OPM_THROW(std::runtime_error, " failed in calculating the shear-multiplier. ");
@@ -307,9 +312,9 @@ namespace Opm {
 
         for (int phaseIdx = 0; phaseIdx < fluid_.numPhases(); ++phaseIdx) {
             const std::vector<PhasePresence>& cond = phaseCondition();
-            const ADB mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
-            const ADB rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
-            computeMassFlux(phaseIdx, transi, kr[canph_[phaseIdx]], mu, rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
+            sd_.rq[phaseIdx].mu = fluidViscosity(canph_[phaseIdx], state.canonical_phase_pressures[canph_[phaseIdx]], state.temperature, state.rs, state.rv, cond);
+            sd_.rq[phaseIdx].rho = fluidDensity(canph_[phaseIdx], sd_.rq[phaseIdx].b, state.rs, state.rv);
+            computeMassFlux(phaseIdx, transi, sd_.rq[phaseIdx].kr, sd_.rq[phaseIdx].mu, sd_.rq[phaseIdx].rho, state.canonical_phase_pressures[canph_[phaseIdx]], state);
 
             residual_.material_balance_eq[ phaseIdx ] =
                 pvdt_ * (sd_.rq[phaseIdx].accum[1] - sd_.rq[phaseIdx].accum[0])
@@ -585,7 +590,7 @@ namespace Opm {
 
     template<class Grid>
     void
-    BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<ADB>& kr,
+    BlackoilPolymerModel<Grid>::computeWaterShearVelocityFaces(const V& transi, const std::vector<typename Base::ReservoirResidualQuant>& rq,
                                                                const std::vector<ADB>& phasePressure, const SolutionState& state,
                                                                std::vector<double>& water_vel, std::vector<double>& visc_mult)
     {
@@ -598,7 +603,7 @@ namespace Opm {
 
         const ADB tr_mult = transMult(state.pressure);
         const ADB mu    = fluidViscosity(canonicalPhaseIdx, phasePressure[canonicalPhaseIdx], state.temperature, state.rs, state.rv, cond);
-        sd_.rq[phase].mob = tr_mult * kr[canonicalPhaseIdx] / mu;
+        sd_.rq[phase].mob = tr_mult * rq[phase].kr / mu;
 
         // compute gravity potensial using the face average as in eclipse and MRST
         const ADB rho   = fluidDensity(canonicalPhaseIdx, sd_.rq[phase].b, state.rs, state.rv);
@@ -617,7 +622,7 @@ namespace Opm {
         const ADB mc = computeMc(state);
         ADB krw_eff = polymer_props_ad_.effectiveRelPerm(state.concentration,
                                                          cmax,
-                                                         kr[canonicalPhaseIdx]);
+                                                         rq[phase].kr);
         ADB inv_wat_eff_visc = polymer_props_ad_.effectiveInvWaterVisc(state.concentration, mu.value());
         sd_.rq[ phase ].mob = tr_mult * krw_eff * inv_wat_eff_visc;