Fix/biomin issue#37

exercises/exercise-biomineralization/README.md

@@ -4,10 +4,10 @@ The aim of this exercise is to get a first glimpse at the _DuMu<sup>x</sup>_ way
 
 ## Problem set-up
 
-The domain has a size of 20 x 15 m and contains a sealing aquitard in the middle. The aquitard is interrupted by a "fault zone" and thereby connects the upper drinking water aquifer and the lower CO<sub>2<sub>-storage aquifer. Initially, the domain is fully saturated with water and biofilm is present in the lower CO<sub>2<sub>-storage aquifer. Calcium and urea are injected in the upper drinking water aquifer by means of a Neumann boundary condition. The remaining parts of the upper boundary and the entire lower boundary are modelled with Neumann no-flow conditions, while on the right-hand side a Dirichlet boundary conditions is assigned, which fixes there the the initial values.
+The domain has a size of 20 x 15 m and contains a sealing aquitard in the middle. The aquitard is interrupted by a "fault zone" and thereby connects the upper drinking water aquifer and the lower CO<sub>2</sub>-storage aquifer. Initially, the domain is fully saturated with water and biofilm is present in the lower CO<sub>2</sub>-storage aquifer. Calcium and urea are injected in the upper drinking water aquifer by means of a Neumann boundary condition. The remaining parts of the upper boundary and the entire lower boundary are modelled with Neumann no-flow conditions, while on the right-hand side a Dirichlet boundary conditions is assigned, which fixes there the the initial values.
 
 
-Disclaimer: Please note, that this is not a realistic scenario. No one would think of storing gaseous CO<sub>2<sub> in this subcritical setting.
+Disclaimer: Please note, that this is not a realistic scenario. No one would think of storing gaseous CO<sub>2</sub> in this subcritical setting.
 
 ![](../extradoc/exercisebiomin_setup.png)
 
@@ -50,15 +50,25 @@ $`\displaystyle mass_{biofilm} = \rho_{biofilm} * \phi_{biofilm}`$
 
 Next, we want to implement the rate of ureolysis. This can be done with the following simplified equation:
 
-$`\displaystyle  r_{urea} = k_{urease} * Z_{urease,biofilm} * m_{urea} / (K_{urea} + m_{urea})`$,
+$`\displaystyle  r_{urea} = k_{urease} * Z_{u,b} * m_{urea} / (K_{urea} + m_{urea})`$,
 
 where $`\displaystyle  r_{urea}`$ is the rate of ureolysis, $`\displaystyle  k_{urease}`$ the urease enzyme activity, $`\displaystyle m_{urea}`$ the molality of urea and $`\displaystyle K_{urea}`$ the half-saturation constant for urea for the ureolysis rate.
 
-Note, that the urease concentration $`\displaystyle Z_{urease,biofilm}`$ has a kinetic term of urease production per biofilm :
+Note, that the urease concentration $`\displaystyle Z_{u,b}`$ has a kinetic term of urease production per biofilm :
 
-$`\displaystyle Z_{urease,biofilm} = k_{urease,biofilm} * mass_{biofilm}`$
+$`\displaystyle Z_{u,b} = k_{u,b} * mass_{biofilm}`$
 
-The last step is defining the source term for each component according to the chemical reaction equations:
+For the sake of simplicity, we consider the rate of calcite precipitation $`\displaystyle r_{prec}`$ to be equal to the rate of ureolysis $`\displaystyle r_{urea}`$:
+
+$`\displaystyle r_{prec} = r_{urea}`$
+
+Note that this assumption is a substantial simplification of the geochemistry, representing a state after sufficient ureolysis has occurred to raise the pH enough to initiate calcite precipitation.
+For a more detailed geochemistry, the exact speciation of inorganic carbon into $`\mathrm{CO_{2}}`$, $`\mathrm{H_{2}CO_{3}}`$, $`\mathrm{HCO_{3}^{-}}`$, and $`\mathrm{CO_{3}^{2-}}`$ would have to be considered.
+This would require to calculate activity coefficients and solve the coupled laws of mass actions of the dissociation reaction of the inorganic carbon species as well as the one for ammonia and ammonium, which are all coupled by the solution pH, or rather the occurrence of $`\mathrm{H^{+}}`$ in all those reactions.
+As such, this would require to solve a coupled system of equations to balance the geochemistry.
+Finally, the product of the activities of carbonate and calcium would need to be compared to the calcite solubility product to determine the saturation state of the solution with respect to calcite and to ultimately calculate the precipitation rate of calcite $`\displaystyle r_{prec}`$
+
+The last step is defining the source term $`\displaystyle q`$ for each component based on the now defined reaction rates $`\displaystyle r_{urea}`$ and $`\displaystyle r_{prec}`$ according to the chemical reaction equations:
 
 $`\displaystyle \mathrm{CO(NH_{2})_{2} + 2 H_{2}O  \rightarrow  2 NH_{3} + H_{2}CO_{3}}`$
 
