Calculators to compute ionic strength and stoichiometric ionic strength. More...

#include <GeochemistryIonicStrength.h>

Public Member Functions
	GeochemistryIonicStrength (Real max_ionic_strength, Real max_stoichiometric_ionic_strength, bool use_only_basis_molality, bool use_only_Cl_molality)

bool	operator== (const GeochemistryIonicStrength &rhs) const

Real	ionicStrength (const ModelGeochemicalDatabase &mgd, const std::vector< Real > &basis_species_molality, const std::vector< Real > &eqm_species_molality, const std::vector< Real > &kin_species_molality) const
	Compute ionic strength. More...

Real	stoichiometricIonicStrength (const ModelGeochemicalDatabase &mgd, const std::vector< Real > &basis_species_molality, const std::vector< Real > &eqm_species_molality, const std::vector< Real > &kin_species_molality) const
	Compute stoichiometric ionic strength. More...

void	setMaxIonicStrength (Real max_ionic_strength)
	Set the maximum ionic strength. More...

Real	getMaxIonicStrength () const
	Return the value of maximum ionic strength. More...

void	setMaxStoichiometricIonicStrength (Real max_stoichiometric_ionic_strength)
	Set the maximum stoichiometric ionic strength. More...

Real	getMaxStoichiometricIonicStrength () const
	Return the value of maximum stoichiometric ionic strength. More...

void	setUseOnlyBasisMolality (bool use_only_basis_molality)
	Set the value of use_only_basis_molality. More...

Real	getUseOnlyBasisMolality () const
	Return the value of use_only_basis_molality. More...

void	setUseOnlyClMolality (bool use_only_Cl_molality)
	Set the value of use_only_Cl_molality. More...

Real	getUseOnlyClMolality () const
	Return the value of use_only_Cl_molality. More...

Private Attributes
Real	_max_ionic_strength
	maximum ionic strength More...

Real	_max_stoichiometric_ionic_strength
	maximum stoichiometric ionic strength More...

bool	_use_only_basis_molality
	use only basis molality in the ionic strength calculations More...

bool	_use_only_Cl_molality
	set the stoichiometric ionic strength to the Cl- molality More...

Detailed Description

Calculators to compute ionic strength and stoichiometric ionic strength.

Definition at line 17 of file GeochemistryIonicStrength.h.

Constructor & Destructor Documentation

◆ GeochemistryIonicStrength()

GeochemistryIonicStrength::GeochemistryIonicStrength	(	Real	max_ionic_strength,
		Real	max_stoichiometric_ionic_strength,
		bool	use_only_basis_molality,
		bool	use_only_Cl_molality
	)

Parameters

max_ionic_strength	maximum value of ionic strength ever returned by this class
max_stoichiometric_ionic_strength	maximum value of stoichiometric ionic strength ever returned by this class
use_only_basis_molality	If true, use only the basis molalities in the ionic strength and stoichiometric ionic strength calculations. Since basis molality is usually much greater than equilibrium molality, and the whole concept of activity coefficients depending on ionic strength is only approximate in practice, setting use_only_basis_molality=true is often a reasonable approximation, and it can help with convergence of geochemistry algorithms since it eliminates problems associated with unphysical huge equilibrium molalities that can occur during Newton-iteration to the equilibrium solution
use_only_Cl_molality	If true, set the stoichiometric ionic strength to the Chlorine- molality. Note that this is the method used by the Geochemists Workbench

Definition at line 13 of file GeochemistryIonicStrength.C.

   : _max_ionic_strength(max_ionic_strength),
     _max_stoichiometric_ionic_strength(max_stoichiometric_ionic_strength),
     _use_only_basis_molality(use_only_basis_molality),
     _use_only_Cl_molality(use_only_Cl_molality)
 {
 }

Member Function Documentation

◆ getMaxIonicStrength()

Real GeochemistryIonicStrength::getMaxIonicStrength ( ) const

Return the value of maximum ionic strength.

Definition at line 129 of file GeochemistryIonicStrength.C.

 {
   return _max_ionic_strength;
 }

◆ getMaxStoichiometricIonicStrength()

Real GeochemistryIonicStrength::getMaxStoichiometricIonicStrength ( ) const

Return the value of maximum stoichiometric ionic strength.

Definition at line 141 of file GeochemistryIonicStrength.C.

 {
   return _max_stoichiometric_ionic_strength;
 }

