diff --git a/opm/simulators/SimulatorIncompTwophase.cpp b/opm/simulators/SimulatorIncompTwophase.cpp
index f51977891..b282b335e 100644
--- a/opm/simulators/SimulatorIncompTwophase.cpp
+++ b/opm/simulators/SimulatorIncompTwophase.cpp
@@ -42,10 +42,9 @@
 #include <opm/core/props/IncompPropertiesInterface.hpp>
 #include <opm/core/props/rock/RockCompressibility.hpp>
 
-#include <opm/core/utility/ColumnExtract.hpp>
 #include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/simulator/WellState.hpp>
-//#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
+#include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
 #include <opm/core/transport/TransportSolverTwophaseImplicit.hpp>
 #include <boost/filesystem.hpp>
 #include <boost/scoped_ptr.hpp>
@@ -87,6 +86,7 @@ namespace Opm
         int max_well_control_iterations_;
         // Parameters for transport solver.
         int num_transport_substeps_;
+        bool use_reorder_;
         bool use_segregation_split_;
         // Observed objects.
         const UnstructuredGrid& grid_;
@@ -98,11 +98,7 @@ namespace Opm
         const FlowBoundaryConditions* bcs_;
         // Solvers
         IncompTpfa psolver_;
-        // this should maybe be packed in a separate file
         boost::scoped_ptr<TransportSolverTwophaseInterface> tsolver_;
-        //ImpliciteTwoPhaseTransportSolver tsolver_;
-        // Needed by column-based gravity segregation solver.
-        std::vector< std::vector<int> > columns_;
         // Misc. data
         std::vector<int> allcells_;
     };
@@ -320,7 +316,9 @@ namespace Opm
                                         const FlowBoundaryConditions* bcs,
                                         LinearSolverInterface& linsolver,
                                         const double* gravity)
-        : grid_(grid),
+        : use_reorder_(param.getDefault("use_reorder", true)),
+          use_segregation_split_(param.getDefault("use_segregation_split", false)),
+          grid_(grid),
           props_(props),
           rock_comp_props_(rock_comp_props),
           wells_manager_(wells_manager),
@@ -333,28 +331,29 @@ namespace Opm
                    param.getDefault("nl_pressure_maxiter", 10),
                    gravity, wells_manager.c_wells(), src, bcs)
     {
-        const bool use_reorder = param.getDefault("use_reorder", true);
-        if (use_reorder) {
-            /*
-            tsolver_.reset(new Opm::TransportModelTwoPhase(grid, props, rock_comp_props, linsolver,
-                                                           param.getDefault("nl_pressure_residual_tolerance", 0.0),
-                                                           param.getDefault("nl_pressure_change_tolerance", 1.0),
-                                                           param.getDefault("nl_pressure_maxiter", 10),
-                                                           gravity, wells_manager.c_wells(), src, bcs));
-                                                           */
+        // Initialize transport solver.
+        if (use_reorder_) {
+            tsolver_.reset(new Opm::TransportSolverTwophaseReorder(grid,
+                                                                   props,
+                                                                   use_segregation_split_ ? gravity : NULL,
+                                                                   param.getDefault("nl_tolerance", 1e-9),
+                                                                   param.getDefault("nl_maxiter", 30)));
+
         } else {
             if (rock_comp_props->isActive()) {
                 THROW("The implicit pressure solver cannot handle rock compressibility.");
             }
+            if (use_segregation_split_) {
+                THROW("The implicit pressure solver is not set up to use segregation splitting.");
+            }
             std::vector<double> porevol;
             computePorevolume(grid, props.porosity(), porevol);
-            tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(
-                               grid,
-                               props,
-                               porevol,
-                               gravity,
-                               psolver_.getHalfTrans(),
-                               param));
+            tsolver_.reset(new Opm::TransportSolverTwophaseImplicit(grid,
+                                                                    props,
+                                                                    porevol,
+                                                                    gravity,
+                                                                    psolver_.getHalfTrans(),
+                                                                    param));
         }
 
         // For output.
@@ -379,7 +378,7 @@ namespace Opm
 
         // Transport related init.
         num_transport_substeps_ = param.getDefault("num_transport_substeps", 1);
-        use_segregation_split_ = param.getDefault("use_segregation_split", false);
+
         // Misc init.
         const int num_cells = grid.number_of_cells;
         allcells_.resize(num_cells);
@@ -445,6 +444,12 @@ namespace Opm
                     outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
                 }
                 outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+                if (use_reorder_) {
+                    // This use of dynamic_cast is not ideal, but should be safe.
+                    outputVectorMatlab(std::string("reorder_it"),
+                                       dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
+                                       timer.currentStepNum(), output_dir_);
+                }
             }
 
             SimulatorReport sreport;
@@ -538,10 +543,7 @@ namespace Opm
             double injected[2] = { 0.0 };
             double produced[2] = { 0.0 };
             for (int tr_substep = 0; tr_substep < num_transport_substeps_; ++tr_substep) {
-                tsolver_->solve(&transport_src[0],
-                                &porevol[0],
-                                stepsize,
-                                state);
+                tsolver_->solve(&transport_src[0], &porevol[0], stepsize, state);
 
                 double substep_injected[2] = { 0.0 };
                 double substep_produced[2] = { 0.0 };
@@ -551,18 +553,20 @@ namespace Opm
                 injected[1] += substep_injected[1];
                 produced[0] += substep_produced[0];
                 produced[1] += substep_produced[1];
-
-
+                if (use_reorder_ && use_segregation_split_) {
+                    // Again, unfortunate but safe use of dynamic_cast.
+                    // Possible solution: refactor gravity solver to its own class.
+                    dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
+                        .solveGravity(&initial_porevol[0], stepsize, state);
+                }
                 watercut.push(timer.currentTime() + timer.currentStepLength(),
                               produced[0]/(produced[0] + produced[1]),
                               tot_produced[0]/tot_porevol_init);
-
                 if (wells_) {
                     wellreport.push(props_, *wells_, state.saturation(),
                                     timer.currentTime() + timer.currentStepLength(),
                                     well_state.bhp(), well_state.perfRates());
                 }
-
             }
             transport_timer.stop();
             double tt = transport_timer.secsSinceStart();
@@ -590,6 +594,12 @@ namespace Opm
                 outputStateVtk(grid_, state, timer.currentStepNum(), output_dir_);
             }
             outputStateMatlab(grid_, state, timer.currentStepNum(), output_dir_);
+            if (use_reorder_) {
+                // This use of dynamic_cast is not ideal, but should be safe.
+                outputVectorMatlab(std::string("reorder_it"),
+                                   dynamic_cast<const TransportSolverTwophaseReorder&>(*tsolver_).getReorderIterations(),
+                                   timer.currentStepNum(), output_dir_);
+                }
             outputWaterCut(watercut, output_dir_);
             if (wells_) {
                 outputWellReport(wellreport, output_dir_);