diff --git a/examples/reactmicp/single/single.cpp b/examples/reactmicp/single/single.cpp
index ff9a22b..126ebd7 100644
--- a/examples/reactmicp/single/single.cpp
+++ b/examples/reactmicp/single/single.cpp
@@ -1,575 +1,701 @@
 /* =============================================================================
 
  Copyright (c)
  2014-2017 F. Georget <fabieng@princeton.edu> Princeton University
  2017-2021 F. Georget <fabien.georget@epfl.ch> EPFL
  All rights reserved.
 
 
 Redistribution and use in source and binary forms, with or without modification,
 are permitted provided that the following conditions are met:
 
     1. Redistributions of source code must retain the above copyright notice,
     this list of conditions and the following disclaimer.
 
     2. Redistributions in binary form must reproduce the above copyright notice,
     this list of conditions and the following disclaimer in the documentation
     and/or other materials provided with the distribution.
 
     3. Neither the name of the copyright holder nor the names of its
     contributors may be used to endorse or promote products derived from this
     software without specific prior written permission.
 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR
 ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
 ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
 
 ============================================================================= */
 
 #include "reactmicp/reactmicp_single.hpp"
 #include "reactmicp/solver/runner.hpp"
 
 #include <string>
 #include <iostream>
 #include <fstream>
 
 #include "specmicp_common/io/safe_config.hpp"
 
 #include <sys/stat.h>
 #include <cmath>
 #include <functional>
 
 using namespace specmicp;
 using namespace specmicp::io;
 using namespace specmicp::reactmicp::systems::single;
 
+
+// ----------------------------------------------------
+
 struct IsothermConf {};
 
+/*   Possible models for the isotherm
+ *
+ * Freundlich
+ * LeachingFreundlich
+ * BindingModelIsotherm
+ *
+ */
+
 struct FreundlichIsothermConf:
     public IsothermConf
 {
     scalar_t freund_coeff;
     scalar_t freund_exp;
 
     FreundlichIsothermConf(scalar_t coeff, scalar_t exp):
         freund_coeff(coeff), freund_exp(exp)
     {}
 };
 
 struct LeachingIsothermConf:
     public IsothermConf
 {
     scalar_t coeff_1;
     scalar_t coeff_2;
     scalar_t freund_exp;
 
 
     LeachingIsothermConf(scalar_t coeff_1, scalar_t coeff_2, scalar_t freund_exp):
         coeff_1(coeff_1), coeff_2(coeff_2), freund_exp(freund_exp)
     {}
 };
 
 struct BindingModelIsothermConf:
     public IsothermConf
 {
     scalar_t factor_csh;
     scalar_t factor_afm;
     scalar_t KCaAd;
     scalar_t KClAd;
     scalar_t KOHAd;
     scalar_t KClAFm;
     scalar_t KOHAFm;
 };
 
