diff --git a/dune/phasefield/porenetwork/pn_1pflow.hh b/dune/phasefield/porenetwork/pn_1pflow.hh
index cc60fce716c9776e5afe4bb4c4928f0f84a13859..6c2e8449a5711e73eea61aff9c1a801e18e26a04 100644
--- a/dune/phasefield/porenetwork/pn_1pflow.hh
+++ b/dune/phasefield/porenetwork/pn_1pflow.hh
@@ -99,6 +99,7 @@ private:
   RF phiSolid = 0;
 
   RF PoreReferenceLength;
+  RF initialPhiSolid = 0;
 
   std::unique_ptr<NS_LOP> nsLop;
   std::unique_ptr<CH_LOP> chLop;
@@ -140,7 +141,7 @@ private:
     RIND rInd(gv,gfs.ch,param.pTree.sub("Refinement"));
     rInd.calculateIndicator(gv,chV,false);
     rInd.markGrid(*gridp,2);
-    rInd.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,1);
+    rInd.refineGrid(*gridp,chV,nsV,gfs.ch,gfs.ns,0);
     // recompute constraints
     gfs.updateConstraintsContainer(bdry);
     interpolateToBdry(gfs,initial.iniNS,initial.iniCH,nsV,chV);
@@ -196,6 +197,21 @@ private:
     phiSolid = integral[0];
   }
 
+  RF calculateGrowthSpeed()
+  {
+    auto productFunction = Dune::PDELab::makeGridFunctionFromCallable (gv, [&](const auto& element, const auto& xlocal){
+      Dune::FieldVector<double,1> phi1;
+      dgf->phiDgf.evaluate(element, xlocal, phi1);
+
+      RF Rspeed = phi1[0] < 0.95 ? param.dw.eval_q(phi1[0]) : 0.0f;
+      RF R1 = -Rspeed/param.eps*param.reaction.eval(1);
+      return Dune::FieldVector<double,1>(-R1);
+    });
+    Dune::FieldVector<double,1> integral;
+    Dune::PDELab::integrateGridFunction(productFunction,integral,2);
+    return integral[0];
+  }
+
 //------------------------------------------------------------------------------------
 
 public:
@@ -256,6 +272,7 @@ public:
     initializeGridOperators();
     applyNs();
     updatePhiSolid();
+    initialPhiSolid = phiSolid;
     updateKf();
   }
 
@@ -264,12 +281,12 @@ public:
     vtkwriter->write(time,Dune::VTK::ascii);
   }
 
-  void grow(RF length)
+  void grow(RF length)  // How much the the solid-fluid interface should grow in normal direction (in m)
   {
     length = length/PoreReferenceLength;
     if(length <= 0)
     {
-      std::cout << "Need positive length" << std::endl;
+      std::cout << "Need positive length (dissolution not yet implemented)" << std::endl;
       return;
     }
     RF dtmax = param.pTree.get("Time.dtMax",0.01);
@@ -287,6 +304,31 @@ public:
     updatePhiSolid();
   }
 
+  void growArea(RF area)  // How much Area the solid phase should grow (in m^2)
+  {
+    RF phiSolid_beforeGrowth = phiSolid;
+    area = area/(PoreReferenceLength*PoreReferenceLength*4);  //4 because of Symmetry
+    if(area <= 0)
+    {
+      std::cout << "Need positive area (dissolution not yet implemented)" << std::endl;
+      return;
+    }
+    RF dtmax = param.pTree.get("Time.dtMax",0.01);
+
+    while(area - (phiSolid - phiSolid_beforeGrowth) > dtmax*calculateGrowthSpeed())
+    {
+      //writeVTK(1-length);
+      param.time.dt = dtmax;
+      applyCh(0);
+      updatePhiSolid();
+    }
+    param.time.dt = (area - (phiSolid - phiSolid_beforeGrowth)) / calculateGrowthSpeed();
+    applyCh(0);
+    applyNs(0);
+    updateKf();
+    updatePhiSolid();
+  }
+
   RF getKf()
   {
     return kf;
@@ -297,6 +339,16 @@ public:
     return phiSolid;
   }
 
+  RF getPoreArea()
+  {
+    return 4*PoreReferenceLength*PoreReferenceLength*(1-phiSolid);  //4 because of Symmetry
+  }
+
+  RF getPoreAreaFraction()
+  {
+    return (1-phiSolid)/(1-initialPhiSolid);
+  }
+
   void setPoreBaseRadius(RF PoreBaseRadius_)
   {
     PoreReferenceLength = PoreBaseRadius_ / (1-param.initial.ptree.get("Radius",0.2) );
diff --git a/src/pntest/pntest.cc b/src/pntest/pntest.cc
index 42fb795384d16f90adeafca7da140faadcf3488d..7c239ff6ee81a17c69f473b2b79be1af573294a5 100644
--- a/src/pntest/pntest.cc
+++ b/src/pntest/pntest.cc
@@ -63,12 +63,20 @@ int main(int argc, char** argv)
 
     //Example of usecase for Pn_1PFlow
     std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
+    std::cout << "Pore Area = " << pn_1pFlow.getPoreArea() << " m^2" << std::endl;
+    std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+
     pn_1pFlow.writeVTK(0.0);
-    pn_1pFlow.grow(2e-6);  //Growth is measured in m
+    //pn_1pFlow.grow(2e-6);  //Growth is measured in m
+    pn_1pFlow.growArea((10e-6)*(10e-6)*0.1);   // Growth of Area in m^2
     pn_1pFlow.writeVTK(1.0);
     std::cout << "Kf = " << pn_1pFlow.getKf() << " phiSolid = " << pn_1pFlow.getPhiSolid() << std::endl;
     std::cout << "Total Flux = - Kf/viscosity * gradient p" << std::endl;
     std::cout << "Kf has unit m^4" << std::endl;
+
+    std::cout << "Pore Area = " << pn_1pFlow.getPoreArea()  << " m^2" << std::endl;
+    std::cout << "Pore Area Fraction (compared to inital values) = " << pn_1pFlow.getPoreAreaFraction() << std::endl;
+
   }
   catch (Dune::Exception &e){
     std::cerr << "Dune reported error: " << e << std::endl;