diff --git a/src/reactmicp/systems/saturated_react/transport_program.cpp b/src/reactmicp/systems/saturated_react/transport_program.cpp
index b63a249..f28b8b8 100644
--- a/src/reactmicp/systems/saturated_react/transport_program.cpp
+++ b/src/reactmicp/systems/saturated_react/transport_program.cpp
@@ -1,278 +1,288 @@
 #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, {}, {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,
                                           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();
     // 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(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)
     {
         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) // 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);
 
+    // diffusion
+    scalar_t flux_diffusion = flux_coeff*(displacement(dof_0) - displacement(dof_1));
 
-    scalar_t flux = flux_coeff*(displacement(dof_0) - displacement(dof_1));
+    element_residual(0) = flux_diffusion;
+    element_residual(1) = - flux_diffusion;
 
-    element_residual(0) = flux;
-    element_residual(1) = - flux;
+    // advection
+
+    if (m_variables->fluid_velocity(element) != 0.0)
+    {
+        scalar_t flux_advection = (m_mesh->get_face_area(element))
+                *m_variables->fluid_velocity(element);
+        if (m_variables->fluid_velocity(element) > 0)
+        {
+            flux_advection *= (displacement(dof_0) - displacement(dof_1));
+            element_residual(1) += flux_advection;
+        }
+        else
+        {
+            flux_advection *= (displacement(dof_1) - displacement(dof_0));
+            element_residual(0) -= flux_advection;
+        }
+    }
 
     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