diff --git a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
index 24899e394..825978d13 100644
--- a/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
+++ b/opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp
@@ -220,7 +220,7 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -237,7 +237,7 @@ public:
     static Evaluation pcnw(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
+        const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
         return OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
     }
 
@@ -316,7 +316,7 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        const Evaluation Sw = 1.0 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -327,7 +327,7 @@ public:
     static Evaluation krw(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
+        const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
         return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
     }
 
@@ -338,27 +338,28 @@ public:
     static Evaluation krn(const Params& params,
                           const FluidState& fluidState)
     {
-        Scalar Swco = params.Swl();
+        const Scalar Swco = params.Swl();
 
-        Evaluation Sw =
+        const Evaluation Sw =
             Opm::max(Evaluation(Swco),
-                         Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
-        Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+                     Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx)));
 
-        Evaluation Sw_ow = Sg + Sw;
-        Evaluation So_go = 1.0 - Sw_ow;
-        const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
-        const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
+        const Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+
+        const Evaluation Sw_ow = Sg + Sw;
+        const Evaluation So_go = 1.0 - Sw_ow;
+        const Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
+        const Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
 
         // avoid the division by zero: chose a regularized kro which is used if Sw - Swco
         // < epsilon/2 and interpolate between the oridinary and the regularized kro between
         // epsilon and epsilon/2
         const Scalar epsilon = 1e-5;
         if (Opm::scalarValue(Sw_ow) - Swco < epsilon) {
-            Evaluation kro2 = (kro_ow + kro_go)/2;;
+            const Evaluation kro2 = (kro_ow + kro_go)/2;;
             if (Opm::scalarValue(Sw_ow) - Swco > epsilon/2) {
-                Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
-                Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
+                const Evaluation kro1 = (Sg*kro_go + (Sw - Swco)*kro_ow)/(Sw_ow - Swco);
+                const Evaluation alpha = (epsilon - (Sw_ow - Swco))/(epsilon/2);
                 return kro2*alpha + kro1*(1 - alpha);
             }
 
@@ -397,12 +398,12 @@ public:
             params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krn_Sw=*/1 - Sg);
         }
         else {
-            Scalar Swco = params.Swl();
+            const Scalar Swco = params.Swl();
             Sw = std::min(Scalar(1.0), std::max(Scalar(0.0), Sw));
             Sg = std::min(Scalar(1.0), std::max(Scalar(0.0), Sg));
 
-            Scalar Sw_ow = Sg + std::max(Swco, Sw);
-            Scalar So_go = 1 + Sw_ow;
+            const Scalar Sw_ow = Sg + std::max(Swco, Sw);
+            const Scalar So_go = 1 + 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);
diff --git a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
index bc0a4e9e1..7d88a791d 100644
--- a/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
+++ b/opm/material/fluidmatrixinteractions/EclEpsTwoPhaseLaw.hpp
@@ -151,16 +151,16 @@ public:
     template <class Evaluation>
     static Evaluation twoPhaseSatPcnw(const Params& params, const Evaluation& SwScaled)
     {
-        const Evaluation& SwUnscaled = scaledToUnscaledSatPc(params, SwScaled);
-        const Evaluation& pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwUnscaled);
+        const Evaluation SwUnscaled = scaledToUnscaledSatPc(params, SwScaled);
+        const Evaluation pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwUnscaled);
         return unscaledToScaledPcnw_(params, pcUnscaled);
     }
 
     template <class Evaluation>
     static Evaluation twoPhaseSatPcnwInv(const Params& params, const Evaluation& pcnwScaled)
     {
-        Evaluation pcnwUnscaled = scaledToUnscaledPcnw_(params, pcnwScaled);
-        Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled);
+        const Evaluation pcnwUnscaled = scaledToUnscaledPcnw_(params, pcnwScaled);
+        const Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled);
         return unscaledToScaledSatPc(params, SwUnscaled);
     }
 
@@ -223,16 +223,16 @@ public:
     template <class Evaluation>
     static Evaluation twoPhaseSatKrw(const Params& params, const Evaluation& SwScaled)
     {
-        const Evaluation& SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled);
-        const Evaluation& krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled);
+        const Evaluation SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled);
+        const Evaluation krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled);
         return unscaledToScaledKrw_(SwScaled, params, krwUnscaled);
     }
 
     template <class Evaluation>
     static Evaluation twoPhaseSatKrwInv(const Params& params, const Evaluation& krwScaled)
     {
-        Evaluation krwUnscaled = scaledToUnscaledKrw_(params, krwScaled);
-        Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled);
+        const Evaluation krwUnscaled = scaledToUnscaledKrw_(params, krwScaled);
+        const Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled);
         return unscaledToScaledSatKrw(params, SwUnscaled);
     }
 
