diff --git a/dumux/common/dimensionlessnumbers.hh b/dumux/common/dimensionlessnumbers.hh
index b13b71f8147f36357db15a81ff0494d672726790..2e072d674643fee93c7f5d190f387548df2407af 100644
--- a/dumux/common/dimensionlessnumbers.hh
+++ b/dumux/common/dimensionlessnumbers.hh
@@ -31,6 +31,7 @@
 #include <iostream>
 
 #include <dune/common/exceptions.hh>
+#include <dune/common/math.hh>
 
 namespace Dumux {
 
@@ -177,6 +178,7 @@ static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
         */
         using std::sqrt;
         using std::pow;
+        using Dune::power;
         Scalar numerator    = 0.037 * pow(reynoldsNumber,0.8) * prandtlNumber ;
         Scalar reToMin01    = pow(reynoldsNumber,-0.1);
         Scalar prTo23       = pow(prandtlNumber, (2./3. ) ) ; // MIND THE pts! :-( otherwise the integer exponent version is chosen
@@ -184,7 +186,7 @@ static Scalar nusseltNumberForced(const Scalar reynoldsNumber,
 
         Scalar nusseltTurbular       = numerator / denominator;
         Scalar nusseltLaminar        = 0.664 * sqrt(reynoldsNumber) * pow(prandtlNumber, (1./3.) );
-        Scalar nusseltSingleSphere   = 2 + sqrt( pow(nusseltLaminar,2.) + pow(nusseltTurbular,2.));
+        Scalar nusseltSingleSphere   = 2 + sqrt( power(nusseltLaminar,2) + power(nusseltTurbular,2));
 
         Scalar funckyFactor           = 1 + 1.5 * (1.-porosity); // for spheres of same size
         Scalar nusseltNumber          = funckyFactor * nusseltSingleSphere  ;
diff --git a/dumux/freeflow/rans/oneeq/volumevariables.hh b/dumux/freeflow/rans/oneeq/volumevariables.hh
index 1ebe6ad28bd837586509bc52c8d32e19774cff90..c3bcbc16141e14f901816daf58a84096a7e4dd9d 100644
--- a/dumux/freeflow/rans/oneeq/volumevariables.hh
+++ b/dumux/freeflow/rans/oneeq/volumevariables.hh
@@ -25,6 +25,7 @@
 #ifndef DUMUX_ONEEQ_VOLUME_VARIABLES_HH
 #define DUMUX_ONEEQ_VOLUME_VARIABLES_HH
 
+#include <dune/common/math.hh>
 #include <dumux/common/parameters.hh>
 #include <dumux/freeflow/rans/volumevariables.hh>
 
@@ -154,14 +155,15 @@ public:
     Scalar fW() const
     {
         using std::pow;
-        return g() * pow(1.0 + pow(cw3(), 6.0) / (pow(g(), 6.0) + pow(cw3(), 6.0)), 1.0/6.0);
+        using Dune::power;
+        return g() * pow(1.0 + power(cw3(), 6) / (power(g(), 6) + power(cw3(), 6)), 1.0/6.0);
     }
 
     //! \brief Returns a model function
     Scalar g() const
     {
-        using std::pow;
-        return r() + cw2() * (pow(r(), 6.0) - r());
+        using Dune::power;
+        return r() + cw2() * (power(r(), 6) - r());
     }
 
     //! \brief Returns a model function
diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
index 6ed4e0eea9c1759935e6c89dd294c1692b60d4f3..0a8a06ba6c3b93d22adce2b12018ea296048c4ae 100644
--- a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
+++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh
@@ -180,7 +180,6 @@ public:
     const Scalar fMu() const
     {
         using std::exp;
-        using std::pow;
         return 1.0 - exp(-0.0115 * RANSParentType::yPlus());
     }
 
diff --git a/dumux/freeflow/rans/zeroeq/problem.hh b/dumux/freeflow/rans/zeroeq/problem.hh
index 3d47f07c53b58dc2e3e2f4322b47abb5e3d3d1f2..84e8577ef1029547df27d5af386c509c69fa18d0 100644
--- a/dumux/freeflow/rans/zeroeq/problem.hh
+++ b/dumux/freeflow/rans/zeroeq/problem.hh
@@ -26,7 +26,7 @@
 
 #include <string>
 
-#include <dune/common/power.hh>
+#include <dune/common/math.hh>
 #include <dumux/common/properties.hh>
 #include <dumux/common/staggeredfvproblem.hh>
 #include <dumux/discretization/localview.hh>
@@ -131,8 +131,8 @@ public:
         using std::abs;
         using std::exp;
         using std::min;
-        using std::pow;
         using std::sqrt;
+        using Dune::power;
         const Scalar aPlus = 26.0;
         const Scalar k = 0.0168;
         const Scalar cCP = 1.6;
@@ -188,7 +188,7 @@ public:
             Scalar yFMax = storedYFMax[asImp_().wallElementIndex(elementIdx)];
             Scalar fMax = storedFMax[asImp_().wallElementIndex(elementIdx)];
             Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax);
-            Scalar fKleb = 1.0 / (1.0 + 5.5 * Dune::power(cKleb * effectiveWallDistance / yFMax, 6));
+            Scalar fKleb = 1.0 / (1.0 + 5.5 * power(cKleb * effectiveWallDistance / yFMax, 6));
             kinematicEddyViscosityOuter[elementIdx] = k * cCP * fWake * fKleb;
 
             kinematicEddyViscosityDifference[elementIdx]
diff --git a/dumux/io/grid/cakegridmanager.hh b/dumux/io/grid/cakegridmanager.hh
index 7dd1c9dcafc8745b000d04c0daac27c96872ca77..e0184eec39d77ce9b1751416b97d58866d5b097f 100644
--- a/dumux/io/grid/cakegridmanager.hh
+++ b/dumux/io/grid/cakegridmanager.hh
@@ -31,6 +31,7 @@
 
 #include <dune/common/dynvector.hh>
 #include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
 #include <dune/grid/common/gridfactory.hh>
 #include <dumux/common/parameters.hh>
 
@@ -159,6 +160,7 @@ public:
         // make the grid
         std::array<std::vector<Scalar>, dim> globalPositions;
         using std::pow;
+        using Dune::power;
         for (int dimIdx = 0; dimIdx < dim; dimIdx++)
         {
             // Each grid direction is subdivided into (numCells + 1) points
@@ -213,13 +215,13 @@ public:
                 // if grading factor is not 1.0, do power law spacing
                 else
                 {
-                    height = (1.0 - gradingFactor) / (1.0 - pow(gradingFactor, numCells));
+                    height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells));
 
                     if (verbose)
                     {
                         std::cout << " -> grading_eff "  << gradingFactor
-                                  << " h_min "  << height * pow(gradingFactor, 0) * length
-                                  << " h_max "  << height * pow(gradingFactor, numCells-1) * length
+                                  << " h_min "  << height * power(gradingFactor, 0) * length
+                                  << " h_max "  << height * power(gradingFactor, numCells-1) * length
                                   << std::endl;
                     }
                 }
