diff --git a/opm/simulators/wells/StandardWellAssemble.cpp b/opm/simulators/wells/StandardWellAssemble.cpp
index 4aa88ef2c..4729657ba 100644
--- a/opm/simulators/wells/StandardWellAssemble.cpp
+++ b/opm/simulators/wells/StandardWellAssemble.cpp
@@ -140,6 +140,30 @@ assembleControlEq(const WellState& well_state,
     }
 }
 
+template<class FluidSystem, class Indices, class Scalar>
+template<class EvalWell>
+void StandardWellAssemble<FluidSystem,Indices,Scalar>::
+assembleInjectivityEq(const EvalWell& eq_pskin,
+                      const EvalWell& eq_wat_vel,
+                      const int pskin_index,
+                      const int wat_vel_index,
+                      const int cell_idx,
+                      const int numWellEq,
+                      StandardWellEquations<Scalar,Indices::numEq>& eqns) const
+{
+    eqns.resWell_[0][pskin_index] = eq_pskin.value();
+    eqns.resWell_[0][wat_vel_index] = eq_wat_vel.value();
+    for (int pvIdx = 0; pvIdx < numWellEq; ++pvIdx) {
+        eqns.duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
+        eqns.duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
+    }
+
+         // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
+    for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
+        eqns.duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
+    }
+}
+
 #define INSTANCE(Dim,...) \
 template class StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>; \
 template void \
@@ -155,7 +179,16 @@ StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA
                       const double, \
                       const int, \
                       StandardWellEquations<double,__VA_ARGS__::numEq>&, \
-                      DeferredLogger&) const;
+                      DeferredLogger&) const; \
+template void \
+StandardWellAssemble<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>,__VA_ARGS__,double>:: \
+    assembleInjectivityEq(const DenseAd::Evaluation<double,-1,Dim>&, \
+                          const DenseAd::Evaluation<double,-1,Dim>&, \
+                          const int, \
+                          const int, \
+                          const int, \
+                          const int, \
+                          StandardWellEquations<double,__VA_ARGS__::numEq>&) const;
 
 // One phase
 INSTANCE(4u, BlackOilOnePhaseIndices<0u,0u,0u,0u,false,false,0u,1u,0u>)
diff --git a/opm/simulators/wells/StandardWellAssemble.hpp b/opm/simulators/wells/StandardWellAssemble.hpp
index 1281add71..e92565630 100644
--- a/opm/simulators/wells/StandardWellAssemble.hpp
+++ b/opm/simulators/wells/StandardWellAssemble.hpp
@@ -61,6 +61,15 @@ public:
                            StandardWellEquations<Scalar,Indices::numEq>& eqns,
                            DeferredLogger& deferred_logger) const;
 
+    //! \brief Assemble injectivity equation.
+    template<class EvalWell>
+    void assembleInjectivityEq(const EvalWell& eq_pskin,
+                               const EvalWell& eq_wat_vel,
+                               const int pskin_index,
+                               const int wat_vel_index,
+                               const int cell_idx,
+                               const int numWellEq,
+                               StandardWellEquations<Scalar,Indices::numEq>& eqns) const;
 
 private:
     const WellInterfaceFluidSystem<FluidSystem>& well_; //!< Reference to well
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index ac6e8d9e8..9b1d24fa9 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -2325,7 +2325,6 @@ namespace Opm
 
         // equation for the water velocity
         const EvalWell eq_wat_vel = this->primary_variables_evaluation_[wat_vel_index] - water_velocity;
-        this->linSys_.resWell_[0][wat_vel_index] = eq_wat_vel.value();
 
         const auto& ws = well_state.well(this->index_of_well_);
         const auto& perf_data = ws.perf_data;
@@ -2340,16 +2339,14 @@ namespace Opm
         const EvalWell eq_pskin = this->primary_variables_evaluation_[pskin_index]
                                   - pskin(throughput, this->primary_variables_evaluation_[wat_vel_index], poly_conc, deferred_logger);
 
-        this->linSys_.resWell_[0][pskin_index] = eq_pskin.value();
-        for (int pvIdx = 0; pvIdx < this->numWellEq_; ++pvIdx) {
-            this->linSys_.duneD_[0][0][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx+Indices::numEq);
-            this->linSys_.duneD_[0][0][pskin_index][pvIdx] = eq_pskin.derivative(pvIdx+Indices::numEq);
-        }
-
-        // the water velocity is impacted by the reservoir primary varaibles. It needs to enter matrix B
-        for (int pvIdx = 0; pvIdx < Indices::numEq; ++pvIdx) {
-            this->linSys_.duneB_[0][cell_idx][wat_vel_index][pvIdx] = eq_wat_vel.derivative(pvIdx);
-        }
+        StandardWellAssemble<FluidSystem,Indices,Scalar>(*this).
+                assembleInjectivityEq(eq_pskin,
+                                      eq_wat_vel,
+                                      pskin_index,
+                                      wat_vel_index,
+                                      cell_idx,
+                                      this->numWellEq_,
+                                      this->linSys_);
     }