first version of micp implementation in flow

opm/models/blackoil/blackoilboundaryratevector.hh

@@ -62,6 +62,7 @@ class BlackOilBoundaryRateVector : public GetPropType<TypeTag, Properties::RateV
     enum { conti0EqIdx = Indices::conti0EqIdx };
     enum { contiEnergyEqIdx = Indices::contiEnergyEqIdx };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
 
     static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
 
@@ -177,6 +178,12 @@ public:
             (*this)[Indices::contiPolymerEqIdx] = extQuants.volumeFlux(FluidSystem::waterPhaseIdx) * insideIntQuants.polymerConcentration();
         }
 
+        if (enableMICP) {
+            (*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();
+        }
+
         // make sure that the right mass conservation quantities are used
         LocalResidual::adaptMassConservationQuantities_(*this, insideIntQuants.pvtRegionIndex());
 

opm/models/blackoil/blackoilextensivequantities.hh

@@ -33,6 +33,7 @@
 #include "blackoilpolymermodules.hh"
 #include "blackoilenergymodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoilmicpmodules.hh"
 #include <opm/models/common/multiphasebaseextensivequantities.hh>
 
 namespace Opm {
@@ -55,6 +56,7 @@ class BlackOilExtensiveQuantities
     , public BlackOilPolymerExtensiveQuantities<TypeTag>
     , public BlackOilEnergyExtensiveQuantities<TypeTag>
     , public BlackOilDiffusionExtensiveQuantities<TypeTag, getPropValue<TypeTag, Properties::EnableDiffusion>()>
+    , public BlackOilMICPExtensiveQuantities<TypeTag>
 {
     using MultiPhaseParent = MultiPhaseBaseExtensiveQuantities<TypeTag>;
 

opm/models/blackoil/blackoilindices.hh

@@ -35,7 +35,7 @@ namespace Opm {
  *
  * \brief The primary variable and equation indices for the black-oil model.
  */
-template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset>
+template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned numMICPsV>
 struct BlackOilIndices
 {
     //! Number of phases active at all times
@@ -58,6 +58,9 @@ struct BlackOilIndices
     //! Shall energy be conserved?
     static const bool enableEnergy = numEnergyV > 0;
 
+    //! Is MICP involved?
+    static const bool enableMICP = numMICPsV > 0;
+
     //! Number of solvent components to be considered
     static const int numSolvents = enableSolvent ? numSolventsV : 0;
 
@@ -76,8 +79,11 @@ struct BlackOilIndices
     //! Number of salt equations to be considered
     static const int numBrine = enableBrine? 1 : 0;
 
+    //! Number of MICP components to be considered
+    static const int numMICPs = enableMICP ? numMICPsV : 0;
+
     //! The number of equations
-    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
+    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs;
 
     //! \brief returns the index of "active" component
     static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx)
@@ -122,17 +128,37 @@ struct BlackOilIndices
     static const int polymerMoleWeightIdx =
         numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000;
 
+    //! Index of the primary variable for the first MICP component
+    static const int microbialConcentrationIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the primary variable for the second MICP component
+    static const int oxygenConcentrationIdx =
+        numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the third MICP component
+    static const int ureaConcentrationIdx =
+        numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fourth MICP component
+    static const int biofilmConcentrationIdx =
+        numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fifth MICP component
+    static const int calciteConcentrationIdx =
+        numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000;
+
     //! Index of the primary variable for the foam
     static const int foamConcentrationIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs : -1000;
 
     //! Index of the primary variable for the brine
     static const int saltConcentrationIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numMICPs + numFoam : -1000;
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
 
 
     ////////
@@ -159,18 +185,38 @@ struct BlackOilIndices
     static const int contiPolymerMWEqIdx =
         numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000;
 
+    //! Index of the continuity equation for the first MICP component
+    static const int contiMicrobialEqIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents + numExtbos : -1000;
+
+    //! Index of the continuity equation for the second MICP component
+    static const int contiOxygenEqIdx =
+        numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the third MICP component
+    static const int contiUreaEqIdx =
+        numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fourth MICP component
+    static const int contiBiofilmEqIdx =
+        numMICPs > 3 ? contiUreaEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fifth MICP component
+    static const int contiCalciteEqIdx =
+        numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000;
+
     //! Index of the continuity equation for the foam component
     static const int contiFoamEqIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs : -1000;
 
     //! Index of the continuity equation for the salt water component
     static const int contiBrineEqIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000;
 
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000;
 };
 
 } // namespace Opm

opm/models/blackoil/blackoilintensivequantities.hh

@@ -36,6 +36,7 @@
 #include "blackoilbrinemodules.hh"
 #include "blackoilenergymodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 #include <opm/material/common/Valgrind.hpp>
 #include <dune/common/fmatrix.hh>
@@ -62,6 +63,7 @@ class BlackOilIntensiveQuantities
     , public BlackOilFoamIntensiveQuantities<TypeTag>
     , public BlackOilBrineIntensiveQuantities<TypeTag>
     , public BlackOilEnergyIntensiveQuantities<TypeTag>
+    , public BlackOilMICPIntensiveQuantities<TypeTag>
 {
     using ParentType = GetPropType<TypeTag, Properties::DiscIntensiveQuantities>;
     using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
@@ -85,6 +87,7 @@ class BlackOilIntensiveQuantities
     enum { enableTemperature = getPropValue<TypeTag, Properties::EnableTemperature>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
     enum { waterCompIdx = FluidSystem::waterCompIdx };
@@ -129,6 +132,7 @@ public:
 
         const auto& problem = elemCtx.problem();
         const auto& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        const auto& linearizationType = elemCtx.linearizationType();
 
         asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
 
@@ -170,7 +174,7 @@ public:
             }
         }
         if (gasEnabled && waterEnabled && !oilEnabled) {
-            Sg = 1.0 - Sw; 
+            Sg = 1.0 - Sw;
         }
 
         Valgrind::CheckDefined(Sg);
@@ -384,11 +388,19 @@ public:
         // deal with water induced rock compaction
         porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
 
+        // the MICP processes change the porosity
+        if (enableMICP){
+          Evaluation biofilm_ = priVars.makeEvaluation(Indices::biofilmConcentrationIdx, timeIdx, linearizationType);
+          Evaluation calcite_ = priVars.makeEvaluation(Indices::calciteConcentrationIdx, timeIdx, linearizationType);
+          porosity_ += - biofilm_ - calcite_;
+        }
+
         asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().zPvtUpdate_();
         asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
         asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
 
         // update the quantities which are required by the chosen
         // velocity model
@@ -474,6 +486,7 @@ private:
     friend BlackOilEnergyIntensiveQuantities<TypeTag>;
     friend BlackOilFoamIntensiveQuantities<TypeTag>;
     friend BlackOilBrineIntensiveQuantities<TypeTag>;
+    friend BlackOilMICPIntensiveQuantities<TypeTag>;
 
 
     Implementation& asImp_()

opm/models/blackoil/blackoillocalresidual.hh

@@ -36,6 +36,7 @@
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 
 namespace Opm {
@@ -88,6 +89,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
     using FoamModule = BlackOilFoamModule<TypeTag>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
 
 public:
     /*!
@@ -161,6 +163,9 @@ public:
 
         // deal with salt (if present)
         BrineModule::addStorage(storage, intQuants);
+
+        // deal with micp (if present)
+        MICPModule::addStorage(storage, intQuants);
     }
 
     /*!
@@ -208,6 +213,9 @@ public:
         // deal with salt (if present)
         BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
+        // deal with micp (if present)
+        MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
         DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
     }
 
@@ -222,6 +230,9 @@ public:
         // retrieve the source term intrinsic to the problem
         elemCtx.problem().source(source, elemCtx, dofIdx, timeIdx);
 
+        // deal with MICP (if present)
+        MICPModule::addSource(source, elemCtx, dofIdx, timeIdx);
+
         // scale the source term of the energy equation
         if (enableEnergy)
             source[Indices::contiEnergyEqIdx] *= getPropValue<TypeTag, Properties::BlackOilEnergyScalingFactor>();

opm/models/blackoil/blackoilmicpmodules.hh

@@ -0,0 +1,675 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \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
+
+#include "blackoilproperties.hh"
+#include <opm/models/io/vtkblackoilmicpmodule.hh>
+#include <opm/models/common/quantitycallbacks.hh>
+
+#if HAVE_ECL_INPUT
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/MICPpara.hpp>
+#endif
+
+#include <opm/material/common/Valgrind.hpp>
+#include <opm/material/common/Unused.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <string>
+
+namespace Opm {
+/*!
+ * \ingroup BlackOil
+ * \brief Contains the high level supplements required to extend the black oil
+ *        model by MICP.
+ */
+template <class TypeTag, bool enableMICPV = getPropValue<TypeTag, Properties::EnableMICP>()>
+class BlackOilMICPModule
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using ExtensiveQuantities = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Model = GetPropType<TypeTag, Properties::Model>;
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using EqVector = GetPropType<TypeTag, Properties::EqVector>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+
+    using Toolbox = MathToolbox<Evaluation>;
+
+    static constexpr unsigned microbialConcentrationIdx = Indices::microbialConcentrationIdx;
+    static constexpr unsigned oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
+    static constexpr unsigned ureaConcentrationIdx = Indices::ureaConcentrationIdx;
+    static constexpr unsigned biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
+    static constexpr unsigned calciteConcentrationIdx = Indices::calciteConcentrationIdx;
+    static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
+    static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
+    static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
+    static constexpr unsigned contiBiofilmEqIdx = Indices::contiBiofilmEqIdx;
+    static constexpr unsigned contiCalciteEqIdx = Indices::contiCalciteEqIdx;
+    static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+    static constexpr unsigned enableMICP = enableMICPV;
+
+    static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
+
+public:
+
+#if HAVE_ECL_INPUT
+    //
+     //* \brief Initialize all internal data structures needed by the MICP module
+     //
+    static void initFromState(const EclipseState& eclState)
+    {
+        // some sanity checks: if MICP is enabled, the MICP keyword must be
+        // present, if MICP is disabled the keyword must not be present.
+        if (enableMICP && !eclState.runspec().micp()) {
+            throw std::runtime_error("Non-trivial MICP treatment requested at compile time, but "
+                                     "the deck does not contain the MICP keyword");
+        }
+        else if (!enableMICP && eclState.runspec().micp()) {
+            throw std::runtime_error("MICP treatment disabled at compile time, but the deck "
+                                     "contains the MICP keyword");
+        }
+
+        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();
+                setMICPpara(MICPpara.getDensityBiofilm(),
+                           MICPpara.getDensityCalcite(),
+                           MICPpara.getDetachmentRate(),
+                           MICPpara.getCriticalPorosity(),
+                           MICPpara.getFittingFactor(),
+                           MICPpara.getHalfVelocityOxygen(),
+                           MICPpara.getHalfVelocityUrea(),
+                           MICPpara.getMaximumGrowthRate(),
+                           MICPpara.getMaximumUreaUtilization(),
+                           MICPpara.getMicrobialAttachmentRate(),
+                           MICPpara.getMicrobialDeathRate(),
+                           MICPpara.getMinimumPermeability(),
+                           MICPpara.getOxygenConsumptionFactor(),
+                           MICPpara.getYieldGrowthCoefficient(),
+                           MICPpara.getMaximumOxygenConcentration(),
+                           MICPpara.getMaximumUreaConcentration(),
+                           MICPpara.getToleranceBeforeClogging());
+        // obtain the porosity for the clamp in the blackoilnewtonmethod
+        phi_ = eclState.fieldProps().get_double("PORO");
+    }
+#endif
+
+    /*!
+     * \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] < MICPparaToleranceBeforeClogging())
+            throw std::logic_error("Clogging has been (almost) reached in at least one cell\n");
+    }
+
+    /*!
+     * \brief Specify the MICP properties a single region.
+     *
+     * The index of specified here must be in range [0, numSatRegions)
+     */
+    static void setMICPpara(const Scalar& MICPparaDensityBiofilm,
+                            const Scalar& MICPparaDensityCalcite,
+                            const Scalar& MICPparaDetachmentRate,
+                            const Scalar& MICPparaCriticalPorosity,
+                            const Scalar& MICPparaFittingFactor,
+                            const Scalar& MICPparaHalfVelocityOxygen,
+                            const Scalar& MICPparaHalfVelocityUrea,
+                            const Scalar& MICPparaMaximumGrowthRate,
+                            const Scalar& MICPparaMaximumUreaUtilization,
+                            const Scalar& MICPparaMicrobialAttachmentRate,
+                            const Scalar& MICPparaMicrobialDeathRate,
+                            const Scalar& MICPparaMinimumPermeability,
+                            const Scalar& MICPparaOxygenConsumptionFactor,
+                            const Scalar& MICPparaYieldGrowthCoefficient,
+                            const Scalar& MICPparaMaximumOxygenConcentration,
+                            const Scalar& MICPparaMaximumUreaConcentration,
+                            const Scalar& MICPparaToleranceBeforeClogging)
+    {
+        MICPparaDensityBiofilm_ = MICPparaDensityBiofilm;
+        MICPparaDensityCalcite_ = MICPparaDensityCalcite;
+        MICPparaDetachmentRate_ = MICPparaDetachmentRate;
+        MICPparaCriticalPorosity_ = MICPparaCriticalPorosity;
+        MICPparaFittingFactor_ = MICPparaFittingFactor;
+        MICPparaHalfVelocityOxygen_ = MICPparaHalfVelocityOxygen;
+        MICPparaHalfVelocityUrea_ = MICPparaHalfVelocityUrea;
+        MICPparaMaximumGrowthRate_ = MICPparaMaximumGrowthRate;
+        MICPparaMaximumUreaUtilization_ = MICPparaMaximumUreaUtilization;
+        MICPparaMicrobialAttachmentRate_ = MICPparaMicrobialAttachmentRate;
+        MICPparaMicrobialDeathRate_ = MICPparaMicrobialDeathRate;
+        MICPparaMinimumPermeability_ = MICPparaMinimumPermeability;
+        MICPparaOxygenConsumptionFactor_ = MICPparaOxygenConsumptionFactor;
+        MICPparaYieldGrowthCoefficient_ = MICPparaYieldGrowthCoefficient;
+        MICPparaMaximumOxygenConcentration_ = MICPparaMaximumOxygenConcentration;
+        MICPparaMaximumUreaConcentration_ = MICPparaMaximumUreaConcentration;
+        MICPparaToleranceBeforeClogging_ = MICPparaToleranceBeforeClogging;
+    }
+
+    /*!
+     * \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();
+    }
+
+    /*!
+     * \brief Register all MICP specific VTK and ECL output modules.
+     */
+    static void registerOutputModules(Model& model,
+                                      Simulator& simulator)
+    {
+        if (!enableMICP)
+            // MICP has been disabled at compile time
+            return;
+
+        model.addOutputModule(new VtkBlackOilMICPModule<TypeTag>(simulator));
+    }
+
+    static bool eqApplies(unsigned eqIdx)
+    {
+        if (!enableMICP)
+            return false;
+
+        // All MICP components are true here
+        return eqIdx == contiMicrobialEqIdx || eqIdx == contiOxygenEqIdx || eqIdx == contiUreaEqIdx || eqIdx == contiBiofilmEqIdx || eqIdx == contiCalciteEqIdx;
+    }
+
+    static Scalar eqWeight(unsigned eqIdx OPM_OPTIM_UNUSED)
+    {
+        assert(eqApplies(eqIdx));
+
+        // TODO: it may be beneficial to chose this differently.
+        return static_cast<Scalar>(1.0);
+    }
+
+    // must be called after water storage is computed
+    template <class LhsEval>
+    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;
+    }
+
+    static void computeFlux(RateVector& flux,
+                            const ElementContext& elemCtx,
+                            unsigned scvfIdx,
+                            unsigned timeIdx)
+
+    {
+        if (!enableMICP)
+            return;
+
+        const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+
+        const unsigned upIdx = extQuants.upstreamIndex(waterPhaseIdx);
+        const unsigned inIdx = extQuants.interiorIndex();
+        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());
+        }
+    }
+
+    // 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)
+
+    {
+        if (!enableMICP)
+            return;
+
+        // 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);
+
+          // compute water velocity from flux
+          Evaluation waterVolumeVelocity = extQuants.volumeFlux(waterPhaseIdx) / (K[0][0] * mobWater);
+          dpW = std::max(dpW, abs(waterVolumeVelocity));
+        }
+
+        // get the model parameters
+        Scalar k_a = MICPparaMicrobialAttachmentRate();
+        Scalar k_d = MICPparaMicrobialDeathRate();
+        Scalar rho_b = MICPparaDensityBiofilm();
+        Scalar rho_c = MICPparaDensityCalcite();
+        Scalar k_str = MICPparaDetachmentRate();
+        Scalar k_o = MICPparaHalfVelocityOxygen();
+        Scalar k_u = MICPparaHalfVelocityUrea() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp)
+        Scalar mu = MICPparaMaximumGrowthRate();
+        Scalar mu_u = MICPparaMaximumUreaUtilization() / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp)
+        Scalar Y_sb = MICPparaYieldGrowthCoefficient();
+        Scalar F = MICPparaOxygenConsumptionFactor();
+        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 MICPparaDensityBiofilm()
+    {
+        return MICPparaDensityBiofilm_;
+    }
+
+    static const Scalar MICPparaDensityCalcite()
+    {
+        return MICPparaDensityCalcite_;
+    }
+
+    static const Scalar MICPparaDetachmentRate()
+    {
+        return MICPparaDetachmentRate_;
+    }
+
+    static const Scalar MICPparaCriticalPorosity()
+    {
+        return MICPparaCriticalPorosity_;
+    }
+
+    static const Scalar MICPparaFittingFactor()
+    {
+        return MICPparaFittingFactor_;
+    }
+
+    static const Scalar MICPparaHalfVelocityOxygen()
+    {
+        return MICPparaHalfVelocityOxygen_;
+    }
+
+    static const Scalar MICPparaHalfVelocityUrea()
+    {
+        return MICPparaHalfVelocityUrea_;
+    }
+
+    static const Scalar MICPparaMaximumGrowthRate()
+    {
+        return MICPparaMaximumGrowthRate_;
+    }
+
+    static const Scalar MICPparaMaximumOxygenConcentration()
+    {
+        return MICPparaMaximumOxygenConcentration_;
+    }
+
+    static const Scalar MICPparaMaximumUreaConcentration()
+    {
+        return MICPparaMaximumUreaConcentration_ / 10.0;//Dividing by scaling factor 10 (see WellInterface_impl.hpp);
+    }
+
+    static const Scalar MICPparaMaximumUreaUtilization()
+    {
+        return MICPparaMaximumUreaUtilization_;
+    }
+
+    static const Scalar MICPparaMicrobialAttachmentRate()
+    {
+        return MICPparaMicrobialAttachmentRate_;
+    }
+
+    static const Scalar MICPparaMicrobialDeathRate()
+    {
+        return MICPparaMicrobialDeathRate_;
+    }
+
+    static const Scalar MICPparaMinimumPermeability()
+    {
+        return MICPparaMinimumPermeability_;
+    }
+
+    static const Scalar MICPparaOxygenConsumptionFactor()
+    {
+        return MICPparaOxygenConsumptionFactor_;
+    }
+
+    static const Scalar MICPparaToleranceBeforeClogging()
+    {
+        return MICPparaToleranceBeforeClogging_;
+    }
+
+    static const Scalar MICPparaYieldGrowthCoefficient()
+    {
+        return MICPparaYieldGrowthCoefficient_;
+    }
+
+    static const std::vector<Scalar> phi()
+    {
+        return phi_;
+    }
+
+private:
+    static Scalar MICPparaDensityBiofilm_;
+    static Scalar MICPparaDensityCalcite_;
+    static Scalar MICPparaDetachmentRate_;
+    static Scalar MICPparaCriticalPorosity_;
+    static Scalar MICPparaFittingFactor_;
+    static Scalar MICPparaHalfVelocityOxygen_;
+    static Scalar MICPparaHalfVelocityUrea_;
+    static Scalar MICPparaMaximumGrowthRate_;
+    static Scalar MICPparaMaximumUreaUtilization_;
+    static Scalar MICPparaMicrobialAttachmentRate_;
+    static Scalar MICPparaMicrobialDeathRate_;
+    static Scalar MICPparaMinimumPermeability_;
+    static Scalar MICPparaOxygenConsumptionFactor_;
+    static Scalar MICPparaYieldGrowthCoefficient_;
+    static Scalar MICPparaMaximumOxygenConcentration_;
+    static Scalar MICPparaMaximumUreaConcentration_;
+    static Scalar MICPparaToleranceBeforeClogging_;
+    static std::vector<Scalar> phi_;
+};
+
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaDensityBiofilm_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaDensityCalcite_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaDetachmentRate_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaCriticalPorosity_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaFittingFactor_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaHalfVelocityOxygen_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaHalfVelocityUrea_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMaximumGrowthRate_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMaximumUreaUtilization_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMicrobialAttachmentRate_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMicrobialDeathRate_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMinimumPermeability_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaOxygenConsumptionFactor_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaYieldGrowthCoefficient_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMaximumOxygenConcentration_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaMaximumUreaConcentration_;
+
+template <class TypeTag, bool enableMICPV>
+typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar
+BlackOilMICPModule<TypeTag, enableMICPV>::MICPparaToleranceBeforeClogging_;
+
+template <class TypeTag, bool enableMICPV>
+std::vector<typename BlackOilMICPModule<TypeTag, enableMICPV>::Scalar>
+BlackOilMICPModule<TypeTag, enableMICPV>::phi_;
+
+
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilMICPIntensiveQuantities
+ *
+ * \brief Provides the volumetric quantities required for the equations needed by the
+ *        MICP extension of the black-oil model.
+ */
+template <class TypeTag, bool enableMICPV = getPropValue<TypeTag, Properties::EnableMICP>()>
+class BlackOilMICPIntensiveQuantities
+{
+    using Implementation = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+
+    using MICPModule = BlackOilMICPModule<TypeTag>;
+
+    static constexpr int microbialConcentrationIdx = Indices::microbialConcentrationIdx;
+    static constexpr int oxygenConcentrationIdx = Indices::oxygenConcentrationIdx;
+    static constexpr int ureaConcentrationIdx = Indices::ureaConcentrationIdx;
+    static constexpr int biofilmConcentrationIdx = Indices::biofilmConcentrationIdx;
+    static constexpr int calciteConcentrationIdx = Indices::calciteConcentrationIdx;
+    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+
+
+public:
+
+    /*!
+     * \brief Update the intensive properties needed to handle MICP from the
+     *        primary variables
+     *
+     */
+    void MICPPropertiesUpdate_(const ElementContext& elemCtx,
+                                  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::MICPparaFittingFactor();
+        Scalar k_min = MICPModule::MICPparaMinimumPermeability();
+        Scalar phi_crit = MICPModule::MICPparaCriticalPorosity();
+
+        microbialConcentration_ = priVars.makeEvaluation(microbialConcentrationIdx, timeIdx, linearizationType);
+        oxygenConcentration_ = priVars.makeEvaluation(oxygenConcentrationIdx, timeIdx, linearizationType);
+        ureaConcentration_ = priVars.makeEvaluation(ureaConcentrationIdx, timeIdx, linearizationType);
+        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& microbialConcentration() const
+    { return microbialConcentration_; }
+
+    const Evaluation& oxygenConcentration() const
+    { return oxygenConcentration_; }
+
+    const Evaluation& ureaConcentration() const
+    { return ureaConcentration_; }
+
+    const Evaluation& biofilmConcentration() const
+    { return biofilmConcentration_; }
+
+    const Evaluation& calciteConcentration() const
+    { return calciteConcentration_; }
+
+
+protected:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+    Evaluation microbialConcentration_;
+    Evaluation oxygenConcentration_;
+    Evaluation ureaConcentration_;
+    Evaluation biofilmConcentration_;
+    Evaluation calciteConcentration_;
+
+};
+
+template <class TypeTag>
+class BlackOilMICPIntensiveQuantities<TypeTag, false>
+{
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+
+public:
+    void MICPPropertiesUpdate_(const ElementContext&,
+                                  unsigned,
+                                  unsigned)
+    { }
+
+    const Evaluation& microbialConcentration() const
+    { throw std::logic_error("microbialConcentration() called but MICP is disabled"); }
+
+    const Evaluation& oxygenConcentration() const
+    { throw std::logic_error("oxygenConcentration() called but MICP is disabled"); }
+
+    const Evaluation& ureaConcentration() const
+    { throw std::logic_error("ureaConcentration() called but MICP is disabled"); }
+
+    const Evaluation& biofilmConcentration() const
+    { throw std::logic_error("biofilmConcentration() called but MICP is disabled"); }
+
+    const Evaluation& calciteConcentration() const
+    { throw std::logic_error("calciteConcentration() called but MICP is disabled"); }
+};
+
+/*!
+ * \ingroup BlackOil
+ * \class Opm::BlackOilMICPExtensiveQuantities
+ *
+ * \brief Provides the MICP specific extensive quantities to the generic black-oil
+ *        module's extensive quantities.
+ */
+template <class TypeTag, bool enableMICPV = getPropValue<TypeTag, Properties::EnableMICP>()>
+class BlackOilMICPExtensiveQuantities
+{
+    using Implementation = GetPropType<TypeTag, Properties::ExtensiveQuantities>;
+
+private:
+    Implementation& asImp_()
+    { return *static_cast<Implementation*>(this); }
+
+};
+
+template <class TypeTag>
+class BlackOilMICPExtensiveQuantities<TypeTag, false>{};
+
+} // namespace Opm
+
+#endif

