Add H2SOL and CO2SOL model. That allows for dissolution of solvent in water

opm/models/blackoil/blackoilintensivequantities.hh

@@ -231,10 +231,12 @@ public:
 
         // deal with solvent
         if constexpr (enableSolvent) {
-            if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
-                So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
-            } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
-                Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+            if(priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
+                if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
+                    So -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+                } else if (FluidSystem::phaseIsActive(gasPhaseIdx)) {
+                    Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
+                }
             }
         }
 

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -310,7 +310,7 @@ protected:
             deltaSg = update[Indices::compositionSwitchIdx];
             deltaSo -= deltaSg;
         }
-        if (enableSolvent) {
+        if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
             deltaSs = update[Indices::solventSaturationIdx];
             deltaSo -= deltaSs;
         }
@@ -365,7 +365,13 @@ protected:
             }
             else if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
                 // solvent saturation updates are also subject to the Appleyard chop
-                delta *= satAlpha;
+                if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
+                    delta *= satAlpha;
+                } else {
+                    //Ensure Rssolw factor does not become negative
+                    if (delta > currentValue[Indices::solventSaturationIdx])
+                        delta = currentValue[Indices::solventSaturationIdx];
+                }
             }
             else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
                 // z fraction updates are also subject to the Appleyard chop
@@ -396,8 +402,10 @@ protected:
             nextValue[pvIdx] = currentValue[pvIdx] - delta;
 
             // keep the solvent saturation between 0 and 1
-            if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
-                nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
+            if (enableSolvent && pvIdx == Indices::solventSaturationIdx) {
+                if (currentValue.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss )
+                    nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
+            }
 
             // keep the z fraction between 0 and 1
             if (enableExtbo && pvIdx == Indices::zFractionIdx)

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -83,6 +83,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     enum { pressureSwitchIdx = Indices::pressureSwitchIdx };
     enum { compositionSwitchIdx = Indices::compositionSwitchIdx };
     enum { saltConcentrationIdx  = Indices::saltConcentrationIdx };
+    enum { solventSaturationIdx  = Indices::solventSaturationIdx };
 
     static constexpr bool compositionSwitchEnabled = Indices::compositionSwitchIdx >= 0;
     static constexpr bool waterEnabled = Indices::waterEnabled;
@@ -150,6 +151,12 @@ public:
         Disabled, // The primary variable is not used
     };
 
+    enum class SolventMeaning {
+        Ss, // solvent saturation
+        Rsolw, // dissolved solvent in water
+        Disabled, // The primary variable is not used
+    };
+
     BlackOilPrimaryVariables()
         : ParentType()
     {
@@ -167,6 +174,7 @@ public:
         Valgrind::SetUndefined(primaryVarsMeaningGas_);
         Valgrind::SetUndefined(primaryVarsMeaningPressure_);
         Valgrind::SetUndefined(primaryVarsMeaningBrine_);
+        Valgrind::SetUndefined(primaryVarsMeaningSolvent_);
 
         pvtRegionIdx_ = 0;
     }
@@ -184,6 +192,7 @@ public:
         result.primaryVarsMeaningGas_ = GasMeaning::Rv;
         result.primaryVarsMeaningPressure_ = PressureMeaning::Pg;
         result.primaryVarsMeaningWater_ = WaterMeaning::Rsw;
+        result.primaryVarsMeaningSolvent_ = SolventMeaning::Ss;
         for (size_t i = 0; i < result.size(); ++i) {
             result[i] = i+1;
         }
@@ -260,6 +269,18 @@ public:
     void setPrimaryVarsMeaningBrine(BrineMeaning newMeaning)
     { primaryVarsMeaningBrine_ = newMeaning; }
 
+
+    SolventMeaning primaryVarsMeaningSolvent() const
+    { return primaryVarsMeaningSolvent_; }
+
+    /*!
+     * \brief Set the interpretation which should be applied to the switching primary
+     *        variables.
+     */
+
+    void setPrimaryVarsMeaningSolvent(SolventMeaning newMeaning)
+    { primaryVarsMeaningSolvent_ = newMeaning; }
+
     /*!
      * \copydoc ImmisciblePrimaryVariables::assignMassConservative
      */