@@ -137,16 +147,16 @@ make exercise_biomin
 
 ### 5. CO<sub>2</sub> injection to test aquitard integrity
 
-Now, the sealed aquitard is tested with a CO<sub>2<sub>-Injection into the lower CO<sub>2<sub>-storage aquifer.
+Now, the sealed aquitard is tested with a CO<sub>2</sub>-Injection into the lower CO<sub>2</sub>-storage aquifer.
 
 __Task:__
 
-Implement a new boundary condition on the left boundary, injecting CO<sub>2<sub> from 2 m to 3 m from the bottom. Make sure, that the injection time for the calcium and urea is finished. You can use the predefined value `gasFlux` directly and divide it by the molar mass of CO<sub>2<sub>.
+Implement a new boundary condition on the left boundary, injecting CO<sub>2</sub> from 2 m to 3 m from the bottom. Make sure, that the injection time for the calcium and urea is finished. You can use the predefined value `gasFlux` directly and divide it by the molar mass of CO<sub>2</sub>.
 Run two simulations and compare them side by side by creating two input files, or overwriting the input file in the command line:
 ```bash
 ./exercise_biomin -Problem.Name biominNoUrea -Injection.ConcUrea 0
 ```
-The result for the biomineralization process during the CO<sub>2<sub> injection should look like this:
+The result for the biomineralization process during the CO<sub>2</sub> injection should look like this:
 
 ![](../extradoc/exercisebiomin_injectionFinal.png)
 
@@ -196,7 +206,7 @@ Note: As both the Kozeny-Carman and the Power-Law relation use the same paramete
     }
 ```
 
-What is the effect of the exchanged permeability calculation on the results, especially the leakage of CO<sub>2<sub>? What if the exponent would be smaller, e.g. $`\displaystyle \eta=2`$, which would mean that the precipitation is less efficient in sealing the leakage?
+What is the effect of the exchanged permeability calculation on the results, especially the leakage of CO<sub>2</sub>? What if the exponent would be smaller, e.g. $`\displaystyle \eta=2`$, which would mean that the precipitation is less efficient in sealing the leakage?
 You can again run two simulations and compare them side by side by creating two input files, or overwriting the input file parameter in the command line:
 ```bash
 ./exercise_biomin -Problem.Name biominPowerLawExponent2 -PowerLaw.Exponent 2.0
@@ -204,11 +214,11 @@ You can again run two simulations and compare them side by side by creating two
 
 ### 7. Use tabulated CO<sub>2</sub> values instead of SimpleCO2
 
-So far we have been using a simplified component for CO<sub>2</sub>, which is based on the ideal gas law. Due to the conditions present in this exercise this is not too inaccurate, but for real applications of CO<sub>2<sub> storage changes to the model are required. We use tabulated data for density and enthalpy of CO<sub>2<sub>, accessed through `GeneratedCO2Tables::CO2Tables` and `Components::CO2` from DuMu<sup>x</sup>.
+So far we have been using a simplified component for CO<sub>2</sub>, which is based on the ideal gas law. Due to the conditions present in this exercise this is not too inaccurate, but for real applications of CO<sub>2</sub> storage changes to the model are required. We use tabulated data for density and enthalpy of CO<sub>2</sub>, accessed through `GeneratedCO2Tables::CO2Tables` and `Components::CO2` from DuMu<sup>x</sup>.
 
 __Task:__
 
