diff --git a/opm/autodiff/BlackoilModel.hpp b/opm/autodiff/BlackoilModel.hpp
index c5006b2174..d4473fed6f 100644
--- a/opm/autodiff/BlackoilModel.hpp
+++ b/opm/autodiff/BlackoilModel.hpp
@@ -26,6 +26,7 @@
 #include <opm/autodiff/BlackoilModelBase.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/autodiff/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/autodiff/StandardWells.hpp>
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 
 namespace Opm {
@@ -40,10 +41,10 @@ namespace Opm {
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
     template<class Grid>
-    class BlackoilModel : public BlackoilModelBase<Grid, BlackoilModel<Grid> >
+    class BlackoilModel : public BlackoilModelBase<Grid, StandardWells, BlackoilModel<Grid> >
     {
     public:
-        typedef BlackoilModelBase<Grid, BlackoilModel<Grid> > Base;
+        typedef BlackoilModelBase<Grid, StandardWells, BlackoilModel<Grid> > Base;
 
         /// Construct the model. It will retain references to the
         /// arguments of this functions, and they are expected to
diff --git a/opm/autodiff/BlackoilModelBase.hpp b/opm/autodiff/BlackoilModelBase.hpp
index 9a59a7f640..d39e4e7196 100644
--- a/opm/autodiff/BlackoilModelBase.hpp
+++ b/opm/autodiff/BlackoilModelBase.hpp
@@ -32,7 +32,6 @@
 #include <opm/autodiff/NewtonIterationBlackoilInterface.hpp>
 #include <opm/autodiff/BlackoilModelEnums.hpp>
 #include <opm/autodiff/VFPProperties.hpp>
-#include <opm/autodiff/StandardWells.hpp>
 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
 
 #include <array>
@@ -105,8 +104,9 @@ namespace Opm {
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
     /// \tparam  Grid            UnstructuredGrid or CpGrid.
+    /// \tparam  WellModel       WellModel employed.
     /// \tparam  Implementation  Provides concrete state types.
-    template<class Grid, class Implementation>
+    template<class Grid, class WellModel, class Implementation>
     class BlackoilModelBase
     {
     public:
@@ -276,7 +276,7 @@ namespace Opm {
         const BlackoilPropsAdInterface& fluid_;
         const DerivedGeology&           geo_;
         const RockCompressibility*      rock_comp_props_;
-        StandardWells                   std_wells_;
+        WellModel                       std_wells_;
         VFPProperties                   vfp_properties_;
         const NewtonIterationBlackoilInterface&    linsolver_;
         // For each canonical phase -> true if active
@@ -328,11 +328,11 @@ namespace Opm {
             return static_cast<const Implementation&>(*this);
         }
 
-        /// return the StandardWells object
-        StandardWells& stdWells() { return std_wells_; }
-        const StandardWells& stdWells() const { return std_wells_; }
+        /// return the WellModel object
+        WellModel& stdWells() { return std_wells_; }
+        const WellModel& stdWells() const { return std_wells_; }
 
-        /// return the Well struct in the StandardWells
+        /// return the Well struct in the WellModel
         const Wells& wells() const { return std_wells_.wells(); }
 
         /// return true if wells are available in the reservoir
diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp
index dffd817698..a9c485b386 100644
--- a/opm/autodiff/BlackoilModelBase_impl.hpp
+++ b/opm/autodiff/BlackoilModelBase_impl.hpp
@@ -143,8 +143,8 @@ namespace detail {
 } // namespace detail
 
 
-    template <class Grid, class Implementation>
-    BlackoilModelBase<Grid, Implementation>::
+    template <class Grid, class WellModel, class Implementation>
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     BlackoilModelBase(const ModelParameters&          param,
                   const Grid&                     grid ,
                   const BlackoilPropsAdInterface& fluid,
@@ -217,9 +217,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     prepareStep(const double dt,
                 ReservoirState& reservoir_state,
                 WellState& /* well_state */)
@@ -234,10 +234,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     template <class NonlinearSolverType>
     IterationReport
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     nonlinearIteration(const int iteration,
                        const double dt,
                        NonlinearSolverType& nonlinear_solver,
@@ -289,9 +289,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     afterStep(const double /* dt */,
               ReservoirState& /* reservoir_state */,
               WellState& /* well_state */)
@@ -303,9 +303,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     int
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     sizeNonLinear() const
     {
         return residual_.sizeNonLinear();
@@ -315,9 +315,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     int
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     linearIterationsLastSolve() const
     {
         return linsolver_.iterations();
@@ -327,9 +327,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     bool
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     terminalOutputEnabled() const
     {
         return terminal_output_;
@@ -339,9 +339,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     int
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     numPhases() const
     {
         return fluid_.numPhases();
@@ -351,9 +351,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     int
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     numMaterials() const
     {
         return material_name_.size();
@@ -363,9 +363,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     const std::string&
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     materialName(int material_index) const
     {
         assert(material_index < numMaterials());
@@ -376,9 +376,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     setThresholdPressures(const std::vector<double>& threshold_pressures)
     {
         const int num_faces = AutoDiffGrid::numFaces(grid_);
@@ -407,8 +407,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    BlackoilModelBase<Grid, Implementation>::ReservoirResidualQuant::ReservoirResidualQuant()
+    template <class Grid, class WellModel, class Implementation>
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    ReservoirResidualQuant::ReservoirResidualQuant()
         : accum(2, ADB::null())
         , mflux(   ADB::null())
         , b    (   ADB::null())
@@ -421,9 +422,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     int
-    BlackoilModelBase<Grid, Implementation>::numWellVars() const
+    BlackoilModelBase<Grid, WellModel, Implementation>::numWellVars() const
     {
         // For each well, we have a bhp variable, and one flux per phase.
         const int nw = stdWells().localWellsActive() ? wells().number_of_wells : 0;
@@ -434,9 +435,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::makeConstantState(SolutionState& state) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    makeConstantState(SolutionState& state) const
     {
         // HACK: throw away the derivatives. this may not be the most
         // performant way to do things, but it will make the state
@@ -463,10 +465,11 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    typename BlackoilModelBase<Grid, Implementation>::SolutionState
-    BlackoilModelBase<Grid, Implementation>::variableState(const ReservoirState& x,
-                                                           const WellState&     xw) const
+    template <class Grid, class WellModel, class Implementation>
+    typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableState(const ReservoirState& x,
+                  const WellState&     xw) const
     {
         std::vector<V> vars0 = asImpl().variableStateInitials(x, xw);
         std::vector<ADB> vars = ADB::variables(vars0);
@@ -477,10 +480,11 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<V>
-    BlackoilModelBase<Grid, Implementation>::variableStateInitials(const ReservoirState& x,
-                                                                   const WellState&     xw) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableStateInitials(const ReservoirState& x,
+                          const WellState&     xw) const
     {
         assert(active_[ Oil ]);
 
@@ -499,9 +503,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::variableReservoirStateInitials(const ReservoirState& x, std::vector<V>& vars0) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableReservoirStateInitials(const ReservoirState& x, std::vector<V>& vars0) const
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -537,9 +542,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::variableWellStateInitials(const WellState&     xw, std::vector<V>& vars0) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableWellStateInitials(const WellState& xw, std::vector<V>& vars0) const
     {
         // Initial well rates.
         if ( stdWells().localWellsActive() )
@@ -571,9 +577,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<int>
-    BlackoilModelBase<Grid, Implementation>::variableStateIndices() const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableStateIndices() const
     {
         assert(active_[Oil]);
         std::vector<int> indices(5, -1);
@@ -594,9 +601,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<int>
-    BlackoilModelBase<Grid, Implementation>::variableWellStateIndices() const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableWellStateIndices() const
     {
         // Black oil model standard is 5 equation.
         // For the pure well solve, only the well equations are picked.
@@ -612,11 +620,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    typename BlackoilModelBase<Grid, Implementation>::SolutionState
-    BlackoilModelBase<Grid, Implementation>::variableStateExtractVars(const ReservoirState& x,
-                                                                      const std::vector<int>& indices,
-                                                                      std::vector<ADB>& vars) const
+    template <class Grid, class WellModel, class Implementation>
+    typename BlackoilModelBase<Grid, WellModel, Implementation>::SolutionState
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableStateExtractVars(const ReservoirState& x,
+                             const std::vector<int>& indices,
+                             std::vector<ADB>& vars) const
     {
         //using namespace Opm::AutoDiffGrid;
         const int nc = Opm::AutoDiffGrid::numCells(grid_);
@@ -684,11 +693,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::variableStateExtractWellsVars(const std::vector<int>& indices,
-                                                                          std::vector<ADB>& vars,
-                                                                          SolutionState& state) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    variableStateExtractWellsVars(const std::vector<int>& indices,
+                                  std::vector<ADB>& vars,
+                                  SolutionState& state) const
     {
         // Qs.
         state.qs = std::move(vars[indices[Qs]]);
@@ -701,10 +711,11 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::computeAccum(const SolutionState& state,
-                                              const int            aix  )
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeAccum(const SolutionState& state,
+                 const int            aix  )
     {
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
 
@@ -763,9 +774,11 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    void BlackoilModelBase<Grid, Implementation>::computeWellConnectionPressures(const SolutionState& state,
-                                                                                 const WellState& xw)
+    template <class Grid, class WellModel, class Implementation>
+    void 
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeWellConnectionPressures(const SolutionState& state,
+                                   const WellState& xw)
     {
         if( ! localWellsActive() ) return ;
 
@@ -780,8 +793,8 @@ namespace detail {
         asImpl().stdWells().computePropertiesForWellConnectionPressures(state, xw, fluid_, active_, phaseCondition_, b_perf, rsmax_perf, rvmax_perf, surf_dens_perf);
 
         // Extract well connection depths.
-        const StandardWells::Vector depth = cellCentroidsZToEigen(grid_);
-        const StandardWells::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells);
+        const typename WellModel::Vector depth = cellCentroidsZToEigen(grid_);
+        const typename WellModel::Vector pdepth = subset(depth, asImpl().stdWells().wellOps().well_cells);
         const int nperf = wells().well_connpos[wells().number_of_wells];
         const std::vector<double> depth_perf(pdepth.data(), pdepth.data() + nperf);
 
@@ -796,9 +809,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     assemble(const ReservoirState& reservoir_state,
              WellState& well_state,
              const bool initial_assembly)
@@ -874,9 +887,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     assembleMassBalanceEq(const SolutionState& state)
     {
         // Compute b_p and the accumulation term b_p*s_p for each phase,
@@ -942,9 +955,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::updateEquationsScaling() {
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    updateEquationsScaling() {
         ADB::V B;
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
         for ( int idx=0; idx<MaxNumPhases; ++idx )
@@ -975,11 +989,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
-                                                                                const SolutionState&,
-                                                                                const WellState&)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    addWellContributionToMassBalanceEq(const std::vector<ADB>& cq_s,
+                                       const SolutionState&,
+                                       const WellState&)
     {
         if ( !asImpl().localWellsActive() )
         {
@@ -1001,11 +1016,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::extractWellPerfProperties(const SolutionState&,
-                                                                       std::vector<ADB>& mob_perfcells,
-                                                                       std::vector<ADB>& b_perfcells) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    extractWellPerfProperties(const SolutionState&,
+                              std::vector<ADB>& mob_perfcells,
+                              std::vector<ADB>& b_perfcells) const
     {
         // If we have wells, extract the mobilities and b-factors for
         // the well-perforated cells.
@@ -1029,9 +1045,11 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    void BlackoilModelBase<Grid, Implementation>::addWellFluxEq(const std::vector<ADB>& cq_s,
-                                                                const SolutionState& state)
+    template <class Grid, class WellModel, class Implementation>
+    void
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    addWellFluxEq(const std::vector<ADB>& cq_s,
+                  const SolutionState& state)
     {
         if( !asImpl().localWellsActive() )
         {
@@ -1055,8 +1073,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    bool BlackoilModelBase<Grid, Implementation>::isVFPActive() const
+    template <class Grid, class WellModel, class Implementation>
+    bool
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    isVFPActive() const
     {
         if( ! localWellsActive() ) {
             return false;
@@ -1089,11 +1109,13 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    bool BlackoilModelBase<Grid, Implementation>::solveWellEq(const std::vector<ADB>& mob_perfcells,
-                                                              const std::vector<ADB>& b_perfcells,
-                                                              SolutionState& state,
-                                                              WellState& well_state)
+    template <class Grid, class WellModel, class Implementation>
+    bool
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    solveWellEq(const std::vector<ADB>& mob_perfcells,
+                const std::vector<ADB>& b_perfcells,
+                SolutionState& state,
+                WellState& well_state)
     {
         V aliveWells;
         const int np = wells().number_of_phases;
@@ -1196,10 +1218,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    void BlackoilModelBase<Grid, Implementation>::addWellControlEq(const SolutionState& state,
-                                                          const WellState& xw,
-                                                          const V& aliveWells)
+    template <class Grid, class WellModel, class Implementation>
+    void
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    addWellControlEq(const SolutionState& state,
+                     const WellState& xw,
+                     const V& aliveWells)
     {
         if( ! localWellsActive() ) return;
 
@@ -1362,8 +1386,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    V BlackoilModelBase<Grid, Implementation>::solveJacobianSystem() const
+    template <class Grid, class WellModel, class Implementation>
+    V
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    solveJacobianSystem() const
     {
         return linsolver_.computeNewtonIncrement(residual_);
     }
@@ -1473,10 +1499,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
-    void BlackoilModelBase<Grid, Implementation>::updateState(const V& dx,
-                                          ReservoirState& reservoir_state,
-                                          WellState& well_state)
+    template <class Grid, class WellModel, class Implementation>
+    void
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    updateState(const V& dx,
+                ReservoirState& reservoir_state,
+                WellState& well_state)
     {
         using namespace Opm::AutoDiffGrid;
         const int np = fluid_.numPhases();
@@ -1695,9 +1723,12 @@ namespace detail {
     }
 
 
-    template <class Grid, class Implementation>
+
+
+
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     computeWellPotentials(const SolutionState& state,
                           const std::vector<ADB>& mob_perfcells,
                           const std::vector<ADB>& b_perfcells,
@@ -1800,9 +1831,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<ADB>
-    BlackoilModelBase<Grid, Implementation>::computeRelPerm(const SolutionState& state) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeRelPerm(const SolutionState& state) const
     {
         using namespace Opm::AutoDiffGrid;
         const int               nc   = numCells(grid_);
@@ -1829,9 +1861,9 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<ADB>
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     computePressures(const ADB& po,
                      const ADB& sw,
                      const ADB& so,
@@ -1866,12 +1898,13 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     V
-    BlackoilModelBase<Grid, Implementation>::computeGasPressure(const V& po,
-                                                       const V& sw,
-                                                       const V& so,
-                                                       const V& sg) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeGasPressure(const V& po,
+                       const V& sw,
+                       const V& so,
+                       const V& sg) const
     {
         assert (active_[Gas]);
         std::vector<ADB> cp = fluid_.capPress(ADB::constant(sw),
@@ -1883,15 +1916,16 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::computeMassFlux(const int               actph ,
-                                                             const V&                transi,
-                                                             const ADB&              kr    ,
-                                                             const ADB&              mu    ,
-                                                             const ADB&              rho   ,
-                                                             const ADB&              phasePressure,
-                                                             const SolutionState&    state)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeMassFlux(const int               actph ,
+                    const V&                transi,
+                    const ADB&              kr    ,
+                    const ADB&              mu    ,
+                    const ADB&              rho   ,
+                    const ADB&              phasePressure,
+                    const SolutionState&    state)
     {
         // Compute and store mobilities.
         const ADB tr_mult = transMult(state.pressure);
@@ -1916,9 +1950,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::applyThresholdPressures(ADB& dp)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    applyThresholdPressures(ADB& dp)
     {
         // We support reversible threshold pressures only.
         // Method: if the potential difference is lower (in absolute
@@ -1948,9 +1983,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     std::vector<double>
-    BlackoilModelBase<Grid, Implementation>::computeResidualNorms() const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    computeResidualNorms() const
     {
         std::vector<double> residualNorms;
 
@@ -1988,9 +2024,9 @@ namespace detail {
     }
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     double
-    BlackoilModelBase<Grid, Implementation>::
+    BlackoilModelBase<Grid, WellModel, Implementation>::
     relativeChange(const SimulationDataContainer& previous,
                    const SimulationDataContainer& current ) const
     {
@@ -2029,16 +2065,17 @@ namespace detail {
         }
     }
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     double
-    BlackoilModelBase<Grid, Implementation>::convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
-                                                                  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
-                                                                  const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
-                                                                  std::vector<double>& R_sum,
-                                                                  std::vector<double>& maxCoeff,
-                                                                  std::vector<double>& B_avg,
-                                                                  std::vector<double>& maxNormWell,
-                                                                  int nc) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    convergenceReduction(const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& B,
+                         const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& tempV,
+                         const Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>& R,
+                         std::vector<double>& R_sum,
+                         std::vector<double>& maxCoeff,
+                         std::vector<double>& B_avg,
+                         std::vector<double>& maxNormWell,
+                         int nc) const
     {
         const int np = asImpl().numPhases();
         const int nm = asImpl().numMaterials();
@@ -2115,9 +2152,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     bool
-    BlackoilModelBase<Grid, Implementation>::getConvergence(const double dt, const int iteration)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    getConvergence(const double dt, const int iteration)
     {
         const double tol_mb    = param_.tolerance_mb_;
         const double tol_cnv   = param_.tolerance_cnv_;
@@ -2239,9 +2277,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     bool
-    BlackoilModelBase<Grid, Implementation>::getWellConvergence(const int iteration)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    getWellConvergence(const int iteration)
     {
         const double tol_wells = param_.tolerance_wells_;
 
@@ -2318,14 +2357,15 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::fluidViscosity(const int               phase,
-                                                            const ADB&              p    ,
-                                                            const ADB&              temp ,
-                                                            const ADB&              rs   ,
-                                                            const ADB&              rv   ,
-                                                            const std::vector<PhasePresence>& cond) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidViscosity(const int               phase,
+                   const ADB&              p    ,
+                   const ADB&              temp ,
+                   const ADB&              rs   ,
+                   const ADB&              rv   ,
+                   const std::vector<PhasePresence>& cond) const
     {
         switch (phase) {
         case Water:
@@ -2343,14 +2383,15 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::fluidReciprocFVF(const int               phase,
-                                                              const ADB&              p    ,
-                                                              const ADB&              temp ,
-                                                              const ADB&              rs   ,
-                                                              const ADB&              rv   ,
-                                                              const std::vector<PhasePresence>& cond) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidReciprocFVF(const int               phase,
+                     const ADB&              p    ,
+                     const ADB&              temp ,
+                     const ADB&              rs   ,
+                     const ADB&              rv   ,
+                     const std::vector<PhasePresence>& cond) const
     {
         switch (phase) {
         case Water:
@@ -2368,12 +2409,13 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::fluidDensity(const int  phase,
-                                                          const ADB& b,
-                                                          const ADB& rs,
-                                                          const ADB& rv) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidDensity(const int  phase,
+                 const ADB& b,
+                 const ADB& rs,
+                 const ADB& rv) const
     {
         const V& rhos = fluid_.surfaceDensity(phase,  cells_);
         const Opm::PhaseUsage& pu = fluid_.phaseUsage();
@@ -2391,11 +2433,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     V
-    BlackoilModelBase<Grid, Implementation>::fluidRsSat(const V&                p,
-                                               const V&                satOil,
-                                               const std::vector<int>& cells) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidRsSat(const V&                p,
+               const V&                satOil,
+               const std::vector<int>& cells) const
     {
         return fluid_.rsSat(ADB::constant(p), ADB::constant(satOil), cells).value();
     }
@@ -2404,11 +2447,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::fluidRsSat(const ADB&              p,
-                                               const ADB&              satOil,
-                                               const std::vector<int>& cells) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidRsSat(const ADB&              p,
+               const ADB&              satOil,
+               const std::vector<int>& cells) const
     {
         return fluid_.rsSat(p, satOil, cells);
     }
@@ -2417,11 +2461,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     V
-    BlackoilModelBase<Grid, Implementation>::fluidRvSat(const V&                p,
-                                               const V&              satOil,
-                                               const std::vector<int>& cells) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidRvSat(const V&                p,
+               const V&              satOil,
+               const std::vector<int>& cells) const
     {
         return fluid_.rvSat(ADB::constant(p), ADB::constant(satOil), cells).value();
     }
@@ -2430,11 +2475,12 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::fluidRvSat(const ADB&              p,
-                                               const ADB&              satOil,
-                                               const std::vector<int>& cells) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    fluidRvSat(const ADB&              p,
+               const ADB&              satOil,
+               const std::vector<int>& cells) const
     {
         return fluid_.rvSat(p, satOil, cells);
     }
@@ -2443,9 +2489,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::poroMult(const ADB& p) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    poroMult(const ADB& p) const
     {
         const int n = p.size();
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -2473,9 +2520,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     ADB
-    BlackoilModelBase<Grid, Implementation>::transMult(const ADB& p) const
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    transMult(const ADB& p) const
     {
         const int n = p.size();
         if (rock_comp_props_ && rock_comp_props_->isActive()) {
@@ -2503,9 +2551,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::classifyCondition(const ReservoirState& state)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    classifyCondition(const ReservoirState& state)
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -2543,9 +2592,10 @@ namespace detail {
 
 
 
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::updatePrimalVariableFromState(const ReservoirState& state)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    updatePrimalVariableFromState(const ReservoirState& state)
     {
         using namespace Opm::AutoDiffGrid;
         const int nc = numCells(grid_);
@@ -2593,9 +2643,10 @@ namespace detail {
 
 
     /// Update the phaseCondition_ member based on the primalVariable_ member.
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::updatePhaseCondFromPrimalVariable()
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    updatePhaseCondFromPrimalVariable()
     {
         if (! active_[Gas]) {
             OPM_THROW(std::logic_error, "updatePhaseCondFromPrimarVariable() logic requires active gas phase.");
@@ -2633,10 +2684,11 @@ namespace detail {
 
     // TODO: only kept for now due to flow_multisegment
     // will be removed soon
-    template <class Grid, class Implementation>
+    template <class Grid, class WellModel, class Implementation>
     void
-    BlackoilModelBase<Grid, Implementation>::updateWellState(const V& dwells,
-                                                             WellState& well_state)
+    BlackoilModelBase<Grid, WellModel, Implementation>::
+    updateWellState(const V& dwells,
+                    WellState& well_state)
     {
         const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
         asImpl().stdWells().updateWellState(dwells, gravity, dpMaxRel(), fluid_.phaseUsage(),
diff --git a/opm/autodiff/BlackoilMultiSegmentModel.hpp b/opm/autodiff/BlackoilMultiSegmentModel.hpp
index 1e2c68c315..3ff7094b29 100644
--- a/opm/autodiff/BlackoilMultiSegmentModel.hpp
+++ b/opm/autodiff/BlackoilMultiSegmentModel.hpp
@@ -25,6 +25,7 @@
 #include <opm/autodiff/BlackoilModelParameters.hpp>
 #include <opm/autodiff/WellStateMultiSegment.hpp>
 #include <opm/autodiff/WellMultiSegment.hpp>
+#include <opm/autodiff/StandardWells.hpp>
 
 namespace Opm {
 
@@ -48,11 +49,11 @@ namespace Opm {
     /// \tparam  Grid            UnstructuredGrid or CpGrid.
     /// \tparam  Implementation  Provides concrete state types.
     template<class Grid>
-    class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid>>
+    class BlackoilMultiSegmentModel : public BlackoilModelBase<Grid, StandardWells, BlackoilMultiSegmentModel<Grid> >
     {
     public:
 
-        typedef BlackoilModelBase<Grid, BlackoilMultiSegmentModel<Grid> > Base; // base class
+        typedef BlackoilModelBase<Grid, StandardWells, BlackoilMultiSegmentModel<Grid> > Base; // base class
         typedef typename Base::ReservoirState ReservoirState;
         typedef typename Base::WellState WellState;
         typedef BlackoilMultiSegmentSolutionState SolutionState;
diff --git a/opm/autodiff/BlackoilSolventModel.hpp b/opm/autodiff/BlackoilSolventModel.hpp
index c7881a3a4f..ddc3aa40ca 100644
--- a/opm/autodiff/BlackoilSolventModel.hpp
+++ b/opm/autodiff/BlackoilSolventModel.hpp
@@ -36,18 +36,18 @@ namespace Opm {
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
     template<class Grid>
-    class BlackoilSolventModel : public BlackoilModelBase<Grid, BlackoilSolventModel<Grid> >
+    class BlackoilSolventModel : public BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> >
     {
     public:
 
         // ---------  Types and enums  ---------
 
-        typedef BlackoilModelBase<Grid, BlackoilSolventModel<Grid> > Base;
+        typedef BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> > Base;
         typedef typename Base::ReservoirState ReservoirState;
         typedef typename Base::WellState WellState;
         // The next line requires C++11 support available in g++ 4.7.
         // friend Base;
-        friend class BlackoilModelBase<Grid, BlackoilSolventModel<Grid> >;
+        friend class BlackoilModelBase<Grid, StandardWellsSolvent, BlackoilSolventModel<Grid> >;
 
         /// Construct the model. It will retain references to the
         /// arguments of this functions, and they are expected to
@@ -104,10 +104,6 @@ namespace Opm {
         const bool is_miscible_;
         std::vector<ADB> mu_eff_;
         std::vector<ADB> b_eff_;
-        StandardWellsSolvent std_wells_;
-        const StandardWellsSolvent& stdWells() const { return std_wells_; }
-        StandardWellsSolvent& stdWells() { return std_wells_; }
-
 
         // Need to declare Base members we want to use here.
         using Base::grid_;
@@ -134,7 +130,7 @@ namespace Opm {
         // ---------  Protected methods  ---------
 
         // Need to declare Base members we want to use here.
-        // using Base::stdWells;
+        using Base::stdWells;
         using Base::wells;
         using Base::variableState;
         using Base::computeGasPressure;
diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp
index bb86bc5785..8dfd0a15e9 100644
--- a/opm/autodiff/BlackoilSolventModel_impl.hpp
+++ b/opm/autodiff/BlackoilSolventModel_impl.hpp
@@ -89,8 +89,7 @@ namespace Opm {
           has_solvent_(has_solvent),
           solvent_pos_(detail::solventPos(fluid.phaseUsage())),
           solvent_props_(solvent_props),
-          is_miscible_(is_miscible),
-          std_wells_(wells_arg, solvent_props, solvent_pos_)
+          is_miscible_(is_miscible)
 
     {
         if (has_solvent_) {
@@ -105,6 +104,8 @@ namespace Opm {
             }
 
             residual_.matbalscale.resize(fluid_.numPhases() + 1, 0.0031); // use the same as gas
+
+            stdWells().initSolvent(&solvent_props_, solvent_pos_);
         }
         if (is_miscible_) {
             mu_eff_.resize(fluid_.numPhases() + 1, ADB::null());
diff --git a/opm/autodiff/StandardWellsSolvent.hpp b/opm/autodiff/StandardWellsSolvent.hpp
index 083b30ab7b..2017824782 100644
--- a/opm/autodiff/StandardWellsSolvent.hpp
+++ b/opm/autodiff/StandardWellsSolvent.hpp
@@ -36,7 +36,10 @@ namespace Opm {
             using Base = StandardWells;
 
             // ---------  Public methods  ---------
-            explicit StandardWellsSolvent(const Wells* wells, const SolventPropsAdFromDeck& solvent_props, const int solvent_pos);
+            explicit StandardWellsSolvent(const Wells* wells);
+
+            // added the Solvent related
+            void initSolvent(const SolventPropsAdFromDeck* solvent_props, const int solvent_pos);
 
             template <class SolutionState, class WellState>
             void computePropertiesForWellConnectionPressures(const SolutionState& state,
@@ -59,8 +62,8 @@ namespace Opm {
                                            std::vector<ADB>& mob_perfcells,
                                            std::vector<ADB>& b_perfcells) const;
         protected:
-            const SolventPropsAdFromDeck& solvent_props_;
-            const int solvent_pos_;
+            const SolventPropsAdFromDeck* solvent_props_;
+            int solvent_pos_;
 
         };
 
diff --git a/opm/autodiff/StandardWellsSolvent_impl.hpp b/opm/autodiff/StandardWellsSolvent_impl.hpp
index ee69a342d3..b55991c353 100644
--- a/opm/autodiff/StandardWellsSolvent_impl.hpp
+++ b/opm/autodiff/StandardWellsSolvent_impl.hpp
@@ -32,12 +32,10 @@ namespace Opm
 
 
 
-    StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg,
-                                               const SolventPropsAdFromDeck& solvent_props,
-                                               const int solvent_pos)
+    StandardWellsSolvent::StandardWellsSolvent(const Wells* wells_arg)
         : Base(wells_arg)
-        , solvent_props_(solvent_props)
-        , solvent_pos_(solvent_pos)
+        , solvent_props_(nullptr)
+        , solvent_pos_(-1)
     {
     }
 
@@ -45,6 +43,16 @@ namespace Opm
 
 
 
+    void StandardWellsSolvent::initSolvent(const SolventPropsAdFromDeck* solvent_props, const int solvent_pos)
+    {
+        solvent_props_ = solvent_props;
+        solvent_pos_ = solvent_pos;
+    }
+
+
+
+
+
     template<class SolutionState, class WellState>
     void
     StandardWellsSolvent::
@@ -129,7 +137,7 @@ namespace Opm
             // to handle solvent related
             {
 
-                const Vector bs = solvent_props_.bSolvent(avg_press_ad,well_cells).value();
+                const Vector bs = solvent_props_->bSolvent(avg_press_ad,well_cells).value();
                 //const V bs_eff = subset(rq_[solvent_pos_].b,well_cells).value();
 
                 // number of cells
@@ -161,7 +169,7 @@ namespace Opm
                 bg = bg * (ones - F_solvent);
                 bg = bg + F_solvent * bs;
 
-                const Vector& rhos = solvent_props_.solventSurfaceDensity(well_cells);
+                const Vector& rhos = solvent_props_->solventSurfaceDensity(well_cells);
                 rhog = ( (ones - F_solvent) * rhog ) + (F_solvent * rhos);
             }
             b.col(pu.phase_pos[BlackoilPhases::Vapour]) = bg;
diff --git a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
index 8c4f00a2f8..71f513552b 100644
--- a/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
+++ b/opm/polymer/fullyimplicit/BlackoilPolymerModel.hpp
@@ -27,6 +27,7 @@
 #include <opm/polymer/fullyimplicit/PolymerPropsAd.hpp>
 #include <opm/polymer/PolymerBlackoilState.hpp>
 #include <opm/polymer/fullyimplicit/WellStateFullyImplicitBlackoilPolymer.hpp>
+#include <opm/autodiff/StandardWells.hpp>
 
 namespace Opm {
 
@@ -40,18 +41,18 @@ namespace Opm {
     /// It uses automatic differentiation via the class AutoDiffBlock
     /// to simplify assembly of the jacobian matrix.
     template<class Grid>
-    class BlackoilPolymerModel : public BlackoilModelBase<Grid, BlackoilPolymerModel<Grid> >
+    class BlackoilPolymerModel : public BlackoilModelBase<Grid, StandardWells, BlackoilPolymerModel<Grid> >
     {
     public:
 
         // ---------  Types and enums  ---------
 
-        typedef BlackoilModelBase<Grid, BlackoilPolymerModel<Grid> > Base;
+        typedef BlackoilModelBase<Grid, StandardWells, BlackoilPolymerModel<Grid> > Base;
         typedef typename Base::ReservoirState ReservoirState;
         typedef typename Base::WellState WellState;
         // The next line requires C++11 support available in g++ 4.7.
         // friend Base;
-        friend class BlackoilModelBase<Grid, BlackoilPolymerModel<Grid> >;
+        friend class BlackoilModelBase<Grid, StandardWells, BlackoilPolymerModel<Grid> >;
 
         /// Construct the model. It will retain references to the
         /// arguments of this functions, and they are expected to