Computes activity coefficients for non-minerals and non-gases (since these species do not have activity coefficients). More...

#include <GeochemistryActivityCoefficientsDebyeHuckel.h>

Inheritance diagram for GeochemistryActivityCoefficientsDebyeHuckel:

[legend]

Public Member Functions
	GeochemistryActivityCoefficientsDebyeHuckel (const GeochemistryIonicStrength &is_calculator, const GeochemicalDatabaseReader &db)

void	setInternalParameters (Real temperature, const ModelGeochemicalDatabase &mgd, const std::vector< Real > &basis_species_molality, const std::vector< Real > &eqm_species_molality, const std::vector< Real > &kin_species_molality) override
	Sets internal parameters, such as the ionic strength and Debye-Huckel parameters, prior to computing activity coefficients and activity of water. More...

Real	waterActivity () const override
	Computes and returns the activity of water. More...

void	buildActivityCoefficients (const ModelGeochemicalDatabase &mgd, std::vector< Real > &basis_activity_coef, std::vector< Real > &eqm_activity_coef) const override
	Compute the activity coefficients and store them in basis_activity_coef and eqm_activity_coef Note: More...

const DebyeHuckelParameters &	getDebyeHuckel () const

Real	getIonicStrength () const
	Return the current value of ionic strength. More...

Real	getStoichiometricIonicStrength () const
	Return the current value of stoichiometric ionic strength. More...

bool	operator== (const GeochemistryActivityCoefficients &) const

Private Attributes
const unsigned	_numT
	number of temperature points in the database file More...

const GeochemistryDebyeHuckel	_database_dh_params
	Debye-Huckel parameters found in the database. More...

const GeochemistryNeutralSpeciesActivity	_database_dh_water
	Debye-Huckel parameters found in the database for computing the water activities. More...

const GeochemistryNeutralSpeciesActivity	_database_dh_neutral
	Debye-Huckel parameters found in the database for computing the neutral (CO2) activities. More...

const GeochemistryIonicStrength &	_is_calculator
	ionic-strength calculator More...

Real	_ionic_strength
	current value of ionic strength More...

Real	_sqrt_ionic_strength
	current value of sqrt(ionic strength) More...

Real	_stoichiometric_ionic_strength
	current value of stoichiometric ionic strength More...

unsigned	_num_basis
	number of basis species More...

unsigned	_num_eqm
	number of equilibrium species More...

DebyeHuckelParameters	_dh
	Debye-Huckel parameters. More...

EquilibriumConstantInterpolator	_interp_A
	Interpolator object for the Debye-Huckel parameter A. More...

EquilibriumConstantInterpolator	_interp_B
	Interpolator object for the Debye-Huckel parameter B. More...

EquilibriumConstantInterpolator	_interp_Bdot
	Interpolator object for the Debye-Huckel parameter Bdot. More...

EquilibriumConstantInterpolator	_interp_a_water
	Interpolator object for the Debye-Huckel parameter a_water. More...

EquilibriumConstantInterpolator	_interp_b_water
	Interpolator object for the Debye-Huckel parameter b_water. More...

EquilibriumConstantInterpolator	_interp_c_water
	Interpolator object for the Debye-Huckel parameter c_water. More...

EquilibriumConstantInterpolator	_interp_d_water
	Interpolator object for the Debye-Huckel parameter d_water. More...

EquilibriumConstantInterpolator	_interp_a_neutral
	Interpolator object for the Debye-Huckel parameter a_neutral. More...

EquilibriumConstantInterpolator	_interp_b_neutral
	Interpolator object for the Debye-Huckel parameter b_neutral. More...

EquilibriumConstantInterpolator	_interp_c_neutral
	Interpolator object for the Debye-Huckel parameter c_neutral. More...

EquilibriumConstantInterpolator	_interp_d_neutral
	Interpolator object for the Debye-Huckel parameter d_neutral. More...

Detailed Description

Computes activity coefficients for non-minerals and non-gases (since these species do not have activity coefficients).

Also computes the activity of water. Uses a Debye-Huckel model

