enable salt precipitation

opm/models/blackoil/blackoilbrinemodules.hh

@@ -35,9 +35,13 @@
 #include <opm/material/common/IntervalTabulated2DFunction.hpp>
 
 #if HAVE_ECL_INPUT
-#include <opm/input/eclipse/Deck/Deck.hpp>
+#include <opm/parser/eclipse/Deck/Deck.hpp>
-#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
-#include <opm/input/eclipse/EclipseState/Tables/PvtwsaltTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/PvtwsaltTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/PermfactTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/SaltSolubilityTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/SimpleTable.hpp>
 #endif
 
 #include <opm/material/common/Valgrind.hpp>
@@ -80,6 +84,7 @@ class BlackOilBrineModule
     static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
     static constexpr unsigned enableBrine = enableBrineV;
+    static constexpr unsigned enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
 
     static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
     static constexpr unsigned numPhases = FluidSystem::numPhases;
@@ -103,7 +108,6 @@ public:
                                      "contains the BRINE keyword");
         }
 
-
         if (!eclState.runspec().phases().active(Phase::BRINE))
             return; // brine treatment is supposed to be disabled
 
@@ -126,6 +130,25 @@ public:
                 bdensityTable_[pvtRegionIdx].setXYContainers(c, bdensityTable);
             }
         }
+
+        if (enableSaltPrecipitation) {
+            const TableContainer& permfactTables = tableManager.getPermfactTables();
+            permfactTable_.resize(numPvtRegions);
+            for (size_t i = 0; i < permfactTables.size(); ++i) {
+                const PermfactTable& permfactTable = permfactTables.getTable<PermfactTable>(i);
+                permfactTable_[i].setXYContainers(permfactTable.getPorosityChangeColumn(), permfactTable.getPermeabilityMultiplierColumn());
+            }
+
+            const TableContainer& saltsolTables = tableManager.getSaltsolTables();
+            if (!saltsolTables.empty()) {
+                saltsolTable_.resize(numPvtRegions);
+                assert(numPvtRegions == saltsolTables.size());
+                for (unsigned pvtRegionIdx = 0; pvtRegionIdx < numPvtRegions; ++ pvtRegionIdx) {
+                    const SaltsolTable& saltsolTable = saltsolTables.getTable<SaltsolTable>(pvtRegionIdx );
+                    saltsolTable_[pvtRegionIdx] = saltsolTable.getSaltsolColumn().front();
+                }
+            }
+        }
     }
 #endif
 
@@ -139,7 +162,6 @@ public:
             return;
     }
 
-
     static bool primaryVarApplies(unsigned pvIdx)
     {
         if (!enableBrine)
@@ -222,7 +244,16 @@ public:
         const LhsEval massBrine = surfaceVolumeWater
                 * Toolbox::template decay<LhsEval>(fs.saltConcentration());
 
-        storage[contiBrineEqIdx] += massBrine;
+        if (enableSaltPrecipitation){
+            double saltDensity = 2170; // Solid salt density kg/m3
+            const LhsEval solidSalt =
+                          Toolbox::template decay<LhsEval>(intQuants.porosity())
+                          * saltDensity
+                          * Toolbox::template decay<LhsEval>(intQuants.saltSaturation());
+
+            storage[contiBrineEqIdx] += massBrine + solidSalt;
+        }
+        else { storage[contiBrineEqIdx] += massBrine;}
     }
 
     static void computeFlux(RateVector& flux,
@@ -234,7 +265,10 @@ public:
         if (!enableBrine)
             return;
 
+        const unsigned contiGasEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);
+        const unsigned contiOilEqIdx = Indices::conti0EqIdx + Indices::canonicalToActiveComponentIndex(FluidSystem::oilCompIdx);
         const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
+        const auto& intQuants = elemCtx.intensiveQuantities(scvfIdx, timeIdx);
 
         const unsigned upIdx = extQuants.upstreamIndex(FluidSystem::waterPhaseIdx);
         const unsigned inIdx = extQuants.interiorIndex();
@@ -245,12 +279,24 @@ public:
                     extQuants.volumeFlux(waterPhaseIdx)
                     *up.fluidState().invB(waterPhaseIdx)
                     *up.fluidState().saltConcentration();