@@ -232,10 +234,10 @@ public:
                     if (useGrading)
                     {
                         if (increasingCellSize)
-                            hI *= pow(gradingFactor, i-1);
+                            hI *= power(gradingFactor, i-1);
 
                         else
-                            hI *= pow(gradingFactor, numCells-i);
+                            hI *= power(gradingFactor, numCells-i);
                     }
                     localPositions[i] = localPositions[i-1] + hI;
                 }
diff --git a/dumux/io/grid/gridmanager_yasp.hh b/dumux/io/grid/gridmanager_yasp.hh
index 5ffbbbc4b3f3cb941c328552b3cbd8ff564ca9b8..6ecb00067e5b62781070e78a5bb0f70abf07374b 100644
--- a/dumux/io/grid/gridmanager_yasp.hh
+++ b/dumux/io/grid/gridmanager_yasp.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_IO_GRID_MANAGER_YASP_HH
 #define DUMUX_IO_GRID_MANAGER_YASP_HH
 
+#include <dune/common/math.hh>
+
 #include <dune/grid/yaspgrid.hh>
 #include <dune/grid/io/file/dgfparser/dgfyasp.hh>
 
@@ -322,6 +324,7 @@ private:
     {
         std::array<std::vector<ctype>, dim> globalPositions;
         using std::pow;
+        using Dune::power;
         for (int dimIdx = 0; dimIdx < dim; dimIdx++)
         {
             for (int zoneIdx = 0; zoneIdx < cells[dimIdx].size(); ++zoneIdx)
@@ -372,13 +375,13 @@ private:
                 // if grading factor is not 1.0, do power law spacing
                 else
                 {
-                    height = (1.0 - gradingFactor) / (1.0 - pow(gradingFactor, numCells));
+                    height = (1.0 - gradingFactor) / (1.0 - power(gradingFactor, numCells));
 
                     if (verbose)
                     {
                         std::cout << " -> grading_eff "  << gradingFactor
-                                  << " h_min "  << height * pow(gradingFactor, 0) * length
-                                  << " h_max "  << height * pow(gradingFactor, numCells-1) * length
+                                  << " h_min "  << height * power(gradingFactor, 0) * length
+                                  << " h_max "  << height * power(gradingFactor, numCells-1) * length
                                   << std::endl;
                     }
                 }
@@ -392,11 +395,11 @@ private:
                     {
                         if (increasingCellSize)
                         {
-                            hI *= pow(gradingFactor, i);
+                            hI *= power(gradingFactor, i);
                         }
                         else
                         {
-                            hI *= pow(gradingFactor, numCells-i-1);
+                            hI *= power(gradingFactor, numCells-i-1);
                         }
                     }
                     localPositions.push_back(localPositions[i] + hI);
diff --git a/dumux/material/binarycoefficients/air_mesitylene.hh b/dumux/material/binarycoefficients/air_mesitylene.hh
index 608a4f2bf0f72706b89899911afb655ede4caf59..14a7254a9d2e4d7bb3b0581e359fa10b5bd86c50 100644
--- a/dumux/material/binarycoefficients/air_mesitylene.hh
+++ b/dumux/material/binarycoefficients/air_mesitylene.hh
@@ -71,6 +71,7 @@ public:
         using std::pow;
         using std::sqrt;
         using std::exp;
+        using Dune::power;
         const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
         const Scalar M_a = 1e3*Air::molarMass(); // [g/mol] molecular weight of air
         const Scalar Tb_m = 437.9;        // [K] boiling temperature of mesitylene
@@ -92,7 +93,7 @@ public:
         const Scalar B_ = 0.00217 - 0.0005*sqrt(1.0/M_a + 1.0/M_m);
         const Scalar Mr = (M_a + M_m)/(M_a*M_m);
         const Scalar D_am = (B_*sqrt(temperature*temperature*temperature*Mr))
-                           /(1e-5*pressure*pow(sigma_am, 2) * Omega); // [cm^2/s]
+                           /(1e-5*pressure*power(sigma_am, 2) * Omega); // [cm^2/s]
 
         return 1e-4*D_am; // [m^2/s]
     }
diff --git a/dumux/material/binarycoefficients/air_xylene.hh b/dumux/material/binarycoefficients/air_xylene.hh
index 1d2e2056fb6232cc80e5bf66c8ee92aaf14543e9..0c48323705072d3b2d87acb85ba94ddcd25f275d 100644
--- a/dumux/material/binarycoefficients/air_xylene.hh
+++ b/dumux/material/binarycoefficients/air_xylene.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_BINARY_COEFF_AIR_XYLENE_HH
 #define DUMUX_BINARY_COEFF_AIR_XYLENE_HH
 
+#include <dune/common/math.hh>
+
 #include <dumux/material/components/air.hh>
 #include <dumux/material/components/xylene.hh>
 
@@ -69,6 +71,7 @@ public:
         pressure = min(pressure, 1e8); // regularization
 
         using std::pow;
+        using Dune::power;
         using std::sqrt;
         using std::exp;
         const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
@@ -90,7 +93,7 @@ public:
         const Scalar B_ = 0.00217 - 0.0005*sqrt(1.0/M_a + 1.0/M_x);
         const Scalar Mr = (M_a + M_x)/(M_a*M_x);
         const Scalar D_ax = (B_*pow(temperature,1.5)*sqrt(Mr))
-                           /(1e-5*pressure*pow(sigma_ax, 2.0)*Omega); // [cm^2/s]
+                           /(1e-5*pressure*power(sigma_ax, 2)*Omega); // [cm^2/s]
 
         return D_ax*1e-4;   //  [m^2/s]
     }
diff --git a/dumux/material/binarycoefficients/brine_co2.hh b/dumux/material/binarycoefficients/brine_co2.hh
index 224e38a0bc18ee0f7eab33bedbd22088a8c0adef..343808b1d58eff2923ac57c1e711aeaad7ab9af5 100644
--- a/dumux/material/binarycoefficients/brine_co2.hh
+++ b/dumux/material/binarycoefficients/brine_co2.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_BINARY_COEFF_BRINE_CO2_HH
 #define DUMUX_BINARY_COEFF_BRINE_CO2_HH
 
+#include <dune/common/math.hh>
+
 #include <dumux/common/parameters.hh>
 #include <dumux/material/components/brine.hh>
 #include <dumux/material/components/h2o.hh>
