diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp
index 38a83be46..269159240 100644
--- a/opm/simulators/wells/BlackoilWellModel_impl.hpp
+++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp
@@ -1972,7 +1972,7 @@ namespace Opm {
 
             auto * perf_pressure = well_state.perfPress().data() + wm.second[1];
             auto * perf_rates = well_state.perfRates().data() + wm.second[1];
-            auto * perf_phase_rates = well_state.perfPhaseRates().data() + wm.second[1]*np;
+            auto * perf_phase_rates = well_state.mutable_perfPhaseRates().data() + wm.second[1]*np;
             const auto& perf_data = this->well_perf_data_[well_index];
 
             for (std::size_t perf_index = 0; perf_index < perf_data.size(); perf_index++) {
diff --git a/opm/simulators/wells/MultisegmentWell_impl.hpp b/opm/simulators/wells/MultisegmentWell_impl.hpp
index d8a124ebf..7bd59cfb2 100644
--- a/opm/simulators/wells/MultisegmentWell_impl.hpp
+++ b/opm/simulators/wells/MultisegmentWell_impl.hpp
@@ -2664,7 +2664,7 @@ namespace Opm
                 // store the perf pressure and rates
                 const int rate_start_offset = (first_perf_ + perf) * number_of_phases_;
                 for (int comp_idx = 0; comp_idx < num_components_; ++comp_idx) {
-                    well_state.perfPhaseRates()[rate_start_offset + ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
+                    well_state.mutable_perfPhaseRates()[rate_start_offset + ebosCompIdxToFlowCompIdx(comp_idx)] = cq_s[comp_idx].value();
                 }
                 well_state.perfPress()[first_perf_ + perf] = perf_press.value();
 
diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp
index 1eae21ac7..36147476d 100644
--- a/opm/simulators/wells/StandardWell_impl.hpp
+++ b/opm/simulators/wells/StandardWell_impl.hpp
@@ -620,7 +620,7 @@ namespace Opm
                 if (has_solvent && componentIdx == contiSolventEqIdx) {
                     well_state.perfRateSolvent()[first_perf_ + perf] = cq_s[componentIdx].value();
                 } else {
-                    well_state.perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
+                    well_state.mutable_perfPhaseRates()[(first_perf_ + perf) * np + ebosCompIdxToFlowCompIdx(componentIdx)] = cq_s[componentIdx].value();
                 }
             }
 
diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
index 34faff6dd..5dcd5667a 100644
--- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
+++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp
@@ -293,13 +293,13 @@ namespace Opm
                             for (int perf_phase_idx = connpos*np;
                                  perf_phase_idx < (connpos + num_perf_this_well)*np; ++perf_phase_idx, ++old_perf_phase_idx )
                             {
-                                perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
+                                mutable_perfPhaseRates()[ perf_phase_idx ] = prevState->perfPhaseRates()[ old_perf_phase_idx ];
                             }
                         } else {
                             const int global_num_perf_this_well = parallel_well_info[w]->communication().sum(num_perf_this_well);
                             for (int perf = connpos; perf < connpos + num_perf_this_well; ++perf) {
                                 for (int p = 0; p < np; ++p) {
-                                    perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(global_num_perf_this_well);
+                                    mutable_perfPhaseRates()[np*perf + p] = wellRates()[np*newIndex + p] / double(global_num_perf_this_well);
                                 }
                             }
                         }
@@ -409,7 +409,7 @@ namespace Opm
         }
 
         /// One rate per phase and well connection.
-        std::vector<double>& perfPhaseRates() { return perfphaserates_; }
+        std::vector<double>& mutable_perfPhaseRates() { return perfphaserates_; }
         const std::vector<double>& perfPhaseRates() const { return perfphaserates_; }
 
         /// One current control per injecting well.
@@ -690,7 +690,7 @@ namespace Opm
                             // maybe the best way is to initialize the fractions first then get the rates
                             for (int perf = 0; perf < n_activeperf; perf++) {
                                 const int perf_pos = start_perf + perf;
-                                perfPhaseRates()[np * perf_pos + gaspos] *= 100.;
+                                mutable_perfPhaseRates()[np * perf_pos + gaspos] *= 100.;
                             }
                         }