Merge pull request #6105 from daavid00/micp

Updating MICP

flow/flow_micp.cpp

@@ -18,12 +18,15 @@
 
 #include <flow/flow_micp.hpp>
 
+#include <opm/grid/CpGrid.hpp>
+
 #include <opm/material/common/ResetLocale.hpp>
 #include <opm/models/blackoil/blackoilonephaseindices.hh>
+#include <opm/models/blackoil/blackoillocalresidualtpfa.hh>
+#include <opm/models/discretization/common/tpfalinearizer.hh>
 
-#include <opm/grid/CpGrid.hpp>
-#include <opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp>
 #include <opm/simulators/flow/Main.hpp>
+#include <opm/simulators/flow/SimulatorFullyImplicitBlackoil.hpp>
 
 namespace Opm {
 namespace Properties {
@@ -36,6 +39,10 @@ template<class TypeTag>
 struct EnableMICP<TypeTag, TTag::FlowMICPProblem> {
     static constexpr bool value = true;
 };
+template<class TypeTag>
+struct Linearizer<TypeTag, TTag::FlowMICPProblem> { using type = TpfaLinearizer<TypeTag>; };
+template<class TypeTag>
+struct LocalResidual<TypeTag, TTag::FlowMICPProblem> { using type = BlackOilLocalResidualTPFA<TypeTag>; };
 //! The indices required by the model
 template<class TypeTag>
 struct Indices<TypeTag, TTag::FlowMICPProblem>
@@ -58,6 +65,13 @@ public:
                                          /*enabledCompIdx=*/FluidSystem::waterCompIdx,
                                          5>; //Five MICP components
 };
+
+template<class TypeTag>
+struct EnableDiffusion<TypeTag, TTag::FlowMICPProblem> { static constexpr bool value = true; };
+
+template<class TypeTag>
+struct EnableDispersion<TypeTag, TTag::FlowMICPProblem> { static constexpr bool value = true; };
+
 }}
 
 namespace Opm {

flowexperimental/BlackOilIntensiveQuantitiesGlobalIndex.hpp

@@ -54,6 +54,7 @@
 
 #include <opm/utility/CopyablePtr.hpp>
 
+#include <stdexcept>
 #include <utility>
 
 #include <fmt/format.h>
@@ -519,6 +520,11 @@ public:
     Scalar referencePorosity() const
     { return referencePorosity_; }
 
+    const Evaluation& permFactor() const
+    {
+        throw std::logic_error("permFactor() is not yet implemented for compositional modeling"); 
+    }
+
 private:
     friend BlackOilSolventIntensiveQuantities<TypeTag>;
     friend BlackOilExtboIntensiveQuantities<TypeTag>;

opm/models/blackoil/blackoilboundaryratevector.hh

@@ -182,6 +182,8 @@ public:
             (*this)[Indices::contiMicrobialEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.microbialConcentration();
             (*this)[Indices::contiOxygenEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.oxygenConcentration();
             (*this)[Indices::contiUreaEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.ureaConcentration();
+            // since the urea concentration can be much larger than 1, then we apply a scaling factor
+            (*this)[Indices::contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
         }
 
         // make sure that the right mass conservation quantities are used

opm/models/blackoil/blackoilbrinemodules.hh

@@ -486,6 +486,6 @@ public:
 
 };
 
-} // namespace Ewoms
+} // namespace Opm
 
 #endif

opm/models/blackoil/blackoildiffusionmodule.hh

@@ -28,6 +28,7 @@
 #ifndef OPM_BLACKOIL_DIFFUSION_MODULE_HH
 #define OPM_BLACKOIL_DIFFUSION_MODULE_HH
 
+#include <opm/models/blackoil/blackoilmicpmodules.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
 
 #include <opm/material/common/Valgrind.hpp>
@@ -100,11 +101,19 @@ class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/true>
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
 
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
     enum { conti0EqIdx = Indices::conti0EqIdx };
 
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+
+    static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
+    static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
+
     using Toolbox = MathToolbox<Evaluation>;
 
 public:
@@ -157,21 +166,48 @@ public:
         if(!FluidSystem::enableDiffusion())
             return;
         const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
-        const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
-        const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
+        const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
+        const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
         const auto& diffusivity = extQuants.diffusivity();
         const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient();
-        addDiffusiveFlux(flux, fluidStateI, fluidStateJ, diffusivity, effectiveDiffusionCoefficient);
+        addDiffusiveFlux(flux, inIq, exIq, diffusivity, effectiveDiffusionCoefficient);
     }
 
-    template<class FluidState,class EvaluationArray>
+    template<class IntensiveQuantities,class EvaluationArray>
     static void addDiffusiveFlux(RateVector& flux,
-                                 const FluidState& fluidStateI,
-                                 const FluidState& fluidStateJ,
+                                 const IntensiveQuantities& inIq,
+                                 const IntensiveQuantities& exIq,
                                  const Evaluation& diffusivity,
                                  const EvaluationArray& effectiveDiffusionCoefficient)
     {
-        unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
+        const auto& inFs = inIq.fluidState();
+        const auto& exFs = exIq.fluidState();
+        Evaluation diffR = 0.0;
+        if constexpr(enableMICP) {
+            // The diffusion coefficients are given for mass concentrations
+            Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
+            diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration());
+            flux[contiMicrobialEqIdx] +=
+                bAvg
+                * diffR
+                * diffusivity
+                * effectiveDiffusionCoefficient[waterPhaseIdx][contiMicrobialEqIdx - 1];
+            diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration());
+            flux[contiOxygenEqIdx] +=
+                bAvg
+                * diffR
+                * diffusivity
+                * effectiveDiffusionCoefficient[waterPhaseIdx][contiOxygenEqIdx - 1];
+            diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration());
+            flux[contiUreaEqIdx] +=
+                bAvg
+                * diffR
+                * diffusivity
+                * effectiveDiffusionCoefficient[waterPhaseIdx][contiUreaEqIdx - 1];
+            return;
+        }
+
+        unsigned pvtRegionIndex = inFs.pvtRegionIndex();
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -188,34 +224,33 @@ public:
             }
 
             // arithmetic mean of the phase's b factor weighed by saturation
-            Evaluation bSAvg = fluidStateI.saturation(phaseIdx) * fluidStateI.invB(phaseIdx);
-            bSAvg += Toolbox::value(fluidStateJ.saturation(phaseIdx)) * Toolbox::value(fluidStateJ.invB(phaseIdx));
+            Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx);
+            bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx));
             bSAvg /= 2;
 
             // phase not present, skip
             if(bSAvg < 1.0e-6)
                 continue;
             Evaluation convFactor = 1.0;
-            Evaluation diffR = 0.0;
             if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
-                Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
+                Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
                 convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg);
-                diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
+                diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
             }
             if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
-                Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
+                Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
                 convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex));
-                diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
+                diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
             }
             if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
-                Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
+                Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
                 convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg);
-                diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
+                diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
             }
             if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
-                Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
+                Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
                 convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex));
-                diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
+                diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
             }
 
             // mass flux of solvent component (oil in oil or gas in gas)
@@ -339,8 +374,13 @@ class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
     using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
+
     enum { numPhases = FluidSystem::numPhases };
     enum { numComponents = FluidSystem::numComponents };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
 public:
     BlackOilDiffusionIntensiveQuantities() = default;
@@ -418,6 +458,14 @@ protected:
                  const IntensiveQuantities& intQuants) {
         using Toolbox = MathToolbox<Evaluation>;
 
+        if constexpr(enableMICP) {
+            unsigned pvtRegionIndex = intQuants.fluidState().pvtRegionIndex();
+            diffusionCoefficient_[waterPhaseIdx][0] = MICPModule::microbialDiffusion(pvtRegionIndex);
+            diffusionCoefficient_[waterPhaseIdx][1] = MICPModule::oxygenDiffusion(pvtRegionIndex);
+            diffusionCoefficient_[waterPhaseIdx][2] = MICPModule::ureaDiffusion(pvtRegionIndex);
+            return;
+        }
+
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -535,6 +583,9 @@ class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
     enum { dimWorld = GridView::dimensionworld };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+    
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
     using DimVector = Dune::FieldVector<Scalar, dimWorld>;
     using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
@@ -569,6 +620,17 @@ public:
                        const IntensiveQuantities& intQuantsInside,
                        const IntensiveQuantities& intQuantsOutside) {
         // opm-models expects per area flux
+        // use the arithmetic average for the effective
+        // diffusion coefficients.
+        if constexpr(enableMICP) {
+            for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
+                effectiveDiffusionCoefficient[waterPhaseIdx][compIdx] = 0.5 *
+                    ( intQuantsInside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) +
+                    intQuantsOutside.effectiveDiffusionCoefficient(waterPhaseIdx, compIdx) );
+            }
+            return;
+        }
+
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -578,8 +640,6 @@ public:
                 continue;
             }
             for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
-                // use the arithmetic average for the effective
-                // diffusion coefficients.
                 effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 *
                     ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) +
                       intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) );
@@ -608,7 +668,7 @@ public:
     { return diffusivity_; }
 
     /*!
-     * \brief The effective diffusion coeffcient of a component in a
+     * \brief The effective diffusion coefficient of a component in a
      *        fluid phase at the face's integration point
      *
      * \copydoc Doxygen::phaseIdxParam

opm/models/blackoil/blackoildispersionmodule.hh

@@ -62,6 +62,7 @@ class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/false>
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
 
     enum { numPhases = FluidSystem::numPhases };
 
@@ -85,10 +86,10 @@ public:
                                   unsigned)
     {}
 
-    template<class FluidState, class Scalar>
+    template<class IntensiveQuantities, class Scalar>
     static void addDispersiveFlux(RateVector&,
-                                  const FluidState&,
-                                  const FluidState&,
+                                  const IntensiveQuantities&,
+                                  const IntensiveQuantities&,
                                   const Evaluation&,
                                   const Scalar&)
     {}
@@ -116,6 +117,12 @@ class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
     enum { numComponents = FluidSystem::numComponents };
     enum { conti0EqIdx = Indices::conti0EqIdx };
     enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+
+    static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
+    static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
 
     using Toolbox = MathToolbox<Evaluation>;
 
@@ -148,11 +155,11 @@ public:
             return;
         }
         const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
-        const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
-        const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
+        const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
+        const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
         const auto& dispersivity = extQuants.dispersivity();
         const auto& normVelocityAvg = extQuants.normVelocityAvg();
-        addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg);
+        addDispersiveFlux(flux, inIq, exIq, dispersivity, normVelocityAvg);
     }
 
     /*!
@@ -175,14 +182,41 @@ public:
      *        \f$D_{ij}\f$ is computed in ecltransmissibility_impl.hh, using the dispersion coefficients \f$E_i\f$
      *        and \f$E_j\f$.
      */