@@ -433,8 +435,8 @@ public:
         const Scalar phiCO2 = fugacityCoeffCO2_(temperature, pgCO2, rhoCO2);
 
         using std::log;
-        using std::pow;
-        const Scalar exponent = A - log(phiCO2) + 2*B*mol_NaCl + C*pow(mol_NaCl,2);
+        using Dune::power;
+        const Scalar exponent = A - log(phiCO2) + 2*B*mol_NaCl + C*power(mol_NaCl,2);
 
         using std::exp;
         const Scalar mol_CO2w = pgCO2 / (1e5 * exp(exponent)); /* paper: equation (6) */
diff --git a/dumux/material/binarycoefficients/h2o_mesitylene.hh b/dumux/material/binarycoefficients/h2o_mesitylene.hh
index 8bd354f004ff9a746affa9a950289de4663733d8..dc051e01af97d3477a7e099dc3324da5a97906cf 100644
--- a/dumux/material/binarycoefficients/h2o_mesitylene.hh
+++ b/dumux/material/binarycoefficients/h2o_mesitylene.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_BINARY_COEFF_H2O_MESITYLENE_HH
 #define DUMUX_BINARY_COEFF_H2O_MESITYLENE_HH
 
+#include <dune/common/math.hh>
+
 #include <dumux/material/components/h2o.hh>
 #include <dumux/material/components/mesitylene.hh>
 
@@ -74,6 +76,7 @@ public:
 
         using std::sqrt;
         using std::pow;
+        using Dune::power;
         using std::exp;
         const Scalar M_m = 1e3*Mesitylene::molarMass(); // [g/mol] molecular weight of mesitylene
         const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
@@ -98,7 +101,7 @@ public:
         const Scalar B_ = 0.00217 - 0.0005*sqrt(1.0/M_w + 1.0/M_m);
         const Scalar Mr = (M_w + M_m)/(M_w*M_m);
         const Scalar D_wm = (B_*pow(temperature, 1.6)*sqrt(Mr))
-                            /(1e-5*pressure*pow(sigma_wm, 2)*Omega); // [cm^2/s]
+                            /(1e-5*pressure*power(sigma_wm, 2)*Omega); // [cm^2/s]
 
         return D_wm*1e-4;   //  [m^2/s]
     }
diff --git a/dumux/material/binarycoefficients/h2o_xylene.hh b/dumux/material/binarycoefficients/h2o_xylene.hh
index 79c51461a0663169bcfe8c675b3e5afdec03f9f1..801f9626a85e820cc91e294860fb3b7454840e3e 100644
--- a/dumux/material/binarycoefficients/h2o_xylene.hh
+++ b/dumux/material/binarycoefficients/h2o_xylene.hh
@@ -24,6 +24,8 @@
 #ifndef DUMUX_BINARY_COEFF_H2O_XYLENE_HH
 #define DUMUX_BINARY_COEFF_H2O_XYLENE_HH
 
+#include <dune/common/math.hh>
+
 #include <dumux/material/components/h2o.hh>
 #include <dumux/material/components/xylene.hh>
 
@@ -78,6 +80,7 @@ public:
         using std::sqrt;
         using std::exp;
         using std::pow;
+        using Dune::power;
         const Scalar M_x = 1e3*Xylene::molarMass(); // [g/mol] molecular weight of xylene
         const Scalar M_w = 1e3*H2O::molarMass(); // [g/mol] molecular weight of water
         const Scalar Tb_x = 412.9;        // [K] boiling temperature of xylene
@@ -99,7 +102,7 @@ public:
         const Scalar  B_ = 0.00217 - 0.0005*sqrt(1.0/M_w + 1.0/M_x);
         const Scalar Mr = (M_w + M_x)/(M_w*M_x);
         const Scalar D_wx = (B_*pow(temperature,1.6)*sqrt(Mr))
-                           /(1e-5*pressure*pow(sigma_wx, 2.0)*Omega); // [cm^2/s]
+                           /(1e-5*pressure*power(sigma_wx, 2)*Omega); // [cm^2/s]
 
         return D_wx*1e-4; // [m^2/s]
     }
diff --git a/dumux/material/components/air.hh b/dumux/material/components/air.hh
index 4a80c820a3225d50645a26ff44bce0852b68e85e..06a76605b5469b0883fac3d52203883576e2093d 100644
--- a/dumux/material/components/air.hh
+++ b/dumux/material/components/air.hh
@@ -24,11 +24,12 @@
 #ifndef DUMUX_AIR_HH
 #define DUMUX_AIR_HH
 
-#include <dumux/common/exceptions.hh>
-#include <dumux/material/idealgas.hh>
+#include <dune/common/math.hh>
 
+#include <dumux/common/exceptions.hh>
 #include <dumux/material/components/base.hh>
 #include <dumux/material/components/gas.hh>