opm/models/blackoil/blackoilmodel.hh

@@ -47,6 +47,7 @@
 #include "blackoilbrinemodules.hh"
 #include "blackoilextbomodules.hh"
 #include "blackoildarcyfluxmodule.hh"
+#include "blackoilmicpmodules.hh"
 
 #include <opm/models/common/multiphasebasemodel.hh>
 #include <opm/models/io/vtkcompositionmodule.hh>
@@ -78,7 +79,8 @@ struct BlackOilModel { using InheritsFrom = std::tuple<VtkComposition,
                                                        VtkBlackOilPolymer,
                                                        VtkBlackOilSolvent,
                                                        VtkBlackOil,
-                                                       MultiPhaseBaseModel>; };
+                                                       MultiPhaseBaseModel,
+                                                       VtkBlackOilMICP>; };
 } // namespace TTag
 
 //! Set the local residual function
@@ -131,7 +133,8 @@ struct Indices<TypeTag, TTag::BlackOilModel>
                                getPropValue<TypeTag, Properties::EnableEnergy>(),
                                getPropValue<TypeTag, Properties::EnableFoam>(),
                                getPropValue<TypeTag, Properties::EnableBrine>(),
-                               /*PVOffset=*/0>; };
+                               /*PVOffset=*/0,
+                               getPropValue<TypeTag, Properties::EnableMICP>()>; };
 
 //! Set the fluid system to the black-oil fluid system by default
 template<class TypeTag>