-    template<class FluidState, class Scalar>
+    template<class IntensiveQuantities, class Scalar>
     static void addDispersiveFlux(RateVector& flux,
-                                  const FluidState& fluidStateI,
-                                  const FluidState& fluidStateJ,
+                                  const IntensiveQuantities& inIq,
+                                  const IntensiveQuantities& exIq,
                                   const Evaluation& dispersivity,
                                   const Scalar& normVelocityAvg)
     {
-        unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
+        const auto& inFs = inIq.fluidState();
+        const auto& exFs = exIq.fluidState();
+        Evaluation diffR = 0.0;
+        if constexpr(enableMICP) {
+            // The dispersion coefficients are given for mass concentrations
+            Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
+            diffR = inIq.microbialConcentration()- Toolbox::value(exIq.microbialConcentration());
+            flux[contiMicrobialEqIdx] +=
+                 bAvg
+                 * normVelocityAvg[waterPhaseIdx]
+                 * dispersivity
+                 * diffR;
+            diffR = inIq.oxygenConcentration()- Toolbox::value(exIq.oxygenConcentration());
+            flux[contiOxygenEqIdx] +=
+                 bAvg
+                 * normVelocityAvg[waterPhaseIdx]
+                 * dispersivity
+                 * diffR;
+            diffR = inIq.ureaConcentration()- Toolbox::value(exIq.ureaConcentration());
+            flux[contiUreaEqIdx] +=
+                 bAvg
+                 * normVelocityAvg[waterPhaseIdx]
+                 * dispersivity
+                 * diffR;
+            return;
+        }
+
+        unsigned pvtRegionIndex = inFs.pvtRegionIndex();
         for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
             if (!FluidSystem::phaseIsActive(phaseIdx)) {
                 continue;
@@ -193,47 +227,38 @@ public:
                 continue;
             }
 
-            // Adding dispersion in the gas phase leads to
+            // adding dispersion in the gas phase leads to
             // convergence issues and unphysical results.
-            // We disable dispersion in the gas phase for now
-            // See comment below
+            // we disable dispersion in the gas phase for now
             if (FluidSystem::gasPhaseIdx == phaseIdx) {
                 continue;
             }
 
-            // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
-            // phase check disabled due to if above, reenable when removing unconditional gas phase disablement
-            if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil())
-                /*&& FluidSystem::gasPhaseIdx == phaseIdx*/) {
-                continue;
-            }
-
             // arithmetic mean of the phase's b factor
-            Evaluation bAvg = fluidStateI.invB(phaseIdx);
-            bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx));
+            Evaluation bAvg = inFs.invB(phaseIdx);
+            bAvg += Toolbox::value(exFs.invB(phaseIdx));
             bAvg /= 2;
 
             Evaluation convFactor = 1.0;
-            Evaluation diffR = 0.0;
             if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
-                Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
+                Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
                 convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
-                diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
+                diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
             }
             if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
-                Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
+                Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
                 convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex));
-                diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
+                diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
             }
             if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
-                Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
+                Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
                 convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
-                diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
+                diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
             }
             if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
-                Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
+                Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
                 convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex));
-                diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
+                diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
             }
 
             // mass flux of solvent component
@@ -498,8 +523,10 @@ public:
             if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
                 continue;
             }
-            // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
-            if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
+            // adding dispersion in the gas phase leads to
+            // convergence issues and unphysical results.
+            // we disable dispersion in the gas phase for now
+            if (FluidSystem::gasPhaseIdx == phaseIdx) {
                 continue;
             }
             // use the arithmetic average for the effective

opm/models/blackoil/blackoilintensivequantities.hh

@@ -132,6 +132,8 @@ class BlackOilIntensiveQuantities
 
     using DirectionalMobilityPtr = Opm::Utility::CopyablePtr<DirectionalMobility<TypeTag, Evaluation>>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
+    using BrineIntQua = BlackOilBrineIntensiveQuantities<TypeTag, enableSaltPrecipitation>;
+    using MICPIntQua = BlackOilMICPIntensiveQuantities<TypeTag, enableMICP>;
 
 
 public:
@@ -509,11 +511,12 @@ public:
         // deal with water induced rock compaction
         porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
 
-        // the MICP processes change the porosity
+        // deal with MICP
         if constexpr (enableMICP){
             Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
             Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
-            porosity_ += - biofilm_ - calcite_;
+            // minimum porosity of 1e-8 to prevent numerical issues 
+            porosity_ -= min(biofilm_ + calcite_, referencePorosity_ - 1e-8);
         }
 
         // deal with salt-precipitation
@@ -704,6 +707,19 @@ public:
     Scalar referencePorosity() const
     { return referencePorosity_; }
 
+    const Evaluation& permFactor() const
+    {
+        if constexpr (enableMICP) {
+            return MICPIntQua::permFactor();
+        }
+        else if constexpr (enableSaltPrecipitation) {
+            return BrineIntQua::permFactor();
+        }
+        else {
+            throw std::logic_error("permFactor() called but salt precipitation or MICP are disabled");
+        }
+    }
+
 private:
     friend BlackOilSolventIntensiveQuantities<TypeTag>;
     friend BlackOilExtboIntensiveQuantities<TypeTag>;

opm/models/blackoil/blackoillocalresidual.hh

@@ -255,7 +255,7 @@ public:
         MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
 
         // scale the source term of the energy equation
-        if (enableEnergy)
+        if constexpr(enableEnergy)
             source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -372,7 +372,10 @@ public:
             const IntensiveQuantities& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
             unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
             // Use arithmetic average (more accurate with harmonic, but that requires recomputing the transmissbility)
-            const Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
+            Evaluation transMult = (intQuantsIn.rockCompTransMultiplier() + Toolbox::value(intQuantsEx.rockCompTransMultiplier()))/2;
+            if constexpr (enableMICP) {
+                transMult *= (intQuantsIn.permFactor() + Toolbox::value(intQuantsEx.permFactor()))/2;
+            }
             Evaluation darcyFlux;
             if (globalUpIndex == globalIndexIn) {
                     darcyFlux = pressureDifference * up.mobility(phaseIdx, facedir) * transMult * (-trans / faceArea);
@@ -396,6 +399,10 @@ public:
                     EnergyModule::template addPhaseEnthalpyFluxes_<Evaluation, Evaluation, FluidState>(
                         flux, phaseIdx, darcyFlux, up.fluidState());
                 }
+                if constexpr (enableMICP) {
+                    MICPModule::template addMICPFluxes_<Evaluation, Evaluation, IntensiveQuantities>(
+                        flux, darcyFlux, intQuantsIn);
+                }
             } else {
                 const auto& invB = getInvB_<FluidSystem, FluidState, Scalar>(up.fluidState(), phaseIdx, pvtRegionIdx);
                 const auto& surfaceVolumeFlux = invB * darcyFlux;
@@ -406,6 +413,11 @@ public:
                         addPhaseEnthalpyFluxes_<Scalar, Evaluation, FluidState>
                         (flux,phaseIdx,darcyFlux, up.fluidState());
                 }
+                if constexpr (enableMICP) {
+                    MICPModule::template
+                        addMICPFluxes_<Scalar, Evaluation, IntensiveQuantities>
+                        (flux, darcyFlux, intQuantsEx);
+                }
             }
 
         }
@@ -469,8 +481,6 @@ public:
         static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling.");
         // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
-
-
         // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
         if constexpr(enableDiffusion){
             typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
@@ -478,12 +488,13 @@ public:
             const Scalar diffusivity = nbInfo.diffusivity;
             const Scalar tmpdiffusivity = diffusivity / faceArea;
             DiffusionModule::addDiffusiveFlux(flux,
-                                              intQuantsIn.fluidState(),
-                                              intQuantsEx.fluidState(),
+                                              intQuantsIn,
+                                              intQuantsEx,
                                               tmpdiffusivity,
                                               effectiveDiffusionCoefficient);
 
         }
+
         // deal with dispersion (if present). opm-models expects per area flux (added in the tmpdispersivity).
         if constexpr(enableDispersion){
             typename DispersionModule::ExtensiveQuantities::ScalarArray normVelocityAvg;
@@ -491,16 +502,17 @@ public:
             const Scalar dispersivity = nbInfo.dispersivity;
             const Scalar tmpdispersivity = dispersivity / faceArea;
             DispersionModule::addDispersiveFlux(flux,
-                                                intQuantsIn.fluidState(),
-                                                intQuantsEx.fluidState(),
+                                                intQuantsIn,
+                                                intQuantsEx,
                                                 tmpdispersivity,
                                                 normVelocityAvg);
 
         }
-        // deal with micp (if present)
-        static_assert(!enableMICP, "Relevant computeFlux() method must be implemented for this module before enabling.");
-        // MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
-
+        
+        // apply the scaling for the urea equation in MICP
+        if constexpr (enableMICP) {
+            MICPModule::applyScaling(flux);
+        }
     }
 
     template <class BoundaryConditionData>
