diff --git a/src/constitutive_relations/eos/eos_evaluator.cc b/src/constitutive_relations/eos/eos_evaluator.cc
index 96e00de78..064b545f2 100644
--- a/src/constitutive_relations/eos/eos_evaluator.cc
+++ b/src/constitutive_relations/eos/eos_evaluator.cc
@@ -348,6 +348,19 @@ EOSEvaluator::EvaluatePartialDerivative_(const State& S,
   }
 }
 
+void
+EOSEvaluator::EnsureCompatibility_ToDeps_(State& S){
+
+  auto akeytag = my_keys_.front();
+  const auto& fac =
+    S.Require<CompositeVector, CompositeVectorSpace>(akeytag.first, akeytag.second);
+
+  for (const auto& dep : dependencies_) {
+    auto& dep_fac = S.Require<CompositeVector, CompositeVectorSpace>(dep.first, dep.second);
+    dep_fac.UpdateSameNumDofs(fac);
+  }
+
+}
 
 } // namespace Relations
 } // namespace Amanzi
diff --git a/src/constitutive_relations/eos/eos_evaluator.hh b/src/constitutive_relations/eos/eos_evaluator.hh
index a6b640a41..d6e4d68d2 100644
--- a/src/constitutive_relations/eos/eos_evaluator.hh
+++ b/src/constitutive_relations/eos/eos_evaluator.hh
@@ -61,6 +61,8 @@ class EOSEvaluator : public EvaluatorSecondaryMonotypeCV {
   }
 
 
+   virtual void EnsureCompatibility_ToDeps_(State& S) override; 
+  
   // Required methods from EvaluatorSecondaryMonotypeCV
   virtual void Evaluate_(const State& S, const std::vector<CompositeVector*>& results) override;
   virtual void EvaluatePartialDerivative_(const State& S,
diff --git a/src/pks/transport/transport_ats.hh b/src/pks/transport/transport_ats.hh
index f42c2380f..7e5218ff0 100644
--- a/src/pks/transport/transport_ats.hh
+++ b/src/pks/transport/transport_ats.hh
@@ -372,7 +372,7 @@ class Transport_ATS : public PK_Physical_Default {
   std::vector<Teuchos::RCP<TransportDomainFunction>> bcs_;
 
   // operators for dispersion/diffusion
-  bool has_diffusion_, has_dispersion_;
+  bool has_diffusion_, has_dispersion_, enforce_bc_;
   Teuchos::RCP<TensorVector> D_; // workspace, disp + diff
   Teuchos::RCP<Operators::BCs> diff_bcs_;
   Teuchos::RCP<Operators::Operator> diff_global_op_;
diff --git a/src/pks/transport/transport_ats_pk.cc b/src/pks/transport/transport_ats_pk.cc
index 057bd2d09..6e0a6d99f 100644
--- a/src/pks/transport/transport_ats_pk.cc
+++ b/src/pks/transport/transport_ats_pk.cc
@@ -63,7 +63,8 @@ Transport_ATS::Transport_ATS(Teuchos::ParameterList& pk_tree,
     dt_stable_(-1.0),
     dt_max_(-1.0),
     has_diffusion_(false),
-    has_dispersion_(false)
+    has_dispersion_(false),
+    enforce_bc_(false)
 {
   // initialize io
   units_.Init(global_plist->sublist("units"));
@@ -290,12 +291,14 @@ Transport_ATS::SetupTransport_()
   if (has_dispersion_ || has_diffusion_) {
     // default boundary conditions (none inside domain and Neumann on its boundary)
     diff_bcs_ = Teuchos::rcp(
-      new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR));
+                             new Operators::BCs(mesh_, AmanziMesh::Entity_kind::FACE, WhetStone::DOF_Type::SCALAR));
+
     PopulateBoundaryData_(-1, *diff_bcs_);
 
     // diffusion operator
     Operators::PDE_DiffusionFactory opfactory;
     Teuchos::ParameterList& op_list = plist_->sublist("diffusion");
+    enforce_bc_ = op_list.get("enforce boundary conditions", false);
     diff_op_ = opfactory.Create(op_list, mesh_, diff_bcs_);
     diff_global_op_ = diff_op_->global_operator();
     diff_acc_op_ =
@@ -306,6 +309,7 @@ Transport_ATS::SetupTransport_()
     diff_sol_ = Teuchos::rcp(new CompositeVector(cvs));
   }
 
+  
   // source term setup
   // --------------------------------------------------------------------------------
   if (plist_->isSublist("source terms")) {
@@ -449,6 +453,7 @@ Transport_ATS::SetupTransport_()
     auto bcs_list = Teuchos::sublist(plist_, "boundary conditions");
     auto conc_bcs_list = Teuchos::sublist(bcs_list, "mole fraction");
 
+    int m = 0;    
     for (const auto& it : *conc_bcs_list) {
       std::string name = it.first;
       if (conc_bcs_list->isSublist(name)) {
@@ -503,13 +508,13 @@ Transport_ATS::SetupTransport_()
 
           // set the component indicies
           for (const auto& n : bc->tcc_names()) {
-            bc->tcc_index().push_back(FindComponentNumber_(n));
+            bc->tcc_index().push_back(FindComponentNumber_(n));            
           }
           bcs_.push_back(bc);
         }
       }
     }
-
+    
 #ifdef ALQUIMIA_ENABLED
     // -- try geochemical conditions
     auto geochem_list = Teuchos::sublist(bcs_list, "geochemical");
@@ -950,7 +955,7 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new)
 {
   if (!has_diffusion_ && !has_dispersion_) return;
   double dt = t_new - t_old;
-
+  
   Epetra_MultiVector& tcc_new =
     *S_->GetW<CompositeVector>(key_, tag_next_, passwd_).ViewComponent("cell", false);
 
@@ -983,6 +988,9 @@ Transport_ATS ::AdvanceDispersionDiffusion_(double t_old, double t_new)
   double md_old(0.0);
 
   for (int i = 0; i != num_aqueous_; ++i) {
+
+    // Set Dirichlet_BC for 
+    if (enforce_bc_) PopulateBoundaryData_(i, *diff_bcs_);
     // add molecular diffusion to the dispersion tensor
     bool changed_tensor(false);
     if (has_diffusion_) {
@@ -1346,8 +1354,7 @@ Transport_ATS::PopulateBoundaryData_(int component, Operators::BCs& bc)
   if (component >= 0) {
     for (int m = 0; m < bcs_.size(); m++) {
       std::vector<int>& tcc_index = bcs_[m]->tcc_index();
-      int ncomp = tcc_index.size();
-
+      int ncomp = tcc_index.size();                                        
       for (auto it = bcs_[m]->begin(); it != bcs_[m]->end(); ++it) {
         int f = it->first;
         std::vector<double>& values = it->second;