Extending the SUMMARY keywords for CO2STORE and H2STORE

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -2070,6 +2070,368 @@ private:
                               }
                   }
             },
+            Entry{ScalarEntry{"BGMIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  Scalar result = getValue(ectx.fs.invB(gasPhaseIdx)) *
+                                                  getValue(ectx.fs.saturation(gasPhaseIdx));
+
+                                  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                                       result += getValue(ectx.fs.Rs()) *
+                                                 getValue(ectx.fs.invB(oilPhaseIdx)) *
+                                                 getValue(ectx.fs.saturation(oilPhaseIdx));
+                                  }
+                                  else {
+                                      result += getValue(ectx.fs.Rsw()) *
+                                                getValue(ectx.fs.invB(waterPhaseIdx)) *
+                                                getValue(ectx.fs.saturation(waterPhaseIdx));
+                                  }
+                                  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
+                                                                                    ectx.intQuants.pvtRegionIndex());
+                                  return result *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         rhoG;
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMGP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
+                                                                                    ectx.intQuants.pvtRegionIndex());
+                                  return getValue(ectx.fs.invB(gasPhaseIdx)) *
+                                         getValue(ectx.fs.saturation(gasPhaseIdx)) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         rhoG;
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMDS",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  Scalar result;
+                                  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                                      result  = getValue(ectx.fs.Rs()) *
+                                                getValue(ectx.fs.invB(oilPhaseIdx)) *
+                                                getValue(ectx.fs.saturation(oilPhaseIdx));
+                                  }
+                                  else {
+                                      result = getValue(ectx.fs.Rsw()) *
+                                               getValue(ectx.fs.invB(waterPhaseIdx)) *
+                                               getValue(ectx.fs.saturation(waterPhaseIdx));
+                                  }
+                                  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
+                                                                                    ectx.intQuants.pvtRegionIndex());
+                                  return result *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         rhoG;
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMST",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
+                                  Scalar strandedGas = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
+                                      strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::min(strandedGas, sg);
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMUS",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));
+                                  Scalar strandedGas = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      const Scalar krg = getValue(ectx.intQuants.relativePermeability(gasPhaseIdx));
+                                      strandedGas = MaterialLaw::strandedGasSaturation(matParams, sg, krg);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::max(Scalar{0.0}, sg - strandedGas);
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMTR",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  Scalar trappedGas = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::min(trappedGas, getValue(ectx.fs.saturation(gasPhaseIdx)));
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGMMO",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  Scalar trappedGas = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      trappedGas = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/true);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - trappedGas);
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGKTR",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));                                
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  if (sg > sgcr) {
+                                      return 0.0;
+                                  }
+                                  else {
+                                      const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                        FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                      : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                      return (1.0 - xgW) *
+                                             model.dofTotalVolume(ectx.globalDofIdx) *
+                                             getValue(ectx.intQuants.porosity()) * 
+                                             getValue(ectx.fs.density(gasPhaseIdx)) *
+                                             getValue(ectx.fs.saturation(gasPhaseIdx));
+                                  }
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGKMO",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));                                
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  if (sgcr >= sg) {
+                                      return 0.0;
+                                  }
+                                  else {
+                                      const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                        FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                      : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                      return (1.0 - xgW) *
+                                             model.dofTotalVolume(ectx.globalDofIdx) *
+                                             getValue(ectx.intQuants.porosity()) * 
+                                             getValue(ectx.fs.density(gasPhaseIdx)) *
+                                             getValue(ectx.fs.saturation(gasPhaseIdx));
+                                  }
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGCDI",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::min(sgcr, getValue(ectx.fs.saturation(gasPhaseIdx))) /
+                                         FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGCDM",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                    FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                  : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                  return (1.0 - xgW) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) * 
+                                         getValue(ectx.fs.density(gasPhaseIdx)) *
+                                         std::max(Scalar{0.0}, getValue(ectx.fs.saturation(gasPhaseIdx)) - sgcr) /
+                                         FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGKDI",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));                                
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  if (sg > sgcr) {
+                                      return 0.0;
+                                  }
+                                  else {
+                                      const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                        FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                      : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                      return (1.0 - xgW) *
+                                             model.dofTotalVolume(ectx.globalDofIdx) *
+                                             getValue(ectx.intQuants.porosity()) * 
+                                             getValue(ectx.fs.density(gasPhaseIdx)) *
+                                             getValue(ectx.fs.saturation(gasPhaseIdx)) /
+                                             FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
+                                  }
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BGKDM",
+                              [&model = this->simulator_.model(),
+                               &problem = this->simulator_.problem()](const Context& ectx)
+                              {
+                                  const auto& scaledDrainageInfo = problem.materialLawManager()
+                                                                   ->oilWaterScaledEpsInfoDrainage(ectx.dofIdx);
+                                  const Scalar sg = getValue(ectx.fs.saturation(gasPhaseIdx));                                
+                                  Scalar sgcr = scaledDrainageInfo.Sgcr;
+                                  if (problem.materialLawManager()->enableHysteresis()) {
+                                      const auto& matParams = problem.materialLawParams(ectx.dofIdx);
+                                      sgcr = MaterialLaw::trappedGasSaturation(matParams, /*maxTrapping*/false);
+                                  }
+                                  if (sgcr >= sg) {
+                                      return 0.0;
+                                  }
+                                  else {
+                                      const Scalar xgW = FluidSystem::phaseIsActive(waterPhaseIdx) ?
+                                        FluidSystem::convertRvwToXgW(getValue(ectx.fs.Rvw()), ectx.intQuants.pvtRegionIndex())
+                                      : FluidSystem::convertRvToXgO(getValue(ectx.fs.Rv()), ectx.intQuants.pvtRegionIndex());
+                                      return (1.0 - xgW) *
+                                             model.dofTotalVolume(ectx.globalDofIdx) *
+                                             getValue(ectx.intQuants.porosity()) * 
+                                             getValue(ectx.fs.density(gasPhaseIdx)) *
+                                             getValue(ectx.fs.saturation(gasPhaseIdx)) /
+                                             FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
+                                  }
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BWCD",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  Scalar result;
+                                  if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                                      result  = getValue(ectx.fs.Rs()) *
+                                                getValue(ectx.fs.invB(oilPhaseIdx)) *
+                                                getValue(ectx.fs.saturation(oilPhaseIdx));
+                                  }
+                                  else {
+                                      result = getValue(ectx.fs.Rsw()) *
+                                               getValue(ectx.fs.invB(waterPhaseIdx)) *
+                                               getValue(ectx.fs.saturation(waterPhaseIdx));
+                                  }
+                                  const Scalar rhoG = FluidSystem::referenceDensity(gasPhaseIdx,
+                                                                                    ectx.intQuants.pvtRegionIndex());
+                                  return result *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         rhoG /
+                                         FluidSystem::molarMass(gasCompIdx, ectx.intQuants.pvtRegionIndex());
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BWIPG",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  Scalar result = 0.0;
+                                  if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+                                      result  = getValue(ectx.fs.Rvw()) *
+                                                getValue(ectx.fs.invB(gasPhaseIdx)) *
+                                                getValue(ectx.fs.saturation(gasPhaseIdx));
+                                  }
+                                  return result *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity());
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BWIPL",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return getValue(ectx.fs.invB(waterPhaseIdx)) *
+                                         getValue(ectx.fs.saturation(waterPhaseIdx)) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.porosity());
+                              }
+                  }
+            },
         };
 
         this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers);