@@ -612,7 +624,7 @@ public:
         // conductive heat flux from boundary
         if constexpr(enableEnergy){
             Evaluation heatFlux;
-            // avoid overload of functions with same numeber of elements in eclproblem
+            // avoid overload of functions with same number of elements in eclproblem
             Scalar alpha = problem.eclTransmissibilities().thermalHalfTransBoundary(globalSpaceIdx, bdyInfo.boundaryFaceIndex);
             unsigned inIdx = 0;//dummy
             // always calculated with derivatives of this cell
@@ -627,7 +639,6 @@ public:
 
         static_assert(!enableSolvent, "Relevant treatment of boundary conditions must be implemented before enabling.");
         static_assert(!enablePolymer, "Relevant treatment of boundary conditions must be implemented before enabling.");
-        static_assert(!enableMICP, "Relevant treatment of boundary conditions must be implemented before enabling.");
 
         // make sure that the right mass conservation quantities are used
         adaptMassConservationQuantities_(bdyFlux, insideIntQuants.pvtRegionIndex());
@@ -677,6 +688,7 @@ public:
 
     static void computeSource(RateVector& source,
                               const Problem& problem,
+                              const IntensiveQuantities& insideIntQuants,
                               unsigned globalSpaceIdex,
                               unsigned timeIdx)
     {
@@ -685,17 +697,16 @@ public:
         problem.source(source, globalSpaceIdex, timeIdx);
 
         // deal with MICP (if present)
-        // deal with micp (if present)
-        static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling.");
-        // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
+        MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
 
         // scale the source term of the energy equation
-        if (enableEnergy)
+        if constexpr(enableEnergy)
             source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     static void computeSourceDense(RateVector& source,
                                    const Problem& problem,
+                                   const IntensiveQuantities& insideIntQuants,
                                    unsigned globalSpaceIdex,
                                    unsigned timeIdx)
     {
@@ -703,12 +714,10 @@ public:
         problem.addToSourceDense(source, globalSpaceIdex, timeIdx);
 
         // deal with MICP (if present)
-        // deal with micp (if present)
-        static_assert(!enableMICP, "Relevant addSource() method must be implemented for this module before enabling.");
-        // MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
+        MICPModule::addSource(source, problem, insideIntQuants, globalSpaceIdex);
 
         // scale the source term of the energy equation
-        if (enableEnergy)
+        if constexpr(enableEnergy)
             source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
@@ -729,7 +738,7 @@ public:
 
         // scale the source term of the energy equation
         if constexpr(enableEnergy)
-                        source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
+            source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();
     }
 
     template <class UpEval, class FluidState>

opm/models/blackoil/blackoilmicpmodules.hh

@@ -25,8 +25,8 @@
  *
  * \brief Contains the classes required to extend the black-oil model by MICP.
  */
-#ifndef EWOMS_BLACK_OIL_MICP_MODULE_HH
-#define EWOMS_BLACK_OIL_MICP_MODULE_HH
+#ifndef OPM_BLACK_OIL_MICP_MODULE_HH
+#define OPM_BLACK_OIL_MICP_MODULE_HH
 
 #include <dune/common/fvector.hh>
 
@@ -59,9 +59,14 @@ class BlackOilMICPModule
     using EqVector = GetPropType<TypeTag, Properties::EqVector>;
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using Problem = GetPropType<TypeTag, Properties::Problem>;
 
     using Toolbox = MathToolbox<Evaluation>;
 
+    using TabulatedFunction = typename BlackOilMICPParams<Scalar>::TabulatedFunction;
+
+    enum { waterCompIdx = FluidSystem::waterCompIdx };
+
     static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx;
     static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
     static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx;
@@ -85,30 +90,13 @@ public:
         params_ = params;
     }
 
-    /*!
-     * \brief The simulator stops if "clogging" has been (almost) reached in any of the cells.
-     *
-     * I.e., porosity - biofilm - calcite < tol_clgg, where tol_clgg is a given tolerance. In the
-     * implemented model a permebaility-porosity relatonship is used where a minimum
-     * permeability value is reached if porosity - biofilm - calcite < phi_crit.
-     */
-    static void checkCloggingMICP(const Model& model, const Scalar phi, unsigned dofIdx)
-    {
-        const PrimaryVariables& priVars = model.solution(/*timeIdx=*/1)[dofIdx];
-        if (phi - priVars[biofilmConcentrationIdx] - priVars[calciteConcentrationIdx] < toleranceBeforeClogging())
-            throw std::logic_error("Clogging has been (almost) reached in at least one cell\n");
-    }
-
     /*!
      * \brief Register all run-time parameters for the black-oil MICP module.
      */
     static void registerParameters()
     {
-        if (!enableMICP)
-            // MICP has been disabled at compile time
-            return;
-
-        VtkBlackOilMICPModule<TypeTag>::registerParameters();
+        if constexpr (enableMICP)
+            VtkBlackOilMICPModule<TypeTag>::registerParameters();
     }
 
     /*!
@@ -117,20 +105,17 @@ public:
     static void registerOutputModules(Model& model,
                                       Simulator& simulator)
     {
-        if (!enableMICP)
-            // MICP has been disabled at compile time
-            return;
-
-        model.addOutputModule(new VtkBlackOilMICPModule<TypeTag>(simulator));
+        if constexpr (enableMICP)
+            model.addOutputModule(new VtkBlackOilMICPModule<TypeTag>(simulator));
     }
 
     static bool eqApplies(unsigned eqIdx)
     {
-        if (!enableMICP)
+        if constexpr (enableMICP)
+            return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx 
+                   || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx;
+        else
             return false;
-
-        // All MICP components are true here
-        return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx;
     }
 
     static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
@@ -146,208 +131,262 @@ public:
     static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
                            const IntensiveQuantities& intQuants)
     {
-        if (!enableMICP)
-            return;
-
-        LhsEval surfaceVolumeWater = Toolbox::template decay<LhsEval>(intQuants.porosity());
-        // avoid singular matrix if no water is present.
-        surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
-        // Suspended microbes in water phase
-        const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.microbialConcentration());
-        LhsEval accumulationMicrobes = massMicrobes;
-        storage[contiMicrobialEqIdx] += accumulationMicrobes;
-        // Oxygen in water phase
-        const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.oxygenConcentration());
-        LhsEval accumulationOxygen = massOxygen;
-        storage[contiOxygenEqIdx] += accumulationOxygen;
-        // Urea in water phase
-        const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.ureaConcentration());
-        LhsEval accumulationUrea = massUrea;
-        storage[contiUreaEqIdx] += accumulationUrea;
-        // Biofilm
-        const LhsEval massBiofilm = Toolbox::template decay<LhsEval>(intQuants.biofilmConcentration());
-        LhsEval accumulationBiofilm = massBiofilm;
-        storage[contiBiofilmEqIdx] += accumulationBiofilm;
-        // Calcite
-        const LhsEval massCalcite = Toolbox::template decay<LhsEval>(intQuants.calciteConcentration());
-        LhsEval accumulationCalcite = massCalcite;
-        storage[contiCalciteEqIdx] += accumulationCalcite;
+        if constexpr (enableMICP) {
+            const auto& fs = intQuants.fluidState();
+            LhsEval surfaceVolumeWater = Toolbox::template decay<LhsEval>(fs.invB(waterPhaseIdx))
+                                         * Toolbox::template decay<LhsEval>(intQuants.porosity());
+            // avoid singular matrix if no water is present
+            surfaceVolumeWater = max(surfaceVolumeWater, 1e-10);
+            // suspended microbes in water phase
+            const LhsEval massMicrobes = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.microbialConcentration());
+            LhsEval accumulationMicrobes = massMicrobes;
+            storage[contiMicrobialEqIdx] += accumulationMicrobes;
+            // oxygen in water phase
+            const LhsEval massOxygen = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.oxygenConcentration());
+            LhsEval accumulationOxygen = massOxygen;
+            storage[contiOxygenEqIdx] += accumulationOxygen;
+            // urea in water phase (applying the scaling factor for the urea equation)
+            const LhsEval massUrea = surfaceVolumeWater * Toolbox::template decay<LhsEval>(intQuants.ureaConcentration());
+            LhsEval accumulationUrea = massUrea;
+            storage[contiUreaEqIdx] += accumulationUrea;
+            storage[contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
+            // biofilm
+            const LhsEval massBiofilm = Toolbox::template decay<LhsEval>(intQuants.biofilmConcentration());
+            LhsEval accumulationBiofilm = massBiofilm;
+            storage[contiBiofilmEqIdx] += accumulationBiofilm;
+            // calcite
+            const LhsEval massCalcite = Toolbox::template decay<LhsEval>(intQuants.calciteConcentration());
+            LhsEval accumulationCalcite = massCalcite;
+            storage[contiCalciteEqIdx] += accumulationCalcite;
+        }
     }
 
-    static void computeFlux(RateVector& flux,
-                            const ElementContext& elemCtx,
-                            unsigned scvfIdx,
-                            unsigned timeIdx)
+    template <class UpEval, class Eval, class IntensiveQuantities>
+    static void addMICPFluxes_(RateVector& flux,
+                               const Eval& volumeFlux,
+                               const IntensiveQuantities& upFs)
+    {
+        if constexpr (enableMICP) {
+            flux[contiMicrobialEqIdx] +=
+                decay<UpEval>(upFs.microbialConcentration())
+                * volumeFlux;
+            flux[contiOxygenEqIdx] +=
+                decay<UpEval>(upFs.oxygenConcentration())
+                * volumeFlux;
+            flux[contiUreaEqIdx] +=
+                decay<UpEval>(upFs.ureaConcentration())
+                * volumeFlux;
+        }
+    }
+
+    // since the urea concentration can be much larger than 1, then we apply a scaling factor
+    static void applyScaling(RateVector& flux)
+    {
+        if constexpr (enableMICP) {
+            flux[contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
+        }
+    }
+
+    static void computeFlux([[maybe_unused]] RateVector& flux,
+                            [[maybe_unused]] const ElementContext& elemCtx,
+                            [[maybe_unused]] unsigned scvfIdx,
+                            [[maybe_unused]] unsigned timeIdx)
 
     {
-        if (!enableMICP)
-            return;
+        if constexpr (enableMICP) {
+            flux[contiMicrobialEqIdx] = 0.0;
+            flux[contiOxygenEqIdx] = 0.0;
+            flux[contiUreaEqIdx] = 0.0;
+            const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+            unsigned focusIdx = elemCtx.focusDofIndex();
+            unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
+            if (upIdx == focusIdx)
+                addMICPFluxes_<Evaluation>(flux, elemCtx, scvfIdx, timeIdx);
+            else
+                addMICPFluxes_<Scalar>(flux, elemCtx, scvfIdx, timeIdx);
+        }
+    }
 
+    template <class UpstreamEval>
+    static void addMICPFluxes_(RateVector& flux,
+                               const ElementContext& elemCtx,
+                               unsigned scvfIdx,
+                               unsigned timeIdx)
+    {
         const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
-
-        const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
-        const unsigned inIdx = extQuants.interiorIndex();
+        unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
         const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
-
-
-        if (upIdx == inIdx) {
-            flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.microbialConcentration();
-            flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.oxygenConcentration();
-            flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * up.ureaConcentration();
-        }
-        else {
-            flux[contiMicrobialEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.microbialConcentration());
-            flux[contiOxygenEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.oxygenConcentration());
-            flux[contiUreaEqIdx] = extQuants.volumeFlux(waterPhaseIdx) * decay<Scalar>(up.ureaConcentration());
-        }
+        const auto& fs = up.fluidState();
+        const auto& volFlux = extQuants.volumeFlux(waterPhaseIdx);
+        addMICPFluxes_<UpstreamEval>(flux, volFlux, fs);
     }
 
-    // See https://doi.org/10.1016/j.ijggc.2021.103256 for the micp processes in the model.
+    // see https://doi.org/10.1016/j.ijggc.2021.103256 for the MICP processes in the model
     static void addSource(RateVector& source,
-                            const ElementContext& elemCtx,
-                            unsigned dofIdx,
-                            unsigned timeIdx)
-
+                          const Problem& problem,
+                          const IntensiveQuantities& intQuants,
+                          unsigned globalSpaceIdex)
     {
-        if (!enableMICP)
-            return;
+        if constexpr (enableMICP) {
+            const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
+            auto velocityInfos = velocityInf[globalSpaceIdex];
+            Scalar normVelocityCell = 0.0;
+            for (const auto& velocityInfo : velocityInfos) {
+                normVelocityCell = max(normVelocityCell, std::abs(velocityInfo.velocity[waterPhaseIdx]));
+            }
 
-        // compute dpW (max norm of the pressure gradient in the cell center)
-        const IntensiveQuantities& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
-        const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, 0);
-        size_t numInteriorFaces = elemCtx.numInteriorFaces(timeIdx);
-        Evaluation dpW = 0;
-        for (unsigned scvfIdx = 0; scvfIdx < numInteriorFaces; scvfIdx++) {
-          const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
-          unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
-          const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
-          const Evaluation& mobWater = up.mobility(waterPhaseIdx);
+            // get the model parameters
+            const auto b = intQuants.fluidState().invB(waterPhaseIdx);
+            unsigned satnumIdx = problem.satnumRegionIndex(globalSpaceIdex);
+            Scalar k_a = microbialAttachmentRate(satnumIdx);
+            Scalar k_d = microbialDeathRate(satnumIdx);
+            Scalar rho_b = densityBiofilm(satnumIdx);
+            Scalar rho_c = densityCalcite(satnumIdx);
+            Scalar k_str = detachmentRate(satnumIdx);
+            Scalar eta = detachmentExponent(satnumIdx);
+            Scalar k_o = halfVelocityOxygen(satnumIdx);
+            Scalar k_u = halfVelocityUrea(satnumIdx);
+            Scalar mu = maximumGrowthRate(satnumIdx);
+            Scalar mu_u = maximumUreaUtilization(satnumIdx);
+            Scalar Y_sb = yieldGrowthCoefficient(satnumIdx);
+            Scalar F = oxygenConsumptionFactor(satnumIdx);
+            Scalar Y_uc = yieldUreaToCalciteCoefficient(satnumIdx);
 
-          // compute water velocity from flux
-          Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / (K[0][0] * mobWater);
-          dpW = std::max(dpW, abs(waterVolumeVelocity));
+            // compute Monod terms (the negative region is replaced by a straight line)
+            // Schäfer et al (1998) https://doi.org/10.1016/S0169-7722(97)00060-0
+            Evaluation k_g = mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration());
+            Evaluation k_c = mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
+            if (intQuants.oxygenConcentration() < 0)
+                k_g = mu * intQuants.oxygenConcentration() / k_o;
+            if (intQuants.ureaConcentration() < 0)
+                k_c = mu_u * intQuants.ureaConcentration() / k_u;
+
+            // compute the processes
+            source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() *
+                                                    b * (Y_sb * k_g - k_d - k_a) +
+                                                    rho_b * intQuants.biofilmConcentration() * k_str * pow(normVelocityCell, eta);
+
+            source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() * 
+                                                  b + rho_b * intQuants.biofilmConcentration()) * F * k_g;
+
+            source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * k_c;
+
+            source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * k_g - k_d -
+                                                  k_str * pow(normVelocityCell, eta) - Y_uc * (rho_b / rho_c) * 
+                                                  intQuants.biofilmConcentration() * k_c / (intQuants.porosity() + 
+                                                  intQuants.biofilmConcentration())) + k_a * intQuants.microbialConcentration() *
+                                                  intQuants.porosity() * b / rho_b;
+
+            source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * k_c;
+            
+            // since the urea concentration can be much larger than 1, then we apply a scaling factor
+            source[Indices::contiUreaEqIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
         }
-
-        // get the model parameters
-        Scalar k_a = microbialAttachmentRate();
-        Scalar k_d = microbialDeathRate();
-        Scalar rho_b = densityBiofilm();
-        Scalar rho_c = densityCalcite();
-        Scalar k_str = detachmentRate();
-        Scalar k_o = halfVelocityOxygen();
-        Scalar k_u = halfVelocityUrea() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp)
-        Scalar mu = maximumGrowthRate();
-        Scalar mu_u = maximumUreaUtilization() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp)
-        Scalar Y_sb = yieldGrowthCoefficient();
-        Scalar F = oxygenConsumptionFactor();
-        Scalar Y_uc = 1.67 * 10; //Multiplying by scaling factor 10 (see WellInterface_impl.hpp)
-
-        // compute the processes
-        source[Indices::contiMicrobialEqIdx] += intQuants.microbialConcentration() * intQuants.porosity() *
-                                                      (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d - k_a)
-                                                + rho_b * intQuants.biofilmConcentration() * k_str * pow(intQuants.porosity() * dpW, 0.58);
-
-        source[Indices::contiOxygenEqIdx] -= (intQuants.microbialConcentration() * intQuants.porosity() + rho_b * intQuants.biofilmConcentration()) *
-                                                    F * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration());
-
-        source[Indices::contiUreaEqIdx] -= rho_b * intQuants.biofilmConcentration() * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
-
-        source[Indices::contiBiofilmEqIdx] += intQuants.biofilmConcentration() * (Y_sb * mu * intQuants.oxygenConcentration() / (k_o + intQuants.oxygenConcentration()) - k_d
-                                                              - k_str * pow(intQuants.porosity() * dpW, 0.58) - Y_uc * (rho_b / rho_c) * intQuants.biofilmConcentration() * mu_u *
-                                                                  (intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration())) / (intQuants.porosity() + intQuants.biofilmConcentration()))
-                                              + k_a * intQuants.microbialConcentration() * intQuants.porosity() / rho_b;
-
-        source[Indices::contiCalciteEqIdx] += (rho_b / rho_c) * intQuants.biofilmConcentration() * Y_uc * mu_u * intQuants.ureaConcentration() / (k_u + intQuants.ureaConcentration());
     }
 
-    static const Scalar densityBiofilm()
+    static void addSource([[maybe_unused]] RateVector& source,
+                          [[maybe_unused]] const ElementContext& elemCtx,
+                          [[maybe_unused]] unsigned dofIdx,
+                          [[maybe_unused]] unsigned timeIdx)
     {
-        return params_.densityBiofilm_;
+        if constexpr (enableMICP) {
+            const auto& problem = elemCtx.problem();
+            const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
+            addSource(source, problem, intQuants, dofIdx);
+        }
     }
 
-    static const Scalar densityCalcite()
+    static const Scalar densityBiofilm(unsigned satnumRegionIdx)
     {
-        return params_.densityCalcite_;
+        return params_.densityBiofilm_[satnumRegionIdx];
     }
 
-    static const Scalar detachmentRate()
+    static const Scalar densityCalcite(unsigned satnumRegionIdx)
     {
-        return params_.detachmentRate_;
+        return params_.densityCalcite_[satnumRegionIdx];
     }
 
-    static const Scalar criticalPorosity()
+    static const Scalar detachmentRate(unsigned satnumRegionIdx)
     {
-        return params_.criticalPorosity_;
+        return params_.detachmentRate_[satnumRegionIdx];
     }
 
-    static const Scalar fittingFactor()
+    static const Scalar detachmentExponent(unsigned satnumRegionIdx)
     {
-        return params_.fittingFactor_;
+        return params_.detachmentExponent_[satnumRegionIdx];
     }
 
-    static const Scalar halfVelocityOxygen()
+    static const Scalar halfVelocityOxygen(unsigned satnumRegionIdx)
     {
-        return params_.halfVelocityOxygen_;
+        return params_.halfVelocityOxygen_[satnumRegionIdx];
     }
 
-    static const Scalar halfVelocityUrea()
+    static const Scalar halfVelocityUrea(unsigned satnumRegionIdx)
     {
-        return params_.halfVelocityUrea_;
+        return params_.halfVelocityUrea_[satnumRegionIdx];
     }
 
-    static const Scalar maximumGrowthRate()
+    static const Scalar maximumGrowthRate(unsigned satnumRegionIdx)
     {
-        return params_.maximumGrowthRate_;
+        return params_.maximumGrowthRate_[satnumRegionIdx];
     }
 
-    static const Scalar maximumOxygenConcentration()
+    static const Scalar maximumUreaUtilization(unsigned satnumRegionIdx)
     {
-        return params_.maximumOxygenConcentration_;
+        return params_.maximumUreaUtilization_[satnumRegionIdx];
     }
 
-    static const Scalar maximumUreaConcentration()
+    static const Scalar microbialAttachmentRate(unsigned satnumRegionIdx)
     {
-        return params_.maximumUreaConcentration_ / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp);
+        return params_.microbialAttachmentRate_[satnumRegionIdx];
     }
 
-    static const Scalar maximumUreaUtilization()
+    static const Scalar microbialDeathRate(unsigned satnumRegionIdx)
     {
-        return params_.maximumUreaUtilization_;
+        return params_.microbialDeathRate_[satnumRegionIdx];
     }
 
-    static const Scalar microbialAttachmentRate()
+    static const Scalar oxygenConsumptionFactor(unsigned satnumRegionIdx)
     {
-        return params_.microbialAttachmentRate_;
+        return params_.oxygenConsumptionFactor_[satnumRegionIdx];
     }
 
-    static const Scalar microbialDeathRate()
+    static const Scalar yieldGrowthCoefficient(unsigned satnumRegionIdx)
     {
-        return params_.microbialDeathRate_;
+        return params_.yieldGrowthCoefficient_[satnumRegionIdx];
     }
 
-    static const Scalar minimumPermeability()
+    static const Scalar yieldUreaToCalciteCoefficient(unsigned satnumRegionIdx)
     {
-        return params_.minimumPermeability_;
+        return params_.yieldUreaToCalciteCoefficient_[satnumRegionIdx];
     }
 
-    static const Scalar oxygenConsumptionFactor()
+    static const Scalar microbialDiffusion(unsigned pvtRegionIdx)
     {
-        return params_.oxygenConsumptionFactor_;
+        return params_.microbialDiffusion_[pvtRegionIdx];
     }
 
-    static const Scalar toleranceBeforeClogging()
+    static const Scalar oxygenDiffusion(unsigned pvtRegionIdx)
     {
-        return params_.toleranceBeforeClogging_;
+        return params_.oxygenDiffusion_[pvtRegionIdx];
     }
 
-    static const Scalar yieldGrowthCoefficient()
+    static const Scalar ureaDiffusion(unsigned pvtRegionIdx)
     {
-        return params_.yieldGrowthCoefficient_;
+        return params_.ureaDiffusion_[pvtRegionIdx];
     }
 
-    static const std::vector<Scalar> phi()
+    static const TabulatedFunction& permfactTable(const ElementContext& elemCtx,
+                                                  unsigned scvIdx,
+                                                  unsigned timeIdx)
     {
-        return params_.phi_;
+        unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
+        return params_.permfactTable_[satnumRegionIdx];
+    }
+
+    static const TabulatedFunction& permfactTable(unsigned satnumRegionIdx)
+    {
+        return params_.permfactTable_[satnumRegionIdx];
     }
 
 private:
@@ -396,17 +435,12 @@ public:
      *
      */
     void MICPPropertiesUpdate_(const ElementContext& elemCtx,
-                                  unsigned dofIdx,
-                                  unsigned timeIdx)
+                               unsigned dofIdx,
+                               unsigned timeIdx)
     {
         const auto linearizationType = elemCtx.linearizationType();
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
-        const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
-        const auto& K = elemCtx.problem().intrinsicPermeability(elemCtx, dofIdx, timeIdx);
-        Scalar referencePorosity_ = elemCtx.problem().porosity(elemCtx, dofIdx, timeIdx);
-        Scalar eta = MICPModule::fittingFactor();
-        Scalar k_min = MICPModule::minimumPermeability();
-        Scalar phi_crit = MICPModule::criticalPorosity();
+        const Scalar referencePorosity_ = elemCtx.problem().referencePorosity(dofIdx, timeIdx);
 
         microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType);
         oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType);
