diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp
index 790dd808..9d96a667 100644
--- a/models/ColorModel.cpp
+++ b/models/ColorModel.cpp
@@ -2329,7 +2329,7 @@ double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_d
 		ComputeGlobalBlobIDs(Nx-2,Ny-2,Nz-2,rank_info,phase,M.Averages->SDs,vF,vS,phase_label,M.Dm->Comm);
 		M.Dm->Comm.barrier();
 
-		// only operate on component "0"
+		// only operate on component "0"ScaLBL_ColorModel &M, 
 		count = 0.0;
 
 		for (int k=0; k<Nz; k++){
@@ -2474,3 +2474,80 @@ double FlowAdaptor::ShellAggregation(ScaLBL_ColorModel &M, const double target_d
 	return delta_volume;
 }
 
+
+double FlowAdaptor::SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil){
+  srand(time(NULL));
+  auto rank = M.rank;
+  auto Np = M.Np; 
+  double mass_loss =0.f;
+  double count =0.f;
+  double *Aq_tmp, *Bq_tmp;
+  
+  Aq_tmp = new double [7*Np];
+  Bq_tmp = new double [7*Np];
+
+  ScaLBL_CopyToHost(Aq_tmp, M.Aq, 7*Np*sizeof(double));
+  ScaLBL_CopyToHost(Bq_tmp, M.Bq, 7*Np*sizeof(double));
+  
+ 
+   for (int n=0; n < M.ScaLBL_Comm->LastExterior(); n++){
+    double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
+    double dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+    double dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+    double phase_id = (dA - dB) / (dA + dB);
+    if (phase_id > 0.0){
+      Aq_tmp[n] -= 0.3333333333333333*random_value;
+      Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
+      
+      Bq_tmp[n] += 0.3333333333333333*random_value;
+      Bq_tmp[n+Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
+    }
+    mass_loss += random_value*seed_water_in_oil;
+  }
+
+  for (int n=M.ScaLBL_Comm->FirstInterior(); n < M.ScaLBL_Comm->LastInterior(); n++){
+    double random_value = seed_water_in_oil*double(rand())/ RAND_MAX;
+    double dA = Aq_tmp[n] + Aq_tmp[n+Np]  + Aq_tmp[n+2*Np] + Aq_tmp[n+3*Np] + Aq_tmp[n+4*Np] + Aq_tmp[n+5*Np] + Aq_tmp[n+6*Np];
+    double dB = Bq_tmp[n] + Bq_tmp[n+Np]  + Bq_tmp[n+2*Np] + Bq_tmp[n+3*Np] + Bq_tmp[n+4*Np] + Bq_tmp[n+5*Np] + Bq_tmp[n+6*Np];
+    double phase_id = (dA - dB) / (dA + dB);
+    if (phase_id > 0.0){
+      Aq_tmp[n] -= 0.3333333333333333*random_value;
+      Aq_tmp[n+Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+2*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+3*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+4*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+5*Np] -= 0.1111111111111111*random_value;
+      Aq_tmp[n+6*Np] -= 0.1111111111111111*random_value;
+      
+      Bq_tmp[n] += 0.3333333333333333*random_value;
+      Bq_tmp[n+Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+2*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+3*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+4*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+5*Np] += 0.1111111111111111*random_value;
+      Bq_tmp[n+6*Np] += 0.1111111111111111*random_value;
+    }
+    mass_loss += random_value*seed_water_in_oil;
+  }
+
+  count= M.Dm->Comm.sumReduce(  count);
+  mass_loss= M.Dm->Comm.sumReduce(  mass_loss);
+  if (rank == 0) printf("Remove mass %f from %f voxels \n",mass_loss,count);
+
+  // Need to initialize Aq, Bq, Den, Phi directly
+  //ScaLBL_CopyToDevice(Phi,phase.data(),7*Np*sizeof(double));
+  ScaLBL_CopyToDevice(M.Aq, Aq_tmp, 7*Np*sizeof(double));
+  ScaLBL_CopyToDevice(M.Bq, Bq_tmp, 7*Np*sizeof(double));
+
+  return(mass_loss);
+}
diff --git a/models/ColorModel.h b/models/ColorModel.h
index 2b76ef2f..d69e3980 100644
--- a/models/ColorModel.h
+++ b/models/ColorModel.h
@@ -101,6 +101,7 @@ public:
 	double ImageInit(ScaLBL_ColorModel &M, std::string Filename);
 	double ShellAggregation(ScaLBL_ColorModel &M, const double delta_volume);
 	double UpdateFractionalFlow(ScaLBL_ColorModel &M);
+	double SeedPhaseField(ScaLBL_ColorModel &M, const double seed_water_in_oil);
 	void Flatten(ScaLBL_ColorModel &M);
 	DoubleArray phi;
 	DoubleArray phi_t;
diff --git a/tests/lbpm_color_simulator.cpp b/tests/lbpm_color_simulator.cpp
index e3450ee9..a82f134f 100644
--- a/tests/lbpm_color_simulator.cpp
+++ b/tests/lbpm_color_simulator.cpp
@@ -83,12 +83,17 @@ int main( int argc, char **argv )
 			int MAX_STEADY_TIME = 1000000;
 			double ENDPOINT_THRESHOLD = 0.1;
 			double FRACTIONAL_FLOW_INCREMENT = 0.0; // this will skip the flow adaptor if not enabled
+			double SEED_WATER = 0.0;
 			if (ColorModel.db->keyExists( "FlowAdaptor" )){
 				auto flow_db = ColorModel.db->getDatabase( "FlowAdaptor" );
 				MAX_STEADY_TIME = flow_db->getWithDefault<int>( "max_steady_timesteps", 1000000 );
 				SKIP_TIMESTEPS = flow_db->getWithDefault<int>( "skip_timesteps", 50000 );
-				FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
 				ENDPOINT_THRESHOLD = flow_db->getWithDefault<double>( "endpoint_threshold", 0.1);
+				/* protocol specific key values */
+				if (PROTOCOL == "fractional flow")
+					FRACTIONAL_FLOW_INCREMENT = flow_db->getWithDefault<double>( "fractional_flow_increment", 0.05);
+				if (PROTOCOL == "seed water")
+					SEED_WATER = flow_db->getWithDefault<double>( "seed_water", 0.01);
 			}
 			/* analysis keys*/
 			int ANALYSIS_INTERVAL = ColorModel.timestepMax;
@@ -121,6 +126,7 @@ int main( int argc, char **argv )
 						ColorModel.timestep =  ColorModel.timestepMax;
 					}
 				}