◆ getUseOnlyBasisMolality()

Real GeochemistryIonicStrength::getUseOnlyBasisMolality ( ) const

Return the value of use_only_basis_molality.

Definition at line 153 of file GeochemistryIonicStrength.C.

Referenced by TEST().

 {
   return _use_only_basis_molality;
 }

◆ getUseOnlyClMolality()

Real GeochemistryIonicStrength::getUseOnlyClMolality ( ) const

Return the value of use_only_Cl_molality.

Definition at line 165 of file GeochemistryIonicStrength.C.

Referenced by TEST().

 {
   return _use_only_Cl_molality;
 }

◆ ionicStrength()

Real GeochemistryIonicStrength::ionicStrength	(	const ModelGeochemicalDatabase &	mgd,
		const std::vector< Real > &	basis_species_molality,
		const std::vector< Real > &	eqm_species_molality,
		const std::vector< Real > &	kin_species_molality
	)		const

Compute ionic strength.

Parameters

mgd	the Model Geochemical database
basis_species_molality	Molalities of the basis species in mgd
eqm_species_molality	Molalities of the equilibrium species in mgd
kin_species_molality	Molalities of the kinetic species

Returns: the ionic strength of the aqueous solution

Definition at line 25 of file GeochemistryIonicStrength.C.

Referenced by GeochemicalSystem::getIonicStrength(), GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters(), GeochemicalSystem::surfacePotPrefactor(), and TEST().

 {
   const unsigned num_basis = mgd.basis_species_charge.size();
   const unsigned num_eqm = mgd.eqm_species_charge.size();
   const unsigned num_kin = mgd.kin_species_charge.size();
   if (num_basis != basis_species_molality.size())
     mooseError("Ionic strength calculation: Number of basis species in mgd not equal to the size "
                "of basis_species_molality");
   if (num_eqm != eqm_species_molality.size())
     mooseError("Ionic strength calculation: Number of equilibrium species in mgd not equal to the "
                "size of eqm_species_molality");
   if (num_kin != kin_species_molality.size())
     mooseError("Ionic strength calculation: Number of kinetic species in mgd not equal to the size "
                "of kin_species_molality");
 
   Real ionic_strength = 0.0;
   for (unsigned i = 0; i < num_basis; ++i)
     ionic_strength += Utility::pow<2>(mgd.basis_species_charge[i]) * basis_species_molality[i];
   if (!_use_only_basis_molality)
   {
     for (unsigned i = 0; i < num_eqm; ++i)
       if (!mgd.surface_sorption_related[i])
         ionic_strength += Utility::pow<2>(mgd.eqm_species_charge[i]) * eqm_species_molality[i];
     for (unsigned i = 0; i < num_kin; ++i)
       ionic_strength += Utility::pow<2>(mgd.kin_species_charge[i]) * kin_species_molality[i] /
                         basis_species_molality[0]; // kin_species_molality is actually the number of
                                                    // moles of the kinetic species
   }
 
   return std::max(0.0, std::min(_max_ionic_strength, 0.5 * ionic_strength));
 }

◆ operator==()

bool GeochemistryIonicStrength::operator== ( const GeochemistryIonicStrength & rhs ) const

inline

Definition at line 39 of file GeochemistryIonicStrength.h.

   {
     return (_max_ionic_strength == rhs._max_ionic_strength) &&
            (_max_stoichiometric_ionic_strength == rhs._max_stoichiometric_ionic_strength) &&
            (_use_only_basis_molality == rhs._use_only_basis_molality) &&
            (_use_only_Cl_molality == rhs._use_only_Cl_molality);
   };

◆ setMaxIonicStrength()

void GeochemistryIonicStrength::setMaxIonicStrength ( Real max_ionic_strength )

Set the maximum ionic strength.

Definition at line 123 of file GeochemistryIonicStrength.C.

Referenced by GeochemicalSolver::solveSystem().

 {
   _max_ionic_strength = max_ionic_strength;
 }

◆ setMaxStoichiometricIonicStrength()

void GeochemistryIonicStrength::setMaxStoichiometricIonicStrength ( Real max_stoichiometric_ionic_strength )

Set the maximum stoichiometric ionic strength.

Definition at line 135 of file GeochemistryIonicStrength.C.

Referenced by GeochemicalSolver::solveSystem().

 {
   _max_stoichiometric_ionic_strength = max_stoichiometric_ionic_strength;
 }