Definition at line 58 of file GeochemistryActivityCoefficientsDebyeHuckel.h.

Constructor & Destructor Documentation

◆ GeochemistryActivityCoefficientsDebyeHuckel()

GeochemistryActivityCoefficientsDebyeHuckel::GeochemistryActivityCoefficientsDebyeHuckel	(	const GeochemistryIonicStrength &	is_calculator,
		const GeochemicalDatabaseReader &	db
	)

Parameters

method	Method used by this class to compute activity coefficients and activity of water
is_calculator	Calculates ionic strengths
db	Original geochemistry database: used to find the temperature interpolation type, Debye-Huckel params, etc

Definition at line 13 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

   : GeochemistryActivityCoefficients(),
     _numT(db.getTemperatures().size()),
     _database_dh_params(db.getDebyeHuckel()),
     _database_dh_water((db.getNeutralSpeciesActivity().count("h2o") == 1)
                            ? db.getNeutralSpeciesActivity().at("h2o")
                            : GeochemistryNeutralSpeciesActivity()),
     _database_dh_neutral((db.getNeutralSpeciesActivity().count("co2") == 1)
                              ? db.getNeutralSpeciesActivity().at("co2")
                              : GeochemistryNeutralSpeciesActivity()),
     _is_calculator(is_calculator),
     _ionic_strength(1.0),
     _sqrt_ionic_strength(1.0),
     _stoichiometric_ionic_strength(1.0),
     _num_basis(0),
     _num_eqm(0),
     _dh(),
     _interp_A(db.getTemperatures(),
               (_database_dh_params.adh.size() == _numT) ? _database_dh_params.adh
                                                         : std::vector<Real>(_numT, 0.0),
               db.getLogKModel()),
     _interp_B(db.getTemperatures(),
               (_database_dh_params.bdh.size() == _numT) ? _database_dh_params.bdh
                                                         : std::vector<Real>(_numT, 0.0),
               db.getLogKModel()),
     _interp_Bdot(db.getTemperatures(),
                  (_database_dh_params.bdot.size() == _numT) ? _database_dh_params.bdot
                                                             : std::vector<Real>(_numT, 0.0),
                  db.getLogKModel()),
     _interp_a_water(db.getTemperatures(),
                     (_database_dh_water.a.size() == _numT) ? _database_dh_water.a
                                                            : std::vector<Real>(_numT, 0.0),
                     db.getLogKModel()),
     _interp_b_water(db.getTemperatures(),
                     (_database_dh_water.b.size() == _numT) ? _database_dh_water.b
                                                            : std::vector<Real>(_numT, 0.0),
                     db.getLogKModel()),
     _interp_c_water(db.getTemperatures(),
                     (_database_dh_water.c.size() == _numT) ? _database_dh_water.c
                                                            : std::vector<Real>(_numT, 0.0),
                     db.getLogKModel()),
     _interp_d_water(db.getTemperatures(),
                     (_database_dh_water.d.size() == _numT) ? _database_dh_water.d
                                                            : std::vector<Real>(_numT, 0.0),
                     db.getLogKModel()),
     _interp_a_neutral(db.getTemperatures(),
                       (_database_dh_neutral.a.size() == _numT) ? _database_dh_neutral.a
                                                                : std::vector<Real>(_numT, 0.0),
                       db.getLogKModel()),
     _interp_b_neutral(db.getTemperatures(),
                       (_database_dh_neutral.b.size() == _numT) ? _database_dh_neutral.b
                                                                : std::vector<Real>(_numT, 0.0),
                       db.getLogKModel()),
     _interp_c_neutral(db.getTemperatures(),
                       (_database_dh_neutral.c.size() == _numT) ? _database_dh_neutral.c
                                                                : std::vector<Real>(_numT, 0.0),
                       db.getLogKModel()),
     _interp_d_neutral(db.getTemperatures(),
                       (_database_dh_neutral.d.size() == _numT) ? _database_dh_neutral.d
                                                                : std::vector<Real>(_numT, 0.0),
                       db.getLogKModel())
 {
   _interp_A.generate();
   _interp_B.generate();
   _interp_Bdot.generate();
   _interp_a_water.generate();
   _interp_b_water.generate();
   _interp_c_water.generate();
   _interp_d_water.generate();
   _interp_a_neutral.generate();
   _interp_b_neutral.generate();
   _interp_c_neutral.generate();
   _interp_d_neutral.generate();
 }