+#include <dumux/material/idealgas.hh>
 
 namespace Dumux {
 namespace Components {
@@ -244,14 +245,15 @@ public:
         const Scalar eta0 = 0.0266958*sqrt(1000.0*molarMass()*temperature)/(sigma*sigma*Omega);
 
         using std::pow;
+        using Dune::power;
         const Scalar tau = criticalTemperature()/temperature;
         const Scalar rhoc = 10.4477; // [mol/m^3]
         const Scalar delta = 0.001*pressure/(temperature*8.3144598)/rhoc;
         const Scalar etaR = 10.72 * pow(tau, 0.2) * delta
-                            + 1.122 * pow(tau, 0.05) * pow(delta, 4)
-                            + 0.002019 * pow(tau, 2.4) * pow(delta, 9)
+                            + 1.122 * pow(tau, 0.05) * power(delta, 4)
+                            + 0.002019 * pow(tau, 2.4) * power(delta, 9)
                             - 8.876 * pow(tau, 0.6) * delta * exp(-delta)
-                            - 0.02916 * pow(tau, 3.6) * pow(delta, 8) * exp(-delta);
+                            - 0.02916 * pow(tau, 3.6) * power(delta, 8) * exp(-delta);
 
         return (eta0 + etaR)*1e-6;
     }
@@ -309,19 +311,20 @@ public:
         Scalar phi = temperature/100;
 
         using std::pow;
+        using Dune::power;
         Scalar c_p = 0.661738E+01
                 -0.105885E+01 * phi
-                +0.201650E+00 * pow(phi,2)
-                -0.196930E-01 * pow(phi,3)
-                +0.106460E-02 * pow(phi,4)
-                -0.303284E-04 * pow(phi,5)
-                +0.355861E-06 * pow(phi,6);
-        c_p +=  -0.549169E+01 * pow(phi,-1)
-                +0.585171E+01 * pow(phi,-2)
-                -0.372865E+01 * pow(phi,-3)
-                +0.133981E+01 * pow(phi,-4)
-                -0.233758E+00 * pow(phi,-5)
-                +0.125718E-01 * pow(phi,-6);
+                +0.201650E+00 * power(phi,2)
+                -0.196930E-01 * power(phi,3)
+                +0.106460E-02 * power(phi,4)
+                -0.303284E-04 * power(phi,5)
+                +0.355861E-06 * power(phi,6);
+        c_p +=  -0.549169E+01 * power(phi,-1)
+                +0.585171E+01 * power(phi,-2)
+                -0.372865E+01 * power(phi,-3)
+                +0.133981E+01 * power(phi,-4)
+                -0.233758E+00 * power(phi,-5)
+                +0.125718E-01 * power(phi,-6);
         c_p *= IdealGas::R / molarMass(); // in J/(mol*K) / (kg/mol)
 
         return  c_p;
diff --git a/dumux/material/components/brine.hh b/dumux/material/components/brine.hh
index a379ae2630e0aa0f467b68bcc075cc7321cc87e6..17f0ad7d2abf40255e4617fb17427d7c5bcb3ede 100644
--- a/dumux/material/components/brine.hh
+++ b/dumux/material/components/brine.hh
@@ -26,12 +26,12 @@
 
 #include <cmath>
 
-#include <dumux/common/parameters.hh>
+#include <dune/common/math.hh>
 
+#include <dumux/common/parameters.hh>
 #include <dumux/material/components/h2o.hh>
 #include <dumux/material/components/nacl.hh>
 #include <dumux/material/components/tabulatedcomponent.hh>
-
 #include <dumux/material/components/base.hh>
 #include <dumux/material/components/liquid.hh>
 #include <dumux/material/components/gas.hh>
@@ -184,11 +184,11 @@ public:
 
         const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
 
-        using std::pow;
+        using Dune::power;
         Scalar d_h = 0;
         for (int i = 0; i<=3; i++) {
             for (int j=0; j<=2; j++) {
-                d_h = d_h + a[i][j] * pow(theta, i) * pow(m, j);
+                d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
             }
         }
 
@@ -415,9 +415,10 @@ public:
         const Scalar salinity = max(0.0, ThisType::salinity());
 
         using std::pow;
+        using Dune::power;
         using std::exp;
         const Scalar T_C = temperature - 273.15;
-        const Scalar A = (0.42*pow((pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
+        const Scalar A = (0.42*power((pow(salinity, 0.8)-0.17), 2) + 0.045)*pow(T_C, 0.8);
         const Scalar mu_brine = 0.1 + 0.333*salinity + (1.65+91.9*salinity*salinity*salinity)*exp(-A);
         assert(mu_brine > 0.0);
         return mu_brine/1000.0;
diff --git a/dumux/material/components/co2.hh b/dumux/material/components/co2.hh
index c26986a9561490607abc3cd13417448a7c480c25..7be6c1446409d63dd3acd176b39ef86e08f1c51e 100644
--- a/dumux/material/components/co2.hh
+++ b/dumux/material/components/co2.hh
@@ -24,12 +24,13 @@
 #ifndef DUMUX_CO2_HH
 #define DUMUX_CO2_HH
 
-#include <dumux/common/exceptions.hh>
-#include <dumux/material/constants.hh>
-
 #include <cmath>
 #include <iostream>
 
+#include <dune/common/math.hh>
+
+#include <dumux/common/exceptions.hh>
+#include <dumux/material/constants.hh>
 #include <dumux/material/components/base.hh>
 #include <dumux/material/components/liquid.hh>
 #include <dumux/material/components/gas.hh>
@@ -366,9 +367,9 @@ public:
         /* dmu : excess viscosity at elevated density */
         rho = gasDensity(temperature, pressure); /* CO2 mass density [kg/m^3] */
 
-        using std::pow;
-        dmu = d11*rho + d21*rho*rho + d64*pow(rho,6)/(TStar*TStar*TStar)
-            + d81*pow(rho,8) + d82*pow(rho,8)/TStar;
+        using Dune::power;
+        dmu = d11*rho + d21*rho*rho + d64*power(rho,6)/(TStar*TStar*TStar)
+            + d81*power(rho,8) + d82*power(rho,8)/TStar;
 
         visco_CO2 = (mu0 + dmu)/1.0E6;   /* conversion to [Pa s] */
 
diff --git a/dumux/material/components/heavyoil.hh b/dumux/material/components/heavyoil.hh
index da8b028e49a272bb5ab8a97e82965421bc97d096..2a7c331039f6d70f3fb34e59964a75c31a67def2 100644
--- a/dumux/material/components/heavyoil.hh
+++ b/dumux/material/components/heavyoil.hh
@@ -24,12 +24,13 @@
 #ifndef DUMUX_HEAVYOIL_HH
 #define DUMUX_HEAVYOIL_HH
 
-#include <dumux/material/idealgas.hh>
-#include <dumux/material/constants.hh>
+#include <dune/common/math.hh>
 
+#include <dumux/material/constants.hh>
 #include <dumux/material/components/base.hh>
-#include <dumux/material/components/liquid.hh>
 #include <dumux/material/components/gas.hh>
+#include <dumux/material/components/liquid.hh>
+#include <dumux/material/idealgas.hh>
 
 namespace Dumux {
 namespace Components {
@@ -122,8 +123,8 @@ public:
         Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
         Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
 
-        using std::pow;
-        return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+        using Dune::power;
+        return A*power(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*power(deltaMolecularWeight,2) + D*deltaMolecularWeight
                 + E*deltaSpecificGravity*deltaMolecularWeight;
     }
 
@@ -139,8 +140,8 @@ public:
         Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
         Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
 
-        using std::pow;
-        return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+        using Dune::power;
+        return A*power(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*power(deltaMolecularWeight,2) + D*deltaMolecularWeight
                 + E*deltaSpecificGravity*deltaMolecularWeight;
     }
 
@@ -156,8 +157,8 @@ public:
         Scalar deltaSpecificGravity = log(refComponentSpecificGravity()/specificGravity());
         Scalar deltaMolecularWeight = log(refComponentMolecularWeight()/molecularWeight());
 
-        using std::pow;
-        return A*pow(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*pow(deltaMolecularWeight,2) + D*deltaMolecularWeight
+        using Dune::power;
+        return A*power(deltaSpecificGravity,2) + B*deltaSpecificGravity + C*power(deltaMolecularWeight,2) + D*deltaMolecularWeight
                 + E*deltaSpecificGravity*deltaMolecularWeight;
     }
 
@@ -199,8 +200,8 @@ public:
     */
     static Scalar boilingTemperature()
     {
-        using std::pow;
-        return refComponentBoilingTemperature() * pow((1 + 2*perbutationFactorBoilingTemperature())/(1 - 2*perbutationFactorBoilingTemperature()),2);
+        using Dune::power;
+        return refComponentBoilingTemperature() * power((1 + 2*perbutationFactorBoilingTemperature())/(1 - 2*perbutationFactorBoilingTemperature()),2);
     }
 
     /*!
@@ -208,8 +209,8 @@ public:
      */
     static Scalar criticalTemperature()
     {
-        using std::pow;
-        return refComponentCriticalTemperature() * pow((1 + 2*perbutationFactorCriticalTemperature())/(1 - 2*perbutationFactorCriticalTemperature()),2);
+        using Dune::power;
+        return refComponentCriticalTemperature() * power((1 + 2*perbutationFactorCriticalTemperature())/(1 - 2*perbutationFactorCriticalTemperature()),2);
     }
 
     /*!
@@ -217,8 +218,8 @@ public:
      */
     static Scalar criticalPressure()
     {
-        using std::pow;
-        return refComponentCriticalPressure() * pow((1 + 2*perbutationFactorCriticalPressure())/(1 - 2*perbutationFactorCriticalPressure()),2);
+        using Dune::power;
+        return refComponentCriticalPressure() * power((1 + 2*perbutationFactorCriticalPressure())/(1 - 2*perbutationFactorCriticalPressure()),2);
     }
 
     /*!
@@ -450,7 +451,8 @@ public:
         Scalar API = 9;
 
         using std::pow;
-        return ((pow(10,0.10231*pow(API,2)-3.9464*API+46.5037))*(pow(temperatureFahrenheit,-0.04542*pow(API,2)+1.70405*API-19.18)))*0.001;
+        using Dune::power;
+        return ((pow(10,0.10231*power(API,2)-3.9464*API+46.5037))*(pow(temperatureFahrenheit,-0.04542*power(API,2)+1.70405*API-19.18)))*0.001;
 
     }
     /*!
diff --git a/dumux/material/components/iapws/common.hh b/dumux/material/components/iapws/common.hh
index 28105284f6f149c6b11bd7a845705bedfad2a983..8938a36f226c6e682cff7e8336d6660cbba666b7 100644
--- a/dumux/material/components/iapws/common.hh
+++ b/dumux/material/components/iapws/common.hh
@@ -30,13 +30,14 @@
 #ifndef DUMUX_IAPWS_COMMON_HH
 #define DUMUX_IAPWS_COMMON_HH
 
-#include <dumux/material/constants.hh>
-
 #include <cmath>
 #include <iostream>
 
-namespace Dumux {
-namespace IAPWS {
+#include <dune/common/math.hh>
+
+#include <dumux/material/constants.hh>
+
+namespace Dumux::IAPWS {
 
 /*!
  * \ingroup IAPWS
@@ -209,6 +210,7 @@ public:
 
         using std::abs;
         using std::pow;
+        using Dune::power;
         Scalar DTbar = abs(Tbar - 1) + thcond_c4;
         Scalar DTbarpow = pow(DTbar, 3./5);
         Scalar Q = 2. + thcond_c5 / DTbarpow;
@@ -224,17 +226,16 @@ public:
         Scalar rhobarQ = pow(rhobar, Q);
 
         lam +=
-            (thcond_d1 / pow(Tbar,10) + thcond_d2) * rhobar18 *
+            (thcond_d1 / power(Tbar,10) + thcond_d2) * rhobar18 *
                 exp(thcond_c1 * (1 - rhobar * rhobar18))
             + thcond_d3 * S * rhobarQ *
                 exp((Q/(1+Q))*(1 - rhobar*rhobarQ))
             + thcond_d4 *
-                exp(thcond_c2 * pow(Troot,3) + thcond_c3 / pow(rhobar,5));
+                exp(thcond_c2 * power(Troot,3) + thcond_c3 / power(rhobar,5));
         return /*thcond_kstar * */ lam;
     }
 };
 
-} // end namespace IAPWS
-} // end namespace Dumux
+} // end namespace Dumux::IAPWS
 
 #endif
diff --git a/dumux/material/components/iapws/region4.hh b/dumux/material/components/iapws/region4.hh
index cccc8f4014cb8efb43358ea66e2e9ee0a83c7ce9..d4db97aff9d50b5a050b78873771bef83da21118 100644
--- a/dumux/material/components/iapws/region4.hh
+++ b/dumux/material/components/iapws/region4.hh
@@ -33,8 +33,9 @@
 #include <cmath>
 #include <iostream>
 
-namespace Dumux {
-namespace IAPWS {
+#include <dune/common/math.hh>
+
+namespace Dumux::IAPWS {
 
 /*!
  * \ingroup IAPWS
@@ -101,21 +102,21 @@ public:
         };
 
         using std::pow;
+        using Dune::power;
         Scalar beta = pow((pressure/1e6 /*from Pa to MPa*/), (1./4.));
-        Scalar beta2 = pow(beta, 2.);
+        Scalar beta2 = power(beta, 2);
         Scalar E = beta2 + n[2] * beta + n[5];
         Scalar F = n[0]*beta2 + n[3]*beta + n[6];
         Scalar G = n[1]*beta2 + n[4]*beta + n[7];
 
         using std::sqrt;
-        Scalar D = ( 2.*G)/(-F -sqrt(pow(F,2.) - 4.*E*G));
-        Scalar temperature = (n[9] + D - sqrt(pow(n[9]+D , 2.) - 4.* (n[8] + n[9]*D)) ) * 0.5;
+        Scalar D = ( 2.*G)/(-F -sqrt(power(F,2) - 4.*E*G));
+        Scalar temperature = (n[9] + D - sqrt(power(n[9]+D , 2) - 4.* (n[8] + n[9]*D)) ) * 0.5;
 
         return temperature;
     }
 };
 