@@ -158,6 +161,8 @@ template<class TypeTag>
 struct EnableFoam<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 template<class TypeTag>
 struct EnableBrine<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
+template<class TypeTag>
+struct EnableMICP<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 
 //! By default, the blackoil model is isothermal and does not conserve energy
 template<class TypeTag>
@@ -287,6 +292,7 @@ class BlackOilModel
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
 
 public:
     BlackOilModel(Simulator& simulator)
@@ -305,6 +311,7 @@ public:
         PolymerModule::registerParameters();
         EnergyModule::registerParameters();
         DiffusionModule::registerParameters();
+        MICPModule::registerParameters();
 
         // register runtime parameters of the VTK output modules
         VtkBlackOilModule<TypeTag>::registerParameters();
@@ -595,6 +602,7 @@ protected:
         SolventModule::registerOutputModules(*this, this->simulator_);
         PolymerModule::registerOutputModules(*this, this->simulator_);
         EnergyModule::registerOutputModules(*this, this->simulator_);
+        MICPModule::registerOutputModules(*this, this->simulator_);
 
         this->addOutputModule(new VtkBlackOilModule<TypeTag>(this->simulator_));
         this->addOutputModule(new VtkCompositionModule<TypeTag>(this->simulator_));

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -32,6 +32,7 @@
 
 #include <opm/models/utils/signum.hh>
 #include <opm/models/nonlinear/newtonmethod.hh>