+				/*********************************************************/
 				/* update the fluid configuration with the flow adapter */
 				int skip_time = 0;
 				timestep = ColorModel.timestep;
@@ -149,6 +155,9 @@ int main( int argc, char **argv )
 							double target_volume_change = FRACTIONAL_FLOW_INCREMENT*initialSaturation - SaturationChange;
 							Adapt.ShellAggregation(ColorModel,target_volume_change);
 						}
+						else if (PROTOCOL == "seed water"){
+							Adapt.SeedPhaseField(ColorModel,SEED_WATER);
+						}
 						/* Run some LBM timesteps to let the system relax a bit */
 						MLUPS = ColorModel.Run(timestep);
 						/* Recompute the volume fraction now that the system has adjusted */
@@ -157,24 +166,13 @@ int main( int argc, char **argv )
 						SaturationChange = volB/(volA + volB) - initialSaturation;
 						skip_time += ANALYSIS_INTERVAL;
 					}
-					if (PROTOCOL == "fractional flow")	{							
-						if (rank==0) printf("  *********************************************************************  \n");
-						if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
-						if (rank==0) printf("  *********************************************************************  \n");
-					}
+					if (rank==0) printf("  *********************************************************************  \n");
+					if (rank==0) printf("   Updated fractional flow with saturation change = %f  \n", SaturationChange);
+					if (rank==0) printf("   Used protocol = %s  \n", PROTOCOL.c_str());
+					if (rank==0) printf("  *********************************************************************  \n");
 				}
-				/* apply timestep skipping algorithm to accelerate steady-state */
-				/* skip_time = 0;
-            	timestep = ColorModel.timestep;
-            	while (skip_time < SKIP_TIMESTEPS){
-            		timestep += ANALYSIS_INTERVAL;
-                	MLUPS = ColorModel.Run(timestep);
-            		Adapt.MoveInterface(ColorModel);
-            		skip_time += ANALYSIS_INTERVAL;
-            	}
-				 */
+				/*********************************************************/
 			}
-			//ColorModel.WriteDebug();
 		}
 
 		else