-} // end namespace IAPWS
-} // end namespace Dumux
+} // end namespace Dumux::IAPWS
 
 #endif
diff --git a/dumux/material/eos/pengrobinson.hh b/dumux/material/eos/pengrobinson.hh
index b8a5164733e02b9d09b571f53a555ddd042e0e32..07561c9fb540125472c217cb4cefbddfbf548a31 100644
--- a/dumux/material/eos/pengrobinson.hh
+++ b/dumux/material/eos/pengrobinson.hh
@@ -31,6 +31,7 @@
 #ifndef DUMUX_PENG_ROBINSON_HH
 #define DUMUX_PENG_ROBINSON_HH
 
+#include <dune/common/math.hh>
 #include <dumux/material/fluidstates/temperatureoverlay.hh>
 #include <dumux/material/idealgas.hh>
 #include <dumux/common/exceptions.hh>
@@ -457,10 +458,10 @@ protected:
         Scalar tau = 1 - Tr;
         Scalar omega = Component::acentricFactor();
         using std::sqrt;
-        using std::pow;
-        Scalar f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*pow(tau, 5))/Tr;
-        Scalar f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*pow(tau, 5))/Tr;
-        Scalar f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*pow(tau, 5))/Tr;
+        using Dune::power;
+        Scalar f0 = (tau*(-5.97616 + sqrt(tau)*(1.29874 - tau*0.60394)) - 1.06841*power(tau, 5))/Tr;
+        Scalar f1 = (tau*(-5.03365 + sqrt(tau)*(1.11505 - tau*5.41217)) - 7.46628*power(tau, 5))/Tr;
+        Scalar f2 = (tau*(-0.64771 + sqrt(tau)*(2.41539 - tau*4.26979)) + 3.25259*power(tau, 5))/Tr;
         using std::exp;
         return Component::criticalPressure()*exp(f0 + omega * (f1 + omega*f2));
     }