@@ -480,7 +501,6 @@ public:
             default:
                 throw std::logic_error("No valid primary variable selected for composision");
         }
-        checkDefined();
     }
 
     /*!
@@ -542,6 +562,31 @@ public:
             }
         }
 
+        // if solvent saturation disappeares:  Ss (Solvent saturation) -> Rsolw (solvent dissolved in water)
+        // if solvent saturation appears: Rsolw (solvent dissolved in water) ->  Ss (Solvent saturation)
+        // Scalar rsolw = 0.0; // not needed at the moment since we dont allow for vapwat in combination with rsolw
+        if constexpr (enableSolvent) {
+            if (SolventModule::isSolubleInWater()) {
+                Scalar p = (*this)[pressureSwitchIdx]; // cap-pressure?
+                Scalar solLimit = SolventModule::solubilityLimit(pvtRegionIndex(), T , p, saltConcentration);
+                if (primaryVarsMeaningSolvent() == SolventMeaning::Ss) {
+                    Scalar solSat = (*this)[solventSaturationIdx];
+                    if (solSat < -eps){ //solvent dissappears
+                        setPrimaryVarsMeaningSolvent(SolventMeaning::Rsolw);
+                        (*this)[solventSaturationIdx] = solLimit; //set rsolw to solubility limit
+
+                    }
+                }
+                else if (primaryVarsMeaningSolvent() == SolventMeaning::Rsolw) {
+                    Scalar rsolw = (*this)[solventSaturationIdx];
+                    if (rsolw > solLimit + eps){ //solvent appears as phase
+                        setPrimaryVarsMeaningSolvent(SolventMeaning::Ss);
+                        (*this)[solventSaturationIdx] = 0.0;
+                    }
+                }
+            }
+        }
+
         // keep track if any meaning has changed
         bool changed = false;
 
@@ -826,7 +871,7 @@ public:
             sg = (*this)[Indices::compositionSwitchIdx];
 
         Scalar ssol = 0.0;
-        if constexpr (enableSolvent)
+        if (primaryVarsMeaningSolvent() == SolventMeaning::Ss)
             ssol =(*this) [Indices::solventSaturationIdx];
 
         Scalar so = 1.0 - sw - sg - ssol;
@@ -843,7 +888,7 @@ public:
             (*this)[Indices::waterSwitchIdx] = sw;
         if (primaryVarsMeaningGas() == GasMeaning::Sg)
             (*this)[Indices::compositionSwitchIdx] = sg;
-        if constexpr (enableSolvent)
+        if (primaryVarsMeaningSolvent() == SolventMeaning::Ss)
             (*this) [Indices::solventSaturationIdx] = ssol;
 
         return !(st==1);
@@ -877,6 +922,7 @@ public:
         Valgrind::CheckDefined(primaryVarsMeaningGas_);
         Valgrind::CheckDefined(primaryVarsMeaningPressure_);
         Valgrind::CheckDefined(primaryVarsMeaningBrine_);
+        Valgrind::CheckDefined(primaryVarsMeaningSolvent_);
 
         Valgrind::CheckDefined(pvtRegionIdx_);
 #endif // NDEBUG
@@ -891,6 +937,7 @@ public:
         serializer(primaryVarsMeaningPressure_);
         serializer(primaryVarsMeaningGas_);
         serializer(primaryVarsMeaningBrine_);
+        serializer(primaryVarsMeaningSolvent_);
         serializer(pvtRegionIdx_);
     }
 
@@ -901,6 +948,7 @@ public:
                this->primaryVarsMeaningPressure_ == rhs.primaryVarsMeaningPressure_ &&
                this->primaryVarsMeaningGas_ == rhs.primaryVarsMeaningGas_ &&
                this->primaryVarsMeaningBrine_ == rhs.primaryVarsMeaningBrine_ &&
+               this->primaryVarsMeaningSolvent_ == rhs.primaryVarsMeaningSolvent_ &&
                this->pvtRegionIdx_ == rhs.pvtRegionIdx_;
     }
 
@@ -913,10 +961,11 @@ private:
 
     Scalar solventSaturation_() const
     {
-        if constexpr (enableSolvent)
-            return (*this)[Indices::solventSaturationIdx];
-        else
-            return 0.0;
+        if constexpr (enableSolvent) {
+            if ( primaryVarsMeaningSolvent() == SolventMeaning::Ss)
+                return (*this)[Indices::solventSaturationIdx];
+        }
+        return 0.0;
     }
 
     Scalar zFraction_() const
@@ -1033,6 +1082,7 @@ private:
     PressureMeaning primaryVarsMeaningPressure_;
     GasMeaning primaryVarsMeaningGas_;
     BrineMeaning primaryVarsMeaningBrine_;
+    SolventMeaning primaryVarsMeaningSolvent_;
     unsigned short pvtRegionIdx_;
 };
 

opm/models/blackoil/blackoilsolventmodules.hh

@@ -78,9 +78,12 @@ class BlackOilSolventModule
     using EqVector = GetPropType<TypeTag, Properties::EqVector>;
     using RateVector = GetPropType<TypeTag, Properties::RateVector>;
     using Indices = GetPropType<TypeTag, Properties::Indices>;
-
     using Toolbox = MathToolbox<Evaluation>;
     using SolventPvt = typename BlackOilSolventParams<Scalar>::SolventPvt;
+    using Co2GasPvt = typename BlackOilSolventParams<Scalar>::Co2GasPvt;
+    using H2GasPvt = typename BlackOilSolventParams<Scalar>::H2GasPvt;
+    using BrineCo2Pvt = typename BlackOilSolventParams<Scalar>::BrineCo2Pvt;
+    using BrineH2Pvt = typename BlackOilSolventParams<Scalar>::BrineH2Pvt;
 
     using TabulatedFunction = typename BlackOilSolventParams<Scalar>::TabulatedFunction;
 
@@ -90,6 +93,7 @@ class BlackOilSolventModule
     static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
     static constexpr unsigned numPhases = FluidSystem::numPhases;
     static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
+    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
 
 
 public:
@@ -111,7 +115,26 @@ public:
         if (!eclState.runspec().phases().active(Phase::SOLVENT))
             return; // solvent treatment is supposed to be disabled
 
-        params_.solventPvt_.initFromState(eclState, schedule);
+        params_.co2sol_ = eclState.runspec().co2Sol();
+        params_.h2sol_ = eclState.runspec().h2Sol();
+
+        if (isCO2Sol() && isH2Sol()) {
+             throw std::runtime_error("CO2SOL and H2SOL can not be used together");
+        }
+
+        if (isCO2Sol() || isH2Sol() ) {
+            if (isCO2Sol()) {
+                params_.co2GasPvt_.initFromState(eclState, schedule);
+                params_.brineCo2Pvt_.initFromState(eclState, schedule);
+            } else { 
+                params_.h2GasPvt_.initFromState(eclState, schedule);
+                params_.brineH2Pvt_.initFromState(eclState, schedule);
+            }
+            if (eclState.getSimulationConfig().hasDISGASW()) {
+                params_.rsSolw_active_ = true;
+            }
+        } else 
+            params_.solventPvt_.initFromState(eclState, schedule);
 
         const auto& tableManager = eclState.getTableManager();
         // initialize the objects which deal with the SSFN keyword
@@ -413,12 +436,25 @@ public:
                         Toolbox::template decay<LhsEval>(intQuants.porosity())
                         * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
                         * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
+                if (isSolubleInWater()) {
+                    storage[contiSolventEqIdx] += Toolbox::template decay<LhsEval>(intQuants.porosity())
+                        * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
+                        * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(waterPhaseIdx)) 
+                        * Toolbox::template decay<LhsEval>(intQuants.rsSolw());
+                }
             }
             else {
                 storage[contiSolventEqIdx] +=
                         Toolbox::template decay<LhsEval>(intQuants.porosity())
                         * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
                         * Toolbox::template decay<LhsEval>(intQuants.solventDensity());
+                if (isSolubleInWater()) {
+                    storage[contiSolventEqIdx] += Toolbox::template decay<LhsEval>(intQuants.porosity())
+                        * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(waterPhaseIdx))
+                        * Toolbox::template decay<LhsEval>(intQuants.fluidState().density(waterPhaseIdx))
+                        * Toolbox::template decay<LhsEval>(intQuants.rsSolw());
+                }
+                         
             }
         }
     }
@@ -437,7 +473,7 @@ public:
             const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
 
             if constexpr (blackoilConserveSurfaceVolume) {
-                if (upIdx == inIdx)
+                if (upIdx == inIdx) 
                     flux[contiSolventEqIdx] =
                             extQuants.solventVolumeFlux()
                             *up.solventInverseFormationVolumeFactor();
@@ -445,6 +481,20 @@ public:
                     flux[contiSolventEqIdx] =
                             extQuants.solventVolumeFlux()
                             *decay<Scalar>(up.solventInverseFormationVolumeFactor());
+
+
+                if (isSolubleInWater()) {
+                    if (upIdx == inIdx) 
+                        flux[contiSolventEqIdx] =
+                                extQuants.volumeFlux(waterPhaseIdx)
+                                * up.fluidState().invB(waterPhaseIdx)
+                                * up.rsSolw();
+                    else
+                        flux[contiSolventEqIdx] =
+                                extQuants.volumeFlux(waterPhaseIdx)
+                                *decay<Scalar>(up.fluidState().invB(waterPhaseIdx))
+                                *decay<Scalar>(up.rsSolw());
+                }
             }
             else {
                 if (upIdx == inIdx)
@@ -455,6 +505,20 @@ public:
                     flux[contiSolventEqIdx] =
                             extQuants.solventVolumeFlux()
                             *decay<Scalar>(up.solventDensity());
+
+
+                if (isSolubleInWater()) {
+                    if (upIdx == inIdx) 
+                        flux[contiSolventEqIdx] =
+                                extQuants.volumeFlux(waterPhaseIdx)
+                                * up.fluidState().density(waterPhaseIdx)
+                                * up.rsSolw();
+                    else
+                        flux[contiSolventEqIdx] =
+                                extQuants.volumeFlux(waterPhaseIdx)
+                                *decay<Scalar>(up.fluidState().density(waterPhaseIdx))
+                                *decay<Scalar>(up.rsSolw());
+                }
             }
         }
     }
@@ -463,10 +527,21 @@ public:
      * \brief Assign the solvent specific primary variables to a PrimaryVariables object
      */
     static void assignPrimaryVars(PrimaryVariables& priVars,
-                                  Scalar solventSaturation)
+                                  Scalar solventSaturation,
+                                  Scalar solventRsw)
     {
-        if constexpr (enableSolvent)
+        if constexpr (!enableSolvent) {
+            priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Disabled);
+            return;
+        }
+        // Determine the meaning of the solvent primary variables
+        if (solventSaturation > 0 || !isSolubleInWater()) {
+            priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Ss);
             priVars[solventSaturationIdx] = solventSaturation;