@@ -414,9 +448,13 @@ public:
         biofilmConcentration_ = priVars.makeEvaluation(biofilmConcentrationIdx, timeIdx, linearizationType);
         calciteConcentration_ = priVars.makeEvaluation(calciteConcentrationIdx, timeIdx, linearizationType);
 
-        // Permeability reduction due to MICP, by adjusting the water mobility
-        asImp_().mobility_[waterPhaseIdx] *= max((pow((intQuants.porosity() - phi_crit) / (referencePorosity_ - phi_crit), eta) + k_min / K[0][0])/(1. + k_min / K[0][0]), k_min / K[0][0]);
+        const Evaluation porosityFactor  = min(1.0 - (biofilmConcentration_+calciteConcentration_)/(referencePorosity_+1e-8), 1.0); //phi/phi_0
+        unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, dofIdx, timeIdx);
+        const auto& permfactTable = MICPModule::permfactTable(satnumRegionIdx);
+        permFactor_ = permfactTable.eval(porosityFactor, /*extrapolation=*/true);
 
+        biofilmMass_ = referencePorosity_*biofilmConcentration_*MICPModule::densityBiofilm(satnumRegionIdx);
+        calciteMass_ = referencePorosity_*calciteConcentration_*MICPModule::densityCalcite(satnumRegionIdx);
     }
 
     const Evaluation& microbialConcentration() const
@@ -434,16 +472,24 @@ public:
     const Evaluation& calciteConcentration() const
     { return calciteConcentration_; }
 
+    const Evaluation biofilmMass() const
+    { return biofilmMass_; }
+
+    const Evaluation calciteMass() const
+    { return calciteMass_; }
+
+    const Evaluation& permFactor() const
+    { return permFactor_; }
 
 protected:
-    Implementation& asImp_()
-    { return *static_cast<Implementation*>(this); }
-
     Evaluation microbialConcentration_;
     Evaluation oxygenConcentration_;
     Evaluation ureaConcentration_;
     Evaluation biofilmConcentration_;
     Evaluation calciteConcentration_;
+    Evaluation biofilmMass_;
+    Evaluation calciteMass_;
+    Evaluation permFactor_;
 
 };
 
@@ -456,8 +502,8 @@ class BlackOilMICPIntensiveQuantities<TypeTag, false>
 
 public:
     void MICPPropertiesUpdate_(const ElementContext&,
-                                  unsigned,
-                                  unsigned)
+                               unsigned,
+                               unsigned)
     { }
 
     const Evaluation& microbialConcentration() const
@@ -474,6 +520,15 @@ public:
 
     const Evaluation& calciteConcentration() const
     { throw std::logic_error("calciteConcentration() called but MICP is disabled"); }