opm/simulators/utils/BlackoilPhases.cpp

@@ -48,12 +48,13 @@ Opm::PhaseUsage::PhaseUsage(std::vector<BlackoilPhases::PhaseIndex> phases)
     if (this->phase_used[BlackoilPhases::Vapour])
         this->num_phases++;
 
-    this->has_solvent   = this->phase_used[BlackoilPhases::Solvent];
-    this->has_polymer   = this->phase_used[BlackoilPhases::Polymer];
-    this->has_energy    = this->phase_used[BlackoilPhases::Energy];
-    this->has_polymermw = this->phase_used[BlackoilPhases::PolymerMW];
-    this->has_foam      = this->phase_used[BlackoilPhases::Foam];
-    this->has_brine     = this->phase_used[BlackoilPhases::Brine];
-    this->has_zFraction = this->phase_used[BlackoilPhases::ZFraction];
-    this->has_micp      = false;
+    this->has_solvent        = this->phase_used[BlackoilPhases::Solvent];
+    this->has_polymer        = this->phase_used[BlackoilPhases::Polymer];
+    this->has_energy         = this->phase_used[BlackoilPhases::Energy];
+    this->has_polymermw      = this->phase_used[BlackoilPhases::PolymerMW];
+    this->has_foam           = this->phase_used[BlackoilPhases::Foam];
+    this->has_brine          = this->phase_used[BlackoilPhases::Brine];
+    this->has_zFraction      = this->phase_used[BlackoilPhases::ZFraction];
+    this->has_micp           = false;
+    this->has_co2_or_h2store = false;
 }

opm/simulators/utils/BlackoilPhases.hpp

@@ -61,6 +61,7 @@ namespace Opm
         bool has_brine{};
         bool has_zFraction{};
         bool has_micp{};
+        bool has_co2_or_h2store{};
 
     };
 

opm/simulators/utils/phaseUsageFromDeck.cpp