+
+            if (enableSaltPrecipitation) {
+                // modify gas and oil flux for mobility change
+                flux[contiGasEqIdx] *= intQuants.permFactor();
+                flux[contiOilEqIdx] *= intQuants.permFactor();
+            }
         }
         else {
             flux[contiBrineEqIdx] =
                     extQuants.volumeFlux(waterPhaseIdx)
                     *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
                     *decay<Scalar>(up.fluidState().saltConcentration());
+
+            if (enableSaltPrecipitation) {
+                // modify gas and oil flux for mobility change
+                flux[contiGasEqIdx] *= decay<Scalar>(intQuants.permFactor());
+                flux[contiOilEqIdx] *= decay<Scalar>(intQuants.permFactor());
+            }
         }
     }
 
@@ -310,13 +356,40 @@ public:
         return bdensityTable_[pvtnumRegionIdx];
     }
 
+    static const TabulatedFunction& permfactTable(const ElementContext& elemCtx,
+                                                  unsigned scvIdx,
+                                                  unsigned timeIdx)
+    {
+        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
+        return permfactTable_[pvtnumRegionIdx];
+    }
+
+    static const Scalar saltsolTable(const ElementContext& elemCtx,
+                                                  unsigned scvIdx,
+                                                  unsigned timeIdx)
+    {
+        unsigned pvtnumRegionIdx = elemCtx.problem().pvtRegionIndex(elemCtx, scvIdx, timeIdx);
+        return saltsolTable_[pvtnumRegionIdx];
+    }
+
     static bool hasBDensityTables()
     {
         return !bdensityTable_.empty();
     }
 
+    static bool hasSaltsolTables()
+    {
+        return !saltsolTable_.empty();
+    }
+
+    static Scalar saltSol(unsigned regionIdx) {
+        return saltsolTable_[regionIdx];
+    }
+
 private:
     static std::vector<TabulatedFunction> bdensityTable_;
+    static std::vector<TabulatedFunction> permfactTable_;
+    static std::vector<Scalar> saltsolTable_;
     static std::vector<Scalar> referencePressure_;
 };
 
@@ -329,6 +402,14 @@ template <class TypeTag, bool enableBrineV>
 std::vector<typename BlackOilBrineModule<TypeTag, enableBrineV>::Scalar>
 BlackOilBrineModule<TypeTag, enableBrineV>::referencePressure_;
 
