diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index c3435e589..78ebb2cde 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -479,14 +479,17 @@ namespace Opm {
             if( isParallel() )
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, true > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this), istlSolver().parallelInformation() );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_,
+                             istlSolver().parallelInformation() );
                 assert( opA.comm() );
                 istlSolver().solve( opA, x, ebosResid, *(opA.comm()) );
             }
             else
             {
                 typedef WellModelMatrixAdapter< Mat, BVector, BVector, ThisType, false > Operator;
-                Operator opA(ebosJac, const_cast< ThisType& > (*this) );
+                Operator opA(ebosJac, const_cast< ThisType& > (*this),
+                             param_.matrix_add_well_contributions_ );
                 istlSolver().solve( opA, x, ebosResid );
             }
 
@@ -532,8 +535,11 @@ namespace Opm {
           };
 
           //! constructor: just store a reference to a matrix
-          WellModelMatrixAdapter (const M& A, WellModel& wellMod, const boost::any& parallelInformation = boost::any() )
-              : A_( A ), wellMod_( wellMod ), comm_()
+          WellModelMatrixAdapter (const M& A, WellModel& wellMod,
+                                  bool matrix_add_well_contributions,
+                                  const boost::any& parallelInformation = boost::any() )
+              : A_( A ), wellMod_( wellMod ), comm_(),
+                matrix_add_well_contributions_(matrix_add_well_contributions)
           {
 #if HAVE_MPI
             if( parallelInformation.type() == typeid(ParallelISTLInformation) )
@@ -548,22 +554,29 @@ namespace Opm {
           virtual void apply( const X& x, Y& y ) const
           {
             A_.mv( x, y );
-            // add well model modification to y
-            wellMod_.applyWellModelAdd(x, y );
+
+            if ( ! matrix_add_well_contributions_ )
+            {
+              // add well model modification to y
+              wellMod_.applyWellModelAdd(x, y );
 
 #if HAVE_MPI
-            if( comm_ )
-              comm_->project( y );
+              if( comm_ )
+                comm_->project( y );
 #endif
+            }
           }
 
           // y += \alpha * A * x
           virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
           {
             A_.usmv(alpha,x,y);
-            // add scaled well model modification to y
-            wellMod_.applyWellModelScaleAdd( alpha, x, y );
 
+            if ( ! matrix_add_well_contributions_ )
+            {
+                // add scaled well model modification to y
+                wellMod_.applyWellModelScaleAdd( alpha, x, y );
+            }
 #if HAVE_MPI
             if( comm_ )
               comm_->project( y );
@@ -581,6 +594,7 @@ namespace Opm {
           const matrix_type& A_ ;
           WellModel& wellMod_;
           std::unique_ptr< communication_type > comm_;
+          bool matrix_add_well_contributions_;
         };
 
         /// Apply an update to the primary variables, chopped if appropriate.
diff --git a/opm/autodiff/StandardWellsDense.hpp b/opm/autodiff/StandardWellsDense.hpp
index 5c2906a49..837e2d52c 100644
--- a/opm/autodiff/StandardWellsDense.hpp
+++ b/opm/autodiff/StandardWellsDense.hpp
@@ -58,6 +58,27 @@
 #include <opm/simulators/WellSwitchingLogger.hpp>
 
 namespace Opm {
+namespace Detail
+{
+    //! calculates ret = A^T * B
+    template< class K, int m, int n, int p >
+    static inline void multMatrixTransposed ( const Dune::FieldMatrix< K, m, n > &A,
+                                              const Dune::FieldMatrix< K, n, p > &B,
+                                              Dune::FieldMatrix< K, m, p > &ret )
+    {
+        typedef typename Dune::FieldMatrix< K, m, p > :: size_type size_type;
+
+        for( size_type i = 0; i < m; ++i )
+        {
+            for( size_type j = 0; j < p; ++j )
+            {
+                ret[ i ][ j ] = K( 0 );
+                for( size_type k = 0; k < n; ++k )
+                    ret[ i ][ j ] += A[ k ][ i ] * B[ k ][ j ];
+            }
+        }
+    }
+}
 
 enum WellVariablePositions {
     XvarWell = 0,
@@ -338,6 +359,12 @@ enum WellVariablePositions {
                                      const WellState& well_state,
                                      const int well_number) const;
 
+            /// \! brief Modifies matrix to include influences of the well perforations.
+            ///
+            /// \param mat The linear system with the assembled mass balance
+            ///            equations
+            void addWellContributions(Mat& mat) const;
+
             using WellMapType = typename WellState::WellMapType;
             using WellMapEntryType = typename WellState::mapentry_t;
 
diff --git a/opm/autodiff/StandardWellsDense_impl.hpp b/opm/autodiff/StandardWellsDense_impl.hpp
index d1af3d5dd..e29964b05 100644
--- a/opm/autodiff/StandardWellsDense_impl.hpp
+++ b/opm/autodiff/StandardWellsDense_impl.hpp
@@ -116,6 +116,36 @@ namespace Opm {
     }
 
 
+    template<class TypeTag>
+    void
+    StandardWellsDense<TypeTag>::
+    addWellContributions(Mat& mat) const
+    {
+        const int nw = wells().number_of_wells;
+
+        for(int w = 0; w < nw; ++w)
+        {
+            for(int perf1 = wells().well_connpos[w] ; perf1 < wells().well_connpos[w+1]; ++perf1) {
+                const auto pi = wells().well_cells[perf1];
+                auto& row = mat[pi];
+
+                for(int perf2 = wells().well_connpos[w] ; perf2 < wells().well_connpos[w+1]; ++perf2) {
+                    const auto pj = wells().well_cells[perf2];
+                    auto col = row.find(pj);
+                    assert(col != row.end());
+                    auto colC = duneC_[w].find(pj);
+                    assert(colC !=  duneC_[w].end());
+                    auto colB = duneB_[w].find(pi);
+                    assert(colB !=  duneB_[w].end());
+
+                    typename Mat::block_type tmp, tmp1;
+                    Dune::FMatrixHelp::multMatrix(invDuneD_[w][w],  (*colC), tmp);
+                    Detail::multMatrixTransposed((*colB), tmp, tmp1);
+                    (*col) -= tmp1;
+                }
+            }
+        }
+    }
 
     template<typename TypeTag>
     SimulatorReport
@@ -243,6 +273,11 @@ namespace Opm {
 
         // do the local inversion of D.
         localInvert( invDuneD_ );
+
+        if ( param_.matrix_add_well_contributions_ )
+        {
+            addWellContributions( ebosJac );
+        }
     }