diff --git a/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh b/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh
index 975d8b3bc1cb738d8292b4fcb0a38d8ddabe4fc6..3434a8f68a4379b7103dee42db2e8797ee4088f0 100644
--- a/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/thermalconductivityjohansen.hh
@@ -128,7 +128,8 @@ public:
 
         const Scalar kappa = 15.6; // fitted to medium quartz sand
         const Scalar rhoBulk = rhoSolid*porosity;
-        const Scalar lSat = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaW, porosity);
+        // lambdaSolid^(1-porosity) * lambdaW^porosity =
+        const Scalar lSat = lambdaSolid * std::pow(lambdaW / lambdaSolid, porosity);
         const Scalar lDry = (0.135*rhoBulk + 64.7)/(rhoSolid - 0.947*rhoBulk);
         const Scalar Ke = (kappa*satW)/(1+(kappa-1)*satW);// Kersten number, equation 13
 
diff --git a/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh b/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh
index 1f27e8e0562415541313e76340cf0ac12b2da440..9bebe3b2596a0455116c86d606370bca01f9732b 100644
--- a/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh
+++ b/dumux/material/fluidmatrixinteractions/2p/thermalconductivitysomerton.hh
@@ -119,9 +119,9 @@ public:
                                                const Scalar rhoSolid = 0.0 /*unused*/)
     {
         const Scalar satW = std::max<Scalar>(0.0, sw);
-        // geometric mean
-        const Scalar lSat = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaW, porosity);
-        const Scalar lDry = std::pow(lambdaSolid, (1.0 - porosity)) * std::pow(lambdaN, porosity);
+        // geometric mean, using ls^(1-p)*l^p = ls*(l/ls)^p
+        const Scalar lSat = lambdaSolid * std::pow(lambdaW / lambdaSolid, porosity);
+        const Scalar lDry = lambdaSolid * std::pow(lambdaN / lambdaSolid, porosity);
 
         return lDry + std::sqrt(satW) * (lSat - lDry);
     }