Member Function Documentation

◆ buildActivityCoefficients()

void GeochemistryActivityCoefficientsDebyeHuckel::buildActivityCoefficients	(	const ModelGeochemicalDatabase &	mgd,
		std::vector< Real > &	basis_activity_coef,
		std::vector< Real > &	eqm_activity_coef
	)		const

overridevirtual

Compute the activity coefficients and store them in basis_activity_coef and eqm_activity_coef Note:

you will probably want to call setInternalParameters prior to calling this method
the activity coefficient for water (basis species = 0) is not computed since it is meaningless: use getWaterActivity() instead
the activity coefficient for any mineral is not computed since minerals do not have activity coefficients
the activity coefficient for any gas is not computed since gases do not have activity coefficients Hence, the elements in basis_activity_coef and eqm_activity_coef corresponding to these species will be undefined after this method returns

Implements GeochemistryActivityCoefficients.

Definition at line 134 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   basis_activity_coef.resize(_num_basis);
   eqm_activity_coef.resize(_num_eqm);
 
   for (unsigned basis_i = 1; basis_i < _num_basis; ++basis_i) // don't loop over water
     if (mgd.basis_species_mineral[basis_i] || mgd.basis_species_gas[basis_i])
       continue; // these should never be used
     else if (mgd.basis_species_radius[basis_i] == -0.5)
     {
       basis_activity_coef[basis_i] = std::pow(
           10.0,
           GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
               _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
     }
     else if (mgd.basis_species_radius[basis_i] == -1.0)
     {
       basis_activity_coef[basis_i] =
           std::pow(10.0,
                    GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(_ionic_strength,
                                                                                    _dh.Bdot));
     }
     else if (mgd.basis_species_radius[basis_i] == -1.5)
     {
       basis_activity_coef[basis_i] = 1.0;
     }
     else if (mgd.basis_species_charge[basis_i] == 0.0)
     {
       basis_activity_coef[basis_i] = 1.0;
     }
     else
     {
       basis_activity_coef[basis_i] = std::pow(
           10.0,
           GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.basis_species_charge[basis_i],
                                                                mgd.basis_species_radius[basis_i],
                                                                _sqrt_ionic_strength,
                                                                _dh.A,
                                                                _dh.B,
                                                                _dh.Bdot));
     }
 
   for (unsigned eqm_j = 0; eqm_j < _num_eqm; ++eqm_j)
     if (mgd.eqm_species_mineral[eqm_j] || mgd.eqm_species_gas[eqm_j])
       continue;
     else if (mgd.eqm_species_radius[eqm_j] == -0.5)
     {
       eqm_activity_coef[eqm_j] = std::pow(
           10.0,
           GeochemistryActivityCalculators::log10ActCoeffDHBdotNeutral(
               _ionic_strength, _dh.a_neutral, _dh.b_neutral, _dh.c_neutral, _dh.d_neutral));
     }
     else if (mgd.eqm_species_radius[eqm_j] == -1.0)
     {
       eqm_activity_coef[eqm_j] =
           std::pow(10.0,
                    GeochemistryActivityCalculators::log10ActCoeffDHBdotAlternative(_ionic_strength,
                                                                                    _dh.Bdot));
     }
     else if (mgd.eqm_species_radius[eqm_j] == -1.5)
     {
       eqm_activity_coef[eqm_j] = 1.0;
     }
     else if (mgd.eqm_species_charge[eqm_j] == 0.0)
     {
       eqm_activity_coef[eqm_j] = 1.0;
     }
     else
     {
       eqm_activity_coef[eqm_j] = std::pow(
           10.0,
           GeochemistryActivityCalculators::log10ActCoeffDHBdot(mgd.eqm_species_charge[eqm_j],
                                                                mgd.eqm_species_radius[eqm_j],
                                                                _sqrt_ionic_strength,
                                                                _dh.A,
                                                                _dh.B,
                                                                _dh.Bdot));
     }
 }