+template <class TypeTag, bool enableBrineV>
+std::vector<typename BlackOilBrineModule<TypeTag, enableBrineV>::Scalar>
+BlackOilBrineModule<TypeTag, enableBrineV>::saltsolTable_;
+
+template <class TypeTag, bool enableBrineV>
+std::vector<typename BlackOilBrineModule<TypeTag, enableBrineV>::TabulatedFunction>
+BlackOilBrineModule<TypeTag, enableBrineV>::permfactTable_;
+
 /*!
  * \ingroup BlackOil
  * \class Ewoms::BlackOilBrineIntensiveQuantities
@@ -354,8 +435,10 @@ class BlackOilBrineIntensiveQuantities
     enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
     static constexpr int saltConcentrationIdx = Indices::saltConcentrationIdx;
     static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
+    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
     static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
     static constexpr unsigned enableBrine = enableBrineV;
+    static constexpr unsigned enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
     static constexpr int contiBrineEqIdx = Indices::contiBrineEqIdx;
 
 public:
@@ -372,9 +455,42 @@ public:
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
 
         auto& fs = asImp_().fluidState_;
-        // set saltconcentration
-        fs.setSaltConcentration(priVars.makeEvaluation(saltConcentrationIdx, timeIdx));
 
+        if (enableSaltPrecipitation) {
+            const auto& saltsolTable = BrineModule::saltsolTable(elemCtx, dofIdx, timeIdx);
+            saltSolubility_ = saltsolTable;
+
+            if (priVars.primaryVarsMeaningBrine() == PrimaryVariables::Sp) {
+                saltSaturation_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
+                fs.setSaltConcentration(saltSolubility_);
+            }
+            else {
+                saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
+                fs.setSaltConcentration(saltConcentration_);
+                saltSaturation_ = 0.0;
+            }
+        }
+        else {
+            saltConcentration_ = priVars.makeEvaluation(saltConcentrationIdx, timeIdx);
+            fs.setSaltConcentration(saltConcentration_);
+        }
+    }
+
+    void saltPropertiesUpdate_(const ElementContext& elemCtx,
+                                  unsigned dofIdx,
+                                  unsigned timeIdx)
+    {
+        if (enableSaltPrecipitation) {
+            const Evaluation porosityFactor  = 1.0 - saltSaturation(); //phi/phi_0
+
+            const auto& permfactTable = BrineModule::permfactTable(elemCtx, dofIdx, timeIdx);
+
+            permFactor_ = permfactTable.eval(scalarValue(porosityFactor));
+            if (permFactor_ < 1 ) {
+                // adjust mobility for changing permeability
+                asImp_().mobility_[waterPhaseIdx] *= permFactor_;
+            }
+         }
     }
 
     const Evaluation& saltConcentration() const
@@ -383,12 +499,24 @@ public:
     const Evaluation& brineRefDensity() const
     { return refDensity_; }
 
+    const Evaluation& saltSaturation() const
+    { return saltSaturation_; }
+
+    Scalar saltSolubility() const
+    { return saltSolubility_; }
+
+    const Evaluation& permFactor() const
+    { return permFactor_; }
+
 protected:
     Implementation& asImp_()
     { return *static_cast<Implementation*>(this); }
 
     Evaluation saltConcentration_;
     Evaluation refDensity_;
+    Evaluation saltSaturation_;
+    Evaluation permFactor_;
+    Scalar saltSolubility_;
 
 };
 
@@ -405,18 +533,28 @@ public:
                                   unsigned)
     { }
 
+    void saltPropertiesUpdate_(const ElementContext&,
+                                  unsigned,
+                                  unsigned)
+    { }
+
     const Evaluation& saltConcentration() const
     { throw std::runtime_error("saltConcentration() called but brine are disabled"); }
 
     const Evaluation& brineRefDensity() const
     { throw std::runtime_error("brineRefDensity() called but brine are disabled"); }
 
+    const Evaluation& saltSaturation() const
+    { throw std::logic_error("saltSaturation() called but salt precipitation is disabled"); }
+
+    const Scalar saltSolubility() const
+    { throw std::logic_error("saltSolubility() called but salt precipitation is disabled"); }
+
+    const Evaluation& permFactor() const
+    { throw std::logic_error("permFactor() called but salt precipitation is disabled"); }
+
 };
 
-
-
-
-
 } // namespace Ewoms
 
 #endif

opm/models/blackoil/blackoilintensivequantities.hh

@@ -402,7 +402,7 @@ public:
         asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
         asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
-        asImp_().MICPPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().saltPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
 
         // update the quantities which are required by the chosen
         // velocity model

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -118,6 +118,7 @@ class BlackOilNewtonMethod : public GetPropType<TypeTag, Properties::DiscNewtonM
     using MICPModule = BlackOilMICPModule<TypeTag>;
 
     static const unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
+    static constexpr bool enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>();
 
 public:
     BlackOilNewtonMethod(Simulator& simulator) : ParentType(simulator)
@@ -369,8 +370,10 @@ protected:
                 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);
 
             // keep the salt concentration above 0
-            if (enableBrine && pvIdx == Indices::saltConcentrationIdx)
+            if (enableBrine && pvIdx == Indices::saltConcentrationIdx) {
+               if (!enableSaltPrecipitation || (enableSaltPrecipitation && currentValue.primaryVarsMeaningBrine() == PrimaryVariables::Cs))
                    nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0); 
+            }
 
             // keep the temperature within given values
             if (enableEnergy && pvIdx == Indices::temperatureIdx)

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -82,6 +82,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { waterSaturationIdx = Indices::waterSaturationIdx };
     enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
     enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
+    enum { saltConcentrationIdx  = Indices::saltConcentrationIdx };
 
     static const bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
     static const bool waterEnabled = Indices::waterEnabled;
@@ -101,6 +102,8 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
+    enum { enableSaltPrecipitation = getPropValue<TypeTag, Properties::EnableSaltPrecipitation>() };
+
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
     enum { gasCompIdx = FluidSystem::gasCompIdx };
@@ -128,6 +131,11 @@ public:
         OnePhase_p, // onephase case
     };
 
+    enum PrimaryVarsMeaningBrine {
+        Cs, // salt concentration
+        Sp, // (precipitated) salt saturation
+    };
+
     BlackOilPrimaryVariables()
         : ParentType()
     {
@@ -181,6 +189,17 @@ public:
     void setPrimaryVarsMeaning(PrimaryVarsMeaning newMeaning)
     { primaryVarsMeaning_ = newMeaning; }
 
+    PrimaryVarsMeaningBrine primaryVarsMeaningBrine() const
+    { return primaryVarsMeaningBrine_; }
+
+    /*!
+     * \brief Set the interpretation which should be applied to the switching primary
+     *        variables.
+     */
+
+    void setPrimaryVarsMeaningBrine(PrimaryVarsMeaningBrine newMeaning)
+    { primaryVarsMeaningBrine_ = newMeaning; }
+
     /*!
      * \copydoc ImmisciblePrimaryVariables::assignMassConservative
      */
