diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp
index dafd31fd..0b7cd3a5 100644
--- a/analysis/SubPhase.cpp
+++ b/analysis/SubPhase.cpp
@@ -430,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 * Porosity * not_water_flow_rate / force_mag;
-        double krw = h * h * nu_w * Porosity* water_flow_rate / force_mag;
+        double krn = h * h * nu_n * Porosity* Porosity * not_water_flow_rate / force_mag;
+        double krw = h * h * nu_w * Porosity* Porosity* water_flow_rate / force_mag;
         /* not counting films */
-        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 krnf = krn - h * h * nu_n * Porosity* Porosity * not_water_film_flow_rate / force_mag;
+        double krwf = krw - h * h * nu_w * Porosity* Porosity * water_film_flow_rate / force_mag;
         double eff_pressure = 1.0 / (krn + krw); // effective pressure drop
 
         fprintf(TIMELOG,
diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 04aa7a4b..f2ec9316 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -912,31 +912,31 @@ double ScaLBL_ColorModel::Run(int returntime) {
             				(vBd_x * dir_x + vBd_y * dir_y + vBd_z * dir_z);
 
             		double kAeff_connected =
-            				h * h * muA * Mask->Porosity()* flow_rate_A_connected / (force_mag);
+            				h * h * muA * Mask->Porosity()*  Mask->Porosity()* flow_rate_A_connected / (force_mag);
             		double kBeff_connected =
-            				h * h * muB * Mask->Porosity()* flow_rate_B_connected / (force_mag);
+            				h * h * muB * Mask->Porosity()* Mask->Porosity()*  flow_rate_B_connected / (force_mag);
 
             		// Saturation normalized effective permeability to account for decoupled phases and
             		// effective porosity.
             		double kAeff_connected_low =
-            				(1.0 - current_saturation) * h * h * muA* Mask->Porosity() *
+            				(1.0 - current_saturation) * h * h * muA* Mask->Porosity() *Mask->Porosity() *
             				flow_rate_A_connected / (force_mag);
             		double kBeff_connected_low = current_saturation * h * h *
-            				muB * Mask->Porosity()* flow_rate_B_connected /
+            				muB * Mask->Porosity()* Mask->Porosity() *flow_rate_B_connected /
             				(force_mag);
 
             		double kAeff_disconnected =
-            				h * h * muA * Mask->Porosity()* flow_rate_A_disconnected / (force_mag);
+            				h * h * muA * Mask->Porosity()*Mask->Porosity() * flow_rate_A_disconnected / (force_mag);
             		double kBeff_disconnected =
-            				h * h * muB * Mask->Porosity()* flow_rate_B_disconnected / (force_mag);
+            				h * h * muB * Mask->Porosity()*Mask->Porosity() * flow_rate_B_disconnected / (force_mag);
 
-            		double kAeff = h * h * muA* Mask->Porosity() * (flow_rate_A) / (force_mag);
-            		double kBeff = h * h * muB* Mask->Porosity() * (flow_rate_B) / (force_mag);
+            		double kAeff = h * h * muA* Mask->Porosity() *Mask->Porosity() * (flow_rate_A) / (force_mag);
+            		double kBeff = h * h * muB* Mask->Porosity() *Mask->Porosity() * (flow_rate_B) / (force_mag);
 
 
             		/* not counting films */
-            		double krnf = kAeff - h * h * muA* Mask->Porosity() * not_water_film_flow_rate / force_mag;
-            		double krwf = kBeff - h * h * muB* Mask->Porosity() * water_film_flow_rate / force_mag;
+            		double krnf = kAeff - h * h * muA* Mask->Porosity()* Mask->Porosity() * not_water_film_flow_rate / force_mag;
+            		double krwf = kBeff - h * h * muB* Mask->Porosity()*  Mask->Porosity() *water_film_flow_rate / force_mag;
 
             		/* not counting films */
             		double krnf_low = (1.0 - current_saturation) * krnf;
@@ -946,7 +946,7 @@ double ScaLBL_ColorModel::Run(int returntime) {
             		// effective porosity.
             		double kAeff_low = (1.0 - current_saturation) * h * h *
             				muA* Mask->Porosity() * (flow_rate_A) / (force_mag);
-            		double kBeff_low = current_saturation * h * h * muB* Mask->Porosity() *
+            		double kBeff_low = current_saturation * h * h * muB* Mask->Porosity() *Mask->Porosity() *
             				(flow_rate_B) / (force_mag);
 
             		double viscous_pressure_drop =
diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp
index b4db467c..b5973c2a 100644
--- a/models/MRTModel.cpp
+++ b/models/MRTModel.cpp
@@ -387,7 +387,7 @@ void ScaLBL_MRTModel::Run() {
 
             double h = Dm->voxel_length;
             double absperm =
-                h * h * mu * Mask->Porosity() * flow_rate / force_mag;
+                h * h * mu * Mask->Porosity() *  Mask->Porosity()* flow_rate / force_mag;
             if (rank == 0) {
                 printf("     %f\n", absperm);
                 FILE *log_file = fopen("Permeability.csv", "a");