@@ -137,8 +137,12 @@ PhaseUsage phaseUsage(const Phases& phases)
 PhaseUsage phaseUsageFromDeck(const EclipseState& eclipseState)
 {
     const auto& phases = eclipseState.runspec().phases();
-
-    if (eclipseState.runspec().micp()) {
+    if (eclipseState.runspec().co2Storage() || eclipseState.runspec().h2Storage()) {
+        PhaseUsage pu = phaseUsage(phases);
+        pu.has_co2_or_h2store = true;
+        return pu;
+    }
+    else if (eclipseState.runspec().micp()) {
         PhaseUsage pu = phaseUsage(phases);
         pu.has_micp = true;
         return pu;

opm/simulators/wells/BlackoilWellModel.hpp

@@ -225,10 +225,10 @@ template<class Scalar> class WellContributions;
 
                 this->assignWellTargets(wsrpt);
                 this->assignShutConnections(wsrpt, this->reportStepIndex());
-                // only used to compute gas injection mass rates for CO2STORE runs
+                // only used to compute gas injection mass rates for CO2STORE/H2STORE runs
                 // The gas reference density is thus the same for all pvt regions
                 // We therefore for simplicity use 0 here 
-                if (eclState().runspec().co2Storage()) {
+                if (eclState().runspec().co2Storage() || eclState().runspec().h2Storage()) {
                     this->assignMassGasRate(wsrpt, FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, 0));
                 }
                 return wsrpt;

opm/simulators/wells/PerfData.cpp

@@ -51,6 +51,7 @@ PerfData<Scalar>::PerfData(const std::size_t num_perf,
     , connection_compaction_tmult(num_perf)
     , satnum_id(num_perf)
     , ecl_index(num_perf)
+    , gas_mass_rates(num_perf)
 {
     if (injector) {
         prepareInjectorContainers();
@@ -95,6 +96,7 @@ PerfData<Scalar> PerfData<Scalar>::serializationTestObject()
     result.water_velocity = {29.0, 30.0};
     result.filtrate_data = ConnFiltrateData<Scalar>::serializationTestObject();
     result.connFracStatistics.assign(3, ConnFracStatistics<Scalar>::serializationTestObject());
+    result.gas_mass_rates = {31.0};
 
     return result;
 }
@@ -139,6 +141,7 @@ bool PerfData<Scalar>::try_assign(const PerfData& other)
     this->water_velocity = other.water_velocity;
     this->filtrate_data = other.filtrate_data;
     this->connFracStatistics = other.connFracStatistics;
+    this->gas_mass_rates = other.gas_mass_rates;
 
     return true;
 }
@@ -170,6 +173,7 @@ bool PerfData<Scalar>::operator==(const PerfData& rhs) const
         && (this->water_velocity == rhs.water_velocity)
         && (this->filtrate_data == rhs.filtrate_data)
         && (this->connFracStatistics == rhs.connFracStatistics)
+        && (this->gas_mass_rates == rhs.gas_mass_rates)
         ;
 }
 

opm/simulators/wells/PerfData.hpp

@@ -80,6 +80,7 @@ public:
         serializer(water_velocity);
         serializer(filtrate_data);
         serializer(connFracStatistics);
+        serializer(gas_mass_rates);
     }
 
     bool operator==(const PerfData&) const;
@@ -108,6 +109,7 @@ public:
     std::vector<Scalar> connection_compaction_tmult{};
     std::vector<int> satnum_id{};
     std::vector<std::size_t> ecl_index{};
+    std::vector<Scalar> gas_mass_rates{};
 
     // The water_throughput, skin_pressure and water_velocity variables are
     // only used for injectors to check the injectivity.

opm/simulators/wells/StandardWell_impl.hpp

@@ -634,6 +634,14 @@ namespace Opm
 
         // Store the perforation pressure for later usage.
         perf_data.pressure[perf] = ws.bhp + this->connections_.pressure_diff(perf);
+
+        // Store the perforation gass mass rate.
+        const auto& pu = well_state.phaseUsage();
+        if (pu.has_co2_or_h2store) {
+            const unsigned gas_comp_idx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+            const Scalar rho = FluidSystem::referenceDensity( FluidSystem::gasPhaseIdx, Base::pvtRegionIdx() );
+            perf_data.gas_mass_rates[perf] = cq_s[gas_comp_idx].value() * rho;
+        }
     }
 
 

opm/simulators/wells/WellState.cpp

@@ -1137,6 +1137,12 @@ reportConnectionPressuresAndRates(const std::size_t well_index,
             connections[i].rates.set(rt::urea, perf_data.urea_rates[i]);
         }
     }
+    
+    if (pu.has_co2_or_h2store) {
+        for (auto i = 0*num_perf_well; i < num_perf_well; ++i) {
+            connections[i].rates.set(rt::mass_gas, perf_data.gas_mass_rates[i]);
+        }
+    }
 }
 
 template <class Scalar>