@@ -273,6 +292,7 @@ public:
         bool oilPresent = FluidSystem::phaseIsActive(oilPhaseIdx)?(fluidState.saturation(oilPhaseIdx) > 0.0):false;
         static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
         bool onlyWater = FluidSystem::phaseIsActive(waterPhaseIdx)?(fluidState.saturation(waterPhaseIdx) > thresholdWaterFilledCell):false;
+        bool precipitatedSaltPresent = enableSaltPrecipitation?(fluidState.saltSaturation() > 0.0):false;
 
         // deal with the primary variables for the energy extension
         EnergyModule::assignPrimaryVars(*this, fluidState);
@@ -304,6 +324,13 @@ public:
                 primaryVarsMeaning_ = Sw_po_Sg;
         }
 
+        if (enableSaltPrecipitation){
+            if (precipitatedSaltPresent)
+                primaryVarsMeaningBrine_ = Sp;
+            else
+                primaryVarsMeaningBrine_ = Cs;
+        }
+
         // assign the actual primary variables
         if (primaryVarsMeaning() == OnePhase_p) {
             if (waterEnabled) {
@@ -311,7 +338,6 @@ public:
             } else {
                 throw std::logic_error("For single-phase runs, only pure water is presently allowed.");
             }
-
         }
         else if (primaryVarsMeaning() == Sw_po_Sg) {
             if (waterEnabled)
@@ -347,6 +373,15 @@ public:
                 (*this)[compositionSwitchIdx] = Rv;
         }
 
+        if (enableSaltPrecipitation) {
+            if (primaryVarsMeaningBrine() == Sp) {
+                (*this)[saltConcentrationIdx] = FsToolbox::value(fluidState.saltSaturation());
+            }
+            else {
+                (*this)[saltConcentrationIdx] = FsToolbox::value(fluidState.saltConcentration());
+            }
+        }
+
         checkDefined();
     }
 
@@ -379,6 +414,24 @@ public:
         if (waterEnabled)
             Sw = (*this)[Indices::waterSaturationIdx];
 
+        if (enableSaltPrecipitation) {
+            Scalar saltSolubility = BrineModule::saltSol(pvtRegionIndex());
+            if (primaryVarsMeaningBrine() == Sp) {
+                Scalar saltSat = (*this)[saltConcentrationIdx];
+                if (saltSat < -eps){ //precipitated salt dissappears
+                    setPrimaryVarsMeaningBrine(Cs);
+                    (*this)[saltConcentrationIdx] = saltSolubility; //set salt concentration to solubility limit
+                }
+            }
+            else if (primaryVarsMeaningBrine() == Cs) {
+                Scalar saltConc = (*this)[saltConcentrationIdx];
+                if (saltConc > saltSolubility + eps){ //salt concentration exceeds solubility limit
+                    setPrimaryVarsMeaningBrine(Sp);
+                    (*this)[saltConcentrationIdx] = 0.0;
+                }
+            }
+         }
+
         if (primaryVarsMeaning() == Sw_po_Sg) {
 
             // special case for cells with almost only water
@@ -816,6 +869,7 @@ private:
     }
 
     PrimaryVarsMeaning primaryVarsMeaning_;
+    PrimaryVarsMeaningBrine primaryVarsMeaningBrine_;
     unsigned short pvtRegionIdx_;
 };