+#include "blackoilmicpmodules.hh"
 
 #include <opm/material/common/Unused.hpp>
 
@@ -111,8 +112,10 @@ class BlackOilNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonM
     using PrimaryVariables = GetPropType<TypeTag, Properties::PrimaryVariables>;
     using EqVector = GetPropType<TypeTag, Properties::EqVector>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
     using Scalar = GetPropType<TypeTag, Properties::Scalar>;
     using Linearizer = GetPropType<TypeTag, Properties::Linearizer>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
 
     static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
 
@@ -248,6 +251,7 @@ protected:
         static constexpr bool enableFoam = Indices::foamConcentrationIdx >= 0;
         static constexpr bool enableBrine = Indices::saltConcentrationIdx >= 0;
         static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
+        static constexpr bool enableMICP = Indices::microbialConcentrationIdx >= 0;
 
         currentValue.checkDefined();
         Valgrind::CheckDefined(update);
@@ -259,7 +263,7 @@ protected:
         Scalar deltaSg = 0.0;
         Scalar deltaSs = 0.0;
 
-        if (Indices::waterEnabled) {
+        if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1) {
             deltaSw = update[Indices::waterSaturationIdx];
             deltaSo = -deltaSw;
         }
@@ -371,6 +375,22 @@ protected:
             // keep the temperature within given values
             if (enableEnergy && pvIdx == Indices::temperatureIdx)
                 nextValue[pvIdx] = std::clamp(nextValue[pvIdx], tempMin_, tempMax_);