@@ -248,16 +248,16 @@ public:
     template <class Evaluation>
     static Evaluation twoPhaseSatKrn(const Params& params, const Evaluation& SwScaled)
     {
-        const Evaluation& SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled);
-        const Evaluation& krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled);
+        const Evaluation SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled);
+        const Evaluation krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled);
         return unscaledToScaledKrn_(params, krnUnscaled);
     }
 
     template <class Evaluation>
     static Evaluation twoPhaseSatKrnInv(const Params& params, const Evaluation& krnScaled)
     {
-        Evaluation krnUnscaled = scaledToUnscaledKrn_(params, krnScaled);
-        Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled);
+        const Evaluation krnUnscaled = scaledToUnscaledKrn_(params, krnScaled);
+        const Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled);
         return unscaledToScaledSatKrn(params, SwUnscaled);
     }
 
diff --git a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp
index 0a22653eb..d893b16a3 100644
--- a/opm/material/fluidmatrixinteractions/EclStone1Material.hpp
+++ b/opm/material/fluidmatrixinteractions/EclStone1Material.hpp
@@ -34,6 +34,8 @@
 
 #include <algorithm>
 #include <cmath>
+#include <stdexcept>
+#include <type_traits>
 
 namespace Opm {
 
@@ -218,7 +220,7 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -235,10 +237,12 @@ public:
     static Evaluation pcnw(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
+        const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
         Valgrind::CheckDefined(Sw);
-        const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
+
+        const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
         Valgrind::CheckDefined(result);
+
         return result;
     }
 
@@ -317,7 +321,7 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -328,7 +332,7 @@ public:
     static Evaluation krw(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
+        const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
         return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
     }
 
@@ -342,16 +346,16 @@ public:
         // the Eclipse docu is inconsistent with naming the variable of connate water: In
         // some places the connate water saturation is represented by "Swl", in others
         // "Swco" is used.
-        Scalar Swco = params.Swl();
+        const Scalar Swco = params.Swl();
 
         // oil relperm at connate water saturations (with Sg=0)
-        Scalar krocw = params.krocw();
+        const Scalar krocw = params.krocw();
 
         const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
         const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
 
-        Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
-        Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco);
+        const Evaluation kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
+        const Evaluation kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Swco);
 
         Evaluation beta;
         if (Sw <= Swco)
@@ -360,9 +364,9 @@ public:
             // there seems to be an error in the ECL documentation: using the approach to
             // the scaled saturations as described there leads to significant deviations
             // from the results produced by Eclipse 100.
-            Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
-            Evaluation SSg = Sg/(1.0 - Swco);
-            Evaluation SSo = 1.0 - SSw - SSg;
+            const Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
+            const Evaluation SSg = Sg/(1.0 - Swco);
+            const Evaluation SSo = 1.0 - SSw - SSg;
 
             if (SSw >= 1.0 || SSg >= 1.0)
                 beta = 1.0;
@@ -383,8 +387,8 @@ public:
     template <class FluidState>
     static void updateHysteresis(Params& params, const FluidState& fluidState)
     {
-        Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
-        Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
+        const Scalar Sw = Opm::scalarValue(fluidState.saturation(waterPhaseIdx));
+        const Scalar Sg = Opm::scalarValue(fluidState.saturation(gasPhaseIdx));
 
         params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
         params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/1 - Sg, /*krnSw=*/1 - Sg);
diff --git a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp
index 6ca08853e..4d0e82742 100644
--- a/opm/material/fluidmatrixinteractions/EclStone2Material.hpp
+++ b/opm/material/fluidmatrixinteractions/EclStone2Material.hpp
@@ -33,6 +33,8 @@
 #include <opm/material/common/MathToolbox.hpp>
 
 #include <algorithm>
+#include <stdexcept>
+#include <type_traits>
 
 namespace Opm {
 
@@ -219,7 +221,7 @@ public:
     static Evaluation pcgn(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
+        const auto Sw = 1.0 - Opm::decay<Evaluation>(fs.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
     }
 
@@ -236,10 +238,12 @@ public:
     static Evaluation pcnw(const Params& params,
                            const FluidState& fs)
     {
-        const auto& Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
+        const auto Sw = Opm::decay<Evaluation>(fs.saturation(waterPhaseIdx));
         Valgrind::CheckDefined(Sw);
-        const auto& result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
+
+        const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
         Valgrind::CheckDefined(result);
+
         return result;
     }
 
@@ -318,7 +322,7 @@ public:
     static Evaluation krg(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        const Evaluation Sw = 1 - Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
         return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
     }
 
@@ -329,7 +333,7 @@ public:
     static Evaluation krw(const Params& params,
                           const FluidState& fluidState)
     {
-        const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
+        const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
         return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
     }
 
@@ -340,17 +344,17 @@ public:
     static Evaluation krn(const Params& params,
                           const FluidState& fluidState)
     {
-        Scalar Swco = params.Swl();
-        const Evaluation& Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
-        const Evaluation& Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
+        const Scalar Swco = params.Swl();
+        const Evaluation Sw = Opm::decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
+        const Evaluation Sg = Opm::decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
 
-        Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
-        Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
-        Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
-        Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg);
-        Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg);
+        const Scalar krocw = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Swco);
+        const Evaluation krow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
+        const Evaluation krw = OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
+        const Evaluation krg = GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), 1 - Sg);
+        const Evaluation krog = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg);
 
-        return krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
+        return krocw * ((krow/krocw + krw) * (krog/krocw + krg) - krw - krg);
     }
 
     /*!