-struct SingleSimpleConf {
-    scalar_t rate_factor;
+// ----------------------------------------------------
+
+/*   Possible models for the isotherm
+ *
+ * Fixed
+ * Leaching
+ * Hydrating
+ *
+ */
+
+struct DeffConf {
     scalar_t diffusion_coeff;
+
+    DeffConf(scalar_t deff):
+        diffusion_coeff(deff)
+    {}
+};
+
+struct LeachedDeffConf:
+    public DeffConf
+{
     scalar_t leached_diffusion_coeff;
+
+    LeachedDeffConf(scalar_t deff, scalar_t leached_deff):
+        leached_diffusion_coeff(leached_deff), DeffConf(deff)
+    {}
+};
+
+struct HydratingDeffConf:
+        public DeffConf
+{
+    scalar_t min_diffusion_coefficient;
+    scalar_t hydration_rate;
+
+    HydratingDeffConf(scalar_t deff,
+                      scalar_t min_deff,
+                      scalar_t hydr_rate):
+        min_diffusion_coefficient(min_deff),
+        hydration_rate(hydr_rate),
+        DeffConf(deff)
+    {}
+};
+
+
+// ----------------------------------------------------
+
+struct SingleSimpleConf {
+    scalar_t rate_factor;
+
     scalar_t porosity;
     scalar_t saturation;
     std::string binding_law;
+    std::string diffusion_model;
 
     std::shared_ptr<IsothermConf> isotherm_conf;
+    std::shared_ptr<DeffConf> deff_conf;
 };
 
+
+// ----------------------------------------------------
+
+
 std::shared_ptr<SingleSimpleConf>
 configure_single_simple(
         YAMLConfigHandle&& configuration
         )
 {
     auto conf = std::make_shared<SingleSimpleConf>();
     conf->rate_factor = configuration.get_required_attribute<scalar_t>("rate_factor");
-    conf->diffusion_coeff = configuration.get_required_attribute<scalar_t>("diffusion_coeff");
-    conf->leached_diffusion_coeff = configuration.get_optional_attribute<scalar_t>("leached_diffusion_coeff", 0);
     conf->porosity = configuration.get_required_attribute<scalar_t>("porosity");
     conf->saturation = configuration.get_required_attribute<scalar_t>("saturation");
 
+    // Binding isotherm law
     conf->binding_law = configuration.get_required_attribute<std::string>("binding_law");
     if (conf->binding_law == "Freundlich") {
         conf->isotherm_conf = std::make_shared<FreundlichIsothermConf>(
                     configuration.get_required_attribute<scalar_t>("freund_coeff"),
                     configuration.get_required_attribute<scalar_t>("freund_exp")
                     );
     } else if (conf->binding_law == "LeachingFreundlich") {
         conf->isotherm_conf = std::make_shared<LeachingIsothermConf>(
                     configuration.get_required_attribute<scalar_t>("freund_coeff_high_pH"),
                     configuration.get_required_attribute<scalar_t>("freund_coeff_low_pH"),
                     configuration.get_required_attribute<scalar_t>("freund_exp")
                     );
     } else if (conf->binding_law == "BindingModelIsotherm") {
         auto bconf = std::make_shared<BindingModelIsothermConf>();
         bconf->factor_csh = configuration.get_required_attribute<scalar_t>("factor_csh");
         bconf->factor_afm = configuration.get_required_attribute<scalar_t>("factor_afm");
         bconf->KCaAd = configuration.get_required_attribute<scalar_t>("KCaAd");
         bconf->KClAd = configuration.get_required_attribute<scalar_t>("KClAd");
         bconf->KOHAd = configuration.get_required_attribute<scalar_t>("KOHAd");
         bconf->KClAFm = configuration.get_required_attribute<scalar_t>("KClAFm");
         bconf->KOHAFm = configuration.get_required_attribute<scalar_t>("KOHAFm");
         conf->isotherm_conf = bconf;
     } else {
         throw std::runtime_error("Unknow binding isotherm law");
     }
+
+    // Diffusion coefficient
+    conf->diffusion_model = configuration.get_required_attribute<std::string>("diffusion_model");
+    auto deff = configuration.get_required_attribute<scalar_t>("diffusion_coeff");
+    if (conf->diffusion_model == "Fixed") {
+        conf->deff_conf = std::make_shared<DeffConf>(deff);
+    } else if (conf->diffusion_model == "Leaching") {
+        conf->deff_conf = std::make_shared<LeachedDeffConf>(
+                    deff,
+                    configuration.get_required_attribute<scalar_t>("leached_diffusion_coeff")
+                    );
+    } else if (conf->diffusion_model == "Hydrating") {
+        conf->deff_conf = std::make_shared<HydratingDeffConf>(
+                    deff,
+                    configuration.get_required_attribute<scalar_t>("min_diffusion_coeff"),
+                    configuration.get_required_attribute<scalar_t>("hydration_rate")
+                    );
+    } else {
+        throw std::runtime_error("Unknow diffusion coefficient model");
+    }
+
     return conf;
 }
 
 scalar_t binding_law(scalar_t conc, scalar_t coeff, scalar_t exp)
 {
     return coeff*std::pow(conc, exp);
 }
 
 scalar_t ph_dependent_binding_law(scalar_t conc, scalar_t coeff1, scalar_t coeff2, scalar_t exp)
 {
     return (2*(coeff2-coeff1)*conc+coeff1)*std::pow(conc, exp);
 }
 scalar_t dph_dependent_binding_law(scalar_t conc, scalar_t coeff1, scalar_t coeff2, scalar_t exp)
 {
     return (2*(coeff2-coeff1)*(1+exp)*std::pow(conc, exp)+exp*coeff1*std::pow(conc, exp-1));
 }
 
 scalar_t competitive_langmuir(scalar_t c1, scalar_t c2, scalar_t A, scalar_t k1, scalar_t k2) {
     return A*k1*c1/(1+k1*c1+k2*c2);
 }
 
 using sfunc_t = std::function<scalar_t (scalar_t)>;
 scalar_t adiff(sfunc_t&& f, scalar_t c) {
     auto h = 1e-8*std::abs(c);
     if (h<1e-16 ) {h=1e-8;}
     auto cp = c+h;
     h = cp-c;
     auto df = (f(cp)-f(c))/h;
     return df;
 }
 
 scalar_t n_sites_csh(scalar_t oh, scalar_t factor, scalar_t KCaAd) {
     return factor*(KCaAd/(std::pow(oh, 2)+KCaAd));
 }
 scalar_t binding_csh(scalar_t cl, scalar_t oh, scalar_t KClAd, scalar_t KOHAd) {
     return (KClAd*cl)/(1+KClAd*cl+KOHAd*oh);
 }
 scalar_t ph_to_oh(scalar_t ph) {
     return std::pow(10.0, ph-14.0);
 }
 scalar_t isotherm_csh(scalar_t cl, scalar_t oh, scalar_t factor, scalar_t KCaAd, scalar_t KClAd, scalar_t KOHAd) {
     return n_sites_csh(oh, factor, KCaAd)*binding_csh(cl, oh, KClAd, KOHAd);
 }
 scalar_t isotherm_afm(scalar_t cl, scalar_t oh, scalar_t factor, scalar_t KClAFm, scalar_t KOHAFm) {
     return factor*(KClAFm*cl)/(
                 1+KClAFm*cl+KOHAFm*oh);
 }
 scalar_t isotherm(scalar_t cl, scalar_t factor_csh, scalar_t factor_afm,
                   scalar_t KCaAd, scalar_t KClAd, scalar_t KOHAd,
                   scalar_t KClAFm, scalar_t KOHAFm)
 {
     auto ph = cl*(12.5-13.5)/0.5 + 13.5;
     auto oh = ph_to_oh(ph);
 
     return isotherm_csh(cl, oh, factor_csh, KCaAd, KClAd, KOHAd)+isotherm_afm(cl, oh, factor_afm, KClAFm, KOHAFm);
 }
 
 scalar_t simple_newton(
         std::function<scalar_t (scalar_t)> f,
         std::function<scalar_t (scalar_t)> df,
         scalar_t b, scalar_t tol, scalar_t x=0, bool print=false)
 {
     auto max_it = 200;
     scalar_t error = f(x)-b;
     auto uerror = tol+1;
     auto it = 0;
     if (print) {
         std::cout << error << " - " << it  << ";" << std::endl;
     }
     while (std::fabs(error) > tol and uerror > tol and it < max_it) {
         scalar_t sk = error/df(x);
         scalar_t x1 = 0.0;
         if ((x-sk) < 0) {
             x1 = 0.1*x;
         } else {
             x1 = x-sk;
         }
         //auto x1 = x-0.01*error/df(x);
         //if (x1 < 0) {x1 = 0.01;}
         //if (print) {
         //    std::cout << x1 << std::endl;
         //}
         uerror =  std::fabs(x1-x);
         error = f(x1)-b;
         x = x1;
         ++it;
     }
     if (print) {
         std::cout << error << " - " << it  << ";" << std::endl;
     }
 
     return x;
 }
 
 class RateFunctorBase
 {
 public:
     RateFunctorBase(scalar_t rate):
         m_rate(rate)
     {}
 
     virtual scalar_t c_eq(index_t node, SingleVariables * const var) = 0;
 
     scalar_t operator()(index_t node, scalar_t dt, SingleVariables * const var) {
         auto c_eq = this->c_eq(node, var);
         auto rate =  m_rate * (var->fluid_concentration(node) - c_eq);
         if (rate < 0) {
             rate = 0;
         }
         auto factor = var->porosity(node)*var->saturation(node);
         auto new_c = var->fluid_concentration(node)-rate/factor*dt;
         if (new_c < 0) {
             rate = factor*var->fluid_concentration(node)/dt;
         }
         return rate;
     }
 private:
     scalar_t m_rate;
 };
 
 class RateFunctor:
         public RateFunctorBase
 {
 public:
     RateFunctor(std::shared_ptr<FreundlichIsothermConf> params, scalar_t rate):
         m_coeff(params->freund_coeff), m_exp(params->freund_exp),
         RateFunctorBase(rate)
     {}
 
      scalar_t c_eq(index_t node, SingleVariables * const var)  {
         return std::pow(var->bound_concentration(node)/m_coeff, 1/m_exp);
     }
 private:
     scalar_t m_coeff;
     scalar_t m_exp;
 };
 
 class HeavysideRateFunctor
 {
 public:
     HeavysideRateFunctor(scalar_t c_eq, scalar_t b_max, scalar_t rate):
         m_ceq(c_eq), m_bmax(b_max), m_rate(rate)
     {}
 
     scalar_t operator()(index_t node, scalar_t dt, SingleVariables * const var) {
         auto rate =  m_rate * (var->fluid_concentration(node) - m_ceq);
         //if (node == 1) {
         //    std::cout << "rate : " << rate << " - " << c_eq << std::endl;
         //}
         if (rate < 0) {
             rate = 0;
         }
         auto factor = var->porosity(node)*var->saturation(node);
         auto new_c = var->fluid_concentration(node)-rate/factor*dt;
         if (new_c < 0) {
             rate = factor*var->fluid_concentration(node)/dt;
         }
         return rate;
     }
 private:
     scalar_t m_ceq;
     scalar_t m_bmax;
     scalar_t m_rate;
 };
 
 class LeachingRateFunctor:
         public RateFunctorBase
 {
 public:
     LeachingRateFunctor(std::shared_ptr<LeachingIsothermConf> params, scalar_t rate):
         m_coeff1(params->coeff_1), m_coeff2(params->coeff_2), m_exp(params->freund_exp),
         RateFunctorBase(rate)
     {
     }
 
     scalar_t c_eq(index_t node, SingleVariables * const var)  {
        return invert(var->bound_concentration(node), false);
    }
 
     scalar_t func(scalar_t conc) {
         return ph_dependent_binding_law(conc, m_coeff1, m_coeff2, m_exp);
     }
     scalar_t dfunc(scalar_t conc) {
         return dph_dependent_binding_law(conc, m_coeff1, m_coeff2, m_exp);
     }
     scalar_t invert(scalar_t bound_conc, bool print) {
         auto val = simple_newton(std::bind(&LeachingRateFunctor::func, this, std::placeholders::_1),
                              std::bind(&LeachingRateFunctor::dfunc, this, std::placeholders::_1),
                              bound_conc, 1e-7, 0.1, print);
         return val;
     }
 
 private:
     scalar_t m_coeff1;
     scalar_t m_coeff2;
     scalar_t m_exp;
 };
 
 class IsothermRateFunctor:
         public RateFunctorBase
 {
 public:
     IsothermRateFunctor(std::shared_ptr<BindingModelIsothermConf> conf,
                         scalar_t rate):
         m_factor_csh(conf->factor_csh), m_factor_afm(conf->factor_afm),
         m_KCaAd(conf->KCaAd), m_KClAd(conf->KClAd), m_KOHAd(conf->KOHAd),
         m_KClAFm(conf->KClAFm), m_KOHAFm(conf->KOHAFm),
         RateFunctorBase(rate)
     {
     }
 
     scalar_t c_eq(index_t node, SingleVariables * const var)  {
        return invert(var->bound_concentration(node), false);
    }
 
     scalar_t func(scalar_t conc) {
         return isotherm(conc, m_factor_csh, m_factor_afm,
                         m_KCaAd, m_KClAd, m_KOHAd,
                         m_KClAFm, m_KOHAFm
                         );
     }
     scalar_t dfunc(scalar_t conc) {
         using namespace std::placeholders;
         return adiff(std::bind(&IsothermRateFunctor::func, this, _1), conc);
     }
     scalar_t invert(scalar_t bound_conc, bool print) {
         using namespace std::placeholders;
         auto val = simple_newton(std::bind(&IsothermRateFunctor::func, this, std::placeholders::_1),
                              std::bind(&IsothermRateFunctor::dfunc, this, std::placeholders::_1),
                              bound_conc, 1e-7, 0.1, print);
         //if (print) {std::cout << val << std::endl;}
         return val;
     }
 
 private:
     scalar_t m_factor_csh;
     scalar_t m_factor_afm;
     scalar_t m_KCaAd;
     scalar_t m_KClAd;
     scalar_t m_KOHAd;
     scalar_t m_KClAFm;
     scalar_t m_KOHAFm;
 };
 
 class VariableUpscalingStagger:
     public specmicp::reactmicp::systems::single::SingleUpscalingStagger
 {
  public:
     VariableUpscalingStagger(std::function<scalar_t (scalar_t)> diff_of_c):
         m_func(diff_of_c)
     {}
 
     void initialize_timestep(scalar_t dt, VariablesBase * const var) override {
         specmicp::reactmicp::systems::single::SingleUpscalingStagger::initialize_timestep(dt, var);
 
         auto truevar = static_cast<SingleVariables * const>(var);
         for (auto node=1; node<truevar->get_mesh()->nb_nodes(); ++node) {
             truevar->diffusion_coefficient(node) = m_func(truevar->fluid_concentration(node));
         }
     }
 private:
     std::function<scalar_t (scalar_t)> m_func;
 };
 
+class HydratingUpscalingStagger:
+    public specmicp::reactmicp::systems::single::SingleUpscalingStagger
+{
+ public:
+    HydratingUpscalingStagger(scalar_t min_deff, scalar_t hydr_rate):
+        m_min_deff(min_deff), m_hydr_rate(hydr_rate)
+    {
+
+    }
+
+    void initialize(VariablesBase * const var) override {
+        auto truevar = static_cast<SingleVariables * const>(var);
+        m_sound_deff = Vector::Zero(truevar->nb_nodes());
+        for (auto node=1; node<truevar->nb_nodes(); ++node) {
+            m_sound_deff(node) = truevar->diffusion_coefficient(node);
+        }
+
+    }
+
+    void initialize_timestep(scalar_t dt, VariablesBase * const var) override {
+        specmicp::reactmicp::systems::single::SingleUpscalingStagger::initialize_timestep(dt, var);
+
+        auto truevar = static_cast<SingleVariables * const>(var);
+        for (auto node=1; node<truevar->nb_nodes(); ++node) {
+            if (truevar->fluid_concentration(node) < 0.01) {
+                auto val = truevar->diffusion_coefficient(node) - m_hydr_rate*dt;
+                if (val < m_min_deff) {
+                    val = m_min_deff;
+                }
+                truevar->diffusion_coefficient(node) = val;
+                m_sound_deff(node) = val;
+            }
+            //truevar->diffusion_coefficient(node) = m_func(truevar->fluid_concentration(node));
+        }
+    }
+private:
+    Vector m_sound_deff;
+    scalar_t m_min_deff;
+    scalar_t m_hydr_rate;
+};
 
 class LeachingDeff
 {
 public:
     LeachingDeff(scalar_t deff, scalar_t leached_deff):
         sound_deff(deff), leached_deff(leached_deff)
     {}
 
     scalar_t operator() (scalar_t conc) {
         return conc*(leached_deff-sound_deff)/0.5+5e-6;
     }
 
 private:
     scalar_t sound_deff;
     scalar_t leached_deff;
 };
 
 
 void output(scalar_t time, specmicp::reactmicp::solver::VariablesBasePtr variables)
 {
     auto vars = static_cast<specmicp::reactmicp::systems::single::SingleVariables*>(variables.get());
 
     auto mesh = vars->get_mesh();
 
     //Eigen::Matrix<scalar_t, Eigen::Dynamic, 2> toplot(mesh->nb_nodes());
     std::ofstream ofile;
     ofile.open("output/"+std::to_string(time)+".csv");
     for (int i=0;i<mesh->nb_nodes();++i) {
         ofile << mesh->get_position(i) << ", " << vars->fluid_concentration(i) << ", " << vars->bound_concentration(i) << "\n";
     }
     ofile.close();
 }
 
 class OutputPolicy
 {
 public:
     OutputPolicy(std::string folder):
         m_folder(folder)
     {
         // make folder
         int status = mkdir(folder.c_str(),0777);
         if ((status < 0) && (errno != EEXIST)) {
             throw std::runtime_error("Could not create folder: " + folder);
         }
 
     }
 
     void operator()(scalar_t time, specmicp::reactmicp::solver::VariablesBasePtr variables)
     {
         auto vars = static_cast<specmicp::reactmicp::systems::single::SingleVariables*>(variables.get());
         auto mesh = vars->get_mesh();
 
         std::ofstream ofile;
         ofile.open(m_folder+"/"+std::to_string(time)+".csv");
         for (int i=0;i<mesh->nb_nodes();++i) {
             ofile << mesh->get_position(i) << ", " << vars->fluid_concentration(i) << ", " << vars->bound_concentration(i) << "\n";
         }
         ofile.close();
     }
 
 private:
     std::string m_folder;
 };
 
 int main(int argc, char *argv[])
 {
 
     if (argc == 1) {
         std::cerr << "Error: expected one argument" << std::endl;
         return 1;
     }
 
     std::string config_filepath(argv[1]);
     auto config = YAMLConfigFile::load(config_filepath);
 
     specmicp::mesh::Mesh1DPtr the_mesh = specmicp::io::configure_mesh(config.get_section("mesh"));
     auto params = configure_single_simple(config.get_section("parameters"));
 
     std::vector<index_t> list_fixed_nodes = {0,};
 
     Vector init_cond = Vector::Zero(the_mesh->nb_nodes());
 
     auto bcs = config.get_section("boundary_conditions");
     init_cond(0) = bcs.get_required_attribute<scalar_t>("concentration");
     Vector init_cond_solid = Vector::Zero(the_mesh->nb_nodes());
 
 
     auto variables = init_variables(the_mesh, list_fixed_nodes, init_cond, init_cond_solid);
     for (int i=1;i<the_mesh->nb_nodes();++i) {
-        variables->diffusion_coefficient(i) = params->diffusion_coeff;
+        variables->diffusion_coefficient(i) = params->deff_conf->diffusion_coeff;
         variables->porosity(i) = params->porosity;
         variables->saturation(i) = params->saturation;
     }
     variables->diffusion_coefficient(0) = 1;
     variables->porosity(0) = 1;
     variables->saturation(0) = 1;
 
 
     // transport stagger
     auto transport_stagger =
             std::make_shared<SingleTransportStagger>(variables, list_fixed_nodes);
     auto& options = transport_stagger->options_solver();
     configure_transport_options(options, config.get_section("transport_options"));
     options.alpha = 1.0;
     //options.residuals_tolerance = 1e-10;
     //options.quasi_newton = 3;
 
     // chemistry stagger
     std::shared_ptr<ChemistryProgram> chem_program = nullptr;
     if (params->binding_law == "Freundlich") {
         auto functor = RateFunctor(std::static_pointer_cast<FreundlichIsothermConf>(params->isotherm_conf), params->rate_factor);
         chem_program = std::make_shared<DefaultChemistryProgram>(functor);
     } else if (params->binding_law == "LeachingFreundlich") {
         auto functor = LeachingRateFunctor(std::static_pointer_cast<LeachingIsothermConf>(params->isotherm_conf), params->rate_factor);
         chem_program = std::make_shared<DefaultChemistryProgram>(functor);
     } else if (params->binding_law == "BindingModelIsotherm") {
         auto functor = IsothermRateFunctor(std::static_pointer_cast<BindingModelIsothermConf>(params->isotherm_conf), params->rate_factor);
         chem_program = std::make_shared<DefaultChemistryProgram>(functor);
     }
 
     auto chemistry_stagger = std::make_shared<ChemistryStagger>(chem_program);
 
 
     // upscaling stagger
-    //auto upscaling_stagger = std::make_shared<SingleUpscalingStagger>();
     std::shared_ptr<SingleUpscalingStagger> upscaling_stagger = nullptr;
-    if (params->leached_diffusion_coeff == 0) {
+    if (params->diffusion_model == "Fixed") {
         upscaling_stagger = std::make_shared<SingleUpscalingStagger>();
-    } else {
-        auto ups_functor = LeachingDeff(params->diffusion_coeff, params->leached_diffusion_coeff);
+    } else if (params->diffusion_model == "Leaching") {
+        auto deff_conf = std::static_pointer_cast<LeachedDeffConf>(params->deff_conf);
+        auto ups_functor = LeachingDeff(deff_conf->diffusion_coeff, deff_conf->leached_diffusion_coeff);
         upscaling_stagger = std::make_shared<VariableUpscalingStagger>(ups_functor);
-
+    } else if (params->diffusion_model == "Hydrating") {
+        auto deff_conf = std::static_pointer_cast<HydratingDeffConf>(params->deff_conf);
+        upscaling_stagger = std::make_shared<HydratingUpscalingStagger>(deff_conf->min_diffusion_coefficient,
+                                                                        deff_conf->hydration_rate);
     }
 
-
+    upscaling_stagger->initialize(variables.get());
     // reactmicp solver
     auto solver = specmicp::reactmicp::solver::ReactiveTransportSolver(transport_stagger, chemistry_stagger, upscaling_stagger);
     auto& opts = solver.get_options();
     specmicp::io::configure_reactmicp_options(opts, config.get_section("reactmicp_options"));
     //opts.maximum_iterations = 51;
     //opts.residuals_tolerance = 1e-5;
     //opts.
 
     // runner
 
     auto conf_runner = config.get_section("runner");
     auto conf_dt = config.get_section("timestepper");
 
 
     auto info = specmicp::reactmicp::solver::SimulationInformation("ex_reactmicp_single", conf_runner.get_required_attribute<scalar_t>("save_time"));
     auto runner = std::make_shared<specmicp::reactmicp::solver::ReactiveTransportRunner>(
                 solver,
                 conf_dt.get_required_attribute<scalar_t>("minimum_dt"),
                 conf_dt.get_required_attribute<scalar_t>("maximum_dt"),
                 info);
-    auto topts = runner->get_timestepper_options();
-    topts.iteration_lower_target = conf_dt.get_optional_attribute<scalar_t>("minimum_dt", 2.1);
-    topts.iteration_upper_target = conf_dt.get_optional_attribute<scalar_t>("maximum_dt", 10.1);
-    topts.increase_factor = conf_dt.get_optional_attribute<scalar_t>("increase_factor", 0.01);
+    auto& topts = runner->get_timestepper_options();
+    //topts.iteration_lower_target = conf_dt.get_optional_attribute<scalar_t>("minimum_dt", 2.1);
+    //topts.iteration_upper_target = conf_dt.get_optional_attribute<scalar_t>("maximum_dt", 10.1);
+    topts.increase_factor = conf_dt.get_optional_attribute<scalar_t>("increase_factor", 1.1);
+    topts.increase_error_minimization = conf_dt.get_optional_attribute<scalar_t>("increase_factor", 1.1);
     auto output_policy = OutputPolicy(conf_runner.get_required_attribute<std::string>("output_folder"));
     runner->set_output_policy(output_policy);
     runner->run_until(conf_runner.get_required_attribute<scalar_t>("run_time"), variables);
 
     for (int i=0; i<the_mesh->nb_nodes(); ++i) {
         std::cout << variables->fluid_concentration(i) << "\n";
     }
     std::cout << "----------\n" << std::endl;
     for (int i=0; i<the_mesh->nb_nodes(); ++i) {
         std::cout << variables->bound_concentration(i) << "\n";
     }
     std::cout << std::endl;
     return 0;
 }