-The CO<sub>2<sub> component is used in the fluidsystem, which is defined in `properties.hh`. Replace the component `SimpleCO2` with `CO2` defined in `dumux/material/components/co2.hh`, with a CO<sub>2<sub> table as an additional template parameter. Use the the table defined in `dumux/material/components/defaultco2table.hh`, noting the different namespace. Take care to include the appropriate headers.
+The CO<sub>2</sub> component is used in the fluidsystem, which is defined in `properties.hh`. Replace the component `SimpleCO2` with `CO2` defined in `dumux/material/components/co2.hh`, with a CO<sub>2</sub> table as an additional template parameter. Use the the table defined in `dumux/material/components/defaultco2table.hh`, noting the different namespace. Take care to include the appropriate headers.
 
 ```c++
 #include <dumux/material/components/co2.hh> //!< CO2 component for use with tabulated values

exercises/exercise-biomineralization/chemistry/simplebiominreactions.hh

@@ -25,7 +25,6 @@
 #define DUMUX_BIOMIN_REACTIONS_HH
 
 #include <dumux/common/parameters.hh>
-// #include <dumux-course/exercises/exercise-biomineralization/components/urea.hh>
 
 namespace Dumux {
 
@@ -45,8 +44,8 @@ public:
     SimpleBiominReactions()
     {   //ureolysis kinetic parameters
         kub_ = getParam<Scalar>("UreolysisCoefficients.Kub");
-        kurease_ = getParam<Scalar>("UreolysisCoefficients.Kurease");
-        ku_ = getParam<Scalar>("UreolysisCoefficients.Ku");
+        kUrease_ = getParam<Scalar>("UreolysisCoefficients.KUrease");
+        kUrea_ = getParam<Scalar>("UreolysisCoefficients.KUrea");
     }
 
     static constexpr int liquidPhaseIdx = FluidSystem::liquidPhaseIdx;
@@ -71,7 +70,7 @@ public:
      */
     static Scalar moleFracToMolality(Scalar moleFracX, Scalar moleFracSalinity, Scalar moleFracCTot)
     {
-        Scalar molalityX = moleFracX / (1 - moleFracSalinity - moleFracCTot) / FluidSystem::molarMass(H2OIdx);
+        const Scalar molalityX = moleFracX / (1 - moleFracSalinity - moleFracCTot) / FluidSystem::molarMass(H2OIdx);
         return molalityX;
     }
 
@@ -84,8 +83,8 @@ public:
     void reactionSource(NumEqVector &q, const VolumeVariables &volVars)
     {
         //   define and compute some parameters for convenience:
-        Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
-        Scalar densityBiofilm = volVars.solidComponentDensity(BiofilmPhaseIdx);
+        const Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
+        const Scalar densityBiofilm = volVars.solidComponentDensity(BiofilmPhaseIdx);
 
         Scalar volFracBiofilm = volVars.solidVolumeFraction(BiofilmPhaseIdx);
 
@@ -95,17 +94,16 @@ public:
         // TODO: dumux-course-task
         // implement mass of biofilm
 
-        Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
+        const Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
                                                  xwCa,
                                                  volVars.moleFraction(liquidPhaseIdx,CO2Idx));  // [mol_urea/kg_H2O]
 
         // TODO: dumux-course-task
         // compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
-        Scalar rurea = 0.0;
+        const Scalar rurea = 0.0;
 
         // compute/set dissolution and precipitation rate of calcite
-        Scalar rprec = 0.0;
-        rprec = rurea;
+        const Scalar rprec = 0.0;
 
         // set source terms
         // TODO: dumux-course-task
@@ -121,8 +119,8 @@ public:
 private:
     // urease parameters
     Scalar kub_;
-    Scalar kurease_;
-    Scalar ku_;
+    Scalar kUrease_;
+    Scalar kUrea_;
 };
 
 } //end namespace Dumux

exercises/exercise-biomineralization/components/urea.hh

-// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
-// vi: set et ts=4 sw=4 sts=4:
-/*****************************************************************************
- *   See the file COPYING for full copying permissions.                      *
- *                                                                           *
- *   This program 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 3 of the License, or       *
- *   (at your option) any later version.                                     *
- *                                                                           *
- *   This program 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 this program.  If not, see <http://www.gnu.org/licenses/>.   *
- *****************************************************************************/
-/*!
- * \file
- * \ingroup Components
- * \brief A class for the urea component properties.
- */
-#ifndef DUMUX_UREA_HH
-#define DUMUX_UREA_HH
-
-#include <dumux/material/components/base.hh>
-
-namespace Dumux {
-namespace Components {
-
-/*!
- * \ingroup Components
- * \brief A class for the urea properties
- */
-template <class Scalar>
-class Urea
-: public Components::Base<Scalar, Urea<Scalar> >
-{
-public:
-    /*!
-     * \brief A human readable name for the urea.
-     */
-    static std::string name()
-    { return "urea"; }
-
-    /*!
-     * \brief The molar mass in \f$\mathrm{[kg/mol]}\f$ of urea.
-     */
-    static constexpr Scalar molarMass()
-    { return 0.0606; } // kg/mol
-
-};
-
-} // end namespace Components
-} // end namespace Dumux
-
-#endif

exercises/exercise-biomineralization/fluidsystems/biomin.hh

@@ -34,7 +34,7 @@
 #include <dumux/material/components/tabulatedcomponent.hh>
 
 #include <dumux/material/components/calciumion.hh>
-#include "../components/urea.hh"
+#include <dumux/material/components/urea.hh>
 
 #include <dumux/material/binarycoefficients/brine_co2.hh>
 

exercises/exercise-biomineralization/params.input

@@ -41,8 +41,8 @@ RhoBiofilm = 6.9 # [kg/m³] density of biofilm
 
 [UreolysisCoefficients]
 Kub = 8.81e-3 # [kg_urease/kg_bio] (max: 0.01)
-Kurease = 1000 # [mol_urea/(kg_urease s)]
-Ku = 0.355 # [mol/kgH2O] Lauchnor et al. 2014
+KUrease = 1000 # [mol_urea/(kg_urease s)]
+KUrea = 0.355 # [mol/kgH2O] Lauchnor et al. 2014
 
 [Brine]
 Salinity = 0.1

exercises/solution/exercise-biomineralization/chemistry/simplebiominreactions.hh

@@ -44,8 +44,8 @@ public:
     SimpleBiominReactions()
     {   //ureolysis kinetic parameters
         kub_ = getParam<Scalar>("UreolysisCoefficients.Kub");
-        kurease_ = getParam<Scalar>("UreolysisCoefficients.Kurease");
-        ku_ = getParam<Scalar>("UreolysisCoefficients.Ku");
+        kUrease_ = getParam<Scalar>("UreolysisCoefficients.KUrease");
+        kUrea_ = getParam<Scalar>("UreolysisCoefficients.KUrea");
     }
 
     static constexpr int liquidPhaseIdx = FluidSystem::liquidPhaseIdx;
@@ -70,7 +70,7 @@ public:
      */
     static Scalar moleFracToMolality(Scalar moleFracX, Scalar moleFracSalinity, Scalar moleFracCTot)
     {
-        Scalar molalityX = moleFracX / (1 - moleFracSalinity - moleFracCTot) / FluidSystem::molarMass(H2OIdx);
+        const Scalar molalityX = moleFracX / (1 - moleFracSalinity - moleFracCTot) / FluidSystem::molarMass(H2OIdx);
         return molalityX;
     }
 
@@ -83,8 +83,8 @@ public:
     void reactionSource(NumEqVector &q, const VolumeVariables &volVars)
     {
         //   define and compute some parameters for convenience:
-        Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
-        Scalar densityBiofilm = volVars.solidComponentDensity(BiofilmPhaseIdx);
+        const Scalar xwCa = volVars.moleFraction(liquidPhaseIdx,CaIdx);
+        const Scalar densityBiofilm = volVars.solidComponentDensity(BiofilmPhaseIdx);
 
         Scalar volFracBiofilm = volVars.solidVolumeFraction(BiofilmPhaseIdx);
 
@@ -93,20 +93,19 @@ public:
 
         // TODO: dumux-course-task
         // implement mass of biofilm
-        Scalar massBiofilm = densityBiofilm * volFracBiofilm;
+        const Scalar massBiofilm = densityBiofilm * volFracBiofilm;
 
-        Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
+        const Scalar molalityUrea = moleFracToMolality(volVars.moleFraction(liquidPhaseIdx,UreaIdx),
                                                  xwCa,
                                                  volVars.moleFraction(liquidPhaseIdx,CO2Idx));  // [mol_urea/kg_H2O]
 
         // TODO: dumux-course-task
         // compute rate of ureolysis by implementing Z_urease,biofilm and r_urea
-        Scalar Zub = kub_ *  massBiofilm; // [kg urease/m³]
-        Scalar rurea = kurease_ * Zub * molalityUrea / (ku_ + molalityUrea); // [mol/m³s]
+        const Scalar Zub = kub_ *  massBiofilm; // [kg urease/m³]
+        const Scalar rurea = kUrease_ * Zub * molalityUrea / (kUrea_ + molalityUrea); // [mol/m³s]
 
         // compute/set dissolution and precipitation rate of calcite
-        Scalar rprec = 0.0;
-        rprec = rurea;
+        const Scalar rprec = rurea;
 
         // set source terms
         // TODO: dumux-course-task
@@ -122,8 +121,8 @@ public:
 private:
     // urease parameters
     Scalar kub_;
-    Scalar kurease_;
-    Scalar ku_;
+    Scalar kUrease_;
+    Scalar kUrea_;
 };
 
 } //end namespace Dumux

exercises/solution/exercise-biomineralization/fluidsystems/biomin.hh

@@ -34,7 +34,7 @@
 #include <dumux/material/components/tabulatedcomponent.hh>
 
 #include <dumux/material/components/calciumion.hh>
-#include "../components/urea.hh"
+#include <dumux/material/components/urea.hh>
 
 #include <dumux/material/binarycoefficients/brine_co2.hh>
 

exercises/solution/exercise-biomineralization/params.input

@@ -40,8 +40,8 @@ RhoBiofilm = 6.9 # [kg/m³] density of biofilm
 
 [UreolysisCoefficients]
 Kub = 8.81e-3 # [kg_urease/kg_bio] (max: 0.01)
-Kurease = 1000 # [mol_urea/(kg_urease s)]
-Ku = 0.355 # [mol/kgH2O] Lauchnor et al. 2014
+KUrease = 1000 # [mol_urea/(kg_urease s)]
+KUrea = 0.355 # [mol/kgH2O] Lauchnor et al. 2014
 
 [Brine]
 Salinity = 0.1