diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index 92f4972d..f2e19638 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -411,6 +411,7 @@ void SubPhase::Basic() {
             dir_z = 1.0;
             force_mag = 1.0;
         }
+        double Porosity = (gwb.V + gnb.V)/Dm->Volume();
         double saturation = gwb.V / (gwb.V + gnb.V);
         double water_flow_rate =
             gwb.V * (gwb.Px * dir_x + gwb.Py * dir_y + gwb.Pz * dir_z) / gwb.M /
@@ -429,11 +430,11 @@ void SubPhase::Basic() {
         //double total_flow_rate = water_flow_rate + not_water_flow_rate;
         //double fractional_flow = water_flow_rate / total_flow_rate;
         double h = Dm->voxel_length;
-        double krn = h * h * nu_n * Mask->Porosity() * not_water_flow_rate / force_mag;
-        double krw = h * h * nu_w * Mask->Porosity()* water_flow_rate / force_mag;
+        double krn = h * h * nu_n * Porosity * not_water_flow_rate / force_mag;
+        double krw = h * h * nu_w * Porosity* water_flow_rate / force_mag;
         /* not counting films */
-        double krnf = krn - h * h * nu_n * Mask->Porosity() * not_water_film_flow_rate / force_mag;
-        double krwf = krw - h * h * nu_w * Mask->Porosity() * water_film_flow_rate / force_mag;
+        double krnf = krn - h * h * nu_n * Porosity * not_water_film_flow_rate / force_mag;
+        double krwf = krw - h * h * nu_w * Porosity * water_film_flow_rate / force_mag;
         double eff_pressure = 1.0 / (krn + krw); // effective pressure drop
 
         fprintf(TIMELOG,