diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp
index 64321231b..6470f8286 100644
--- a/opm/simulators/wells/BlackoilWellModel.hpp
+++ b/opm/simulators/wells/BlackoilWellModel.hpp
@@ -283,7 +283,12 @@ template<class Scalar> class WellContributions;
 
                 this->assignWellTargets(wsrpt);
                 this->assignShutConnections(wsrpt, this->reportStepIndex());
-
+                // only used to compute gas injection mass rates for CO2STORE runs
+                // The gas reference density is thus the same for all pvt regions
+                // We therefore for simplicity use 0 here 
+                if (eclState().runspec().co2Storage()) {
+                    this->assignMassGasRate(wsrpt, FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, 0));
+                }
                 return wsrpt;
             }
 
diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
index 32b9ea5ef..d9b030724 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp
@@ -1678,6 +1678,20 @@ getWellsForTesting(const int timeStepIdx,
       return {};
 }
 
+template<class Scalar>
+void BlackoilWellModelGeneric<Scalar>::
+assignMassGasRate(data::Wells& wsrpt,
+                  const Scalar& gasDensity) const
+{
+    using rt = data::Rates::opt;
+    for (auto& wrpt : wsrpt) {
+        auto& well_rates = wrpt.second.rates;
+        const auto w_mass_rate = well_rates.get(rt::gas) * gasDensity;
+        well_rates.set(rt::mass_gas, w_mass_rate);
+    }
+}
+
+
 template<class Scalar>
 void BlackoilWellModelGeneric<Scalar>::
 assignWellTracerRates(data::Wells& wsrpt,
diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.hpp b/opm/simulators/wells/BlackoilWellModelGeneric.hpp
index 98f512622..ab518b2d5 100644
--- a/opm/simulators/wells/BlackoilWellModelGeneric.hpp
+++ b/opm/simulators/wells/BlackoilWellModelGeneric.hpp
@@ -440,6 +440,8 @@ protected:
     using MswTracerRates = std::map<std::tuple<std::string, std::string, std::size_t>, Scalar>;
     void assignMswTracerRates(data::Wells& wsrpt,
                               const MswTracerRates& mswTracerRates) const;
+    void assignMassGasRate(data::Wells& wsrpt,
+                           const Scalar& gasDensity) const;
 
     Schedule& schedule_;
     const SummaryState& summaryState_;