+
+    const Evaluation& biofilmMass() const
+    { throw std::logic_error("biofilmMass() called but MICP is disabled"); }
+
+    const Evaluation& calciteMass() const
+    { throw std::logic_error("calciteMass() called but MICP is disabled"); }
+
+    const Evaluation& permFactor() const
+    { throw std::logic_error("permFactor() called but MICP is disabled"); }
 };
 
 /*!
@@ -487,11 +542,6 @@ template <class TypeTag, bool enableMICPV = getPropValue<TypeTag, Properties::En
 class BlackOilMICPExtensiveQuantities
 {
     using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
-
-private:
-    Implementation& asImp_()
-    { return *static_cast<Implementation*>(this); }
-
 };
 
 template <class TypeTag>

opm/models/blackoil/blackoilmicpparams.cpp

@@ -26,7 +26,11 @@
 
 #if HAVE_ECL_INPUT
 #include <opm/input/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/EclipseState/MICPpara.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/BiofilmTable.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/DiffMICPTable.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/PermfactTable.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/SimpleTable.hpp>
+#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
 #endif
 
 #include <algorithm>
@@ -59,33 +63,71 @@ initFromState(const EclipseState& eclState)
     if (!eclState.runspec().micp())
         return; // MICP treatment is supposed to be disabled*/
 
-    // initialize the objects which deal with the MICPpara keyword
-    const auto& MICPpara = eclState.getMICPpara();
-    densityBiofilm_ = MICPpara.getDensityBiofilm();
-    densityCalcite_ = MICPpara.getDensityCalcite();
-    detachmentRate_ = MICPpara.getDetachmentRate();
-    criticalPorosity_ = MICPpara.getCriticalPorosity();
-    fittingFactor_ = MICPpara.getFittingFactor();
-    halfVelocityOxygen_ = MICPpara.getHalfVelocityOxygen();
-    halfVelocityUrea_ = MICPpara.getHalfVelocityUrea();
-    maximumGrowthRate_ = MICPpara.getMaximumGrowthRate();
-    maximumUreaUtilization_ = MICPpara.getMaximumUreaUtilization();
-    microbialAttachmentRate_ = MICPpara.getMicrobialAttachmentRate();
-    microbialDeathRate_ = MICPpara.getMicrobialDeathRate();
-    minimumPermeability_ = MICPpara.getMinimumPermeability();
-    oxygenConsumptionFactor_ = MICPpara.getOxygenConsumptionFactor();
-    yieldGrowthCoefficient_ = MICPpara.getYieldGrowthCoefficient();
-    maximumOxygenConcentration_ = MICPpara.getMaximumOxygenConcentration();
-    maximumUreaConcentration_ = MICPpara.getMaximumUreaConcentration();
-    toleranceBeforeClogging_ = MICPpara.getToleranceBeforeClogging();
+    const auto& tableManager = eclState.getTableManager();
+    unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
+    unsigned numPvtRegions = tableManager.getTabdims().getNumPVTTables();
 
-    // obtain the porosity for the clamp in the blackoilnewtonmethod
-    if constexpr (std::is_same_v<Scalar, float>) {
-        const auto phi = eclState.fieldProps().get_double("PORO");
-        phi_.resize(phi.size());
-        std::copy(phi.begin(), phi.end(), phi_.begin());
-    } else {
-        phi_ = eclState.fieldProps().get_double("PORO");
+    // initialize the objects which deal with the MICP parameters
+    const TableContainer& biofilmTables = tableManager.getBiofilmTables();
+    if (biofilmTables.empty()) {
+        throw std::runtime_error("MICP requires the BIOFPARA keyword");
+    }
+    densityBiofilm_.resize(numSatRegions);
+    densityCalcite_.resize(numSatRegions);
+    detachmentRate_.resize(numSatRegions);
+    detachmentExponent_.resize(numSatRegions);
+    halfVelocityOxygen_.resize(numSatRegions);
+    halfVelocityUrea_.resize(numSatRegions);
+    maximumGrowthRate_.resize(numSatRegions);
+    maximumUreaUtilization_.resize(numSatRegions);
+    microbialAttachmentRate_.resize(numSatRegions);
+    microbialDeathRate_.resize(numSatRegions);
+    oxygenConsumptionFactor_.resize(numSatRegions);
+    yieldGrowthCoefficient_.resize(numSatRegions);
+    yieldUreaToCalciteCoefficient_.resize(numSatRegions);
+    for (unsigned stnRegionIdx = 0; stnRegionIdx < numSatRegions; ++stnRegionIdx) {
+        const BiofilmTable& biofilmTable = biofilmTables.getTable<BiofilmTable>(stnRegionIdx);
+        densityBiofilm_[stnRegionIdx] = biofilmTable.getDensityBiofilm().front();
+        densityCalcite_[stnRegionIdx] = biofilmTable.getDensityCalcite().front();
+        detachmentRate_[stnRegionIdx] = biofilmTable.getDetachmentRate().front();
+        detachmentExponent_[stnRegionIdx] = biofilmTable.getDetachmentExponent().front();
+        halfVelocityOxygen_[stnRegionIdx] = biofilmTable.getHalfVelocityOxygen().front();
+        halfVelocityUrea_[stnRegionIdx] = biofilmTable.getHalfVelocityUrea().front();
+        maximumGrowthRate_[stnRegionIdx] = biofilmTable.getMaximumGrowthRate().front();
+        maximumUreaUtilization_[stnRegionIdx] = biofilmTable.getMaximumUreaUtilization().front();
+        microbialAttachmentRate_[stnRegionIdx] = biofilmTable.getMicrobialAttachmentRate().front();
+        microbialDeathRate_[stnRegionIdx] = biofilmTable.getMicrobialDeathRate().front();
+        oxygenConsumptionFactor_[stnRegionIdx] = biofilmTable.getOxygenConsumptionFactor().front();
+        yieldGrowthCoefficient_[stnRegionIdx] = biofilmTable.getYieldGrowthCoefficient().front();
+        yieldUreaToCalciteCoefficient_[stnRegionIdx] = biofilmTable.getYieldUreaToCalciteCoefficient().front();
+    }
+
+    const TableContainer& diffMICPTables = tableManager.getDiffMICPTables();
+    microbialDiffusion_.resize(numPvtRegions);
+    oxygenDiffusion_.resize(numPvtRegions);
+    ureaDiffusion_.resize(numPvtRegions);
+    for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++pvtRegionIdx) {
+        if (!diffMICPTables.empty()) {
+            const DiffMICPTable& diffMICPTable = diffMICPTables.getTable<DiffMICPTable>(pvtRegionIdx);
+            microbialDiffusion_[pvtRegionIdx] = diffMICPTable.getMicrobialDiffusion().front();
+            oxygenDiffusion_[pvtRegionIdx] = diffMICPTable.getOxygenDiffusion().front();
+            ureaDiffusion_[pvtRegionIdx] = diffMICPTable.getUreaDiffusion().front();
+        }
+        else {
+            microbialDiffusion_[pvtRegionIdx] = 0.0;
+            oxygenDiffusion_[pvtRegionIdx] = 0.0;
+            ureaDiffusion_[pvtRegionIdx] = 0.0;
+        }
+    }
+
+    const TableContainer& permfactTables = tableManager.getPermfactTables();
+    if (permfactTables.empty()) {
+        throw std::runtime_error("MICP requires the PERMFACT keyword");
+    }
+    permfactTable_.resize(numSatRegions);
+    for (std::size_t i = 0; i < permfactTables.size(); ++i) {
+        const PermfactTable& permfactTable = permfactTables.getTable<PermfactTable>(i);
+        permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn());
     }
 }
 #endif

opm/models/blackoil/blackoilmicpparams.hpp

@@ -28,6 +28,8 @@
 #ifndef OPM_BLACK_OIL_MICP_PARAMS_HPP
 #define OPM_BLACK_OIL_MICP_PARAMS_HPP
 
+#include <opm/material/common/Tabulated1DFunction.hpp>
+
 #include <vector>
 
 namespace Opm {
@@ -45,24 +47,26 @@ struct BlackOilMICPParams
     void initFromState(const EclipseState& eclState);
 #endif
 
-    Scalar densityBiofilm_;
-    Scalar densityCalcite_;
-    Scalar detachmentRate_;
-    Scalar criticalPorosity_;
-    Scalar fittingFactor_;
-    Scalar halfVelocityOxygen_;
-    Scalar halfVelocityUrea_;
-    Scalar maximumGrowthRate_;
-    Scalar maximumUreaUtilization_;
-    Scalar microbialAttachmentRate_;
-    Scalar microbialDeathRate_;
-    Scalar minimumPermeability_;
-    Scalar oxygenConsumptionFactor_;
-    Scalar yieldGrowthCoefficient_;
-    Scalar maximumOxygenConcentration_;
-    Scalar maximumUreaConcentration_;
-    Scalar toleranceBeforeClogging_;
-    std::vector<Scalar> phi_;
+    using TabulatedFunction = Tabulated1DFunction<Scalar>;
+
+    std::vector<TabulatedFunction> permfactTable_;
+
+    std::vector<Scalar> densityBiofilm_;
+    std::vector<Scalar> densityCalcite_;
+    std::vector<Scalar> detachmentRate_;
+    std::vector<Scalar> detachmentExponent_;
+    std::vector<Scalar> halfVelocityOxygen_;
+    std::vector<Scalar> halfVelocityUrea_;
+    std::vector<Scalar> maximumGrowthRate_;
+    std::vector<Scalar> maximumUreaUtilization_;
+    std::vector<Scalar> microbialAttachmentRate_;
+    std::vector<Scalar> microbialDeathRate_;
+    std::vector<Scalar> oxygenConsumptionFactor_;
+    std::vector<Scalar> yieldGrowthCoefficient_;
+    std::vector<Scalar> yieldUreaToCalciteCoefficient_;
+    std::vector<Scalar> microbialDiffusion_;
+    std::vector<Scalar> oxygenDiffusion_;
+    std::vector<Scalar> ureaDiffusion_;
 };
 
 } // namespace Opm

opm/models/blackoil/blackoilmodel.hh

@@ -197,6 +197,19 @@ public:
     static constexpr Scalar value = 1.0/(30.0*4184.0*alpha);
 };
 
+//! similarly to the energy equation, the urea equation in MICP is scaled by a factor of 0.1
+template<class TypeTag>
+struct BlackOilUreaScalingFactor<TypeTag, TTag::BlackOilModel>
+{
+private:
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    static constexpr Scalar alpha = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>() ? 1000.0 : 1.0;
+
+public:
+    using type = Scalar;
+    static constexpr Scalar value = 1.0/(10.0*alpha);
+};
+
 // by default, ebos formulates the conservation equations in terms of mass not surface
 // volumes
 template<class TypeTag>

opm/models/blackoil/blackoilnewtonmethod.hpp

@@ -372,26 +372,29 @@ protected:
                 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], bparams_.pressMin_, bparams_.pressMax_);
             }
 
