diff --git a/src/reactmicp/systems/saturated_react/transport_program.cpp b/src/reactmicp/systems/saturated_react/transport_program.cpp
index 5e4ce0d..b63a249 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.cpp
@@ -1,241 +1,278 @@
 #include "transport_program.hpp"
 #include "dfpm/meshes/mesh1d.hpp"
 #include "variables.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 //class SaturatedDiffusion::
 
 
 SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
                    std::vector<index_t> list_fixed_nodes):
     m_mesh(variables->get_mesh()),
     m_variables(variables),
     m_is_in_residual_computation(false)
 {
-    number_equations(list_fixed_nodes);
+    number_equations(list_fixed_nodes, {}, {0,});
+}
+
+
+SaturatedDiffusion::SaturatedDiffusion(SaturatedVariablesPtr variables,
+                   std::vector<index_t> list_fixed_nodes,
+                   std::map<index_t, index_t>  list_slave_nodes,
+                   std::vector<index_t> list_immobile_components):
+    m_mesh(variables->get_mesh()),
+    m_variables(variables),
+    m_is_in_residual_computation(false)
+{
+    number_equations(list_fixed_nodes, list_slave_nodes, list_immobile_components);
 }
 
 //! \brief Return the number of degrees of freedom per node
 index_t SaturatedDiffusion::get_ndf() const {return 2*m_variables->nb_component();}
 
 //! \brief Return the total number of degrees of freedom
 index_t SaturatedDiffusion::get_tot_ndf() const {return 2*m_variables->nb_component()*m_variables->get_mesh()->nb_nodes();}
 
 
-void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes)
+void SaturatedDiffusion::number_equations(std::vector<index_t> list_fixed_nodes,
+                                          std::map<index_t, index_t>  list_slave_nodes,
+                                          std::vector<index_t> list_immobile_components)
 {
     m_ideq.resizeLike(m_variables->displacement());
     m_ideq.setZero();
-    for (auto it=list_fixed_nodes.begin(); it<list_fixed_nodes.end(); ++it)
+    // flag fixed nodes
+    for (index_t node: list_fixed_nodes)
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
-            m_ideq(m_variables->dof_aqueous_concentration(*it, component)) = no_equation;
+            m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
         }
     }
-
+    // flag slaves nodes
+    // we flag them by making their ideq more negative than no_equation
+    for (auto slave_pair: list_slave_nodes)
+    {
+        for (index_t component=1; component<m_variables->nb_component(); ++component)
+        {
+            m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) = no_equation-1;
+        }
+    }
+    // set equation numbers
     index_t neq = 0;
     for (index_t node=0; node<m_mesh->nb_nodes(); ++node)
     {
-        m_ideq(m_variables->dof_aqueous_concentration(node, 0)) = no_equation;
-        for (index_t component=1; component<m_variables->nb_component(); ++component)
+        for (index_t component: list_immobile_components)
+        {
+            m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
+            m_ideq(m_variables->dof_aqueous_concentration(node, component)) = no_equation;
+        }
+        for (index_t component=0; component<m_variables->nb_component(); ++component)
         {
             index_t dof = m_variables->dof_aqueous_concentration(node, component);
-            if (m_ideq(dof) != no_equation)
+            if (m_ideq(dof) > no_equation) // don't attribute an equation number if it is a slave or a fixed node
             {
                 m_ideq(dof) = neq;
                 ++neq;
             }
             dof = m_variables->dof_solid_concentration(node, component);
             m_ideq(dof) = no_equation;
         }
     }
+    // slave nodes
+    // attribute the correct equation number
+    for (auto slave_pair: list_slave_nodes)
+    {
+        for (index_t component=1; component<m_variables->nb_component(); ++component)
+        {
+            m_ideq(m_variables->dof_aqueous_concentration(slave_pair.first, component)) =
+                    m_ideq(m_variables->dof_aqueous_concentration(slave_pair.second, component));
+        }
+    }
     m_neq = neq;