+        } else {
+            priVars.setPrimaryVarsMeaningSolvent(PrimaryVariables::SolventMeaning::Rsolw);
+            priVars[solventSaturationIdx] = solventRsw;
+        } 
     }
 
     /*!
@@ -530,7 +605,30 @@ public:
     }
 
     static const SolventPvt& solventPvt()
-    { return params_.solventPvt_; }
+    { 
+        return params_.solventPvt_;    
+    }
+
+
+    static const Co2GasPvt& co2GasPvt()
+    { 
+            return params_.co2GasPvt_;
+    }
+
+    static const H2GasPvt& h2GasPvt()
+    { 
+            return params_.h2GasPvt_;
+    }
+
+    static const BrineCo2Pvt& brineCo2Pvt()
+    { 
+            return params_.brineCo2Pvt_;
+    }
+
+    static const BrineH2Pvt& brineH2Pvt()
+    { 
+            return params_.brineH2Pvt_;
+    }
 
     static const TabulatedFunction& ssfnKrg(const ElementContext& elemCtx,
                                             unsigned scvIdx,
@@ -633,6 +731,34 @@ public:
         return params_.isMiscible_;
     }
 
+    template <class Value>
+    static const Value solubilityLimit(unsigned pvtIdx, const Value& temperature, const Value& pressure, const Value& saltConcentration) 
+    {
+        if (!isSolubleInWater())
+            return 0.0;
+
+        assert(isCO2Sol() || isH2Sol());
+        if (isCO2Sol())
+            return brineCo2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
+        else
+            return brineH2Pvt().saturatedGasDissolutionFactor(pvtIdx, temperature, pressure, saltConcentration);
+    }
+
+    static bool isSolubleInWater() 
+    {
+        return params_.rsSolw_active_;
+    } 
+
+    static bool isCO2Sol() 
+    {
+        return params_.co2sol_;
+    }
+
+    static bool isH2Sol() 
+    {
+        return params_.h2sol_;
+    }  
+
 private:
     static BlackOilSolventParams<Scalar> params_; // the krg(Fs) column of the SSFN table
 };
@@ -683,8 +809,13 @@ public:
                                   unsigned timeIdx)
     {
         const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+
         auto& fs = asImp_().fluidState_;
-        solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
+        solventSaturation_ = 0.0;
+        if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
+            solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
+        }
+
         hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
 
         // apply a cut-off. Don't waste calculations if no solvent
@@ -712,6 +843,16 @@ public:
         auto& fs = asImp_().fluidState_;
         fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
 
+
+        // update rsSolw. This needs to be done after the pressure is defined in the fluid state. 
+        rsSolw_ = 0.0;
+        const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
+        if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Ss) {
+            rsSolw_ = SolventModule::solubilityLimit(asImp_().pvtRegionIndex(), fs.temperature(waterPhaseIdx), fs.pressure(waterPhaseIdx), fs.saltConcentration());
+        } else if (priVars.primaryVarsMeaningSolvent() == PrimaryVariables::SolventMeaning::Rsolw) {
+            rsSolw_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx, elemCtx.linearizationType());
+        }
+
         solventMobility_ = 0.0;
 
         // apply a cut-off. Don't waste calculations if no solvent
@@ -726,7 +867,6 @@ public:
 
             // compute capillary pressure for miscible fluid
             const auto& problem = elemCtx.problem();
-            const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
             Evaluation pgMisc = 0.0;
             std::array<Evaluation, numPhases> pC;
             const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
@@ -824,18 +964,58 @@ public:
                            unsigned timeIdx)
     {
         const auto& iq = asImp_();
-        const auto& fs = iq.fluidState();
-        const auto& solventPvt = SolventModule::solventPvt();
-
         unsigned pvtRegionIdx = iq.pvtRegionIndex();
-        solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
-        const Evaluation& T = fs.temperature(gasPhaseIdx);
-        const Evaluation& p = fs.pressure(gasPhaseIdx);
-        solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+        const Evaluation& T = iq.fluidState().temperature(gasPhaseIdx);
+        const Evaluation& p = iq.fluidState().pressure(gasPhaseIdx);
+
+        const Evaluation rv = 0.0;
+        const Evaluation rvw = 0.0;
+        if (SolventModule::isCO2Sol() || SolventModule::isH2Sol() ){
+            if (SolventModule::isCO2Sol()) {
+                const auto& co2gasPvt = SolventModule::co2GasPvt();
+                solventInvFormationVolumeFactor_ = co2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
+                solventRefDensity_ = co2gasPvt.gasReferenceDensity(pvtRegionIdx);
+                solventViscosity_ = co2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw);
+
+                const auto& brineCo2Pvt = SolventModule::brineCo2Pvt();
+                auto& fs = asImp_().fluidState_;
+                const auto& bw = brineCo2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
+
+                const auto denw = bw*brineCo2Pvt.waterReferenceDensity(pvtRegionIdx)
+                        + rsSolw()*bw*brineCo2Pvt.gasReferenceDensity(pvtRegionIdx);
+                fs.setDensity(waterPhaseIdx, denw);
+                fs.setInvB(waterPhaseIdx, bw);
+                Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
+                const auto& muWat = fs.viscosity(waterPhaseIdx);
+                const auto& muWatEff = brineCo2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
+                mobw *= muWat / muWatEff;
+            } else { 
+                const auto& h2gasPvt = SolventModule::h2GasPvt();
+                solventInvFormationVolumeFactor_ = h2gasPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rv, rvw);
+                solventRefDensity_ = h2gasPvt.gasReferenceDensity(pvtRegionIdx);
+                solventViscosity_ = h2gasPvt.viscosity(pvtRegionIdx, T, p, rv, rvw); 
+            
+                const auto& brineH2Pvt = SolventModule::brineH2Pvt();
+                auto& fs = asImp_().fluidState_;
+                const auto& bw = brineH2Pvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p, rsSolw());
+
+                const auto denw = bw*brineH2Pvt.waterReferenceDensity(pvtRegionIdx)
+                        + rsSolw()*bw*brineH2Pvt.gasReferenceDensity(pvtRegionIdx);
+                fs.setDensity(waterPhaseIdx, denw);
+                fs.setInvB(waterPhaseIdx, bw);
+                Evaluation& mobw = asImp_().mobility_[waterPhaseIdx];
+                const auto& muWat = fs.viscosity(waterPhaseIdx);
+                const auto& muWatEff = brineH2Pvt.viscosity(pvtRegionIdx, T, p, rsSolw());
+                mobw *= muWat / muWatEff;
+            }
+        } else {
+            const auto& solventPvt = SolventModule::solventPvt();
+            solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
+            solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
+            solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
+        }   
 
         solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
-        solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
-
         effectiveProperties(elemCtx, scvIdx, timeIdx);
 
         solventMobility_ /= solventViscosity_;
@@ -846,6 +1026,9 @@ public:
     const Evaluation& solventSaturation() const
     { return solventSaturation_; }
 
+    const Evaluation& rsSolw() const
+    { return rsSolw_; }
+
     const Evaluation& solventDensity() const
     { return solventDensity_; }
 
@@ -1036,6 +1219,7 @@ protected:
 
     Evaluation hydrocarbonSaturation_;
     Evaluation solventSaturation_;
+    Evaluation rsSolw_;
     Evaluation solventDensity_;
     Evaluation solventViscosity_;
     Evaluation solventMobility_;
@@ -1071,6 +1255,9 @@ public:
     const Evaluation& solventSaturation() const
     { throw std::runtime_error("solventSaturation() called but solvents are disabled"); }
 
+    const Evaluation& rsSolw() const
+    { throw std::runtime_error("rsSolw() called but solvents are disabled"); }
+
     const Evaluation& solventDensity() const
     { throw std::runtime_error("solventDensity() called but solvents are disabled"); }
 

opm/models/blackoil/blackoilsolventparams.hh

@@ -29,6 +29,11 @@
 #define EWOMS_BLACK_OIL_SOLVENT_PARAMS_HH
 
 #include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
+#include <opm/material/fluidsystems/blackoilpvt/Co2GasPvt.hpp>
+#include <opm/material/fluidsystems/blackoilpvt/H2GasPvt.hpp>
+#include <opm/material/fluidsystems/blackoilpvt/BrineCo2Pvt.hpp>
+#include <opm/material/fluidsystems/blackoilpvt/BrineH2Pvt.hpp>
+
 #include <opm/material/common/Tabulated1DFunction.hpp>
 
 namespace Opm {
@@ -40,6 +45,19 @@ struct BlackOilSolventParams {
 
     using SolventPvt = ::Opm::SolventPvt<Scalar>;
     SolventPvt solventPvt_;
+
+    using Co2GasPvt = ::Opm::Co2GasPvt<Scalar>;
+    Co2GasPvt co2GasPvt_;
+
+    using H2GasPvt = ::Opm::H2GasPvt<Scalar>;
+    H2GasPvt h2GasPvt_;
+    
+    using BrineCo2Pvt = ::Opm::BrineCo2Pvt<Scalar>;
+    BrineCo2Pvt brineCo2Pvt_;
+    
+    using BrineH2Pvt = ::Opm::BrineH2Pvt<Scalar>;
+    BrineH2Pvt brineH2Pvt_;
+
     std::vector<TabulatedFunction> ssfnKrg_; // the krg(Fs) column of the SSFN table
     std::vector<TabulatedFunction> ssfnKrs_; // the krs(Fs) column of the SSFN table
     std::vector<TabulatedFunction> sof2Krn_; // the krn(Sn) column of the SOF2 table
@@ -55,6 +73,9 @@ struct BlackOilSolventParams {
     std::vector<TabulatedFunction> tlPMixTable_; // the tlpmixpa(Po) column of the TLPMIXPA table
 
     bool isMiscible_;
+    bool rsSolw_active_ = false;
+    bool co2sol_;
+    bool h2sol_;
 
     /*!
      * \brief Specify the number of satuation regions.

opm/models/io/vtkblackoilsolventmodule.hh

@@ -53,6 +53,8 @@ struct VtkBlackOilSolvent {};
 template<class TypeTag, class MyTypeTag>
 struct VtkWriteSolventSaturation { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
+struct VtkWriteSolventRsw { using type = UndefinedProperty; };
+template<class TypeTag, class MyTypeTag>
 struct VtkWriteSolventDensity { using type = UndefinedProperty; };
 template<class TypeTag, class MyTypeTag>
 struct VtkWriteSolventViscosity { using type = UndefinedProperty; };
@@ -63,6 +65,8 @@ struct VtkWriteSolventMobility { using type = UndefinedProperty; };
 template<class TypeTag>
 struct VtkWriteSolventSaturation<TypeTag, TTag::VtkBlackOilSolvent> { static constexpr bool value = true; };
 template<class TypeTag>
+struct VtkWriteSolventRsw<TypeTag, TTag::VtkBlackOilSolvent> { static constexpr bool value = true; };
+template<class TypeTag>
 struct VtkWriteSolventDensity<TypeTag, TTag::VtkBlackOilSolvent> { static constexpr bool value = true; };
 template<class TypeTag>
 struct VtkWriteSolventViscosity<TypeTag, TTag::VtkBlackOilSolvent> { static constexpr bool value = true; };
@@ -112,6 +116,9 @@ public:
         EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventSaturation,
                              "Include the \"saturation\" of the solvent component "
                              "in the VTK output files");
+        EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventRsw,
+                             "Include the \"dissolved volume in water\" of the solvent component "
+                             "in the VTK output files");
         EWOMS_REGISTER_PARAM(TypeTag, bool, VtkWriteSolventDensity,
                              "Include the \"density\" of the solvent component "
                              "in the VTK output files");
@@ -137,6 +144,8 @@ public:
 
         if (solventSaturationOutput_())
             this->resizeScalarBuffer_(solventSaturation_);
+        if (solventRswOutput_())
+            this->resizeScalarBuffer_(solventRsw_);
         if (solventDensityOutput_())
             this->resizeScalarBuffer_(solventDensity_);
         if (solventViscosityOutput_())
@@ -166,6 +175,10 @@ public:
                 solventSaturation_[globalDofIdx] =
                     Toolbox::scalarValue(intQuants.solventSaturation());
 
+            if (solventRswOutput_())
+                solventRsw_[globalDofIdx] =
+                    Toolbox::scalarValue(intQuants.rsSolw());
+
             if (solventDensityOutput_())
                 solventDensity_[globalDofIdx] =
                     Toolbox::scalarValue(intQuants.solventDensity());
@@ -195,6 +208,9 @@ public:
         if (solventSaturationOutput_())
             this->commitScalarBuffer_(baseWriter, "saturation_solvent", solventSaturation_);
 
+        if (solventRswOutput_())
+            this->commitScalarBuffer_(baseWriter, "dissolved_solvent", solventRsw_);
+
         if (solventDensityOutput_())
             this->commitScalarBuffer_(baseWriter, "density_solvent", solventDensity_);
 
@@ -212,6 +228,12 @@ private:
         return val;
     }
 
+    static bool solventRswOutput_()
+    {
+        static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventRsw);
+        return val;
+    }
+
     static bool solventDensityOutput_()
     {
         static bool val = EWOMS_GET_PARAM(TypeTag, bool, VtkWriteSolventDensity);
@@ -231,6 +253,7 @@ private:
     }
 
     ScalarBuffer solventSaturation_;
+    ScalarBuffer solventRsw_;
     ScalarBuffer solventDensity_;
     ScalarBuffer solventViscosity_;
     ScalarBuffer solventMobility_;