+
+            // 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], 0.0, MICPModule::MICPparaDensityBiofilm());
+            if (enableMICP && pvIdx == Indices::oxygenConcentrationIdx)
+                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::MICPparaMaximumOxygenConcentration());
+            if (enableMICP && pvIdx == Indices::ureaConcentrationIdx)
+                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::MICPparaMaximumUreaConcentration());
+            if (enableMICP && pvIdx == Indices::biofilmConcentrationIdx)
+                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::MICPparaToleranceBeforeClogging());
+            if (enableMICP && pvIdx == Indices::calciteConcentrationIdx)
+                nextValue[pvIdx] = std::clamp(nextValue[pvIdx], 0.0, MICPModule::phi()[globalDofIdx] - MICPModule::MICPparaToleranceBeforeClogging());
         }
 
         // switch the new primary variables to something which is physically meaningful.

opm/models/blackoil/blackoilonephaseindices.hh

@@ -23,7 +23,7 @@
 /*!
  * \file
  *
- * \copydoc Ewoms::BlackOilTwoPhaseIndices
+ * \copydoc Ewoms::BlackOilOnePhaseIndices
  */
 #ifndef EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH
 #define EWOMS_BLACK_OIL_ONE_PHASE_INDICES_HH
@@ -37,7 +37,7 @@ namespace Opm {
  *
  * \brief The primary variable and equation indices for the black-oil model.
  */
-template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned canonicalCompIdx>
+template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned canonicalCompIdx, unsigned numMICPsV>
 struct BlackOilOnePhaseIndices
 {
     //! Is phase enabled or not
@@ -57,6 +57,9 @@ struct BlackOilOnePhaseIndices
     //! Shall energy be conserved?
     static const bool enableEnergy = numEnergyV > 0;
 
+    //! Is MICP involved?
+    static const bool enableMICP = numMICPsV > 0;
+
     //! Number of solvent components to be considered
     static const int numSolvents = enableSolvent ? numSolventsV : 0;
 
@@ -78,8 +81,11 @@ struct BlackOilOnePhaseIndices
     //! The number of fluid phases
     static const int numPhases = 1;
 
+    //! Number of MICP components to be considered
+    static const int numMICPs = enableMICP ? numMICPsV : 0;
+
     //! The number of equations
-    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
+    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs;
 
     //////////////////////////////
     // Primary variable indices
@@ -115,17 +121,37 @@ struct BlackOilOnePhaseIndices
     static const int polymerMoleWeightIdx =
         numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000;
 
+    //! Index of the primary variable for the first MICP component
+    static const int microbialConcentrationIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the primary variable for the second MICP component
+    static const int oxygenConcentrationIdx =
+        numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the third MICP component
+    static const int ureaConcentrationIdx =
+        numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fourth MICP component
+    static const int biofilmConcentrationIdx =
+        numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fifth MICP component
+    static const int calciteConcentrationIdx =
+        numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000;
+
     //! Index of the primary variable for the foam
     static const int foamConcentrationIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000;
 
     //! Index of the primary variable for the salt
     static const int saltConcentrationIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000;
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: - 1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: - 1000;
 
     //////////////////////
     // Equation indices
@@ -170,17 +196,37 @@ struct BlackOilOnePhaseIndices
     static const int contiPolymerMWEqIdx =
         numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000;
 
+    //! Index of the continuity equation for the first MICP component
+    static const int contiMicrobialEqIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the continuity equation for the second MICP component
+    static const int contiOxygenEqIdx =
+        numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the third MICP component
+    static const int contiUreaEqIdx =
+        numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fourth MICP component
+    static const int contiBiofilmEqIdx =
+        numMICPs > 3 ? contiUreaEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fifth MICP component
+    static const int contiCalciteEqIdx =
+        numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000;
+
     //! Index of the continuity equation for the foam component
     static const int contiFoamEqIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000;
 
     //! Index of the continuity equation for the salt component
     static const int contiBrineEqIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam : -1000;
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000;
 };
 
 } // namespace Opm

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -2,20 +2,16 @@
 // vi: set et ts=4 sw=4 sts=4:
 /*
   This file is part of the Open Porous Media project (OPM).
-
   OPM is free software: you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 2 of the License, or
   (at your option) any later version.
-
   OPM is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
-
   You should have received a copy of the GNU General Public License
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
-
   Consult the COPYING file in the top-level source directory of this
   module for the precise wording of the license and the list of
   copyright holders.
@@ -35,6 +31,7 @@
 #include "blackoilenergymodules.hh"
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
+#include "blackoilmicpmodules.hh"
 
 #include <opm/models/discretization/common/fvbaseprimaryvariables.hh>
 
@@ -105,6 +102,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { gasCompIdx = FluidSystem::gasCompIdx };
     enum { waterCompIdx = FluidSystem::waterCompIdx };
     enum { oilCompIdx = FluidSystem::oilCompIdx };
@@ -117,6 +115,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     using EnergyModule = BlackOilEnergyModule<TypeTag, enableEnergy>;
     using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
     using BrineModule = BlackOilBrineModule<TypeTag, enableBrine>;
+    using MICPModule = BlackOilMICPModule<TypeTag, enableMICP>;
 
     static_assert(numPhases == 3, "The black-oil model assumes three phases!");
     static_assert(numComponents == 3, "The black-oil model assumes three components!");
@@ -312,13 +311,13 @@ public:
             } else {
                 throw std::logic_error("For single-phase runs, only pure water is presently allowed.");
             }
-            
+
         }
         else if (primaryVarsMeaning() == Sw_po_Sg) {
             if (waterEnabled)
                 (*this)[waterSaturationIdx] = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
-            if (gasEnabled && waterEnabled && !oilEnabled) { 
+            if (gasEnabled && waterEnabled && !oilEnabled) {
-                //-> water-gas system   
+                //-> water-gas system
                 (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
             }
             else if (oilEnabled) {
@@ -603,13 +602,13 @@ public:
             Scalar Ssol = 0.0;
             if (enableSolvent)
                 Ssol =(*this) [Indices::solventSaturationIdx];
-            
+
             Scalar So = 1.0 - Sw - Sg - Ssol;
             Sw = std::min(std::max(Sw,0.0),1.0);
             So = std::min(std::max(So,0.0),1.0);
             Sg = std::min(std::max(Sg,0.0),1.0);
             Ssol = std::min(std::max(Ssol,0.0),1.0);
-            
+
             Scalar St = Sw + So + Sg + Ssol;
             Sw = Sw/St;
             Sg = Sg/St;
@@ -622,8 +621,8 @@ public:
             if (enableSolvent)
                 (*this)[Indices::solventSaturationIdx] = Ssol;
             return not(St==1);
-            
+
-        }else if (primaryVarsMeaning() == Sw_po_Rs) { 
+        }else if (primaryVarsMeaning() == Sw_po_Rs) {
             Scalar Ssol = 0.0;
             if (enableSolvent)
                 Ssol = (*this)[Indices::solventSaturationIdx];
@@ -639,7 +638,7 @@ public:
             if (waterEnabled)
                 (*this)[Indices::waterSaturationIdx]= Sw;
             if (enableSolvent)
-                (*this)[Indices::solventSaturationIdx] = Ssol;            
+                (*this)[Indices::solventSaturationIdx] = Ssol;
             return not(St==1);
         }else{
             assert(primaryVarsMeaning() == Sw_pg_Rv);
@@ -660,11 +659,11 @@ public:
                 (*this)[Indices::waterSaturationIdx]= Sw;
             if (enableSolvent)
                 (*this)[Indices::solventSaturationIdx] = Ssol;
-            return not(St==1);            
+            return not(St==1);
         }
-        assert(false); 
+        assert(false);
     }
-    
+
     BlackOilPrimaryVariables& operator=(const BlackOilPrimaryVariables& other) = default;
     BlackOilPrimaryVariables& operator=(Scalar value)
     {
@@ -749,6 +748,46 @@ private:
         return (*this)[Indices::temperatureIdx];
     }
 
+    Scalar microbialConcentration_() const
+    {
+        if (!enableMICP)
+            return 0.0;
+
+        return (*this)[Indices::microbialConcentrationIdx];
+    }
+
+    Scalar oxygenConcentration_() const
+    {
+        if (!enableMICP)
+            return 0.0;
+
+        return (*this)[Indices::oxygenConcentrationIdx];
+    }
+
+    Scalar ureaConcentration_() const
+    {
+        if (!enableMICP)
+            return 0.0;
+
+        return (*this)[Indices::ureaConcentrationIdx];
+    }
+
+    Scalar biofilmConcentration_() const
+    {
+        if (!enableMICP)
+            return 0.0;
+
+        return (*this)[Indices::biofilmConcentrationIdx];
+    }
+
+    Scalar calciteConcentration_() const
+    {
+        if (!enableMICP)
+            return 0.0;
+
+        return (*this)[Indices::calciteConcentrationIdx];
+    }
+
     template <class Container>
     void computeCapillaryPressures_(Container& result,
                                     Scalar So,

opm/models/blackoil/blackoilproperties.hh

@@ -58,6 +58,9 @@ struct EnableFoam { using type = UndefinedProperty; };
 //! Enable the ECL-blackoil extension for salt
 template<class TypeTag, class MyTypeTag>
 struct EnableBrine { using type = UndefinedProperty; };
+//! Enable the ECL-blackoil extension for MICP.
+template<class TypeTag, class MyTypeTag>
+struct EnableMICP { using type = UndefinedProperty; };
 
 
 //! Allow the spatial and temporal domains to exhibit non-constant temperature

opm/models/blackoil/blackoilratevector.hh

@@ -60,6 +60,7 @@ class BlackOilRateVector
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using FoamModule = BlackOilFoamModule<TypeTag>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
+    using MICPModule = BlackOilMICPModule<TypeTag>;
 
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
@@ -71,6 +72,7 @@ class BlackOilRateVector
     enum { enablePolymerMolarWeight = getPropValue<TypeTag, Properties::EnablePolymerMW>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     using Toolbox = MathToolbox<Evaluation>;
     using ParentType = Dune::FieldVector<Evaluation, numEq>;
 
@@ -144,6 +146,10 @@ public:
             throw std::logic_error("setMolarRate() not implemented for salt water");
         }
 
+        if ( enableMICP ) {
+            throw std::logic_error("setMolarRate() not implemented for MICP");
+        }
+
         // convert to "surface volume" if requested
         if (getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>()) {
             if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx)) {

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -37,7 +37,7 @@ namespace Opm {
  *
  * \brief The primary variable and equation indices for the black-oil model.
  */
