diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index 20879065a..4304635d3 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -321,7 +321,11 @@ namespace Opm {
                 // Apply the update, with considering model-dependent limitations and
                 // chopping of the update.
                 updateState(x,iteration);
-                wellModel().updateWellState(xw, well_state);
+
+                if( nw > 0 )
+                {
+                    wellModel().applySolutionWellState(x, well_state);
+                }
                 report.update_time += perfTimer.stop();
             }
 
diff --git a/opm/autodiff/StandardWell.hpp b/opm/autodiff/StandardWell.hpp
index 5de51cd2c..356599af0 100644
--- a/opm/autodiff/StandardWell.hpp
+++ b/opm/autodiff/StandardWell.hpp
@@ -141,6 +141,11 @@ namespace Opm
         // r = r - C D^-1 Rw
         virtual void apply(BVector& r) const;
 
+        // using the solution x to recover the solution xw for wells and applying
+        // xw to update Well State
+        virtual void applySolutionWellState(const BVector& x, const ModelParameters& param,
+                                            WellState& well_state) const;
+
         using WellInterface<TypeTag>::phaseUsage;
         using WellInterface<TypeTag>::active;
         using WellInterface<TypeTag>::numberOfPerforations;
@@ -167,6 +172,9 @@ namespace Opm
         // TODO: maybe this function can go to some helper file.
         void localInvert(Mat& istlA) const;
 
+        // xw = inv(D)*(rw - C*x)
+        void recoverSolutionWell(const BVector& x, BVector& xw) const;
+
         // TODO: decide wether to use member function to refer to private member later
         using WellInterface<TypeTag>::vfp_properties_;
         using WellInterface<TypeTag>::gravity_;
diff --git a/opm/autodiff/StandardWell_impl.hpp b/opm/autodiff/StandardWell_impl.hpp
index 274eaa1f3..a9973877e 100644
--- a/opm/autodiff/StandardWell_impl.hpp
+++ b/opm/autodiff/StandardWell_impl.hpp
@@ -1904,4 +1904,36 @@ namespace Opm
         duneC_.mmtv(invDrw_, r);
     }
 
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    recoverSolutionWell(const BVector& x, BVector& xw) const
+    {
+        BVector resWell = resWell_;
+        // resWell = resWell - B * x
+        duneB_.mmv(x, resWell);
+        // xw = D^-1 * resWell
+        invDuneD_.mv(resWell, xw);
+    }
+
+
+
+
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    applySolutionWellState(const BVector& x,
+                           const ModelParameters& param,
+                           WellState& well_state) const
+    {
+        BVector xw(1);
+        recoverSolutionWell(x, xw);
+        updateWellState(xw, param, well_state);
+    }
+
 }
diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp
index 9c62c10d7..c81835c39 100644
--- a/opm/autodiff/StandardWellsDense.hpp
+++ b/opm/autodiff/StandardWellsDense.hpp
@@ -98,6 +98,8 @@ enum WellVariablePositions {
             static const int solventSaturationIdx = BlackoilIndices::solventSaturationIdx;
             static const int polymerConcentrationIdx = BlackoilIndices::polymerConcentrationIdx;
 
+            // TODO: where we should put these types, WellInterface or Well Model?
+            // or there is some other strategy, like TypeTag
             typedef Dune::FieldVector<Scalar, numEq    > VectorBlockType;
             typedef Dune::FieldMatrix<Scalar, numEq, numEq > MatrixBlockType;
             typedef Dune::BCRSMatrix <MatrixBlockType> Mat;
@@ -175,8 +177,9 @@ enum WellVariablePositions {
             // apply well model with scaling of alpha
             void applyScaleAdd(const Scalar alpha, const BVector& x, BVector& Ax) const;
 
-            // xw = inv(D)*(rw - C*x)
-            void recoverVariable(const BVector& x, BVector& xw) const;
+            // using the solution x to recover the solution xw for wells and applying
+            // xw to update Well State
+            void applySolutionWellState(const BVector& x, WellState& well_state) const;
 
             int flowPhaseToEbosCompIdx( const int phaseIdx ) const;
 
diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp
index 1577e6e08..64e81ae0e 100644
--- a/opm/autodiff/StandardWellsDense_impl.hpp
+++ b/opm/autodiff/StandardWellsDense_impl.hpp
@@ -431,6 +431,7 @@ namespace Opm {
     StandardWellsDense<TypeTag>::
     apply(const BVector& x, BVector& Ax) const
     {
+        // TODO: do we still need localWellsActive()?
         if ( ! localWellsActive() ) {
             return;
         }
@@ -486,28 +487,20 @@ namespace Opm {
 
 
 
-    // xw = D^-1(resWell - B * x)
-    // TODO: this function should be moved to StandardWell
-    // xw should not appear in StandardWellsDense to avoid
-    // the data type to hold different types of xw
     template<typename TypeTag>
     void
     StandardWellsDense<TypeTag>::
-    recoverVariable(const BVector& x, BVector& xw) const
+    applySolutionWellState(const BVector& x, WellState& well_state) const
     {
-        if ( ! localWellsActive() ) {
-             return;
+        for (auto& well : well_container_) {
+            well->applySolutionWellState(x, param_, well_state);
         }
-        BVector resWell = resWell_;
-        // resWell = resWell - B * x
-        duneB_.mmv(x, resWell);
-        // xw = D^-1 * resWell
-        invDuneD_.mv(resWell, xw);
     }
 
 
 
 
+
     template<typename TypeTag>
     int
     StandardWellsDense<TypeTag>::
diff --git a/opm/autodiff/WellInterface.hpp b/opm/autodiff/WellInterface.hpp
index 6423ef519..c9fb75e99 100644
--- a/opm/autodiff/WellInterface.hpp
+++ b/opm/autodiff/WellInterface.hpp
@@ -192,6 +192,11 @@ namespace Opm
         // r = r - C D^-1 Rw
         virtual void apply(BVector& r) const = 0;
 
+        // using the solution x to recover the solution xw for wells and applying
+        // xw to update Well State
+        virtual void applySolutionWellState(const BVector& x, const ModelParameters& param,
+                                            WellState& well_state) const = 0;
+
     protected:
         // TODO: some variables shared by all the wells should be made static
         // well name