diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index a1619b9..a71c349 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -1,142 +1,162 @@
 ###### Installation #######################
 
 set(EXAMPLES_INSTALL_DIR ${SHARE_INSTALL_DIR}/examples)
 
 # SpecMiCP
 # ========
 
 set (SPECMICP_EXAMPLES
     specmicp/adimensional/thermocarbo.cpp
     specmicp/adimensional/equilibrium_curve.cpp
     specmicp/adimensional/momas_thermo.cpp
     specmicp/adimensional/carbofe.cpp
 )
 install(FILES ${SPECMICP_EXAMPLES}
     DESTINATION ${EXAMPLES_INSTALL_DIR}/specmicp
 )
 
 # ReactMiCP
 # =========
 
 set (REACTMICP_SATURATED_EXAMPLES
     reactmicp/equilibrium_curve.cpp
+
     reactmicp/saturated_react/react_leaching.cpp
     reactmicp/saturated_react/carbonation.cpp
     reactmicp/saturated_react/carbonation.yaml # the configuration file
     reactmicp/saturated_react/carbonationfe.cpp
     reactmicp/saturated_react/leaching_kinetic.cpp
+
+    reactmicp/unsaturated/drying.cpp
 )
 
 install(FILES ${REACTMICP_SATURATED_EXAMPLES}
     DESTINATION ${EXAMPLES_INSTALL_DIR}/reactmicp/saturated
 )
 
 file(COPY reactmicp/saturated_react/carbonation.yaml
      DESTINATION ${CMAKE_CURRENT_BINARY_DIR}
 )
 
 ###### Build (optional) ####################
 
 if (NOT SPECMICP_BUILD_EXAMPLE)
     set(EXAMPLE_IS_OPTIONAL EXCLUDE_FROM_ALL)
 else()
     set(EXAMPLE_IS_OPTIONAL )
 endif()
 
 # SpecMiCP
 # ========
 
 #            thermocarbo
 
 add_executable(ex_adim_thermocarbo
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/thermocarbo.cpp
 )
 target_link_libraries(ex_adim_thermocarbo ${SPECMICP_LIBS})
 
 
 #            carbofe
 
 add_executable(ex_adim_carbofe
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/carbofe.cpp
 )
 target_link_libraries(ex_adim_carbofe ${SPECMICP_LIBS})
 
 #           Momas thermodynamic example
 
 add_executable(ex_adim_momas
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/momas_thermo.cpp
 )
 target_link_libraries(ex_adim_momas ${SPECMICP_LIBS})
 
 
 #            equilibrium curve
 
 add_executable(ex_adim_equilibriumcurve
     ${EXAMPLE_IS_OPTIONAL}
     specmicp/adimensional/equilibrium_curve.cpp
 )
 target_link_libraries(ex_adim_equilibriumcurve ${SPECMICP_LIBS})
 
 
 
 # ReactMiCP
 # =========
 
 # EquilibriumCurve
 
 add_executable(ex_reactmicp_equilibriumcurve
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/equilibrium_curve.cpp
 )
 target_link_libraries(ex_reactmicp_equilibriumcurve ${REACTMICP_LIBS})
 
 # Leaching
 
 add_executable(ex_reactmicp_leaching
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/saturated_react/react_leaching.cpp
 )
 target_link_libraries(ex_reactmicp_leaching ${REACTMICP_LIBS})
 
 # Momas benchmark
 
 add_executable(ex_reactmicp_momas
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/saturated_react/momas_benchmark.cpp
 )
 target_link_libraries(ex_reactmicp_momas ${REACTMICP_LIBS})
 
 # Carbonation
 
 add_executable(ex_reactmicp_carbo
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/saturated_react/carbonation.cpp
 )
 target_link_libraries(ex_reactmicp_carbo ${REACTMICP_LIBS})
 
 add_executable(ex_reactmicp_carbo_static
     EXCLUDE_FROM_ALL
     reactmicp/saturated_react/carbonation.cpp
 )
 set_source_files_properties(reactmicp/saturated_react/carbonation.cpp PROPERTIES COMPILE_FLAGS -flto)
 target_link_libraries(ex_reactmicp_carbo_static ${REACTMICP_STATIC_LIBS})
 
 
 # Carbonation with Fe
 
 add_executable(ex_reactmicp_carbofe
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/saturated_react/carbonationfe.cpp
 )
 target_link_libraries(ex_reactmicp_carbofe ${REACTMICP_LIBS})
 
 # Leaching with hydration
 
 add_executable(ex_reactmicp_leaching_kinetic
     ${EXAMPLE_IS_OPTIONAL}
     reactmicp/saturated_react/leaching_kinetic.cpp
 )
 target_link_libraries(ex_reactmicp_leaching_kinetic ${REACTMICP_STATIC_LIBS})