-template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
+template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx, unsigned numMICPsV>
 struct BlackOilTwoPhaseIndices
 {
     //! Is phase enabled or not
@@ -57,6 +57,9 @@ struct BlackOilTwoPhaseIndices
     //! Shall energy be conserved?
     static const bool enableEnergy = numEnergyV > 0;
 
+    //! Is MICP involved?
+    static const bool enableMICP = numMICPsV > 0;
+
     //! Number of solvent components to be considered
     static const int numSolvents = enableSolvent ? numSolventsV : 0;
 
@@ -78,8 +81,11 @@ struct BlackOilTwoPhaseIndices
     //! The number of fluid phases
     static const int numPhases = 2;
 
+    //! Number of MICP components to be considered
+    static const int numMICPs = enableMICP ? numMICPsV : 0;
+
     //! The number of equations
-    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
+    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine + numMICPs;
 
     //////////////////////////////
     // Primary variable indices
@@ -115,17 +121,37 @@ struct BlackOilTwoPhaseIndices
     static const int polymerMoleWeightIdx =
         numPolymers > 1 ? polymerConcentrationIdx + 1 : -1000;
 
+    //! Index of the primary variable for the first MICP component
+    static const int microbialConcentrationIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the primary variable for the second MICP component
+    static const int oxygenConcentrationIdx =
+        numMICPs > 1 ? microbialConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the third MICP component
+    static const int ureaConcentrationIdx =
+        numMICPs > 2 ? oxygenConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fourth MICP component
+    static const int biofilmConcentrationIdx =
+        numMICPs > 3 ? ureaConcentrationIdx + 1 : -1000;
+
+    //! Index of the primary variable for the fifth MICP component
+    static const int calciteConcentrationIdx =
+        numMICPs > 4 ? biofilmConcentrationIdx + 1 : -1000;
+
     //! Index of the primary variable for the foam
     static const int foamConcentrationIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000;
 
     //! Index of the primary variable for the salt
     static const int saltConcentrationIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs + numFoam : -1000;
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine : - 1000;
 
     //////////////////////
     // Equation indices
@@ -188,17 +214,37 @@ struct BlackOilTwoPhaseIndices
     static const int contiPolymerMWEqIdx =
         numPolymers > 1 ? contiPolymerEqIdx + 1 : -1000;
 
+    //! Index of the continuity equation for the first MICP component
+    static const int contiMicrobialEqIdx =
+        enableMICP ? PVOffset + numPhases + numSolvents : -1000;
+
+    //! Index of the continuity equation for the second MICP component
+    static const int contiOxygenEqIdx =
+        numMICPs > 1 ? contiMicrobialEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the third MICP component
+    static const int contiUreaEqIdx =
+        numMICPs > 2 ? contiOxygenEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fourth MICP component
+    static const int contiBiofilmEqIdx =
+        numMICPs > 3 ? contiUreaEqIdx + 1 : -1000;
+
+    //! Index of the continuity equation for the fifth MICP component
+    static const int contiCalciteEqIdx =
+        numMICPs > 4 ? contiBiofilmEqIdx + 1 : -1000;
+
     //! Index of the continuity equation for the foam component
     static const int contiFoamEqIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs : -1000;
 
     //! Index of the continuity equation for the salt component
     static const int contiBrineEqIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numMICPs + numFoam : -1000;
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : -1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numMICPs + numFoam + numBrine: -1000;
 };
 
 } // namespace Opm

