diff --git a/opm/autodiff/BlackoilModelEbos.hpp b/opm/autodiff/BlackoilModelEbos.hpp
index 9c4558fa1..92f562672 100644
--- a/opm/autodiff/BlackoilModelEbos.hpp
+++ b/opm/autodiff/BlackoilModelEbos.hpp
@@ -246,7 +246,7 @@ namespace Opm {
             report.total_linearizations = 1;
 
             try {
-                report += assemble(timer, iteration);
+                report += assembleReservoir(timer, iteration);
                 report.assemble_time += perfTimer.stop();
             }
             catch (...) {
@@ -293,6 +293,11 @@ namespace Opm {
                 const int nc = UgGridHelpers::numCells(grid_);
                 BVector x(nc);
 
+                // apply the Schur compliment of the well model to the reservoir linearized
+                // equations
+                wellModel().linearize(ebosSimulator().model().linearizer().jacobian(),
+                                      ebosSimulator().model().linearizer().residual());
+
                 try {
                     solveJacobianSystem(x);
                     report.linear_solve_time += perfTimer.stop();
@@ -360,13 +365,13 @@ namespace Opm {
         /// \param[in]      reservoir_state   reservoir state variables
         /// \param[in, out] well_state        well state variables
         /// \param[in]      initial_assembly  pass true if this is the first call to assemble() in this timestep
-        SimulatorReport assemble(const SimulatorTimerInterface& timer,
-                                 const int iterationIdx)
+        SimulatorReport assembleReservoir(const SimulatorTimerInterface& timer,
+                                          const int iterationIdx)
         {
             // -------- Mass balance equations --------
             ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx);
             ebosSimulator_.problem().beginIteration();
-            ebosSimulator_.model().linearizer().linearize();
+            ebosSimulator_.model().linearizer().linearizeDomain();
             ebosSimulator_.problem().endIteration();
 
             return wellModel().lastReport();
@@ -463,19 +468,6 @@ namespace Opm {
 
             auto& ebosJac = ebosSimulator_.model().linearizer().jacobian();
             auto& ebosResid = ebosSimulator_.model().linearizer().residual();
-            // J = [A, B; C, D], where A is the reservoir equations, B and C the interaction of well
-            // with the reservoir and D is the wells itself.
-            // The full system is reduced to a number of cells X number of cells system via Schur complement
-            // A -= B^T D^-1 C
-            // If matrix_add_well_contribution is false, the Ax operator is modified. i.e Ax -= B^T D^-1 C x in the WellModelMatrixAdapter
-            // instead of A.
-            // The residual is modified similarly.
-            // r = [r, r_well], where r is the residual and r_well the well residual.
-            // r -= B^T * D^-1 r_well
-            wellModel().apply(ebosResid);
-            if (param_.matrix_add_well_contributions_) {
-                wellModel().addWellContributions(ebosJac.istlMatrix());
-            }
 
             // set initial guess
             x = 0.0;
diff --git a/opm/autodiff/BlackoilWellModel_impl.hpp b/opm/autodiff/BlackoilWellModel_impl.hpp
index fff1a6218..aca7821c2 100644
--- a/opm/autodiff/BlackoilWellModel_impl.hpp
+++ b/opm/autodiff/BlackoilWellModel_impl.hpp
@@ -101,18 +101,23 @@ namespace Opm {
         if (!localWellsActive())
             return;
 
-        // we don't what to add the schur complement
-        // here since it affects the getConvergence method
-        /*
+        if (!param_.matrix_add_well_contributions_) {
+            // if the well contributions are not supposed to be included explicitly in
+            // the matrix, we only apply the vector part of the Schur complement here.
+            for (const auto& well: well_container_) {
+                // r = r - duneC_^T * invDuneD_ * resWell_
+                well->apply(res);
+            }
+            return;
+        }
+
         for (const auto& well: well_container_) {
-            if (param_.matrix_add_well_contributions_)
-                well->addWellContributions(mat);
+            well->addWellContributions(mat.istlMatrix());
 
             // applying the well residual to reservoir residuals
             // r = r - duneC_^T * invDuneD_ * resWell_
             well->apply(res);
         }
-        */
     }