-
-            // Limit the variables to [0, cmax] values to improve the convergence.
-            // For the microorganisms we set this value equal to the biomass density value.
-            // For the oxygen and urea we set this value to the maximum injected
-            // concentration (the urea concentration has been scaled by 10). For
-            // the biofilm and calcite, we set this value equal to the porosity minus the clogging tolerance.
-            if (enableMICP && pvIdx == Indices::microbialConcentrationIdx) {
-                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::densityBiofilm());
-            }
-            if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx) {
-                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumOxygenConcentration());
-            }
-            if (enableMICP && pvIdx == Indices::ureaConcentrationIdx) {
-                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::maximumUreaConcentration());
-            }
-            if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx) {
-                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
-            }
-            if (enableMICP && pvIdx == Indices::calciteConcentrationIdx) {
-                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0}, MICPModule::phi()[globalDofIdx] - MICPModule::toleranceBeforeClogging());
+            // keep the values above 0
+            // for the biofilm and calcite, we set an upper limit equal to the initial porosity
+            // minus 1e-8. This prevents singularities (e.g., one of the calcite source term is 
+            // evaluated at 1/(iniPoro - calcite)). The value 1e-8 is taken from the salt precipitation
+            // clapping above. 
+            if constexpr (enableMICP) {
+                if (pvIdx == Indices::microbialConcentrationIdx) {
+                    nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
+                }
+                if (pvIdx == Indices::oxygenConcentrationIdx) {
+                    nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
+                }
+                if (pvIdx == Indices::ureaConcentrationIdx) {
+                    nextValue[pvIdx] = std::max(nextValue[pvIdx], Scalar{0.0});
+                }
+                if (pvIdx == Indices::biofilmConcentrationIdx) {
+                    nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0},
+                                                                    this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
+                }
+                if (pvIdx == Indices::calciteConcentrationIdx) {
+                    nextValue[pvIdx] = std::clamp(nextValue[pvIdx], Scalar{0.0},
+                                                                    this->problem().referencePorosity(globalDofIdx, 0) - 1e-8);
+                }
             }
         }
 

opm/models/blackoil/blackoilproperties.hh

@@ -90,6 +90,9 @@ struct EnableMech { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct BlackOilEnergyScalingFactor { using type = UndefinedProperty; };
 
+//! Similarly to the energy equation, a scaling is applied to the urea equation in MICP
+template<class TypeTag, class MyTypeTag>
+struct BlackOilUreaScalingFactor { using type = UndefinedProperty; };
 
 } // namespace Opm::Properties
 

opm/models/discretization/common/tpfalinearizer.hh

@@ -108,6 +108,7 @@ class TpfaLinearizer
     static const bool linearizeNonLocalElements = getPropValue<TypeTag, Properties::LinearizeNonLocalElements>();
     static const bool enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>();
     static const bool enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>();
+    static const bool enableMICP = getPropValue<TypeTag, Properties::EnableMICP>();
 
     // copying the linearizer is not a good idea
     TpfaLinearizer(const TpfaLinearizer&) = delete;
@@ -543,7 +544,7 @@ private:
                               simulator_().problem().eclWriter().outputModule().getFlows().hasBlockFlows();
         const bool anyFlores = simulator_().problem().eclWriter().outputModule().getFlows().anyFlores();
         const bool enableDispersion = simulator_().vanguard().eclState().getSimulationConfig().rock_config().dispersion();
-        if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && !enableDispersion) {
+        if (((!anyFlows || !flowsInfo_.empty()) && (!anyFlores || !floresInfo_.empty())) && (!enableDispersion && !enableMICP)) {
             return;
         }
         const auto& model = model_();
@@ -569,7 +570,7 @@ private:
         if (anyFlores) {
             floresInfo_.reserve(numCells, 6 * numCells);
         }
-        if (enableDispersion) {
+        if (enableDispersion || enableMICP) {
             velocityInfo_.reserve(numCells, 6 * numCells);
         }
 
@@ -616,7 +617,7 @@ private:
                 if (anyFlores) {
                     floresInfo_.appendRow(loc_flinfo.begin(), loc_flinfo.end());
                 }
-                if (enableDispersion) {
+                if (enableDispersion || enableMICP) {
                     velocityInfo_.appendRow(loc_vlinfo.begin(), loc_vlinfo.end());
                 }
             }
@@ -757,7 +758,7 @@ private:
                 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
                 LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo,  problem_().moduleParams());
                 adres *= nbInfo.res_nbinfo.faceArea;
-                if (enableDispersion) {
+                if (enableDispersion || enableMICP) {
                     for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {
                         velocityInfo_[globI][loc].velocity[phaseIdx] = darcyFlux[phaseIdx].value() / nbInfo.res_nbinfo.faceArea;
                     }
@@ -831,9 +832,9 @@ private:
             bMat = 0.0;
             adres = 0.0;
             if (separateSparseSourceTerms_) {
-                LocalResidual::computeSourceDense(adres, problem_(), globI, 0);
+                LocalResidual::computeSourceDense(adres, problem_(), intQuantsIn, globI, 0);
             } else {
-                LocalResidual::computeSource(adres, problem_(), globI, 0);
+                LocalResidual::computeSource(adres, problem_(), intQuantsIn, globI, 0);
             }
             adres *= -volume;
             setResAndJacobi(res, bMat, adres);

opm/models/io/vtkblackoilmicpmodule.hpp

@@ -142,7 +142,7 @@ public:
 
                 if (params_.ureaConcentrationOutput_) {
                     ureaConcentration_[globalDofIdx] =
-                        10 * scalarValue(intQuants.ureaConcentration());//Multypliging by scaling factor 10 (see WellInterface_impl.hpp)
+                        scalarValue(intQuants.ureaConcentration());
                 }
 
                 if (params_.biofilmConcentrationOutput_) {

opm/simulators/flow/FIPContainer.cpp

@@ -337,6 +337,91 @@ assignCo2InWater(const unsigned globalDofIdx,
     }
 }
 
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasMicrobialMass() const
+{
+    return has(Inplace::Phase::MicrobialMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignMicrobialMass(const unsigned globalDofIdx,
+                    const Scalar   microbialMass)
+{
+    this->fip_[Inplace::Phase::MicrobialMass][globalDofIdx] = microbialMass;
+}
+
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasOxygenMass() const
+{
+    return has(Inplace::Phase::OxygenMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignOxygenMass(const unsigned globalDofIdx,
+                 const Scalar   oxygenMass)
+{
+    this->fip_[Inplace::Phase::OxygenMass][globalDofIdx] = oxygenMass;
+}
+
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasUreaMass() const
+{
+    return has(Inplace::Phase::UreaMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignUreaMass(const unsigned globalDofIdx,
+               const Scalar   ureaMass)
+{
+    this->fip_[Inplace::Phase::UreaMass][globalDofIdx] = ureaMass;
+}
+
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasBiofilmMass() const
+{
+    return has(Inplace::Phase::BiofilmMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignBiofilmMass(const unsigned globalDofIdx,
+                  const Scalar   biofilmMass)
+{
+    this->fip_[Inplace::Phase::BiofilmMass][globalDofIdx] = biofilmMass;
+}
+
+template<class FluidSystem>
+bool
+FIPContainer<FluidSystem>::
+hasCalciteMass() const
+{
+    return has(Inplace::Phase::CalciteMass);
+}
+
+template<class FluidSystem>
+void
+FIPContainer<FluidSystem>::
+assignCalciteMass(const unsigned globalDofIdx,
+                  const Scalar   calciteMass)
+{
+    this->fip_[Inplace::Phase::CalciteMass][globalDofIdx] = calciteMass;
+}
+
 template<class FluidSystem>
 void
 FIPContainer<FluidSystem>::

opm/simulators/flow/FIPContainer.hpp

@@ -86,6 +86,26 @@ public:
                           const Scalar   co2InWater,
                           const Scalar   mM);
 
+    bool hasMicrobialMass() const;
+    void assignMicrobialMass(const unsigned globalDofIdx,
+                             const Scalar   microbialMass);
+
+    bool hasOxygenMass() const;
+    void assignOxygenMass(const unsigned globalDofIdx,
+                          const Scalar   oxygenMass);
+
+    bool hasUreaMass() const;
+    void assignUreaMass(const unsigned globalDofIdx,
+                        const Scalar   ureaMass);
+
+    bool hasBiofilmMass() const;
+    void assignBiofilmMass(const unsigned globalDofIdx,
+                           const Scalar   biofilmMass);
+
+    bool hasCalciteMass() const;
+    void assignCalciteMass(const unsigned globalDofIdx,
+                           const Scalar   calciteMass);
+
     void assignGasWater(const unsigned  globalDofIdx,
                         const std::array<Scalar, numPhases>& fip,
                         const Scalar    gasInPlaceWater,

opm/simulators/flow/FlowGenericProblem.hpp

@@ -168,7 +168,7 @@ public:
     /*!
      * \brief Returns the initial polymer concentration for a given a cell index
      */
-    Scalar  polymerConcentration(unsigned elemIdx) const;
+    Scalar polymerConcentration(unsigned elemIdx) const;
 
     /*!
     * \brief Returns the polymer molecule weight for a given cell index
@@ -179,27 +179,27 @@ public:
     /*!
      * \brief Returns the initial microbial concentration for a given a cell index
      */
-    Scalar  microbialConcentration(unsigned elemIdx) const;
+    Scalar microbialConcentration(unsigned elemIdx) const;
 
     /*!
      * \brief Returns the initial oxygen concentration for a given a cell index
      */
-    Scalar  oxygenConcentration(unsigned elemIdx) const;
+    Scalar oxygenConcentration(unsigned elemIdx) const;
 
     /*!
      * \brief Returns the initial urea concentration for a given a cell index
      */
-    Scalar  ureaConcentration(unsigned elemIdx) const;
+    Scalar ureaConcentration(unsigned elemIdx) const;
 
     /*!
      * \brief Returns the initial biofilm concentration for a given a cell index
      */
-    Scalar  biofilmConcentration(unsigned elemIdx) const;
+    Scalar biofilmConcentration(unsigned elemIdx) const;
 
     /*!
      * \brief Returns the initial calcite concentration for a given a cell index
      */
-    Scalar  calciteConcentration(unsigned elemIdx) const;
+    Scalar calciteConcentration(unsigned elemIdx) const;
 
     /*!
      * \brief Returns the index the relevant PVT region given a cell index

opm/simulators/flow/FlowProblem.hpp

@@ -1102,7 +1102,7 @@ public:
 
 
     /*!
-     * \brief Return the well transmissibility multiplier due to rock changues.
+     * \brief Return the well transmissibility multiplier due to rock changes.
      */
     template <class LhsEval>
     LhsEval wellTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
@@ -1153,6 +1153,15 @@ public:
         case BCComponent::POLYMER:
             this->handlePolymerBC(bc, rate);
             break;
+        case BCComponent::MICR:
+            this->handleMicrBC(bc, rate);
+            break;
+        case BCComponent::OXYG:
+            this->handleOxygBC(bc, rate);
+            break;
+        case BCComponent::UREA:
+            this->handleUreaBC(bc, rate);
+            break;
         case BCComponent::NONE:
             throw std::logic_error("you need to specify the component when RATE type is set in BC");
             break;
@@ -1701,6 +1710,12 @@ protected:
 
     virtual void handlePolymerBC(const BCProp::BCFace&, RateVector&) const = 0;
 
+    virtual void handleMicrBC(const BCProp::BCFace&, RateVector&) const = 0;
+
+    virtual void handleOxygBC(const BCProp::BCFace&, RateVector&) const = 0;
+
+    virtual void handleUreaBC(const BCProp::BCFace&, RateVector&) const = 0;
+
     BCData<int> bcindex_;
     bool nonTrivialBoundaryConditions_ = false;
     bool first_step_ = true;

opm/simulators/flow/FlowProblemBlackoil.hpp

@@ -486,17 +486,6 @@ public:
                     });
             });
         }
-
-        // Deal with "clogging" for the MICP model
-        if constexpr (enableMICP) {
-            auto& model = this->model();
-            const auto& residual = model.linearizer().residual();
-
-            for (unsigned globalDofIdx = 0; globalDofIdx < residual.size(); ++globalDofIdx) {
-                auto& phi = this->referencePorosity_[/*timeIdx=*/1][globalDofIdx];
-                MICPModule::checkCloggingMICP(model, phi, globalDofIdx);
-            }
-        }
     }
 
     /*!
@@ -629,6 +618,11 @@ public:
             if constexpr (enablePolymer) {
                 rate[Indices::polymerConcentrationIdx] += source.rate({ijk, SourceComponent::POLYMER}) / this->model().dofTotalVolume(globalDofIdx);
             }
+            if constexpr (enableMICP) {
+                rate[Indices::microbialConcentrationIdx] += source.rate({ijk, SourceComponent::MICR}) / this->model().dofTotalVolume(globalDofIdx);
+                rate[Indices::oxygenConcentrationIdx] += source.rate({ijk, SourceComponent::OXYG}) / this->model().dofTotalVolume(globalDofIdx);
+                rate[Indices::ureaConcentrationIdx] += source.rate({ijk, SourceComponent::UREA}) / (this->model().dofTotalVolume(globalDofIdx));
+            }
             if constexpr (enableEnergy) {
                 for (unsigned i = 0; i < phidx_map.size(); ++i) {
                     const auto phaseIdx = phidx_map[i];
@@ -692,7 +686,7 @@ public:
     LhsEval permFactTransMultiplier(const IntensiveQuantities& intQuants, unsigned elementIdx) const
     {
         OPM_TIMEBLOCK_LOCAL(permFactTransMultiplier);
-        if (!enableSaltPrecipitation)
+        if (!enableSaltPrecipitation && !enableMICP)
             return 1.0;
 
         const auto& fs = intQuants.fluidState();
@@ -700,8 +694,15 @@ public:
         
         LhsEval porosityFactor = decay<LhsEval>(1. - fs.saltSaturation());
         porosityFactor = min(porosityFactor, 1.0);
-        const auto& permfactTable = BrineModule::permfactTable(tableIdx);
-        return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
+        if constexpr (enableSaltPrecipitation) {
+            const auto& permfactTable = BrineModule::permfactTable(tableIdx);
+            return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
+        }
+        if constexpr (enableMICP){
+            const auto& permfactTable = MICPModule::permfactTable(tableIdx);
+            return permfactTable.eval(porosityFactor, /*extrapolation=*/true);
+        }
+        return 1.0;
     }
 
     // temporary solution to facilitate output of initial state from flow
@@ -763,6 +764,9 @@ public:
                         break;
                     case BCComponent::SOLVENT:
                     case BCComponent::POLYMER:
+                    case BCComponent::MICR:
+                    case BCComponent::OXYG:
+                    case BCComponent::UREA:
                     case BCComponent::NONE:
                         throw std::logic_error("you need to specify a valid component (OIL, WATER or GAS) when DIRICHLET type is set in BC");
                 }
@@ -1521,6 +1525,32 @@ protected:
         rate[Indices::polymerConcentrationIdx] = bc.rate;
     }
 