◆ setUseOnlyBasisMolality()

void GeochemistryIonicStrength::setUseOnlyBasisMolality ( bool use_only_basis_molality )

Set the value of use_only_basis_molality.

Definition at line 147 of file GeochemistryIonicStrength.C.

 {
   _use_only_basis_molality = use_only_basis_molality;
 }

◆ setUseOnlyClMolality()

void GeochemistryIonicStrength::setUseOnlyClMolality ( bool use_only_Cl_molality )

Set the value of use_only_Cl_molality.

Definition at line 159 of file GeochemistryIonicStrength.C.

 {
   _use_only_Cl_molality = use_only_Cl_molality;
 }

◆ stoichiometricIonicStrength()

Real GeochemistryIonicStrength::stoichiometricIonicStrength	(	const ModelGeochemicalDatabase &	mgd,
		const std::vector< Real > &	basis_species_molality,
		const std::vector< Real > &	eqm_species_molality,
		const std::vector< Real > &	kin_species_molality
	)		const

Compute stoichiometric ionic strength.

Parameters

mgd	the Model Geochemical database
basis_species_molality	Molalities of the basis species in mgd
eqm_species_molality	Molalities of the equilibrium species in mgd
kin_species_molality	Molalities of the kinetic species

Returns: the stoichiometric ionic strength of the aqueous solution

Definition at line 61 of file GeochemistryIonicStrength.C.

Referenced by GeochemicalSystem::getStoichiometricIonicStrength(), GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters(), and TEST().

 {
   const unsigned num_basis = mgd.basis_species_charge.size();
   const unsigned num_eqm = mgd.eqm_species_charge.size();
   const unsigned num_kin = mgd.kin_species_charge.size();
   if (num_basis != basis_species_molality.size())
     mooseError("Stoichiometric ionic strength calculation: Number of basis species in mgd not "
                "equal to the size of basis_species_molality");
   if (num_eqm != eqm_species_molality.size())
     mooseError("Stoichiometric ionic strength calculation: Number of equilibrium species in mgd "
                "not equal to the size of eqm_species_molality");
   if (num_kin != kin_species_molality.size())
     mooseError("Stoichiometric ionic strength calculation: Number of kinetic species in mgd not "
                "equal to the size of kin_species_molality");
 
   Real ionic_strength = 0.0;
   if (_use_only_Cl_molality)
   {
     if (mgd.basis_species_index.count("Cl-"))
       ionic_strength = basis_species_molality[mgd.basis_species_index.at("Cl-")];
     else if (mgd.eqm_species_index.count("Cl-"))
       ionic_strength = eqm_species_molality[mgd.eqm_species_index.at("Cl-")];
     else if (mgd.kin_species_index.count("Cl-"))
       ionic_strength = kin_species_molality[mgd.kin_species_index.at("Cl-")];
     else
       mooseError("GeochemistryIonicStrength: attempting to compute stoichiometric ionic strength "
                  "using only the Cl- molality, but Cl- does not appear in the geochemical system");
     return std::max(0.0, std::min(_max_stoichiometric_ionic_strength, ionic_strength));
   }
 
   for (unsigned i = 0; i < num_basis; ++i)
     ionic_strength += Utility::pow<2>(mgd.basis_species_charge[i]) * basis_species_molality[i];
   if (!_use_only_basis_molality)
   {
     for (unsigned i = 0; i < num_eqm; ++i)
       if (!mgd.surface_sorption_related[i])
       {
 
         if (mgd.eqm_species_charge[i] != 0.0)
           ionic_strength += Utility::pow<2>(mgd.eqm_species_charge[i]) * eqm_species_molality[i];
         else
         {
           for (unsigned j = 0; j < num_basis; ++j)
             ionic_strength += Utility::pow<2>(mgd.basis_species_charge[j]) *
                               eqm_species_molality[i] * mgd.eqm_stoichiometry(i, j);
         }
       }
     for (unsigned i = 0; i < num_kin;
          ++i) // kin_species_molality is actually the number of moles, not molality
       if (mgd.kin_species_charge[i] != 0.0)
         ionic_strength += Utility::pow<2>(mgd.kin_species_charge[i]) * kin_species_molality[i] /
                           basis_species_molality[0];
   }
 
   return std::max(0.0, std::min(_max_stoichiometric_ionic_strength, 0.5 * ionic_strength));
 }