diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 971e803d33..b2bf94bf3f 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -317,6 +317,60 @@ namespace Opm {
                 }
             }
 
+            std::vector<std::vector<int>> getMaxWellConnections() const
+            {
+                std::vector<std::vector<int>> wells;
+                // Create cartesian to compressed mapping
+                const auto& globalCell = grid().globalCell();
+                const auto& cartesianSize = grid().logicalCartesianSize();
+
+                auto size = cartesianSize[0]*cartesianSize[1]*cartesianSize[2];
+
+                std::vector<int> cartesianToCompressed(size, -1);
+                auto begin = globalCell.begin();
+
+                for ( auto cell = begin, end= globalCell.end(); cell != end; ++cell )
+                {
+                    cartesianToCompressed[ *cell ] = cell - begin;
+                }
+
+                auto schedule_wells = schedule().getWellsatEnd();
+                schedule_wells.erase(std::remove_if(schedule_wells.begin(), schedule_wells.end(), not_on_process_), schedule_wells.end());
+                wells.reserve(schedule_wells.size());
+
+                // initialize the additional cell connections introduced by wells.
+                for ( const auto& well : schedule_wells )
+                {
+                    std::vector<int> compressed_well_perforations;
+                    // All possible completions of the well
+                    const auto& completionSet = well.getConnections();
+                    compressed_well_perforations.reserve(completionSet.size());
+
+                    for ( size_t c=0; c < completionSet.size(); c++ )
+                    {
+                        const auto& completion = completionSet.get(c);
+                        int i = completion.getI();
+                        int j = completion.getJ();
+                        int k = completion.getK();
+                        int cart_grid_idx = i + cartesianSize[0]*(j + cartesianSize[1]*k);
+                        int compressed_idx = cartesianToCompressed[cart_grid_idx];
+
+                        if ( compressed_idx >= 0 ) // Ignore completions in inactive/remote cells.
+                        {
+                            compressed_well_perforations.push_back(compressed_idx);
+                        }
+                    }
+
+                    if ( ! compressed_well_perforations.empty() )
+                    {
+                        std::sort(compressed_well_perforations.begin(),
+                                  compressed_well_perforations.end());
+
+                        wells.push_back(compressed_well_perforations);
+                    }
+                }
+                return wells;
+            }
             
             
             void addWellPressureEquationsStruct(PressureMatrix& jacobian) const
@@ -327,11 +381,21 @@ namespace Opm {
                     int wdof = rdofs + i; 
                     jacobian.entry(wdof,wdof) = 1.0;// better scaling ?
                 }
+                std::vector<std::vector<int>> wellconnections = getMaxWellConnections();
+                for(int i=0; i < nw; i++){
+                    const auto& perfcells = wellconnections[i];
+                    for(int perfcell : perfcells){
+                        int wdof = rdofs + i; 
+                        jacobian.entry(wdof,perfcell) = 0.0;
+                        jacobian.entry(perfcell, wdof) = 0.0;
+                    }
+                }
                 for (const auto& well : well_container_) {
                     well->addWellPressureEquationsStruct(jacobian);
                 }
             }
 
+
             void initGliftEclWellMap(GLiftEclWells &ecl_well_map);
 
             /// \brief Get list of local nonshut wells