+    void handleMicrBC(const BCProp::BCFace& bc, RateVector& rate) const override
+    {
+        if constexpr (!enableMICP)
+            throw std::logic_error("MICP is disabled and you're trying to add microbes to BC");
+
+        rate[Indices::microbialConcentrationIdx] = bc.rate;
+    }
+
+    void handleOxygBC(const BCProp::BCFace& bc, RateVector& rate) const override
+    {
+        if constexpr (!enableMICP)
+            throw std::logic_error("MICP is disabled and you're trying to add oxygen to BC");
+
+        rate[Indices::oxygenConcentrationIdx] = bc.rate;
+    }
+
+    void handleUreaBC(const BCProp::BCFace& bc, RateVector& rate) const override
+    {
+        if constexpr (!enableMICP)
+            throw std::logic_error("MICP is disabled and you're trying to add urea to BC");
+
+        rate[Indices::ureaConcentrationIdx] = bc.rate;
+        // since the urea concentration can be much larger than 1, then we apply a scaling factor
+        rate[Indices::ureaConcentrationIdx] *= getPropValue<TypeTag, Properties::BlackOilUreaScalingFactor>();
+    }
+
     void updateExplicitQuantities_(const bool first_step_after_restart)
     {
         OPM_TIMEBLOCK(updateExplicitQuantities);

opm/simulators/flow/FlowProblemComp.hpp

@@ -607,6 +607,21 @@ private:
         throw std::logic_error("polymer is disabled for compositional modeling and you're trying to add polymer to BC");
     }
 
+    void handleMicrBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
+    {
+        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add microbes to BC");
+    }
+
+    void handleOxygBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
+    {
+        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add oxygen to BC");
+    }
+
+    void handleUreaBC(const BCProp::BCFace& /* bc */, RateVector& /* rate */) const override
+    {
+        throw std::logic_error("MICP is disabled for compositional modeling and you're trying to add urea to BC");
+    }
+
     FlowThresholdPressure<TypeTag> thresholdPressures_;
 
     std::vector<InitialFluidState> initialFluidStates_;

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -760,7 +760,7 @@ doAllocBuffers(const unsigned bufferSize,
        Entry{&cFoam_,                             "", enableFoam_},
        Entry{&cSalt_,                             "", enableBrine_},
        Entry{&pSalt_,                             "", enableSaltPrecipitation_},
-       Entry{&permFact_,                          "", enableSaltPrecipitation_},
+       Entry{&permFact_,                          "", enableSaltPrecipitation_ || enableMICP_},
        Entry{&soMax_,                             "", oilvap.getType() == OilVapP::VAPPARS},
        Entry{&soMax_,                             "", hysteresisConfig &&
                                                       hysteresisConfig->enableNonWettingHysteresis() &&

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -114,6 +114,7 @@ class OutputBlackOilModule : public GenericOutputBlackoilModule<GetPropType<Type
     static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
     static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
 
     template<int idx, class VectorType>
     static Scalar value_or_zero(const VectorType& v)
@@ -873,6 +874,25 @@ private:
         {
             this->updateCO2InWater(globalDofIdx, pv, fs);
         }
+
+        if constexpr(enableMICP) {
+            const auto surfVolWat = pv * getValue(fs.invB(waterPhaseIdx));
+            if (this->fipC_.hasMicrobialMass()) {
+                this->updateMicrobialMass(globalDofIdx, intQuants, surfVolWat);
+            }
+            if (this->fipC_.hasOxygenMass()) {
+                this->updateOxygenMass(globalDofIdx, intQuants, surfVolWat);
+            }
+            if (this->fipC_.hasUreaMass()) {
+                this->updateUreaMass(globalDofIdx, intQuants, surfVolWat);
+            }
+            if (this->fipC_.hasBiofilmMass()) {
+                this->updateBiofilmMass(globalDofIdx, intQuants, totVolume);
+            }
+            if (this->fipC_.hasCalciteMass()) {
+                this->updateCalciteMass(globalDofIdx, intQuants, totVolume);
+            }
+        }
     }
 
     template <typename FluidState, typename FIPArray>
@@ -991,6 +1011,56 @@ private:
         return xoG * pv * rhoo * so / mM;
     }
 
+    template <typename IntensiveQuantities>
+    void updateMicrobialMass(const unsigned             globalDofIdx,
+                             const IntensiveQuantities& intQuants,
+                             const double               surfVolWat)
+    {
+        const Scalar mass = surfVolWat * intQuants.microbialConcentration().value();
+
+        this->fipC_.assignMicrobialMass(globalDofIdx, mass);
+    }
+
+    template <typename IntensiveQuantities>
+    void updateOxygenMass(const unsigned             globalDofIdx,
+                          const IntensiveQuantities& intQuants,
+                          const double               surfVolWat)
+    {
+        const Scalar mass = surfVolWat * intQuants.oxygenConcentration().value();
+
+        this->fipC_.assignOxygenMass(globalDofIdx, mass);
+    }
+
+    template <typename IntensiveQuantities>
+    void updateUreaMass(const unsigned             globalDofIdx,
+                        const IntensiveQuantities& intQuants,
+                        const double               surfVolWat)
+    {
+        const Scalar mass = surfVolWat * intQuants.ureaConcentration().value();
+
+        this->fipC_.assignUreaMass(globalDofIdx, mass);
+    }
+
+    template <typename IntensiveQuantities>
+    void updateBiofilmMass(const unsigned             globalDofIdx,
+                           const IntensiveQuantities& intQuants,
+                           const double               totVolume)
+    {
+        const Scalar mass = totVolume * intQuants.biofilmMass().value();
+
+        this->fipC_.assignBiofilmMass(globalDofIdx, mass);
+    }
+
+    template <typename IntensiveQuantities>
+    void updateCalciteMass(const unsigned             globalDofIdx,
+                           const IntensiveQuantities& intQuants,
+                           const double               totVolume)
+    {
+        const Scalar mass = totVolume * intQuants.calciteMass().value();
+
+        this->fipC_.assignCalciteMass(globalDofIdx, mass);
+    }
+
     //! \brief Setup extractors for element-level data.
     void setupElementExtractors_()
     {
@@ -1441,8 +1511,7 @@ private:
                       micpC.assign(ectx.globalDofIdx,
                                    ectx.intQuants.microbialConcentration().value(),
                                    ectx.intQuants.oxygenConcentration().value(),
-                                   // Rescaling back the urea concentration (see WellInterface_impl.hpp)
-                                   10 * ectx.intQuants.ureaConcentration().value(),
+                                   ectx.intQuants.ureaConcentration().value(),
                                    ectx.intQuants.biofilmConcentration().value(),
                                    ectx.intQuants.calciteConcentration().value());
                   }, this->micpC_.allocated()
@@ -1958,6 +2027,49 @@ private:
                             }
                   }
             },
+            Entry{ScalarEntry{"BMMIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return getValue(ectx.intQuants.microbialConcentration()) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         model.dofTotalVolume(ectx.globalDofIdx);
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BMOIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return getValue(ectx.intQuants.oxygenConcentration()) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         model.dofTotalVolume(ectx.globalDofIdx);
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BMUIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return getValue(ectx.intQuants.ureaConcentration()) *
+                                         getValue(ectx.intQuants.porosity()) *
+                                         model.dofTotalVolume(ectx.globalDofIdx) * 1;
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BMBIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.biofilmMass());
+                              }
+                  }
+            },
+            Entry{ScalarEntry{"BMCIP",
+                              [&model = this->simulator_.model()](const Context& ectx)
+                              {
+                                  return model.dofTotalVolume(ectx.globalDofIdx) *
+                                         getValue(ectx.intQuants.calciteMass());
+                              }
+                  }
+            },
         };
 
         this->blockExtractors_ = BlockExtractor::setupExecMap(this->blockData_, handlers);

opm/simulators/utils/BlackoilPhases.cpp

@@ -55,4 +55,5 @@ Opm::PhaseUsage::PhaseUsage(std::vector<BlackoilPhases::PhaseIndex> phases)
     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;
 }

opm/simulators/utils/BlackoilPhases.hpp

@@ -60,9 +60,10 @@ namespace Opm
         bool has_foam{};
         bool has_brine{};
         bool has_zFraction{};
+        bool has_micp{};
 
     };
 
 } // namespace Opm
 
-#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED
+#endif // OPM_BLACKOILPHASES_HEADER_INCLUDED

opm/simulators/utils/phaseUsageFromDeck.cpp