opm/models/io/vtkblackoilmicpmodule.hh

@@ -0,0 +1,264 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ * \copydoc Opm::VtkBlackOilMICPModule
+ */
+#ifndef EWOMS_VTK_BLACK_OIL_MICP_MODULE_HH
+#define EWOMS_VTK_BLACK_OIL_MICP_MODULE_HH
+
+#include <opm/material/densead/Math.hpp>
+
+#include "vtkmultiwriter.hh"
+#include "baseoutputmodule.hh"
+
+#include <opm/models/utils/propertysystem.hh>
+#include <opm/models/utils/parametersystem.hh>
+#include <opm/models/blackoil/blackoilproperties.hh>
+
+#include <dune/common/fvector.hh>
+
+#include <cstdio>
+
+namespace Opm::Properties {
+
+namespace TTag {
+
+// create new type tag for the VTK multi-phase output
+struct VtkBlackOilMICP {};
+
+} // namespace TTag
+
+// create the property tags needed for the MICP output module
+template<class TypeTag, class MyTypeTag>
+struct VtkWriteMicrobialConcentration { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
+struct VtkWriteOxygenConcentration { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
+struct VtkWriteUreaConcentration { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
+struct VtkWriteBiofilmConcentration { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
+struct VtkWriteCalciteConcentration { using type = UndefinedProperty; };
+
+// set default values for what quantities to output
+template<class TypeTag>
+struct VtkWriteMicrobialConcentration<TypeTag, TTag::VtkBlackOilMICP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct VtkWriteOxygenConcentration<TypeTag, TTag::VtkBlackOilMICP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct VtkWriteUreaConcentration<TypeTag, TTag::VtkBlackOilMICP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct VtkWriteBiofilmConcentration<TypeTag, TTag::VtkBlackOilMICP> { static constexpr bool value = true; };
+template<class TypeTag>
+struct VtkWriteCalciteConcentration<TypeTag, TTag::VtkBlackOilMICP> { static constexpr bool value = true; };
+
+} // namespace Opm::Properties
+
+namespace Opm {
+/*!
+ * \ingroup Vtk
+ *
+ * \brief VTK output module for the MICP model's related quantities.
+ */
+template <class TypeTag>
+class VtkBlackOilMICPModule : public BaseOutputModule<TypeTag>
+{
+    using ParentType = BaseOutputModule<TypeTag>;
+
+    using Simulator = GetPropType<TypeTag, Properties::Simulator>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using ElementContext = GetPropType<TypeTag, Properties::ElementContext>;
+
+    static const int vtkFormat = getPropValue<TypeTag, Properties::VtkOutputFormat>();
+    using VtkMultiWriter = ::Opm::VtkMultiWriter<GridView, vtkFormat>;
+
+    enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
+
+    using ScalarBuffer = typename ParentType::ScalarBuffer;
+
+public:
+    VtkBlackOilMICPModule(const Simulator& simulator)
+        : ParentType(simulator)
+    { }
+
+    /*!
+     * \brief Register all run-time parameters for the multi-phase VTK output
+     * module.
+     */
+    static void registerParameters()
+    {
+        if (!enableMICP)
+            return;
+
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteMicrobialConcentration,
+                             "Include the concentration of the microbial component in the water phase "
+                             "in the VTK output files");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteOxygenConcentration,
+                             "Include the concentration of the oxygen component in the water phase "
+                             "in the VTK output files");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteUreaConcentration,
+                             "Include the concentration of the urea component in the water phase "
+                             "in the VTK output files");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteBiofilmConcentration,
+                             "Include the biofilm volume fraction "
+                             "in the VTK output files");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteCalciteConcentration,
+                             "Include the calcite volume fraction "
+                             "in the VTK output files");
+    }
+
+    /*!
+     * \brief Allocate memory for the scalar fields we would like to
+     *        write to the VTK file.
+     */
+    void allocBuffers()
+    {
+        if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
+            return;
+
+        if (!enableMICP)
+            return;
+
+        if (microbialConcentrationOutput_())
+            this->resizeScalarBuffer_(microbialConcentration_);
+        if (oxygenConcentrationOutput_())
+            this->resizeScalarBuffer_(oxygenConcentration_);
+        if (ureaConcentrationOutput_())
+            this->resizeScalarBuffer_(ureaConcentration_);
+        if (biofilmConcentrationOutput_())
+            this->resizeScalarBuffer_(biofilmConcentration_);
+        if (calciteConcentrationOutput_())
+            this->resizeScalarBuffer_(calciteConcentration_);
+    }
+
+    /*!
+     * \brief Modify the internal buffers according to the intensive quantities relevant for
+     *        an element
+     */
+    void processElement(const ElementContext& elemCtx)
+    {
+        if (!EWOMS_GET_PARAM(TypeTag, bool, EnableVtkOutput))
+            return;
+
+        if (!enableMICP)
+            return;
+
+        for (unsigned dofIdx = 0; dofIdx < elemCtx.numPrimaryDof(/*timeIdx=*/0); ++dofIdx) {
+            const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);
+            unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, /*timeIdx=*/0);
+
+            if (microbialConcentrationOutput_())
+                microbialConcentration_[globalDofIdx] =
+                    scalarValue(intQuants.microbialConcentration());
+
+            if (oxygenConcentrationOutput_())
+                oxygenConcentration_[globalDofIdx] =
+                    scalarValue(intQuants.oxygenConcentration());
+
+            if (ureaConcentrationOutput_())
+                ureaConcentration_[globalDofIdx] =
+                    10 * scalarValue(intQuants.ureaConcentration());//Multypliging by scaling factor 10 (see WellInterface_impl.hpp)
+
+            if (biofilmConcentrationOutput_())
+                biofilmConcentration_[globalDofIdx] =
+                    scalarValue(intQuants.biofilmConcentration());
+
+            if (calciteConcentrationOutput_())
+                calciteConcentration_[globalDofIdx] =
+                    scalarValue(intQuants.calciteConcentration());
+
+        }
+    }
+
+    /*!
+     * \brief Add all buffers to the VTK output writer.
+     */
+    void commitBuffers(BaseOutputWriter& baseWriter)
+    {
+        VtkMultiWriter *vtkWriter = dynamic_cast<VtkMultiWriter*>(&baseWriter);
+        if (!vtkWriter)
+            return;
+
+        if (!enableMICP)
+            return;
+
+        if (microbialConcentrationOutput_())
+            this->commitScalarBuffer_(baseWriter, "microbial concentration", microbialConcentration_);
+
+        if (oxygenConcentrationOutput_())
+            this->commitScalarBuffer_(baseWriter, "oxygen concentration", oxygenConcentration_);
+
+        if (ureaConcentrationOutput_())
+            this->commitScalarBuffer_(baseWriter, "urea concentration", ureaConcentration_);
+
+        if (biofilmConcentrationOutput_())
+            this->commitScalarBuffer_(baseWriter, "biofilm fraction", biofilmConcentration_);
+
+        if (calciteConcentrationOutput_())
+            this->commitScalarBuffer_(baseWriter, "calcite fraction", calciteConcentration_);
+
+    }
+
+private:
+    static bool microbialConcentrationOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteMicrobialConcentration);
+        return val;
+    }
+
+    static bool oxygenConcentrationOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteOxygenConcentration);
+        return val;
+    }
+
+    static bool ureaConcentrationOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteUreaConcentration);
+        return val;
+    }
+
+    static bool biofilmConcentrationOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteBiofilmConcentration);
+        return val;
+    }
+
+    static bool calciteConcentrationOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteCalciteConcentration);
+        return val;
+    }
+
+    ScalarBuffer microbialConcentration_;
+    ScalarBuffer oxygenConcentration_;
+    ScalarBuffer ureaConcentration_;
+    ScalarBuffer biofilmConcentration_;
+    ScalarBuffer calciteConcentration_;
+};
+} // namespace Opm
+
+#endif

opm/models/io/vtkblackoilmodule.hh

@@ -225,8 +225,10 @@ public:
             const auto& primaryVars = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0);
 
             unsigned pvtRegionIdx = elemCtx.primaryVars(dofIdx, /*timeIdx=*/0).pvtRegionIndex();
-            Scalar SoMax = std::max(getValue(fs.saturation(oilPhaseIdx)),
+            Scalar SoMax = 0.0;
-                                    elemCtx.problem().maxOilSaturation(globalDofIdx));
+            if (FluidSystem::phaseIsActive(oilPhaseIdx))
+                SoMax = std::max(getValue(fs.saturation(oilPhaseIdx)),
+                               elemCtx.problem().maxOilSaturation(globalDofIdx));
 
             if (FluidSystem::phaseIsActive(gasPhaseIdx) && FluidSystem::phaseIsActive(oilPhaseIdx)) {
                 Scalar x_oG = getValue(fs.moleFraction(oilPhaseIdx, gasCompIdx));