◆ getDebyeHuckel()

const DebyeHuckelParameters & GeochemistryActivityCoefficientsDebyeHuckel::getDebyeHuckel ( ) const

Returns: the Debye-Huckel parameters

Definition at line 121 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   return _dh;
 }

◆ getIonicStrength()

Real GeochemistryActivityCoefficientsDebyeHuckel::getIonicStrength ( ) const

Return the current value of ionic strength.

Definition at line 218 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   return _ionic_strength;
 }

◆ getStoichiometricIonicStrength()

Real GeochemistryActivityCoefficientsDebyeHuckel::getStoichiometricIonicStrength ( ) const

Return the current value of stoichiometric ionic strength.

Definition at line 224 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   return _stoichiometric_ionic_strength;
 }

◆ operator==()

bool GeochemistryActivityCoefficients::operator== ( const GeochemistryActivityCoefficients & ) const

inlineinherited

Definition at line 23 of file GeochemistryActivityCoefficients.h.

23 { return true; };

◆ setInternalParameters()

void GeochemistryActivityCoefficientsDebyeHuckel::setInternalParameters	(	Real	temperature,
		const ModelGeochemicalDatabase &	mgd,
		const std::vector< Real > &	basis_species_molality,
		const std::vector< Real > &	eqm_species_molality,
		const std::vector< Real > &	kin_species_molality
	)

overridevirtual

Sets internal parameters, such as the ionic strength and Debye-Huckel parameters, prior to computing activity coefficients and activity of water.

If using a Debye-Huckel activity model, you must ensure the ionic strength calculator (eg, maxIonicStrength) is set appropriately before calling this function.

Parameters

temperature	the temperature in degC
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

Implements GeochemistryActivityCoefficients.

Definition at line 90 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   _ionic_strength = _is_calculator.ionicStrength(
       mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
   _sqrt_ionic_strength = std::sqrt(_ionic_strength);
   _stoichiometric_ionic_strength = _is_calculator.stoichiometricIonicStrength(
       mgd, basis_species_molality, eqm_species_molality, kin_species_molality);
   _num_basis = mgd.basis_species_index.size();
   _num_eqm = mgd.eqm_species_index.size();
   // Debye-Huckel base parameters
   _dh.A = _interp_A.sample(temperature);
   _dh.B = _interp_B.sample(temperature);
   _dh.Bdot = _interp_Bdot.sample(temperature);
   // water Debye-Huckel
   _dh.a_water = _interp_a_water.sample(temperature);
   _dh.b_water = _interp_b_water.sample(temperature);
   _dh.c_water = _interp_c_water.sample(temperature);
   _dh.d_water = _interp_d_water.sample(temperature);
   // neutral species Debye-Huckel
   _dh.a_neutral = _interp_a_neutral.sample(temperature);
   _dh.b_neutral = _interp_b_neutral.sample(temperature);
   _dh.c_neutral = _interp_c_neutral.sample(temperature);
   _dh.d_neutral = _interp_d_neutral.sample(temperature);
 }

◆ waterActivity()

Real GeochemistryActivityCoefficientsDebyeHuckel::waterActivity ( ) const

overridevirtual

Computes and returns the activity of water.

Note that you will probably want to call setInternalParameters prior to calling this method

Returns: the activity of water

Implements GeochemistryActivityCoefficients.

Definition at line 127 of file GeochemistryActivityCoefficientsDebyeHuckel.C.

Referenced by TEST().

 {
   return std::exp(GeochemistryActivityCalculators::lnActivityDHBdotWater(
       _stoichiometric_ionic_strength, _dh.A, _dh.a_water, _dh.b_water, _dh.c_water, _dh.d_water));
 }