@@ -138,7 +138,14 @@ PhaseUsage phaseUsageFromDeck(const EclipseState& eclipseState)
 {
     const auto& phases = eclipseState.runspec().phases();
 
-    return phaseUsage(phases);
+    if (eclipseState.runspec().micp()) {
+        PhaseUsage pu = phaseUsage(phases);
+        pu.has_micp = true;
+        return pu;
+    }
+    else {
+        return phaseUsage(phases);
+    }
 }
 
 /// Looks at presence of WATER, OIL and GAS keywords in deck

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -98,6 +98,10 @@ namespace Opm
             OPM_THROW(std::runtime_error, "water evaporation is not supported by multisegment well yet");
         }
 
+        if constexpr (Base::has_micp) {
+            OPM_THROW(std::runtime_error, "MICP is not supported by multisegment well yet");
+        }
+
         if(this->rsRvInj() > 0) {
             OPM_THROW(std::runtime_error,
                       "dissolved gas/ vapporized oil in injected oil/gas not supported by multisegment well yet."

opm/simulators/wells/PerfData.cpp

@@ -42,7 +42,9 @@ PerfData<Scalar>::PerfData(const std::size_t num_perf,
     , polymer_rates(num_perf)
     , brine_rates(num_perf)
     , prod_index(num_perf * num_phases)
-    , micp_rates(num_perf)
+    , microbial_rates(num_perf)
+    , oxygen_rates(num_perf)
+    , urea_rates(num_perf)
     , cell_index(num_perf)
     , connection_transmissibility_factor(num_perf)
     , connection_d_factor(num_perf)
@@ -79,7 +81,9 @@ PerfData<Scalar> PerfData<Scalar>::serializationTestObject()
     result.polymer_rates = {10.0, 11.0, 12.0};
     result.brine_rates = {13.0};
     result.prod_index = {14.0, 15.0};
-    result.micp_rates = {16.0};
+    result.microbial_rates = {16.0};
+    result.oxygen_rates = {16.0};
+    result.urea_rates = {16.0};
     result.cell_index = {17, 18, 19, 20};
     result.connection_transmissibility_factor = {21.0};
     result.connection_d_factor = {21.5};
@@ -127,7 +131,9 @@ bool PerfData<Scalar>::try_assign(const PerfData& other)
     this->polymer_rates = other.polymer_rates;
     this->brine_rates = other.brine_rates;
     this->prod_index = other.prod_index;
-    this->micp_rates = other.micp_rates;
+    this->microbial_rates = other.microbial_rates;
+    this->oxygen_rates = other.oxygen_rates;
+    this->urea_rates = other.urea_rates;
     this->water_throughput = other.water_throughput;
     this->skin_pressure = other.skin_pressure;
     this->water_velocity = other.water_velocity;
@@ -150,7 +156,9 @@ bool PerfData<Scalar>::operator==(const PerfData& rhs) const
         && (this->polymer_rates == rhs.polymer_rates)
         && (this->brine_rates == rhs.brine_rates)
         && (this->prod_index == rhs.prod_index)
-        && (this->micp_rates == rhs.micp_rates)
+        && (this->microbial_rates == rhs.microbial_rates)
+        && (this->oxygen_rates == rhs.oxygen_rates)
+        && (this->urea_rates == rhs.urea_rates)
         && (this->cell_index == rhs.cell_index)
         && (this->connection_transmissibility_factor == rhs.connection_transmissibility_factor)
         && (this->connection_d_factor == rhs.connection_d_factor)

opm/simulators/wells/PerfData.hpp

@@ -66,7 +66,9 @@ public:
         serializer(polymer_rates);
         serializer(brine_rates);
         serializer(prod_index);
-        serializer(micp_rates);
+        serializer(microbial_rates);
+        serializer(oxygen_rates);
+        serializer(urea_rates);
         serializer(cell_index);
         serializer(connection_transmissibility_factor);
         serializer(connection_d_factor);
@@ -97,7 +99,9 @@ public:
     std::vector<Scalar> polymer_rates{};
     std::vector<Scalar> brine_rates{};
     std::vector<Scalar> prod_index{};
-    std::vector<Scalar> micp_rates{};
+    std::vector<Scalar> microbial_rates{};
+    std::vector<Scalar> oxygen_rates{};
+    std::vector<Scalar> urea_rates{};
     std::vector<std::size_t> cell_index{};
     std::vector<Scalar> connection_transmissibility_factor{};
     std::vector<Scalar> connection_d_factor{};

opm/simulators/wells/SingleWellState.cpp

@@ -193,6 +193,24 @@ Scalar SingleWellState<Scalar>::sum_solvent_rates() const
     return this->sum_connection_rates(this->perf_data.solvent_rates);
 }
 
+template<class Scalar>
+Scalar SingleWellState<Scalar>::sum_microbial_rates() const
+{
+    return this->sum_connection_rates(this->perf_data.microbial_rates);
+}
+
+template<class Scalar>
+Scalar SingleWellState<Scalar>::sum_oxygen_rates() const
+{
+    return this->sum_connection_rates(this->perf_data.oxygen_rates);
+}
+
+template<class Scalar>
+Scalar SingleWellState<Scalar>::sum_urea_rates() const
+{
+    return this->sum_connection_rates(this->perf_data.urea_rates);
+}
+
 template<class Scalar>
 Scalar SingleWellState<Scalar>::sum_filtrate_rate() const
 {

opm/simulators/wells/SingleWellState.hpp

@@ -144,6 +144,9 @@ public:
     Scalar sum_solvent_rates() const;
     Scalar sum_polymer_rates() const;
     Scalar sum_brine_rates() const;
+    Scalar sum_microbial_rates() const;
+    Scalar sum_oxygen_rates() const;
+    Scalar sum_urea_rates() const;
 
     Scalar sum_filtrate_rate() const;
     Scalar sum_filtrate_total() const;

opm/simulators/wells/StandardWellConnections.cpp

@@ -756,7 +756,10 @@ std::tuple<typename StandardWellConnections<FluidSystem,Indices>::Eval,
            typename StandardWellConnections<FluidSystem,Indices>::Eval,
            typename StandardWellConnections<FluidSystem,Indices>::Eval>
 StandardWellConnections<FluidSystem,Indices>::
-connectionRatesMICP(const std::vector<EvalWell>& cq_s,
+connectionRatesMICP(Scalar& rate_m,
+                    Scalar& rate_o,
+                    Scalar& rate_u,
+                    const std::vector<EvalWell>& cq_s,
                     const std::variant<Scalar,EvalWell>& microbialConcentration,
                     const std::variant<Scalar,EvalWell>& oxygenConcentration,
                     const std::variant<Scalar,EvalWell>& ureaConcentration) const
@@ -769,6 +772,8 @@ connectionRatesMICP(const std::vector<EvalWell>& cq_s,
         cq_s_microbe *= std::get<EvalWell>(microbialConcentration);
     }
 
+    rate_m = cq_s_microbe.value();
+
     EvalWell cq_s_oxygen = cq_s[waterCompIdx];
     if (well_.isInjector()) {
         cq_s_oxygen *= std::get<Scalar>(oxygenConcentration);
@@ -776,6 +781,8 @@ connectionRatesMICP(const std::vector<EvalWell>& cq_s,
         cq_s_oxygen *= std::get<EvalWell>(oxygenConcentration);
     }
 
+    rate_o = cq_s_oxygen.value();
+
     EvalWell cq_s_urea = cq_s[waterCompIdx];
     if (well_.isInjector()) {
         cq_s_urea *= std::get<Scalar>(ureaConcentration);
@@ -783,6 +790,8 @@ connectionRatesMICP(const std::vector<EvalWell>& cq_s,
         cq_s_urea *= std::get<EvalWell>(ureaConcentration);
     }
 
+    rate_u = cq_s_urea.value();
+
     return {well_.restrictEval(cq_s_microbe),
             well_.restrictEval(cq_s_oxygen),
             well_.restrictEval(cq_s_urea)};

opm/simulators/wells/StandardWellConnections.hpp

@@ -124,7 +124,10 @@ public:
                           const std::variant<Scalar,EvalWell>& polymerConcentration) const;
 
     std::tuple<Eval,Eval,Eval>
-    connectionRatesMICP(const std::vector<EvalWell>& cq_s,
+    connectionRatesMICP(Scalar& rate_m,
+                        Scalar& rate_o,
+                        Scalar& rate_u,
+                        const std::vector<EvalWell>& cq_s,
                         const std::variant<Scalar,EvalWell>& microbialConcentration,
                         const std::variant<Scalar,EvalWell>& oxygenConcentration,
                         const std::variant<Scalar,EvalWell>& ureaConcentration) const;

opm/simulators/wells/StandardWell_impl.hpp

@@ -623,7 +623,10 @@ namespace Opm
             std::tie(connectionRates[perf][Indices::contiMicrobialEqIdx],
                      connectionRates[perf][Indices::contiOxygenEqIdx],
                      connectionRates[perf][Indices::contiUreaEqIdx]) =
-                this->connections_.connectionRatesMICP(cq_s,
+                this->connections_.connectionRatesMICP(perf_data.microbial_rates[perf],
+                                                       perf_data.oxygen_rates[perf],
+                                                       perf_data.urea_rates[perf],
+                                                       cq_s,
                                                        microbialConcentration,
                                                        oxygenConcentration,
                                                        ureaConcentration);

opm/simulators/wells/WellInterfaceGeneric.cpp

@@ -792,7 +792,7 @@ Scalar WellInterfaceGeneric<Scalar>::wurea_() const
 
     if (injectorType == InjectorType::WATER) {
         WellMICPProperties urea = this->well_ecl_.getMICPProperties();
-        const Scalar urea_injection_concentration = urea.m_ureaConcentration / 10.; //Dividing by scaling factor 10
+        const Scalar urea_injection_concentration = urea.m_ureaConcentration;
         return urea_injection_concentration;
     } else {
         // Not a water injection well => no urea.

opm/simulators/wells/WellInterface_impl.hpp

@@ -173,12 +173,6 @@ namespace Opm
       return 0.0;
     }
 
-    // The urea injection concentration is scaled down by a factor of 10, since its value
-    // can be much bigger than 1 (not doing this slows the simulations). The
-    // corresponding values are scaled accordingly in blackoilmicpmodules.hh when computing
-    // the reactions and also when writing the output files (vtk and eclipse format, i.e.,
-    // vtkblackoilmicpmodule.hh and ecloutputblackoilmodel.hh respectively).
-
     template<typename TypeTag>
     typename WellInterface<TypeTag>::Scalar
     WellInterface<TypeTag>::

opm/simulators/wells/WellState.cpp

@@ -577,6 +577,12 @@ WellState<Scalar>::report(const int* globalCellIdxMap,
             well.rates.set(rt::brine, ws.sum_brine_rates());
         }
 
+        if (pu.has_micp) {
+            well.rates.set(rt::microbial, ws.sum_microbial_rates());
+            well.rates.set(rt::oxygen, ws.sum_oxygen_rates());
+            well.rates.set(rt::urea, ws.sum_urea_rates());
+        }
+
         if (ws.producer) {
             well.rates.set(rt::alq, getALQ(wname));
         }
@@ -1123,6 +1129,14 @@ reportConnectionPressuresAndRates(const std::size_t well_index,
             connections[i].rates.set(rt::solvent, perf_data.solvent_rates[i]);
         }
     }
+
+    if (pu.has_micp) {
+        for (auto i = 0*num_perf_well; i < num_perf_well; ++i) {
+            connections[i].rates.set(rt::microbial, perf_data.microbial_rates[i]);
+            connections[i].rates.set(rt::oxygen, perf_data.oxygen_rates[i]);
+            connections[i].rates.set(rt::urea, perf_data.urea_rates[i]);
+        }
+    }
 }
 
 template <class Scalar>