+
+
+
+# Unsaturated system
+# ==================
+
+add_executable(ex_reactmicp_drying
+    ${EXAMPLE_IS_OPTIONAL}
+    reactmicp/unsaturated/drying.cpp
+)
+
+set(CEMDATA_PATH \"../data/cemdata.yaml\")
+set_source_files_properties(reactmicp/unsaturated/drying.cpp PROPERTIES COMPILE_DEFINITIONS
+"CEMDATA_PATH=${CEMDATA_PATH}"
+)
+
+target_link_libraries(ex_reactmicp_drying ${REACTMICP_STATIC_LIBS})
diff --git a/examples/reactmicp/unsaturated/drying.cpp b/examples/reactmicp/unsaturated/drying.cpp
new file mode 100644
index 0000000..4f0adf3
--- /dev/null
+++ b/examples/reactmicp/unsaturated/drying.cpp
@@ -0,0 +1,564 @@
+/*-------------------------------------------------------------------------------
+
+Copyright (c) 2015 F. Georget <fabieng@princeton.edu>, Princeton University
+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.
+
+-----------------------------------------------------------------------------*/
+
+
+// ====================================================== //
+//                                                        //
+//                     Drying example                     //
+//                     ==============                     //
+//                                                        //
+// ====================================================== //
+//
+//
+// This is a simple example to show how the unsaturated
+// system of ReactMiCP is able to solve a complex drying
+// equation.
+//
+// It implements a special chemistry stagger (DryingStagger)
+// which only computes the vapor pressure.
+
+
+// Macro to define the problem
+// ===========================
+
+#define NB_NODES 50
+#define DX 0.05*1e-2                 // cm
+#define CROSS_SECTION 5.0*1e-4      // cm^2
+
+#define GAMMA_WATER_GLASS 0.1  // N/m
+#define GAMMA_WATER_AIR 71.97e-3  // N/m = kg/s^2
+#define R_SMALL_BEAD 0.015*1e-2     // cm
+#define R_BIG_BEAD 0.04*1e-2        // cm
+
+#define M_v 18.015e-3          // kg/mol
+
+// Includes
+// ========
+
+#include "database/database.hpp"
+
+#include "reactmicp/systems/unsaturated/variables.hpp"
+#include "reactmicp/systems/unsaturated/variables_sub.hpp"
+#include "reactmicp/systems/unsaturated/transport_stagger.hpp"
+#include "reactmicp/systems/unsaturated/transport_constraints.hpp"
+
+#include "reactmicp/solver/staggers_base/chemistry_stagger_base.hpp"
+#include "reactmicp/solver/staggers_base/upscaling_stagger_base.hpp"
+#include "reactmicp/solver/staggers_base/stagger_structs.hpp"
+
+#include "reactmicp/solver/runner.hpp"
+
+#include "dfpm/meshes/generic_mesh1d.hpp"
+
+#include "physics/laws.hpp"
+#include "physics/units.hpp"
+
+#include "utils/log.hpp"
+#include "utils/io/csv_formatter.hpp"
+
+#include <functional>
+#include <iostream>
+
+// Namespace
+// =========
+using namespace specmicp;
+using namespace specmicp::reactmicp;
+using namespace specmicp::reactmicp::systems;
+using namespace specmicp::reactmicp::systems::unsaturated;
+
+// Declarations
+// ============
+
+// Constants
+// ==========
+static const scalar_t rho_l_si = constants::water_density_25;
+static const scalar_t rho_l = 1e-6*constants::water_density_25;
+static const scalar_t cw_tilde_si = rho_l_si / M_v;
+static const scalar_t cw_tilde = 1e-6*cw_tilde_si;
+
+int main();
+
+database::RawDatabasePtr get_database();
+mesh::Mesh1DPtr get_mesh();
+
+scalar_t leverett_function(scalar_t s_eff);
+
+scalar_t big_bead_cap_pressure(scalar_t saturation);
+scalar_t small_bead_cap_pressure(scalar_t saturation);
+scalar_t vapor_pressure(scalar_t pc);
+
+// Drying Stagger
+// ==============
+//
+// This stagger computes the exchange between water vapor and liquid water
+class DryingStagger: public solver::ChemistryStaggerBase
+{
+public:
+    DryingStagger() {}
+
+    void initialize(solver::VariablesBasePtr var) {}
+
+     //! \brief Initialize the stagger at the beginning of an iteration
+    //!
+    //! This is where the predictor can be saved, the first trivial iteration done, ...
+    //!
+    //! \param dt the duration of the timestep
+    //! \param var a shared_ptr to the variables
+    void initialize_timestep(scalar_t dt, solver::VariablesBasePtr var) {
+        m_dt = dt;
+    }
+
+    //! \brief Solve the equation for the timestep
+    //!
+    //! \param var a shared_ptr to the variables
+    solver::StaggerReturnCode restart_timestep(solver::VariablesBasePtr var);
+
+    void compute_one_node(index_t node, UnsaturatedVariables* vars);
+
+private:
+    scalar_t m_dt {-1.0};
+};
+
+// Upscaling Stagger
+// =================
+//
+// Define how we compute capillary pressure, vapor pressure, ...
+class UpscalingDryingStagger: public solver::UpscalingStaggerBase
+{
+public:
+    UpscalingDryingStagger(std::vector<bool>& is_big_bead):
+        m_is_big_bead(is_big_bead)
+    {}
+
+    //! \brief Initialize the stagger at the beginning of the computation
+    //!
+    //! \param var a shared_ptr to the variables
+    void initialize(solver::VariablesBasePtr var);
+
+    //! \brief Initialize the stagger at the beginning of an iteration
+    //!
+    //! This is where the predictor can be saved, the first trivial iteration done, ...
+    //!
+    //! \param dt the duration of the timestep
+    //! \param var a shared_ptr to the variables
+    void initialize_timestep(scalar_t dt, solver::VariablesBasePtr var) {}
+
+    //! \brief Solve the equation for the timestep
+    //!
+    //! \param var a shared_ptr to the variables
+    solver::StaggerReturnCode restart_timestep(solver::VariablesBasePtr var) {
+        return solver::StaggerReturnCode::ResidualMinimized;
+    }
+
+
+    scalar_t capillary_pressure_model(index_t node, scalar_t saturation);
+
+    scalar_t vapor_pressure_model(index_t node, scalar_t saturation);
+
+    scalar_t relative_liquid_permeability_model(index_t node, scalar_t saturation);
+
+    scalar_t relative_liquid_diffusivity_model(index_t node, scalar_t saturation);
+
+    scalar_t relative_gas_diffusivity_model(index_t node, scalar_t saturation);
+
+private:
+    std::vector<bool> m_is_big_bead;
+};
+
+
+
+// Implementation
+// ==============
+
+scalar_t UpscalingDryingStagger::capillary_pressure_model(index_t node, scalar_t saturation)
+{
+    scalar_t pc = 0.0;
+    if (m_is_big_bead[node])
+    {
+        pc = big_bead_cap_pressure(saturation);
+    }
+    else
+        pc = small_bead_cap_pressure(saturation);
+    return pc;
+}
+
+scalar_t UpscalingDryingStagger::vapor_pressure_model(index_t node, scalar_t saturation)
+{
+    return vapor_pressure(capillary_pressure_model(node, saturation));
+}
+
+
+scalar_t UpscalingDryingStagger::relative_liquid_permeability_model(index_t node, scalar_t saturation)
+{
+    return std::pow(saturation, 3);
+}
+
+scalar_t UpscalingDryingStagger::relative_liquid_diffusivity_model(index_t node, scalar_t saturation)
+{
+    return saturation;
+}
+
+scalar_t UpscalingDryingStagger::relative_gas_diffusivity_model(index_t node, scalar_t saturation)
+{
+    return (1.0 - saturation);
+}
+
+void UpscalingDryingStagger::initialize(solver::VariablesBasePtr var)
+{
+    UnsaturatedVariables* true_vars = static_cast<UnsaturatedVariables*>(var.get());
+
+    true_vars->set_capillary_pressure_model(std::bind(std::mem_fn(
+                            &UpscalingDryingStagger::capillary_pressure_model
+                            ), this, std::placeholders::_1, std::placeholders::_2));
+
+    true_vars->set_vapor_pressure_model(std::bind(std::mem_fn(
+                            &UpscalingDryingStagger::vapor_pressure_model
+                            ), this, std::placeholders::_1, std::placeholders::_2));
+
+    true_vars->set_relative_liquid_permeability_model(std::bind(std::mem_fn(
+                            &UpscalingDryingStagger::relative_liquid_permeability_model
+                            ), this, std::placeholders::_1, std::placeholders::_2));
+
+    true_vars->set_relative_liquid_diffusivity_model(std::bind(std::mem_fn(
+                            &UpscalingDryingStagger::relative_liquid_diffusivity_model
+                            ), this, std::placeholders::_1, std::placeholders::_2));
+
+    true_vars->set_relative_gas_diffusivity_model(std::bind(std::mem_fn(
+                            &UpscalingDryingStagger::relative_gas_diffusivity_model
+                            ), this, std::placeholders::_1, std::placeholders::_2));
+
+    MainVariable& vapor_pressure = true_vars->get_pressure_main_variables(0);
+    MainVariable& saturation = true_vars->get_liquid_saturation();
+
+    SecondaryTransientVariable& porosity   = true_vars->get_porosity();
+    SecondaryVariable& liquid_diffusivity  = true_vars->get_liquid_diffusivity();
+    SecondaryVariable& liquid_permeability = true_vars->get_liquid_permeability();
+    SecondaryVariable& gas_diffusivity     = true_vars->get_gas_diffusivity();
+
+
+    porosity(0) = 1;
+    gas_diffusivity(0) = 0.282*1e-4;
+
+
+    for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node)
+    {
+        vapor_pressure(node) = vapor_pressure_model(node, saturation(node));
+
+        if (m_is_big_bead[node]) {
+            porosity(node) = 0.371;
+            liquid_diffusivity(node) = 0;
+            liquid_permeability(node) = 3.52e-11;
+            gas_diffusivity(node) = 2*0.371/(3-0.371)*0.282*1e-4;
+        } else {
+            porosity(node) = 0.387;
+            liquid_diffusivity(node) = 0;
+            liquid_permeability(node) = 8.41e-12;
+            gas_diffusivity(node) = 2*0.387/(3-0.387)*0.282*1e-4;
+        }
+    }
+
+    porosity.predictor = porosity.variable;
+
+    true_vars->set_relative_variables();
+}
+
+
+
+solver::StaggerReturnCode DryingStagger::restart_timestep(solver::VariablesBasePtr var)
+{
+    UnsaturatedVariables* true_vars = static_cast<UnsaturatedVariables*>(var.get());
+
+    for (index_t node=1; node<true_vars->get_mesh()->nb_nodes(); ++node) {
+        compute_one_node(node, true_vars);
+    }
+
+    return solver::StaggerReturnCode::ResidualMinimized;
+}
+
+void DryingStagger::compute_one_node(index_t node, UnsaturatedVariables *vars)
+{
+    MainVariable& satvars = vars->get_liquid_saturation();
+    MainVariable& presvars = vars->get_pressure_main_variables(0);
+
+
+    scalar_t rt = vars->get_rt();
+    scalar_t saturation = satvars(node);
+    scalar_t pressure = presvars(node);
+    scalar_t porosity = vars->get_porosity()(node);
+    auto vapor_pressure_f = vars->get_vapor_pressure_model();
+
+    const scalar_t tot_conc = cw_tilde_si*saturation + (1.0-saturation)*pressure/rt;
+
+    pressure = vapor_pressure_f(node, saturation);
+    scalar_t cw_hat = pressure/rt;
+    saturation = (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);
+
+    pressure = vapor_pressure_f(node, saturation);
+    cw_hat = pressure/rt;
+
+    scalar_t res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat;
+
+    while (std::abs(res)/tot_conc > 1e-12) {
+        saturation =  (tot_conc - cw_hat) / (cw_tilde_si-cw_hat);
+        pressure = vapor_pressure_f(node, saturation);
+        cw_hat = pressure/rt;
+
+        res = tot_conc - cw_tilde_si*saturation - (1.0-saturation)*cw_hat;
+    }
+
+    satvars(node) = saturation;
+    satvars.velocity(node) = (saturation - satvars.predictor(node))/m_dt;
+
+    presvars(node) = pressure;
+    presvars.velocity(node) = (pressure - presvars.predictor(node))/m_dt;
+    presvars.chemistry_rate(node) =  - porosity/(rt*m_dt) * (
+                (pressure*(1.0-saturation)) -
+                (presvars.predictor(node)*(1.0-satvars.predictor(node)))
+                )
+            + presvars.transport_fluxes(node);
+}
+
+
+database::RawDatabasePtr get_database()
+{
+
+    specmicp::database::Database thedatabase(CEMDATA_PATH);
+    std::vector<std::string> to_keep = {"H2O", "H[+]"};
+    thedatabase.keep_only_components(to_keep);
+    thedatabase.remove_gas_phases();
+    database::RawDatabasePtr raw_data = thedatabase.get_database();
+    raw_data->freeze_db();
+    return raw_data;
+}
+
+mesh::Mesh1DPtr get_mesh()
+{
+    mesh::Uniform1DMeshGeometry geom;
+    geom.dx = DX;
+    geom.nb_nodes = NB_NODES;
+    geom.section = CROSS_SECTION;
+    return mesh::uniform_mesh1d(geom);
+}
+
+scalar_t big_bead_cap_pressure(scalar_t saturation) {
+    return GAMMA_WATER_AIR/std::sqrt(3.52e-11/0.371)*leverett_function(saturation);
+}
+
+scalar_t small_bead_cap_pressure(scalar_t saturation) {
+    return GAMMA_WATER_AIR/std::sqrt(8.41e-12/0.387)*leverett_function(saturation);
+}
+
+scalar_t vapor_pressure(scalar_t pc) {
+    return 3e3*std::exp(-M_v * pc / (rho_l_si * (8.314*(25+273.15))));
+}
+
+
+scalar_t leverett_function(scalar_t s_eff) {
+    return 0.325*std::pow((1.0/s_eff - 1), 0.217);
+}
+
+
+class Printer
+{
+public:
+    Printer(UnsaturatedVariables* vars, std::size_t reserve_size);
+    
+    
+    void register_timestep(scalar_t time, solver::VariablesBasePtr base_vars);
+    
+    void print(const std::string& name);
+    
+private:
+    mesh::Mesh1DPtr m_mesh;
+
+    std::vector<std::vector<scalar_t>> m_pressure;
+    std::vector<std::vector<scalar_t>> m_saturation;
+    std::vector<scalar_t> m_timestep;
+};
+
+Printer::Printer(UnsaturatedVariables *vars, std::size_t reserve_size):
+    m_mesh(vars->get_mesh()),
+    m_pressure(m_mesh->nb_nodes()),
+    m_saturation(m_mesh->nb_nodes())
+{
+    m_timestep.reserve(reserve_size);
+    for (auto node: m_mesh->range_nodes())
+    {
+        m_pressure[node].reserve(reserve_size);
+        m_saturation[node].reserve(reserve_size);
+    }
+}
+
+void Printer::register_timestep(scalar_t time, solver::VariablesBasePtr base_vars)
+{
+    UnsaturatedVariables* vars = static_cast<UnsaturatedVariables*>(base_vars.get());
+
+    m_timestep.push_back(time);
+    MainVariable& saturation = vars->get_liquid_saturation();
+    MainVariable& pressure = vars->get_pressure_main_variables(0);
+
+    for (auto node: m_mesh->range_nodes())
+    {
+        m_pressure[node].push_back(pressure(node));
+        m_saturation[node].push_back(saturation(node));
+    }
+}
+
+void Printer::print(const std::string& name)
+{
+
+    io::CSVFile pressure_file(name+"_pressure.dat");
+    io::CSVFile saturation_file(name+"_saturation.dat");
+
+    pressure_file << "Node";
+    saturation_file << "Node";
+    for (auto& time: m_timestep)
+    {
+        pressure_file.separator();
+        pressure_file << time;
+        saturation_file.separator();
+        saturation_file << time;
+    }
+    pressure_file.eol();
+    saturation_file.eol();
+
+    for (auto node: m_mesh->range_nodes())
+    {
+        pressure_file << m_mesh->get_position(node);
+        saturation_file << m_mesh->get_position(node);
+        for (std::size_t ind=0; ind<m_timestep.size(); ++ind)
+        {
+            pressure_file.separator();
+            pressure_file << m_pressure[node][ind];
+            saturation_file.separator();
+            saturation_file << m_saturation[node][ind];
+        }
+        pressure_file.eol();
+        saturation_file.eol();
+    }
+
+    pressure_file.flush();
+    saturation_file.flush();
+}
+
+
+
+
+int main()
+{
+    init_logger(&std::cerr, logger::Warning);
+
+    mesh::Mesh1DPtr the_mesh = get_mesh();
+    database::RawDatabasePtr raw_data = get_database();
+    std::vector<bool> has_gas = {true, false, false};
+
+    std::vector<bool> is_big_bead(the_mesh->nb_nodes(), true);
+    //for (int node=25; node<NB_NODES; ++node)
+    //    is_big_bead[node] = true;
+
+    UnsaturatedVariablesPtr vars = std::make_shared<UnsaturatedVariables>(
+                the_mesh, raw_data, has_gas, units::LengthUnit::meter);
+
+
+    TransportConstraints constraints;
+    constraints.add_fixed_node(0);
+    constraints.add_gas_node(0);
+
+
+
+    // variables initialisation
+    MainVariable& saturation = vars->get_liquid_saturation();
+    MainVariable& vapor_pressure = vars->get_pressure_main_variables(0);
+    SecondaryTransientVariable& water_aqueous_concentration = vars->get_water_aqueous_concentration();
+    water_aqueous_concentration.predictor = water_aqueous_concentration.variable;
+
+    SecondaryVariable& rel_gas_d = vars->get_relative_gas_diffusivity();
+    rel_gas_d(0) = 1.0;
+
+    vapor_pressure(0) = 100;
+
+    for (index_t node = 1; node < the_mesh->nb_nodes(); ++node) {
+        saturation(node) = 0.9;
+    }
+
+    water_aqueous_concentration.set_constant(cw_tilde_si);
+
+
+
+    std::shared_ptr<UpscalingDryingStagger> upscaling_stagger = std::make_shared<UpscalingDryingStagger>(is_big_bead);
+    upscaling_stagger->initialize(vars);
+
+    std::shared_ptr<DryingStagger> drying_stagger = std::make_shared<DryingStagger>();
+    // init vapor pressure
+    drying_stagger->initialize_timestep(1.0, vars);
+    drying_stagger->restart_timestep(vars);
+
+
+    std::shared_ptr<UnsaturatedTransportStagger> transport_stagger =
+            std::make_shared<UnsaturatedTransportStagger>(vars, constraints);
+
+    solver::ReactiveTransportSolver react_solver(transport_stagger, drying_stagger, upscaling_stagger);
+    solver::ReactiveTransportOptions& opts = react_solver.get_options();
+    opts.residuals_tolerance = 1e-2;
+    opts.good_enough_tolerance = 0.99;
+
+    vapor_pressure.chemistry_rate.setZero();
+
+    //for (auto it=0; it<100; ++it)
+    //{
+//    auto retcode = react_solver.solve_timestep(0.1, vars);
+//    std::cout << (int) retcode << std::endl;
+    //}
+
+    //std::cout << saturation.variable << std::endl;
+    //std::cout << vapor_pressure.variable << std::endl;
+    //std::cout << vars->get_gas_diffusivity().variable << std::endl;
+    //std::cout << vars->get_relative_gas_diffusivity().variable << std::endl;
+
+    int nb_hours = 9;
+    Printer printer(vars.get(), 4*nb_hours);
+    solver::SimulationInformation simul_info("drying", 900);
+    solver::ReactiveTransportRunner runner(react_solver, 0.01, 10, simul_info);
+
+    solver::output_f out_func = std::bind(&Printer::register_timestep,
+                             &printer,
+                             std::placeholders::_1, std::placeholders::_2);
+
+    runner.set_output_policy(out_func);
+
+    runner.run_until(nb_hours*3600, vars);
+
+    printer.print("out_drying");
+
+    std::cout << saturation.variable << std::endl;
+    std::cout << vapor_pressure.variable << std::endl;
+}