diff --git a/dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh b/dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh
index d4a60bb2ba911ceeb556bc39040add0af7575935..2e9cbe25c75ad88b1112fdd07f495e93c399085c 100644
--- a/dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh
+++ b/dumux/material/fluidmatrixinteractions/2pia/awnsurfacepolynomial2ndorder.hh
@@ -30,6 +30,7 @@
 #include <algorithm>
 #include <cmath>
 #include <assert.h>
+#include <dune/common/math.hh>
 
 #warning "This header is deprecated. Removal after 3.3. Use new material laws."
 
@@ -67,8 +68,8 @@ public:
         const Scalar a01 = params.a01();
         const Scalar a02 = params.a02();
 
-        using std::pow;
-        const Scalar aAlphaBeta = a00 + a10 * Sw + a20 * pow(Sw,2) + a11*Sw*pc +  a01*pc + a02*pow(pc,2);
+        using Dune::power;
+        const Scalar aAlphaBeta = a00 + a10 * Sw + a20 * power(Sw,2) + a11*Sw*pc +  a01*pc + a02*power(pc,2);
         return aAlphaBeta;
     }
 
diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
index 538a739e6da0ffe4ef2fa8b595afd62141ed3772..a0b8a369ecbe6d83612103047c1bc0dc19e83e7b 100644
--- a/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
+++ b/dumux/material/fluidmatrixinteractions/frictionlaws/manning.hh
@@ -26,6 +26,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <dune/common/math.hh>
 #include "frictionlaw.hh"
 
 namespace Dumux {
@@ -64,11 +65,12 @@ public:
     Dune::FieldVector<Scalar, 2> shearStress(const VolumeVariables& volVars) const final
     {
         using std::pow;
+        using Dune::power;
         using std::hypot;
 
         Dune::FieldVector<Scalar, 2> shearStress(0.0);
 
-        Scalar roughnessHeight = pow(25.68/(1.0/manningN_),6.0);
+        Scalar roughnessHeight = power(25.68/(1.0/manningN_),6);
         roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth());
         const Scalar c = pow((volVars.waterDepth() + roughnessHeight),1.0/6.0) * 1.0/(manningN_);
         const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1));
diff --git a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
index 7a85827458dc260b9f2b1296f44511927945a975..10a8d15cc46e6f8ab91767ebebad1c37bc4891f4 100644
--- a/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
+++ b/dumux/material/fluidmatrixinteractions/frictionlaws/nikuradse.hh
@@ -26,6 +26,7 @@
 
 #include <algorithm>
 #include <cmath>