-
 }
 
 void SaturatedDiffusion::compute_residuals(
         const Vector& displacement,
         const Vector& velocity,
         Vector& residual)
 {
     residual = Vector::Zero(get_neq());
     m_is_in_residual_computation = true;
     for (index_t element: m_mesh->range_elements())
     {
         for (index_t component=1; component<m_variables->nb_component(); ++component)
         {
             Eigen::Vector2d element_residual;
             element_residual.setZero();
             residuals_element_component(element, component, displacement, velocity, element_residual);
             for (index_t en=0; en<2; ++en)
             {
                 const index_t node = m_mesh->get_node(element, en);
                 const index_t id = m_ideq(m_variables->dof_aqueous_concentration(node, component));
                 if (id != no_equation) {residual(id) += element_residual(en);}
             }
         }
     }
     m_is_in_residual_computation = false;
 }
 
 void SaturatedDiffusion::residuals_element_component(
         index_t element,
         index_t component,
         const Vector& displacement,
         const Vector& velocity,
         Eigen::Vector2d& element_residual
         )
 {
     const scalar_t mass_coeff_0 = -m_mesh->get_volume_cell_element(element, 0);
     const scalar_t mass_coeff_1 = -m_mesh->get_volume_cell_element(element, 1);
 
     const index_t node_0 = m_mesh->get_node(element, 0);
     const index_t node_1 = m_mesh->get_node(element, 1);
 
     const scalar_t diff_coeff = 1.0/(0.5/m_variables->diffusion_coefficient(node_0) +
                                      0.5/m_variables->diffusion_coefficient(node_1));
 
 
     scalar_t flux_coeff = -( m_mesh->get_face_area(element) / m_mesh->get_dx(element)
                            * diff_coeff
                            ) ;
 
     const index_t dof_0 = m_variables->dof_aqueous_concentration(node_0, component);
     const index_t dof_1 = m_variables->dof_aqueous_concentration(node_1, component);
 
 
     scalar_t flux = flux_coeff*(displacement(dof_0) - displacement(dof_1));
 
     element_residual(0) = flux;
     element_residual(1) = - flux;
 
     if (m_is_in_residual_computation)
     {
         m_variables->aqueous_concentration(node_0, component,
                            m_variables->transport_rate()) += element_residual(0);
         m_variables->aqueous_concentration(node_1, component,
                            m_variables->transport_rate()) += element_residual(1);
     }
     // velocity
     element_residual(0) += mass_coeff_0*(velocity(dof_0)*m_variables->porosity(node_0)
                                          +m_variables->vel_porosity(node_0)*displacement(dof_0));
     element_residual(1) += mass_coeff_1*(velocity(dof_1)*m_variables->porosity(node_1)
                                          + m_variables->vel_porosity(node_1)*displacement(dof_1));
     // external rate
     element_residual(0) +=  mass_coeff_0*m_variables->solid_concentration(node_0, component,
                 m_variables->chemistry_rate());
     element_residual(1) +=  mass_coeff_1*m_variables->solid_concentration(node_1, component,
                 m_variables->chemistry_rate());
 
 
 //    for (index_t en=0; en<2; ++en)
 //    {
 //        element_residual(en) +=  mass_coeff_0*m_variables->solid_concentration(
 //                    m_mesh->get_node(element, en), component,
 //                    m_variables->chemistry_rate());
 //    }
 }
 
 
 void SaturatedDiffusion::compute_jacobian(
         Vector& displacement,
         Vector& velocity,
         Eigen::SparseMatrix<scalar_t>& jacobian,
         scalar_t alphadt
         )
 {
     dfpm::list_triplet_t jacob;
     const index_t ncomp = m_variables->nb_component();
     const index_t estimation = m_mesh->nb_nodes()*(ncomp*m_mesh->nen);
     jacob.reserve(estimation);
     for (index_t element: m_mesh->range_elements())
     {
         jacobian_element(element, displacement, velocity, jacob, alphadt);
     }
     jacobian = Eigen::SparseMatrix<scalar_t>(get_neq(), get_neq());
     jacobian.setFromTriplets(jacob.begin(), jacob.end());
 }
 
 void SaturatedDiffusion::jacobian_element(
         index_t element,
         Vector& displacement,
         Vector& velocity,
         dfpm::list_triplet_t& jacobian,
         scalar_t alphadt)
 {
     for (index_t component=1; component<m_variables->nb_component(); ++component)
     {
         Eigen::Vector2d element_residual_orig(Eigen::Vector2d::Zero());
         residuals_element_component(element, component, displacement, velocity, element_residual_orig);
 
         for (index_t en=0; en<2; ++en)
         {
             Eigen::Vector2d element_residual(Eigen::Vector2d::Zero());
 
             const index_t node = m_mesh->get_node(element, en);
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t idc = m_ideq(dof);
 
             if (idc == no_equation) continue;
 
             const scalar_t tmp_v = velocity(dof);
             const scalar_t tmp_d = displacement(dof);
 
             scalar_t h = eps_jacobian*std::abs(tmp_v);
             if (h < 1e-4*eps_jacobian) h = eps_jacobian;
             velocity(dof) = tmp_v + h;
             h = velocity(dof) - tmp_v;
 
             displacement(dof) = tmp_d + alphadt*h;
 
             residuals_element_component(element, component, displacement, velocity, element_residual);
             velocity(dof) = tmp_v;
             displacement(dof) = tmp_d;
 
             for (index_t enr=0; enr<2; ++enr)
             {
                 const index_t noder = m_mesh->get_node(element, enr);
                 const index_t idr = m_ideq(m_variables->dof_aqueous_concentration(noder, component));
 
                 if (idr == no_equation) continue;
                 jacobian.push_back(dfpm::triplet_t(
                               idr,
                               idc,
                               (element_residual(enr) - element_residual_orig(enr))/h
                               ));
             }
         }
     }
 }
 
 //! \brief Update the solutions
 void SaturatedDiffusion::update_solution(const Vector& update,
                      scalar_t lambda,
                      scalar_t alpha_dt,
                      Vector& predictor,
                      Vector& displacement,
                      Vector& velocity)
 {
     for (index_t node: m_mesh->range_nodes())
     {
         for (index_t component=1; component< m_variables->nb_component(); ++component)
         {
             const index_t dof = m_variables->dof_aqueous_concentration(node, component);
             const index_t id = m_ideq(dof);
             if (id == no_equation) continue;
             velocity(dof) += lambda*update(id);
         }
     }
     //displacement = m_variables->predictor() + alpha_dt*velocity;
     displacement = predictor + alpha_dt*velocity;
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/transport_program.hpp b/src/reactmicp/systems/saturated_react/transport_program.hpp
index 45cd2a0..166eb6e 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.hpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.hpp
@@ -1,102 +1,108 @@
 #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP
 
 #include "dfpmsolver/parabolic_program.hpp"
 
 #include "variablesfwd.hpp"
 #include "dfpm/meshes/mesh1dfwd.hpp"
 #include "dfpm/types.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 class SaturatedDiffusion:  public dfpmsolver::ParabolicProgram<SaturatedDiffusion>
 {
 public:
 
+    SaturatedDiffusion(SaturatedVariablesPtr variables,
+                       std::vector<index_t> list_fixed_nodes,
+                       std::map<index_t, index_t>  list_slave_nodes,
+                       std::vector<index_t> list_immobile_components);
+
     SaturatedDiffusion(SaturatedVariablesPtr variables,
                        std::vector<index_t> list_fixed_nodes);
 
     //! \brief Return the number of equations
     index_t get_neq() const {return m_neq;}
     //! \brief Return the number of degrees of freedom per node
     index_t get_ndf() const;
     //! \brief Return the total number of degrees of freedom
     index_t get_tot_ndf() const;
 
     //! \brief Method to update the variables
     void set_variables(SaturatedVariablesPtr variables) {m_variables = variables;}
 
     //! \brief Return the id of the equation corresponding to the degree of freedom 'id_dof'
     //!
     //! Return 'no_equation' if no equation exist
     index_t id_equation(index_t id_dof) const {return m_ideq(id_dof);}
 
     //! \brief Compute the residuals
     void compute_residuals(const Vector& displacement,
                            const Vector& velocity,
                            Vector& residual
                            );
     //! Compute the residuals inside 'element' for 'component'
     void residuals_element_component(
             index_t element,
             index_t component,
             const Vector& displacement,
             const Vector& velocity,
             Eigen::Vector2d& element_residual
             );
 
     //! \brief Compute the jacobian
     void compute_jacobian(Vector& displacement,
                           Vector& velocity,
                           Eigen::SparseMatrix<scalar_t>& jacobian,
                           scalar_t alphadt
                           );
     //! \brief Compute the contribution of 'element' in the jacobian
     void jacobian_element(
             index_t element,
             Vector& displacement,
             Vector& velocity,
             dfpm::list_triplet_t& jacobian,
             scalar_t alphadt);
 
     //! \brief Update the solutions
     void update_solution(const Vector& update,
                          scalar_t lambda,
                          scalar_t alpha_dt,
                          Vector& predictor,
                          Vector& displacement,
                          Vector& velocity);
 
     //! \brief Apply boundary conditions to the velocity vector
     //!
     //! by default do nothing.
     void apply_bc(scalar_t dt,
                   const Vector& displacement,
                   Vector& velocity)
     {}
 
 
 private:
     // number the equations
-    void number_equations(std::vector<index_t> list_fixed_nodes);
+    void number_equations(std::vector<index_t> list_fixed_nodes,
+                          std::map<index_t, index_t>  list_slave_nodes, std::vector<index_t> list_immobile_components);
 
 
     index_t m_neq;
     Eigen::Matrix<index_t, Eigen::Dynamic, 1> m_ideq;
 
     mesh::Mesh1DPtr m_mesh;
     SaturatedVariablesPtr m_variables;
 
     bool m_is_in_residual_computation;
 
 };
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTPROGRAM_HPP
diff --git a/src/reactmicp/systems/saturated_react/transport_stagger.cpp b/src/reactmicp/systems/saturated_react/transport_stagger.cpp
index 29c9af8..1569d67 100644
--- a/src/reactmicp/systems/saturated_react/transport_stagger.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_stagger.cpp
@@ -1,114 +1,125 @@
 #include "transport_stagger.hpp"
 #include "variables.hpp"
 #include "reactmicp/solver/staggers_base/stagger_structs.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 // class SaturatedTransportStagger::
 
 namespace internal {
 
 //! \brief Translate a ParabolicDriverReturnCode to a StaggerReturnCode
 solver::StaggerReturnCode Parabolic2StaggerReturnCode(dfpmsolver::ParabolicDriverReturnCode retcode)
 {
     using ParRC = dfpmsolver::ParabolicDriverReturnCode;
     using StaRC = solver::StaggerReturnCode;
 
     switch (retcode)
     {
     case ParRC::ResidualMinimized:
             return StaRC::ResidualMinimized;
             break;
     case ParRC::ErrorMinimized:
             return StaRC::ErrorMinimized;
             break;
     default: // error
         switch(retcode)
         {
         case ParRC::MaxIterations:
             return StaRC::MaximumIterationsReached;
             break;
         case ParRC::StationaryPoint:
             return StaRC::StationaryPoint;
             break;
         case ParRC::NotConvergedYet:
             return StaRC::NotConvergedYet;
             break;
         default:
             return StaRC::UnknownError;
         }
     }
 }
 
 } // end namespace internal
 
 SaturatedTransportStagger::SaturatedTransportStagger(SaturatedVariablesPtr variables,
                                                      std::vector<index_t> list_fixed_nodes):
     m_program(variables, list_fixed_nodes),
     m_solver(m_program)
 {
 
 }
 
+SaturatedTransportStagger::SaturatedTransportStagger(
+        SaturatedVariablesPtr variables,
+        std::vector<index_t> list_fixed_nodes,
+        std::map<index_t, index_t>  list_slave_nodes,
+        std::vector<index_t> list_immobile_components):
+    m_program(variables, list_fixed_nodes, list_slave_nodes, list_immobile_components),
+    m_solver(m_program)
+{
+
+}
+
 //! \brief Initialize the stagger at the beginning of an iteration
 void SaturatedTransportStagger::initialize_timestep(scalar_t dt, VariablesBasePtr var)
 {
     m_dt = dt;
 
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
     true_var->predictor() = true_var->displacement();
 
     true_var->velocity().setZero();
     true_var->transport_rate().setZero();
     m_solver.initialize_timestep(dt, true_var->displacement());
 
     Vector tmp = true_var->chemistry_rate();
     true_var->chemistry_rate().setZero();
     Eigen::VectorXd residuals;
     m_program.compute_residuals(true_var->displacement(), true_var->velocity(), residuals);
     m_residual_0 = residuals.norm();
 
     true_var->chemistry_rate() = tmp;
 
     m_program.set_variables(true_var);
 }
 
 
 //! \brief Solve the equation for the timestep
 solver::StaggerReturnCode SaturatedTransportStagger::restart_timestep(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
     m_solver.velocity() = true_var->velocity();
     dfpmsolver::ParabolicDriverReturnCode retcode = m_solver.restart_timestep(true_var->displacement());
     // copy variables if successful
     if (retcode > dfpmsolver::ParabolicDriverReturnCode::NotConvergedYet)
     {
         true_var->velocity() = m_solver.velocity();
     }
 
     return internal::Parabolic2StaggerReturnCode(retcode);
 }
 
 scalar_t SaturatedTransportStagger::get_update(VariablesBasePtr var)
 {
     return cast_var_from_base(var)->velocity().norm();
 }
 
 
 //! \brief Compute the residuals norm
 scalar_t SaturatedTransportStagger::get_residual(VariablesBasePtr var)
 {
     SaturatedVariablesPtr true_var = cast_var_from_base(var);
 
     Eigen::VectorXd residuals;
     m_program.compute_residuals(true_var->displacement(), true_var->velocity(), residuals);
     return residuals.norm();
 }
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
diff --git a/src/reactmicp/systems/saturated_react/transport_stagger.hpp b/src/reactmicp/systems/saturated_react/transport_stagger.hpp
index 4b3e2c8..52219ad 100644
--- a/src/reactmicp/systems/saturated_react/transport_stagger.hpp
+++ b/src/reactmicp/systems/saturated_react/transport_stagger.hpp
@@ -1,57 +1,62 @@
 #ifndef SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP
 #define SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP
 
 #include "variablesfwd.hpp"
 #include "reactmicp/solver/staggers_base/transport_stagger_base.hpp"
 
 #include "transport_program.hpp"
 #include "dfpmsolver/parabolic_driver.hpp"
 
 namespace specmicp {
 namespace reactmicp {
 namespace systems {
 namespace satdiff {
 
 using VariablesBasePtr = solver::VariablesBasePtr;
 
 class SaturatedTransportStagger: public solver::TransportStaggerBase
 {
 public:
 
     SaturatedTransportStagger(SaturatedVariablesPtr variables,
                               std::vector<index_t> list_fixed_nodes);
 
+    SaturatedTransportStagger(SaturatedVariablesPtr variables,
+                              std::vector<index_t> list_fixed_nodes,
+                              std::map<index_t, index_t>  list_slave_nodes,
+                              std::vector<index_t> list_immobile_components);
+
     //! \brief Return the options of the solver
     dfpmsolver::ParabolicDriverOptions& options_solver() {return m_solver.get_options();}
 
     //! \brief Initialize the stagger at the beginning of an iteration
     void initialize_timestep(scalar_t dt, VariablesBasePtr var) override;
 
     //! \brief Solve the equation for the timestep
     solver::StaggerReturnCode restart_timestep(VariablesBasePtr var) override;
 
     //! \brief Compute the residuals norm
     scalar_t get_residual(VariablesBasePtr var) override;
     //! \brief Compute the residuals norm
     scalar_t get_residual_0(VariablesBasePtr var) override {
         // TODO : a slightly more optimized version would use the value from the solver
         return m_residual_0;
     }
 
     //! \brief Obtain the update
     scalar_t get_update(VariablesBasePtr var) override;
 
 private:
     scalar_t m_dt;
     scalar_t m_residual_0;
     SaturatedDiffusion m_program;
     dfpmsolver::ParabolicDriver<SaturatedDiffusion> m_solver;
 
 };
 
 } // end namespace satdiff
 } // end namespace systems
 } // end namespace reactmicp
 } // end namespace specmicp
 
 #endif // SPECMICP_REACTMICP_SYSTEMS_SATURATED_TRANSPORTSTAGGER_HPP