+#include <dune/common/math.hh>
 #include "frictionlaw.hh"
 
 namespace Dumux {
@@ -62,7 +63,7 @@ public:
      */
     Dune::FieldVector<Scalar, 2> shearStress(const VolumeVariables& volVars) const final
     {
-        using std::pow;
+        using Dune::power;
         using std::log;
         using std::hypot;
 
@@ -70,7 +71,7 @@ public:
 
         Scalar roughnessHeight = ks_;
         roughnessHeight = this->limitRoughH(roughnessHeight, volVars.waterDepth());
-        const Scalar ustarH = pow(0.41,2.0)/pow(log((12*(volVars.waterDepth() + roughnessHeight))/ks_),2.0);
+        const Scalar ustarH = power(0.41,2)/power(log((12*(volVars.waterDepth() + roughnessHeight))/ks_),2);
         const Scalar uv = hypot(volVars.velocity(0),volVars.velocity(1));
 
         shearStress[0] = -ustarH * volVars.velocity(0) * uv;
diff --git a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
index f519cef75fbb410b56a386a75b361b320eb23980..6c2f73924408b7fe825d3ab004047f33a2eaa6db 100644
--- a/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
+++ b/dumux/material/fluidmatrixinteractions/permeabilitykozenycarman.hh
@@ -26,6 +26,7 @@
 
 #include <cmath>
 #include <dune/common/fmatrix.hh>
+#include <dune/common/math.hh>
 
 namespace Dumux {
 
@@ -50,8 +51,8 @@ public:
     template<class Scalar>
     PermeabilityType evaluatePermeability(PermeabilityType refPerm, Scalar refPoro, Scalar poro) const
     {
-        using std::pow;
-        auto factor = pow((1.0 - refPoro)/(1.0 - poro), 2) * pow(poro/refPoro, 3);
+        using Dune::power;
+        auto factor = power((1.0 - refPoro)/(1.0 - poro), 2) * power(poro/refPoro, 3);
         refPerm *= factor;
         return refPerm;
     }
diff --git a/dumux/material/fluidsystems/brine.hh b/dumux/material/fluidsystems/brine.hh
index b83a4f45b3e5520cec0823a412065685aab8c451..edc282da01db37f4d9cacf559de641585837ec4e 100644
--- a/dumux/material/fluidsystems/brine.hh
+++ b/dumux/material/fluidsystems/brine.hh
@@ -24,17 +24,16 @@
 #ifndef DUMUX_BRINE_FLUID_SYSTEM_HH
 #define DUMUX_BRINE_FLUID_SYSTEM_HH
 
-#include <dumux/material/fluidsystems/base.hh>
+#include <dune/common/math.hh>
 
+#include <dumux/common/exceptions.hh>
+#include <dumux/io/name.hh>
+#include <dumux/material/fluidsystems/base.hh>
 #include <dumux/material/constants.hh>
 #include <dumux/material/components/h2o.hh>
 #include <dumux/material/components/nacl.hh>
 #include <dumux/material/components/tabulatedcomponent.hh>
 
-#include <dumux/common/exceptions.hh>
-
-#include <dumux/io/name.hh>
-
 namespace Dumux {
 namespace FluidSystems {
 
@@ -286,13 +285,14 @@ public:
         const Scalar xNaCl = fluidState.massFraction(phaseIdx, NaClIdx);
 
         using std::pow;
+        using Dune::power;
         using std::exp;
         using std::max;
         const Scalar T = max(temperature, 275.0);
         const Scalar salinity = max(0.0, xNaCl);
 
         const Scalar T_C = T - 273.15;
-        const Scalar A = ((0.42*pow((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
+        const Scalar A = ((0.42*power((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
         const Scalar mu_brine = 0.1 + (0.333*salinity) + (1.65+(91.9*salinity*salinity*salinity))*exp(-A); // [cP]
         assert(mu_brine > 0.0);
         return mu_brine/1000.0; // [Pa·s]
@@ -364,11 +364,11 @@ public:
 
         const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
 
-        using std::pow;
+        using Dune::power;
         Scalar d_h = 0;
         for (int i = 0; i<=3; i++)
             for (int j=0; j<=2; j++)
-                d_h = d_h + a[i][j] * pow(theta, i) * pow(m, j);
+                d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
 
         /* heat of dissolution for halite according to Michaelides 1971 */
         const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
diff --git a/examples/biomineralization/doc/fluidmaterial.md b/examples/biomineralization/doc/fluidmaterial.md
index c4fdb397eb0cb9ba617f80a06e349eaea9462ce8..6e353b999f22897b9ddcd8496b2b4bfb05ec3846 100644
--- a/examples/biomineralization/doc/fluidmaterial.md
+++ b/examples/biomineralization/doc/fluidmaterial.md
@@ -558,6 +558,8 @@ and the "brine" fluidsystem assuming water and a single salt determining the bri
 we include all necessary components
 
 ```cpp
+#include <dune/common/math.hh>
+
 #include <dumux/material/fluidsystems/base.hh>
 #include <dumux/material/constants.hh>
 #include <dumux/material/components/h2o.hh>
@@ -757,13 +759,14 @@ of each phase depending on temperature, pressure, and phase composition
                            + fluidState.massFraction(phaseIdx, CaIdx);
 
         using std::pow;
+        using Dune::power;
         using std::exp;
         using std::max;
         const Scalar T = max(temperature, 275.0);
         const Scalar salinity = max(0.0, xNaCl);
 
         const Scalar T_C = T - 273.15;
-        const Scalar A = ((0.42*pow((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
+        const Scalar A = ((0.42*power((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
         const Scalar mu_brine = 0.1 + (0.333*salinity) + (1.65+(91.9*salinity*salinity*salinity))*exp(-A); // [cP]
         assert(mu_brine > 0.0);
         return mu_brine/1000.0; // [Pa·s]
diff --git a/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh b/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh
index 2e0186627280c721e7961503dd8213b80c010a51..6faa11403db233053c59f7ac913b9bcf2b81f215 100644
--- a/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh
+++ b/examples/biomineralization/material/fluidsystems/icpcomplexsalinitybrine.hh
@@ -32,6 +32,8 @@
 // ### Include files
 // [[details]] includes
 // we include all necessary components
+#include <dune/common/math.hh>
+
 #include <dumux/material/fluidsystems/base.hh>
 #include <dumux/material/constants.hh>
 #include <dumux/material/components/h2o.hh>
@@ -273,13 +275,14 @@ public:
                            + fluidState.massFraction(phaseIdx, CaIdx);
 
         using std::pow;
+        using Dune::power;
         using std::exp;
         using std::max;
         const Scalar T = max(temperature, 275.0);
         const Scalar salinity = max(0.0, xNaCl);
 
         const Scalar T_C = T - 273.15;
-        const Scalar A = ((0.42*pow((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
+        const Scalar A = ((0.42*power((pow(salinity, 0.8)-0.17), 2)) + 0.045)*pow(T_C, 0.8);
         const Scalar mu_brine = 0.1 + (0.333*salinity) + (1.65+(91.9*salinity*salinity*salinity))*exp(-A); // [cP]
         assert(mu_brine > 0.0);
         return mu_brine/1000.0; // [Pa·s]
@@ -340,11 +343,11 @@ public:
 
         const Scalar m = (1E3/58.44)*(salinity/(1-salinity));
 
-        using std::pow;
+        using Dune::power;
         Scalar d_h = 0;
         for (int i = 0; i<=3; i++)
             for (int j=0; j<=2; j++)
-                d_h = d_h + a[i][j] * pow(theta, i) * pow(m, j);
+                d_h = d_h + a[i][j] * power(theta, i) * power(m, j);
 
         /* heat of dissolution for halite according to Michaelides 1971 */
         const Scalar delta_h = (4.184/(1E3 + (58.44 * m)))*d_h;
diff --git a/test/common/integrate/test_integratescalar.cc b/test/common/integrate/test_integratescalar.cc
index acc5407e28bcab696e0bbaf9014ddf383815c4af..21be225b9487d262d4f4a70efe21c729fa73a0cd 100644
--- a/test/common/integrate/test_integratescalar.cc
+++ b/test/common/integrate/test_integratescalar.cc
@@ -6,6 +6,7 @@
 
 #include <dune/common/exceptions.hh>
 #include <dune/common/float_cmp.hh>
+#include <dune/common/math.hh>
 #include <dumux/common/integrate.hh>
 #include <dumux/nonlinear/findscalarroot.hh>
 
@@ -48,7 +49,7 @@ int main(int argc, char* argv[])
 
     {
         // test some integral where we know the analytical expression
-        const auto func = [](const double x){ return pow(1.0 - x, 5.0)*pow(x, -1.0/3.0); };
+        const auto func = [](const double x){ return Dune::power(1.0 - x, 5)*pow(x, -1.0/3.0); };
         const auto exactIntegral = 2187.0/5236.0;
         const auto numIntegral = integrateScalarFunction(func, 0.0, 1.0, 1e-8);
         if (Dune::FloatCmp::ne<double, Dune::FloatCmp::CmpStyle::absolute>(exactIntegral, numIntegral, 1e-8))
diff --git a/test/porousmediumflow/2p/sequential/mcwhorteranalyticsolution.hh b/test/porousmediumflow/2p/sequential/mcwhorteranalyticsolution.hh
index 311dc2479c12e389ec3b3aa2470730a82d8aac37..79bf0974a9eb06235cc6d1e4631bd92f1718845e 100644
--- a/test/porousmediumflow/2p/sequential/mcwhorteranalyticsolution.hh
+++ b/test/porousmediumflow/2p/sequential/mcwhorteranalyticsolution.hh
@@ -19,6 +19,7 @@
 #ifndef DUMUX_MCWHORTER_ANALYTIC_HH
 #define DUMUX_MCWHORTER_ANALYTIC_HH
 
+#include <dune/common/math.hh>
 #include <dumux/porousmediumflow/2p/sequential/properties.hh>
 
 /**
@@ -238,6 +239,7 @@ private:
         int k = 0;
 
         using std::pow;
+        using Dune::power;
         using std::abs;
 
         while (diff> tolAnalytic_)
@@ -254,7 +256,7 @@ private:
             I0=0;
             for (int i=0; i<intervalNum_; i++)
             {
-                a_[i] = 0.5 * h_ * sInit_ *(gk_[i] + gk_[i+1])+ pow(h_, 2) / 6* ((3* i + 1) * gk_[i]
+                a_[i] = 0.5 * h_ * sInit_ *(gk_[i] + gk_[i+1])+ power(h_, 2) / 6* ((3* i + 1) * gk_[i]
                         + (3 * i + 2) * gk_[i+1]);
                 b_[i] = 0.5 * h_ * (gk_[i] + gk_[i+1]);
                 I0 += (a_[i] - sInit_ * b_[i]);
@@ -272,7 +274,7 @@ private:
             }
 
             // with f(sInit) = 0: relationship between A and sInit
-            Ak = pow((0.5*porosity_/pow((1 - fInit_), 2)*I0), 0.5);
+            Ak = pow((0.5*porosity_/power((1 - fInit_), 2)*I0), 0.5);
             diff=abs(Ak - Akm1);
             // std::cout<<"diff = "<<diff<<std::endl;
         }
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
index 0ad703c3037902a7e77592d3ed065cb5f755658b..e1102bea270492c06eacc618bfd1d8121290763c 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/combustionlocalresidual.hh
@@ -25,6 +25,7 @@
 #define DUMUX_ENERGY_COMBUSTION_LOCAL_RESIDUAL_HH
 
 #include <cmath>
+#include <dune/common/math.hh>
 #include <dumux/common/spline.hh>
 #include <dumux/common/exceptions.hh>
 #include <dumux/common/properties.hh>
@@ -186,6 +187,7 @@ public:
         const Scalar gamma(0.0589);
         const Scalar TSolid = volVars.temperatureSolid();
         using std::pow;
+        using Dune::power;
         const Scalar as = volVars.fluidSolidInterfacialArea();
         const Scalar mul = fs.viscosity(0);
         const Scalar deltahv = fs.enthalpy(1) - fs.enthalpy(0);
@@ -196,7 +198,7 @@ public:
         // If a different state is to be simulated, please use the actual fluid temperature instead.
         const Scalar Tsat = FluidSystem::vaporTemperature(fs, 1 ) ;
         const Scalar deltaT = TSolid - Tsat;
-        const Scalar secondBracket = pow( (cp *deltaT / (0.006 * deltahv)  ) , 3.0 );
+        const Scalar secondBracket = power( (cp *deltaT / (0.006 * deltahv)  ) , 3);
         const Scalar Prl = volVars.prandtlNumber(0);
         const Scalar thirdBracket = pow( 1/Prl , (1.7/0.33));
         const Scalar QBoil = satW * as * mul * deltahv * firstBracket * secondBracket * thirdBracket;
diff --git a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
index f273a4e84fa35c445bed63342c03e2b034358822..460d1bbf1afa727ef2343fb7e13889acd4251a1f 100644
--- a/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
+++ b/test/porousmediumflow/mpnc/implicit/thermalnonequilibrium/spatialparams.hh
@@ -28,6 +28,7 @@
 #define DUMUX_COMBUSTION_SPATIALPARAMS_HH
 
 #include <dune/common/parametertreeparser.hh>
+#include <dune/common/math.hh>
 
 #include <dumux/material/spatialparams/fvnonequilibrium.hh>
 #include <dumux/material/fluidmatrixinteractions/2p/heatpipelaw.hh>
@@ -75,8 +76,8 @@ public:
 
         characteristicLength_ =getParam<Scalar>("SpatialParams.PorousMedium.meanPoreSize");
 
-        using std::pow;
-        intrinsicPermeability_  =  (pow(characteristicLength_,2.0)  * pow(porosity_, 3.0)) / (150.0 * pow((1.0-porosity_),2.0)); // 1.69e-10 ; //
+        using Dune::power;
+        intrinsicPermeability_  =  (power(characteristicLength_,2)  * power(porosity_, 3)) / (150.0 * power((1.0-porosity_),2)); // 1.69e-10 ; //
 
         factorEnergyTransfer_ = getParam<Scalar>("SpatialParams.PorousMedium.factorEnergyTransfer");
         lengthPM_ = getParam<Scalar>("Grid.lengthPM");
diff --git a/test/porousmediumflow/richards/implicit/analytical/problem.hh b/test/porousmediumflow/richards/implicit/analytical/problem.hh
index 98d5df3b448ace185543d308b371f9e47b76eb17..346582759e5702d9388032d4766807282c10ec5d 100644
--- a/test/porousmediumflow/richards/implicit/analytical/problem.hh
+++ b/test/porousmediumflow/richards/implicit/analytical/problem.hh
@@ -30,6 +30,7 @@
 #define DUMUX_RICHARDS_ANALYTICALPROBLEM_HH
 
 #include <cmath>
+#include <dune/common/math.hh>
 #include <dune/geometry/quadraturerules.hh>
 #include <dune/grid/yaspgrid.hh>
 
@@ -188,15 +189,15 @@ public:
 
         // linear model with complex solution
         // calculated with Matlab script "Richards.m"
-        using std::pow;
+        using Dune::power;
         using std::tanh;
 
-        values = (pow(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2.0)*(1.0/1.0E1)
-            -1.0/1.0E1)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))*4.0E-8-((pow(tanh(globalPos[1]
-            *5.0+time*(1.0/1.0E1)-1.5E1),2.0)*5.0-5.0)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))-1.0E3)
-            *(pow(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2.0)*5.0-5.0)*(pwBottom
+        values = (power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*(1.0/1.0E1)
+            -1.0/1.0E1)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))*4.0E-8-((power(tanh(globalPos[1]
+            *5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))-1.0E3)
+            *(power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom
             *(1.0/2.0)-pwTop*(1.0/2.0))*5.0E-16+tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1)
-            *(pow(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2.0)*5.0-5.0)*(pwBottom
+            *(power(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)-1.5E1),2)*5.0-5.0)*(pwBottom
             *(1.0/2.0)-pwTop*(1.0/2.0))*(pwBottom*5.0E-16-(tanh(globalPos[1]*5.0+time*(1.0/1.0E1)
             -1.5E1)+1.0)*(pwBottom*(1.0/2.0)-pwTop*(1.0/2.0))*5.0E-16+4.99995E-6)*1.0E1;
         return values;