Module pychnosz.geochemistry
Geochemistry package for CHNOSZ.
This package provides specialized geochemical calculations including mineral equilibria, redox reactions, and environmental applications.
Sub-modules
pychnosz.geochemistry.minerals-
Mineral equilibria and phase diagram calculations for CHNOSZ …
pychnosz.geochemistry.redox-
Redox reactions and Eh-pH diagram calculations for CHNOSZ …
Functions
def eh(pe_value: float, T: float = 298.15) ‑> float-
Expand source code
def eh(pe_value: float, T: float = 298.15) -> float: """ Convert pe to Eh. Parameters ---------- pe_value : float Electron activity T : float, default 298.15 Temperature in Kelvin Returns ------- float Eh in Volts """ return _redox_calculator.eh_from_pe(pe_value, T)Convert pe to Eh.
Parameters
pe_value:float- Electron activity
T:float, default298.15- Temperature in Kelvin
Returns
float- Eh in Volts
def eh_ph(element: str, **kwargs) ‑> Dict[str, Any]-
Expand source code
def eh_ph(element: str, **kwargs) -> Dict[str, Any]: """ Create Eh-pH diagram for an element. Parameters ---------- element : str Element symbol **kwargs Additional parameters for eh_ph_diagram() Returns ------- dict Eh-pH diagram data """ return _redox_calculator.eh_ph_diagram(element, **kwargs)Create Eh-pH diagram for an element.
Parameters
element:str- Element symbol
**kwargs- Additional parameters for eh_ph_diagram()
Returns
dict- Eh-pH diagram data
def logfO2(pe_value: float, pH: float = 7.0, **kwargs) ‑> float-
Expand source code
def logfO2(pe_value: float, pH: float = 7.0, **kwargs) -> float: """ Calculate log oxygen fugacity. Parameters ---------- pe_value : float Electron activity pH : float, default 7.0 Solution pH **kwargs Additional parameters Returns ------- float log fO2 """ return _redox_calculator.oxygen_fugacity(pe_value, pH, **kwargs)Calculate log oxygen fugacity.
Parameters
pe_value:float- Electron activity
pH:float, default7.0- Solution pH
**kwargs- Additional parameters
Returns
float- log fO2
def mineral_solubility(mineral: str,
pH_range: Tuple[float, float] = (0, 14),
T: float = 298.15,
P: float = 1.0) ‑> Dict[str, Any]-
Expand source code
def mineral_solubility(mineral: str, pH_range: Tuple[float, float] = (0, 14), T: float = 298.15, P: float = 1.0) -> Dict[str, Any]: """ Calculate mineral solubility vs pH. Parameters ---------- mineral : str Mineral name pH_range : tuple, default (0, 14) pH range for calculation T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- dict Solubility data """ return _mineral_calculator.mineral_solubility(mineral, T, P, pH_range)Calculate mineral solubility vs pH.
Parameters
mineral:str- Mineral name
pH_range:tuple, default(0, 14)- pH range for calculation
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
Returns
dict- Solubility data
def pe(redox_couple: str | Dict[str, float],
concentrations: Dict[str, float],
**kwargs) ‑> float-
Expand source code
def pe(redox_couple: Union[str, Dict[str, float]], concentrations: Dict[str, float], **kwargs) -> float: """ Calculate pe for a redox couple. Parameters ---------- redox_couple : str or dict Redox couple or reaction concentrations : dict Species concentrations **kwargs Additional parameters Returns ------- float pe value """ return _redox_calculator.pe_calculation(redox_couple, concentrations, **kwargs)Calculate pe for a redox couple.
Parameters
redox_couple:strordict- Redox couple or reaction
concentrations:dict- Species concentrations
**kwargs- Additional parameters
Returns
float- pe value
def phase_boundary(reaction: Dict[str, float], T_range: Tuple[float, float] = (273.15, 673.15)) ‑> Dict[str, Any]-
Expand source code
def phase_boundary(reaction: Dict[str, float], T_range: Tuple[float, float] = (273.15, 673.15)) -> Dict[str, Any]: """ Calculate phase boundary for a reaction. Parameters ---------- reaction : dict Reaction stoichiometry {species: coefficient} T_range : tuple, default (273.15, 673.15) Temperature range in Kelvin Returns ------- dict Phase boundary data """ return _mineral_calculator.phase_equilibrium(reaction, T_range=T_range)Calculate phase boundary for a reaction.
Parameters
reaction:dict- Reaction stoichiometry {species: coefficient}
T_range:tuple, default(273.15, 673.15)- Temperature range in Kelvin
Returns
dict- Phase boundary data
def stability_field(minerals: List[str],
variables: Tuple[str, str],
ranges: Tuple[Tuple[float, float], Tuple[float, float]],
T: float = 298.15,
P: float = 1.0) ‑> Dict[str, Any]-
Expand source code
def stability_field(minerals: List[str], variables: Tuple[str, str], ranges: Tuple[Tuple[float, float], Tuple[float, float]], T: float = 298.15, P: float = 1.0) -> Dict[str, Any]: """ Calculate mineral stability fields. Parameters ---------- minerals : list of str Minerals to consider variables : tuple of str (x_variable, y_variable) for diagram ranges : tuple of tuples ((x_min, x_max), (y_min, y_max)) T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- dict Stability diagram data """ x_var, y_var = variables x_range, y_range = ranges return _mineral_calculator.stability_diagram( minerals, x_var, x_range, y_var, y_range, T, P)Calculate mineral stability fields.
Parameters
minerals:listofstr- Minerals to consider
variables:tupleofstr- (x_variable, y_variable) for diagram
ranges:tupleoftuples- ((x_min, x_max), (y_min, y_max))
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
Returns
dict- Stability diagram data
Classes
class MineralEquilibria-
Expand source code
class MineralEquilibria: """ Mineral equilibria calculator for geological systems. This class handles mineral stability calculations, including: - Mineral solubility equilibria - Phase stability diagrams - Mineral-water reactions - Solid solution equilibria """ def __init__(self): """Initialize the mineral equilibria calculator.""" self.equilibrium_solver = EquilibriumSolver() self.speciation_engine = SpeciationEngine() # Common mineral dissolution reactions (simplified) self.dissolution_reactions = { 'quartz': {'SiO2': 1, 'H2O': 0, 'H4SiO4': -1}, 'calcite': {'CaCO3': 1, 'H+': 2, 'Ca+2': -1, 'HCO3-': -1, 'H2O': -1}, 'pyrite': {'FeS2': 1, 'H+': 8, 'SO4-2': -2, 'Fe+2': -1, 'H2O': -4}, 'hematite': {'Fe2O3': 1, 'H+': 6, 'Fe+3': -2, 'H2O': -3}, 'magnetite': {'Fe3O4': 1, 'H+': 8, 'Fe+2': -1, 'Fe+3': -2, 'H2O': -4}, 'albite': {'NaAlSi3O8': 1, 'H+': 1, 'H2O': 4, 'Na+': -1, 'Al+3': -1, 'H4SiO4': -3}, 'anorthite': {'CaAl2Si2O8': 1, 'H+': 8, 'H2O': 0, 'Ca+2': -1, 'Al+3': -2, 'H4SiO4': -2}, 'kaolinite': {'Al2Si2O5(OH)4': 1, 'H+': 6, 'Al+3': -2, 'H4SiO4': -2, 'H2O': -1} } def mineral_solubility(self, mineral: str, T: float = 298.15, P: float = 1.0, pH_range: Optional[Tuple[float, float]] = None, ionic_strength: float = 0.0) -> Dict[str, Any]: """ Calculate mineral solubility as a function of pH. Parameters ---------- mineral : str Mineral name T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar pH_range : tuple, optional pH range (min, max). Default: (0, 14) ionic_strength : float, default 0.0 Solution ionic strength Returns ------- dict Solubility data vs pH """ if pH_range is None: pH_range = (0, 14) pH_values = np.linspace(pH_range[0], pH_range[1], 100) if mineral.lower() not in self.dissolution_reactions: warnings.warn(f"Dissolution reaction not defined for {mineral}") return {'pH': pH_values, 'solubility': np.zeros_like(pH_values)} reaction = self.dissolution_reactions[mineral.lower()] try: # Get equilibrium constant for dissolution species_names = list(reaction.keys()) coefficients = list(reaction.values()) # Calculate logK at T, P result = subcrt(species_names, coefficients, T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 0.0 warnings.warn(f"Could not calculate logK for {mineral} dissolution") except Exception as e: warnings.warn(f"Error calculating equilibrium constant: {e}") logK = 0.0 # Calculate solubility at each pH solubility = np.zeros_like(pH_values) for i, pH in enumerate(pH_values): try: # Simplified solubility calculation # This would be more complex in a full implementation if 'H+' in reaction: h_coeff = reaction['H+'] # Basic pH dependence log_solubility = logK - h_coeff * pH else: log_solubility = logK # Convert to molality solubility[i] = 10**log_solubility # Apply ionic strength corrections if needed if ionic_strength > 0: # Simple activity coefficient correction gamma = self._estimate_activity_coefficient(ionic_strength) solubility[i] *= gamma except: solubility[i] = 1e-20 # Very low solubility return { 'pH': pH_values, 'solubility': solubility, 'logK': logK, 'mineral': mineral, 'T': T, 'P': P } def stability_diagram(self, minerals: List[str], x_axis: str, x_range: Tuple[float, float], y_axis: str, y_range: Tuple[float, float], T: float = 298.15, P: float = 1.0, resolution: int = 50) -> Dict[str, Any]: """ Create mineral stability diagram. Parameters ---------- minerals : list of str Mineral names to consider x_axis : str X-axis variable ('pH', 'log_activity_SiO2', etc.) x_range : tuple X-axis range (min, max) y_axis : str Y-axis variable y_range : tuple Y-axis range (min, max) T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar resolution : int, default 50 Grid resolution Returns ------- dict Stability diagram data """ x_values = np.linspace(x_range[0], x_range[1], resolution) y_values = np.linspace(y_range[0], y_range[1], resolution) X, Y = np.meshgrid(x_values, y_values) # Initialize stability field array stable_mineral = np.zeros_like(X, dtype=int) # Calculate stability at each grid point for i in range(len(y_values)): for j in range(len(x_values)): x_val, y_val = X[i, j], Y[i, j] # Find most stable mineral at this point min_energy = np.inf stable_idx = 0 for k, mineral in enumerate(minerals): try: # Calculate relative stability # This would involve full equilibrium calculations # For now, use simplified energy estimates energy = self._calculate_stability_energy( mineral, x_axis, x_val, y_axis, y_val, T, P) if energy < min_energy: min_energy = energy stable_idx = k except Exception as e: warnings.warn(f"Error calculating stability for {mineral}: {e}") continue stable_mineral[i, j] = stable_idx return { x_axis: x_values, y_axis: y_values, 'stable_mineral': stable_mineral, 'mineral_names': minerals, 'T': T, 'P': P } def phase_equilibrium(self, reaction: Dict[str, float], T_range: Optional[Tuple[float, float]] = None, P_range: Optional[Tuple[float, float]] = None) -> Dict[str, Any]: """ Calculate phase equilibrium curve. Parameters ---------- reaction : dict Reaction stoichiometry {species: coefficient} T_range : tuple, optional Temperature range (min, max) in K P_range : tuple, optional Pressure range (min, max) in bar Returns ------- dict Equilibrium curve data """ if T_range is None and P_range is None: T_range = (273.15, 673.15) # 0-400°C if T_range is not None: # T-P curve at equilibrium T_values = np.linspace(T_range[0], T_range[1], 50) P_equilibrium = np.zeros_like(T_values) for i, T in enumerate(T_values): try: # Calculate equilibrium pressure # This would involve iterative solution # For now, use simplified correlation P_equilibrium[i] = self._calculate_equilibrium_pressure(reaction, T) except: P_equilibrium[i] = 1.0 # Default return { 'T': T_values, 'P': P_equilibrium, 'reaction': reaction } else: # P-T curve - similar logic P_values = np.linspace(P_range[0], P_range[1], 50) T_equilibrium = np.zeros_like(P_values) for i, P in enumerate(P_values): try: T_equilibrium[i] = self._calculate_equilibrium_temperature(reaction, P) except: T_equilibrium[i] = 298.15 return { 'P': P_values, 'T': T_equilibrium, 'reaction': reaction } def solid_solution(self, end_members: List[str], compositions: List[float], T: float = 298.15, P: float = 1.0, model: str = 'ideal') -> Dict[str, Any]: """ Calculate solid solution thermodynamics. Parameters ---------- end_members : list of str End-member mineral names compositions : list of float Mole fractions of end-members (must sum to 1) T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar model : str, default 'ideal' Mixing model ('ideal', 'regular', 'margules') Returns ------- dict Solid solution properties """ if len(end_members) != len(compositions): raise ValueError("Number of end-members must match number of compositions") if abs(sum(compositions) - 1.0) > 1e-6: raise ValueError("Compositions must sum to 1.0") # Calculate end-member properties end_member_props = [] for mineral in end_members: try: # This would get mineral properties from database props = self._get_mineral_properties(mineral, T, P) end_member_props.append(props) except: # Default properties props = {'G': -1000000, 'H': -1200000, 'S': 100, 'V': 50, 'Cp': 100} end_member_props.append(props) # Calculate mixing properties if model == 'ideal': # Ideal mixing G_mix = 8.314 * T * sum(x * np.log(x) for x in compositions if x > 0) H_mix = 0.0 S_mix = -8.314 * sum(x * np.log(x) for x in compositions if x > 0) elif model == 'regular': # Regular solution (symmetric) W = 10000 # Interaction parameter (J/mol) - would be mineral-specific G_mix = W * compositions[0] * compositions[1] # Binary case H_mix = G_mix S_mix = 0.0 else: # Default to ideal G_mix = 8.314 * T * sum(x * np.log(x) for x in compositions if x > 0) H_mix = 0.0 S_mix = -8.314 * sum(x * np.log(x) for x in compositions if x > 0) # Total properties total_props = {} for prop in ['G', 'H', 'S', 'V', 'Cp']: total_props[prop] = sum(compositions[i] * end_member_props[i][prop] for i in range(len(end_members))) # Add mixing contributions total_props['G'] += G_mix total_props['H'] += H_mix total_props['S'] += S_mix return { 'properties': total_props, 'mixing': {'G': G_mix, 'H': H_mix, 'S': S_mix}, 'end_members': end_members, 'compositions': compositions, 'model': model } def _estimate_activity_coefficient(self, ionic_strength: float) -> float: """Estimate activity coefficient from ionic strength.""" if ionic_strength <= 0: return 1.0 else: # Simple Debye-Hückel approximation A = 0.509 return 10**(-A * np.sqrt(ionic_strength)) def _calculate_stability_energy(self, mineral: str, x_axis: str, x_val: float, y_axis: str, y_val: float, T: float, P: float) -> float: """Calculate relative stability energy for mineral.""" # Simplified stability calculation # Would involve full thermodynamic analysis base_energy = np.random.random() * 1000 # Placeholder # Add dependencies on axes if 'pH' in x_axis: base_energy += abs(x_val - 7) * 100 # pH dependence if 'log' in y_axis: base_energy += abs(y_val + 3) * 50 # Activity dependence return base_energy def _calculate_equilibrium_pressure(self, reaction: Dict[str, float], T: float) -> float: """Calculate equilibrium pressure for reaction at given T.""" # Simplified using Clapeyron equation approximation # dP/dT = ΔS/ΔV # Estimate volume and entropy changes delta_V = 5.0e-6 # m³/mol (typical for mineral reactions) delta_S = 20.0 # J/(mol·K) (typical) # Integrate from reference conditions T0, P0 = 298.15, 1.0 dP_dT = delta_S / delta_V * 1e-5 # Convert to bar/K P_eq = P0 + dP_dT * (T - T0) return max(P_eq, 0.1) # Minimum 0.1 bar def _calculate_equilibrium_temperature(self, reaction: Dict[str, float], P: float) -> float: """Calculate equilibrium temperature for reaction at given P.""" # Inverse of pressure calculation delta_V = 5.0e-6 delta_S = 20.0 T0, P0 = 298.15, 1.0 dT_dP = delta_V / delta_S * 1e5 # K/bar T_eq = T0 + dT_dP * (P - P0) return max(T_eq, 273.15) # Minimum 0°C def _get_mineral_properties(self, mineral: str, T: float, P: float) -> Dict[str, float]: """Get thermodynamic properties of a mineral.""" # This would look up mineral in database # For now, return typical values typical_props = { 'quartz': {'G': -856300, 'H': -910700, 'S': 41.5, 'V': 22.7, 'Cp': 44.6}, 'calcite': {'G': -1128800, 'H': -1207400, 'S': 91.7, 'V': 36.9, 'Cp': 83.5}, 'pyrite': {'G': -160200, 'H': -171500, 'S': 52.9, 'V': 23.9, 'Cp': 62.2}, 'hematite': {'G': -744400, 'H': -824200, 'S': 87.4, 'V': 30.3, 'Cp': 103.9} } return typical_props.get(mineral.lower(), { 'G': -500000, 'H': -600000, 'S': 50, 'V': 30, 'Cp': 80 })Mineral equilibria calculator for geological systems.
This class handles mineral stability calculations, including: - Mineral solubility equilibria - Phase stability diagrams - Mineral-water reactions - Solid solution equilibria
Initialize the mineral equilibria calculator.
Methods
def mineral_solubility(self,
mineral: str,
T: float = 298.15,
P: float = 1.0,
pH_range: Tuple[float, float] | None = None,
ionic_strength: float = 0.0) ‑> Dict[str, Any]-
Expand source code
def mineral_solubility(self, mineral: str, T: float = 298.15, P: float = 1.0, pH_range: Optional[Tuple[float, float]] = None, ionic_strength: float = 0.0) -> Dict[str, Any]: """ Calculate mineral solubility as a function of pH. Parameters ---------- mineral : str Mineral name T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar pH_range : tuple, optional pH range (min, max). Default: (0, 14) ionic_strength : float, default 0.0 Solution ionic strength Returns ------- dict Solubility data vs pH """ if pH_range is None: pH_range = (0, 14) pH_values = np.linspace(pH_range[0], pH_range[1], 100) if mineral.lower() not in self.dissolution_reactions: warnings.warn(f"Dissolution reaction not defined for {mineral}") return {'pH': pH_values, 'solubility': np.zeros_like(pH_values)} reaction = self.dissolution_reactions[mineral.lower()] try: # Get equilibrium constant for dissolution species_names = list(reaction.keys()) coefficients = list(reaction.values()) # Calculate logK at T, P result = subcrt(species_names, coefficients, T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 0.0 warnings.warn(f"Could not calculate logK for {mineral} dissolution") except Exception as e: warnings.warn(f"Error calculating equilibrium constant: {e}") logK = 0.0 # Calculate solubility at each pH solubility = np.zeros_like(pH_values) for i, pH in enumerate(pH_values): try: # Simplified solubility calculation # This would be more complex in a full implementation if 'H+' in reaction: h_coeff = reaction['H+'] # Basic pH dependence log_solubility = logK - h_coeff * pH else: log_solubility = logK # Convert to molality solubility[i] = 10**log_solubility # Apply ionic strength corrections if needed if ionic_strength > 0: # Simple activity coefficient correction gamma = self._estimate_activity_coefficient(ionic_strength) solubility[i] *= gamma except: solubility[i] = 1e-20 # Very low solubility return { 'pH': pH_values, 'solubility': solubility, 'logK': logK, 'mineral': mineral, 'T': T, 'P': P }Calculate mineral solubility as a function of pH.
Parameters
mineral:str- Mineral name
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
pH_range:tuple, optional- pH range (min, max). Default: (0, 14)
ionic_strength:float, default0.0- Solution ionic strength
Returns
dict- Solubility data vs pH
def phase_equilibrium(self,
reaction: Dict[str, float],
T_range: Tuple[float, float] | None = None,
P_range: Tuple[float, float] | None = None) ‑> Dict[str, Any]-
Expand source code
def phase_equilibrium(self, reaction: Dict[str, float], T_range: Optional[Tuple[float, float]] = None, P_range: Optional[Tuple[float, float]] = None) -> Dict[str, Any]: """ Calculate phase equilibrium curve. Parameters ---------- reaction : dict Reaction stoichiometry {species: coefficient} T_range : tuple, optional Temperature range (min, max) in K P_range : tuple, optional Pressure range (min, max) in bar Returns ------- dict Equilibrium curve data """ if T_range is None and P_range is None: T_range = (273.15, 673.15) # 0-400°C if T_range is not None: # T-P curve at equilibrium T_values = np.linspace(T_range[0], T_range[1], 50) P_equilibrium = np.zeros_like(T_values) for i, T in enumerate(T_values): try: # Calculate equilibrium pressure # This would involve iterative solution # For now, use simplified correlation P_equilibrium[i] = self._calculate_equilibrium_pressure(reaction, T) except: P_equilibrium[i] = 1.0 # Default return { 'T': T_values, 'P': P_equilibrium, 'reaction': reaction } else: # P-T curve - similar logic P_values = np.linspace(P_range[0], P_range[1], 50) T_equilibrium = np.zeros_like(P_values) for i, P in enumerate(P_values): try: T_equilibrium[i] = self._calculate_equilibrium_temperature(reaction, P) except: T_equilibrium[i] = 298.15 return { 'P': P_values, 'T': T_equilibrium, 'reaction': reaction }Calculate phase equilibrium curve.
Parameters
reaction:dict- Reaction stoichiometry {species: coefficient}
T_range:tuple, optional- Temperature range (min, max) in K
P_range:tuple, optional- Pressure range (min, max) in bar
Returns
dict- Equilibrium curve data
def solid_solution(self,
end_members: List[str],
compositions: List[float],
T: float = 298.15,
P: float = 1.0,
model: str = 'ideal') ‑> Dict[str, Any]-
Expand source code
def solid_solution(self, end_members: List[str], compositions: List[float], T: float = 298.15, P: float = 1.0, model: str = 'ideal') -> Dict[str, Any]: """ Calculate solid solution thermodynamics. Parameters ---------- end_members : list of str End-member mineral names compositions : list of float Mole fractions of end-members (must sum to 1) T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar model : str, default 'ideal' Mixing model ('ideal', 'regular', 'margules') Returns ------- dict Solid solution properties """ if len(end_members) != len(compositions): raise ValueError("Number of end-members must match number of compositions") if abs(sum(compositions) - 1.0) > 1e-6: raise ValueError("Compositions must sum to 1.0") # Calculate end-member properties end_member_props = [] for mineral in end_members: try: # This would get mineral properties from database props = self._get_mineral_properties(mineral, T, P) end_member_props.append(props) except: # Default properties props = {'G': -1000000, 'H': -1200000, 'S': 100, 'V': 50, 'Cp': 100} end_member_props.append(props) # Calculate mixing properties if model == 'ideal': # Ideal mixing G_mix = 8.314 * T * sum(x * np.log(x) for x in compositions if x > 0) H_mix = 0.0 S_mix = -8.314 * sum(x * np.log(x) for x in compositions if x > 0) elif model == 'regular': # Regular solution (symmetric) W = 10000 # Interaction parameter (J/mol) - would be mineral-specific G_mix = W * compositions[0] * compositions[1] # Binary case H_mix = G_mix S_mix = 0.0 else: # Default to ideal G_mix = 8.314 * T * sum(x * np.log(x) for x in compositions if x > 0) H_mix = 0.0 S_mix = -8.314 * sum(x * np.log(x) for x in compositions if x > 0) # Total properties total_props = {} for prop in ['G', 'H', 'S', 'V', 'Cp']: total_props[prop] = sum(compositions[i] * end_member_props[i][prop] for i in range(len(end_members))) # Add mixing contributions total_props['G'] += G_mix total_props['H'] += H_mix total_props['S'] += S_mix return { 'properties': total_props, 'mixing': {'G': G_mix, 'H': H_mix, 'S': S_mix}, 'end_members': end_members, 'compositions': compositions, 'model': model }Calculate solid solution thermodynamics.
Parameters
end_members:listofstr- End-member mineral names
compositions:listoffloat- Mole fractions of end-members (must sum to 1)
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
model:str, default'ideal'- Mixing model ('ideal', 'regular', 'margules')
Returns
dict- Solid solution properties
def stability_diagram(self,
minerals: List[str],
x_axis: str,
x_range: Tuple[float, float],
y_axis: str,
y_range: Tuple[float, float],
T: float = 298.15,
P: float = 1.0,
resolution: int = 50) ‑> Dict[str, Any]-
Expand source code
def stability_diagram(self, minerals: List[str], x_axis: str, x_range: Tuple[float, float], y_axis: str, y_range: Tuple[float, float], T: float = 298.15, P: float = 1.0, resolution: int = 50) -> Dict[str, Any]: """ Create mineral stability diagram. Parameters ---------- minerals : list of str Mineral names to consider x_axis : str X-axis variable ('pH', 'log_activity_SiO2', etc.) x_range : tuple X-axis range (min, max) y_axis : str Y-axis variable y_range : tuple Y-axis range (min, max) T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar resolution : int, default 50 Grid resolution Returns ------- dict Stability diagram data """ x_values = np.linspace(x_range[0], x_range[1], resolution) y_values = np.linspace(y_range[0], y_range[1], resolution) X, Y = np.meshgrid(x_values, y_values) # Initialize stability field array stable_mineral = np.zeros_like(X, dtype=int) # Calculate stability at each grid point for i in range(len(y_values)): for j in range(len(x_values)): x_val, y_val = X[i, j], Y[i, j] # Find most stable mineral at this point min_energy = np.inf stable_idx = 0 for k, mineral in enumerate(minerals): try: # Calculate relative stability # This would involve full equilibrium calculations # For now, use simplified energy estimates energy = self._calculate_stability_energy( mineral, x_axis, x_val, y_axis, y_val, T, P) if energy < min_energy: min_energy = energy stable_idx = k except Exception as e: warnings.warn(f"Error calculating stability for {mineral}: {e}") continue stable_mineral[i, j] = stable_idx return { x_axis: x_values, y_axis: y_values, 'stable_mineral': stable_mineral, 'mineral_names': minerals, 'T': T, 'P': P }Create mineral stability diagram.
Parameters
minerals:listofstr- Mineral names to consider
x_axis:str- X-axis variable ('pH', 'log_activity_SiO2', etc.)
x_range:tuple- X-axis range (min, max)
y_axis:str- Y-axis variable
y_range:tuple- Y-axis range (min, max)
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
resolution:int, default50- Grid resolution
Returns
dict- Stability diagram data
class RedoxCalculator-
Expand source code
class RedoxCalculator: """ Redox equilibria calculator for geochemical systems. This class handles redox reactions, pe-pH diagrams, and electron activity calculations. """ def __init__(self): """Initialize the redox calculator.""" self.equilibrium_solver = EquilibriumSolver() # Standard electrode potentials (V) at 25°C self.standard_potentials = { 'O2/H2O': 1.229, 'H+/H2': 0.000, 'Fe+3/Fe+2': 0.771, 'NO3-/NO2-': 0.835, 'SO4-2/HS-': -0.217, 'CO2/CH4': -0.244, 'N2/NH4+': -0.277, 'Fe+2/Fe': -0.447, 'S/HS-': -0.065, } # Common redox couples and their reactions self.redox_reactions = { # Oxygen reduction 'O2_H2O': {'O2': 1, 'H+': 4, 'e-': 4, 'H2O': -2}, 'H2O_H2': {'H2O': 2, 'e-': 2, 'H2': -1, 'OH-': -2}, # Iron redox 'Fe3_Fe2': {'Fe+3': 1, 'e-': 1, 'Fe+2': -1}, 'Fe2_Fe': {'Fe+2': 1, 'e-': 2, 'Fe': -1}, 'Fe2O3_Fe2': {'Fe2O3': 1, 'H+': 6, 'e-': 2, 'Fe+2': -2, 'H2O': -3}, # Nitrogen redox 'NO3_NO2': {'NO3-': 1, 'H+': 2, 'e-': 2, 'NO2-': -1, 'H2O': -1}, 'NO2_NH4': {'NO2-': 1, 'H+': 8, 'e-': 6, 'NH4+': -1, 'H2O': -2}, # Sulfur redox 'SO4_HS': {'SO4-2': 1, 'H+': 9, 'e-': 8, 'HS-': -1, 'H2O': -4}, 'S_HS': {'S': 1, 'H+': 1, 'e-': 2, 'HS-': -1}, # Carbon redox 'CO2_CH4': {'CO2': 1, 'H+': 8, 'e-': 8, 'CH4': -1, 'H2O': -2}, 'HCO3_CH4': {'HCO3-': 1, 'H+': 9, 'e-': 8, 'CH4': -1, 'H2O': -3} } def eh_ph_diagram(self, element: str, pH_range: Tuple[float, float] = (0, 14), pe_range: Tuple[float, float] = (-10, 15), T: float = 298.15, P: float = 1.0, total_concentration: float = 1e-6, resolution: int = 100) -> Dict[str, Any]: """ Create Eh-pH (pe-pH) diagram for an element. Parameters ---------- element : str Element symbol (e.g., 'Fe', 'S', 'N', 'C') pH_range : tuple, default (0, 14) pH range for diagram pe_range : tuple, default (-10, 15) pe (electron activity) range T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar total_concentration : float, default 1e-6 Total element concentration (molal) resolution : int, default 100 Grid resolution Returns ------- dict Eh-pH diagram data """ pH_grid = np.linspace(pH_range[0], pH_range[1], resolution) pe_grid = np.linspace(pe_range[0], pe_range[1], resolution) pH_mesh, pe_mesh = np.meshgrid(pH_grid, pe_grid) # Get species for this element element_species = self._get_element_species(element) if not element_species: raise ValueError(f"No species found for element {element}") # Calculate predominance at each point predominant = np.zeros_like(pH_mesh, dtype=int) activities = {species: np.zeros_like(pH_mesh) for species in element_species} for i in range(len(pe_grid)): for j in range(len(pH_grid)): pH, pe = pH_mesh[i, j], pe_mesh[i, j] # Calculate speciation at this point spec_result = self._calculate_redox_speciation( element_species, pH, pe, T, P, total_concentration) # Find predominant species max_activity = -np.inf max_idx = 0 for k, species in enumerate(element_species): activity = spec_result.get(species, 1e-20) activities[species][i, j] = activity if np.log10(activity) > max_activity: max_activity = np.log10(activity) max_idx = k predominant[i, j] = max_idx # Add water stability limits water_limits = self._water_stability_lines(pH_grid, T, P) return { 'pH': pH_grid, 'pe': pe_grid, 'pH_mesh': pH_mesh, 'pe_mesh': pe_mesh, 'predominant': predominant, 'activities': activities, 'species_names': element_species, 'water_limits': water_limits, 'element': element, 'T': T, 'P': P, 'total_concentration': total_concentration } def pe_calculation(self, redox_couple: Union[str, Dict[str, float]], concentrations: Dict[str, float], pH: float = 7.0, T: float = 298.15, P: float = 1.0) -> float: """ Calculate pe (electron activity) for a redox couple. Parameters ---------- redox_couple : str or dict Redox couple name or reaction dictionary concentrations : dict Species concentrations {species: concentration} pH : float, default 7.0 Solution pH T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- float pe value """ if isinstance(redox_couple, str): if redox_couple not in self.redox_reactions: raise ValueError(f"Unknown redox couple: {redox_couple}") reaction = self.redox_reactions[redox_couple] else: reaction = redox_couple # Calculate equilibrium constant try: species_names = [sp for sp in reaction.keys() if sp != 'e-'] coefficients = [reaction[sp] for sp in species_names] result = subcrt(species_names, coefficients, T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 0.0 except Exception as e: warnings.warn(f"Could not calculate logK: {e}") logK = 0.0 # Apply Nernst equation n_electrons = abs(reaction.get('e-', 1)) # Number of electrons # Calculate activity quotient log_Q = 0.0 for species, coeff in reaction.items(): if species == 'e-': continue elif species == 'H+': activity = 10**(-pH) elif species in concentrations: activity = concentrations[species] # Apply activity coefficients if needed gamma = self._get_activity_coefficient(species, concentrations, T) activity *= gamma else: activity = 1.0 # Default for species not specified if activity > 0: log_Q += coeff * np.log10(activity) # Nernst equation: pe = pe° + (1/n) * log(Q) # where pe° = logK/n for the half-reaction pe = logK / n_electrons + log_Q / n_electrons return pe def eh_from_pe(self, pe: float, T: float = 298.15) -> float: """ Convert pe to Eh (redox potential). Parameters ---------- pe : float Electron activity (pe) T : float, default 298.15 Temperature in Kelvin Returns ------- float Redox potential (Eh) in Volts """ # Eh = (RT/F) * ln(10) * pe # where R = 8.314 J/(mol·K), F = 96485 C/mol RT_F = 8.314 * T / 96485 Eh = RT_F * 2.302585 * pe # ln(10) = 2.302585 return Eh def pe_from_eh(self, eh: float, T: float = 298.15) -> float: """ Convert Eh (redox potential) to pe. Parameters ---------- eh : float Redox potential (Eh) in Volts T : float, default 298.15 Temperature in Kelvin Returns ------- float Electron activity (pe) """ # pe = F * Eh / (RT * ln(10)) RT_F = 8.314 * T / 96485 pe = eh / (RT_F * 2.302585) return pe def oxygen_fugacity(self, pe: float, pH: float = 7.0, T: float = 298.15, P: float = 1.0) -> float: """ Calculate oxygen fugacity from pe and pH. Parameters ---------- pe : float Electron activity pH : float, default 7.0 Solution pH T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- float log fO2 (log oxygen fugacity) """ # O2 + 4H+ + 4e- = 2H2O # At equilibrium: log fO2 = 4*pe + 4*pH - logK try: # Calculate logK for oxygen-water reaction result = subcrt(['O2', 'H+', 'H2O'], [1, 4, -2], T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 83.1 # Approximate value at 25°C except: logK = 83.1 log_fO2 = 4 * pe + 4 * pH - logK return log_fO2 def _get_element_species(self, element: str) -> List[str]: """Get list of species containing the specified element.""" # Simplified species lists for common elements element_species = { 'Fe': ['Fe+2', 'Fe+3', 'Fe2O3', 'FeOH+', 'Fe(OH)2', 'Fe(OH)3'], 'S': ['SO4-2', 'SO3-2', 'S2O3-2', 'HS-', 'S0', 'S-2'], 'N': ['NO3-', 'NO2-', 'NH4+', 'NH3', 'N2O', 'N2'], 'C': ['CO2', 'HCO3-', 'CO3-2', 'CH4', 'HCOOH', 'CH3COO-'], 'Mn': ['Mn+2', 'MnO4-', 'MnO2', 'Mn+3', 'MnOH+'], 'As': ['AsO4-3', 'AsO3-3', 'H3AsO4', 'H3AsO3', 'As0'] } return element_species.get(element, [f'{element}+2', f'{element}+3']) def _calculate_redox_speciation(self, species: List[str], pH: float, pe: float, T: float, P: float, total_conc: float) -> Dict[str, float]: """Calculate speciation at given pH and pe.""" # Simplified speciation calculation # Full implementation would solve equilibrium system activities = {} for sp in species: # Simplified pe-pH dependence if '+' in sp: # Cations (more stable at low pH, high pe) charge = self._extract_charge(sp) log_activity = -6 + charge * (14 - pH) / 14 + pe / 10 elif '-' in sp: # Anions (more stable at high pH, high pe) charge = abs(self._extract_charge(sp)) log_activity = -6 + charge * pH / 14 + pe / 10 else: # Neutral (less pH/pe dependence) log_activity = -6 + (pe - pH) / 20 activities[sp] = 10**log_activity # Normalize to total concentration total_calculated = sum(activities.values()) if total_calculated > 0: factor = total_conc / total_calculated activities = {sp: act * factor for sp, act in activities.items()} return activities def _extract_charge(self, species: str) -> int: """Extract charge from species name.""" if '+' in species: parts = species.split('+') if len(parts) > 1 and parts[-1].isdigit(): return int(parts[-1]) else: return 1 elif '-' in species: parts = species.split('-') if len(parts) > 1 and parts[-1].isdigit(): return -int(parts[-1]) else: return -1 else: return 0 def _water_stability_lines(self, pH_values: np.ndarray, T: float, P: float) -> Dict[str, np.ndarray]: """Calculate water stability limits (H2O/H2 and O2/H2O).""" # Upper limit: O2/H2O # O2 + 4H+ + 4e- = 2H2O # pe = 20.75 - pH (at 25°C, 1 atm O2) pe_upper = 20.75 - pH_values # Lower limit: H2O/H2 # 2H2O + 2e- = H2 + 2OH- # pe = -pH (at 25°C, 1 atm H2) pe_lower = -pH_values # Temperature corrections (simplified) if T != 298.15: dT = T - 298.15 pe_upper += dT * 0.001 # Small temperature dependence pe_lower -= dT * 0.001 return { 'upper': pe_upper, # O2/H2O line 'lower': pe_lower # H2O/H2 line } def _get_activity_coefficient(self, species: str, composition: Dict[str, float], T: float) -> float: """Get activity coefficient for species.""" # Simplified - would use proper activity models charge = abs(self._extract_charge(species)) if charge == 0: return 1.0 else: # Simple ionic strength correction I = 0.001 # Assume low ionic strength return 10**(-0.509 * charge**2 * np.sqrt(I))Redox equilibria calculator for geochemical systems.
This class handles redox reactions, pe-pH diagrams, and electron activity calculations.
Initialize the redox calculator.
Methods
def eh_from_pe(self, pe: float, T: float = 298.15) ‑> float-
Expand source code
def eh_from_pe(self, pe: float, T: float = 298.15) -> float: """ Convert pe to Eh (redox potential). Parameters ---------- pe : float Electron activity (pe) T : float, default 298.15 Temperature in Kelvin Returns ------- float Redox potential (Eh) in Volts """ # Eh = (RT/F) * ln(10) * pe # where R = 8.314 J/(mol·K), F = 96485 C/mol RT_F = 8.314 * T / 96485 Eh = RT_F * 2.302585 * pe # ln(10) = 2.302585 return EhConvert pe to Eh (redox potential).
Parameters
pe:float- Electron activity (pe)
T:float, default298.15- Temperature in Kelvin
Returns
float- Redox potential (Eh) in Volts
def eh_ph_diagram(self,
element: str,
pH_range: Tuple[float, float] = (0, 14),
pe_range: Tuple[float, float] = (-10, 15),
T: float = 298.15,
P: float = 1.0,
total_concentration: float = 1e-06,
resolution: int = 100) ‑> Dict[str, Any]-
Expand source code
def eh_ph_diagram(self, element: str, pH_range: Tuple[float, float] = (0, 14), pe_range: Tuple[float, float] = (-10, 15), T: float = 298.15, P: float = 1.0, total_concentration: float = 1e-6, resolution: int = 100) -> Dict[str, Any]: """ Create Eh-pH (pe-pH) diagram for an element. Parameters ---------- element : str Element symbol (e.g., 'Fe', 'S', 'N', 'C') pH_range : tuple, default (0, 14) pH range for diagram pe_range : tuple, default (-10, 15) pe (electron activity) range T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar total_concentration : float, default 1e-6 Total element concentration (molal) resolution : int, default 100 Grid resolution Returns ------- dict Eh-pH diagram data """ pH_grid = np.linspace(pH_range[0], pH_range[1], resolution) pe_grid = np.linspace(pe_range[0], pe_range[1], resolution) pH_mesh, pe_mesh = np.meshgrid(pH_grid, pe_grid) # Get species for this element element_species = self._get_element_species(element) if not element_species: raise ValueError(f"No species found for element {element}") # Calculate predominance at each point predominant = np.zeros_like(pH_mesh, dtype=int) activities = {species: np.zeros_like(pH_mesh) for species in element_species} for i in range(len(pe_grid)): for j in range(len(pH_grid)): pH, pe = pH_mesh[i, j], pe_mesh[i, j] # Calculate speciation at this point spec_result = self._calculate_redox_speciation( element_species, pH, pe, T, P, total_concentration) # Find predominant species max_activity = -np.inf max_idx = 0 for k, species in enumerate(element_species): activity = spec_result.get(species, 1e-20) activities[species][i, j] = activity if np.log10(activity) > max_activity: max_activity = np.log10(activity) max_idx = k predominant[i, j] = max_idx # Add water stability limits water_limits = self._water_stability_lines(pH_grid, T, P) return { 'pH': pH_grid, 'pe': pe_grid, 'pH_mesh': pH_mesh, 'pe_mesh': pe_mesh, 'predominant': predominant, 'activities': activities, 'species_names': element_species, 'water_limits': water_limits, 'element': element, 'T': T, 'P': P, 'total_concentration': total_concentration }Create Eh-pH (pe-pH) diagram for an element.
Parameters
element:str- Element symbol (e.g., 'Fe', 'S', 'N', 'C')
pH_range:tuple, default(0, 14)- pH range for diagram
pe_range:tuple, default(-10, 15)- pe (electron activity) range
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
total_concentration:float, default1e-6- Total element concentration (molal)
resolution:int, default100- Grid resolution
Returns
dict- Eh-pH diagram data
def oxygen_fugacity(self, pe: float, pH: float = 7.0, T: float = 298.15, P: float = 1.0) ‑> float-
Expand source code
def oxygen_fugacity(self, pe: float, pH: float = 7.0, T: float = 298.15, P: float = 1.0) -> float: """ Calculate oxygen fugacity from pe and pH. Parameters ---------- pe : float Electron activity pH : float, default 7.0 Solution pH T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- float log fO2 (log oxygen fugacity) """ # O2 + 4H+ + 4e- = 2H2O # At equilibrium: log fO2 = 4*pe + 4*pH - logK try: # Calculate logK for oxygen-water reaction result = subcrt(['O2', 'H+', 'H2O'], [1, 4, -2], T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 83.1 # Approximate value at 25°C except: logK = 83.1 log_fO2 = 4 * pe + 4 * pH - logK return log_fO2Calculate oxygen fugacity from pe and pH.
Parameters
pe:float- Electron activity
pH:float, default7.0- Solution pH
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
Returns
float- log fO2 (log oxygen fugacity)
def pe_calculation(self,
redox_couple: str | Dict[str, float],
concentrations: Dict[str, float],
pH: float = 7.0,
T: float = 298.15,
P: float = 1.0) ‑> float-
Expand source code
def pe_calculation(self, redox_couple: Union[str, Dict[str, float]], concentrations: Dict[str, float], pH: float = 7.0, T: float = 298.15, P: float = 1.0) -> float: """ Calculate pe (electron activity) for a redox couple. Parameters ---------- redox_couple : str or dict Redox couple name or reaction dictionary concentrations : dict Species concentrations {species: concentration} pH : float, default 7.0 Solution pH T : float, default 298.15 Temperature in Kelvin P : float, default 1.0 Pressure in bar Returns ------- float pe value """ if isinstance(redox_couple, str): if redox_couple not in self.redox_reactions: raise ValueError(f"Unknown redox couple: {redox_couple}") reaction = self.redox_reactions[redox_couple] else: reaction = redox_couple # Calculate equilibrium constant try: species_names = [sp for sp in reaction.keys() if sp != 'e-'] coefficients = [reaction[sp] for sp in species_names] result = subcrt(species_names, coefficients, T=T, P=P, show=False) if result.out is not None and 'logK' in result.out.columns: logK = result.out['logK'].iloc[0] else: logK = 0.0 except Exception as e: warnings.warn(f"Could not calculate logK: {e}") logK = 0.0 # Apply Nernst equation n_electrons = abs(reaction.get('e-', 1)) # Number of electrons # Calculate activity quotient log_Q = 0.0 for species, coeff in reaction.items(): if species == 'e-': continue elif species == 'H+': activity = 10**(-pH) elif species in concentrations: activity = concentrations[species] # Apply activity coefficients if needed gamma = self._get_activity_coefficient(species, concentrations, T) activity *= gamma else: activity = 1.0 # Default for species not specified if activity > 0: log_Q += coeff * np.log10(activity) # Nernst equation: pe = pe° + (1/n) * log(Q) # where pe° = logK/n for the half-reaction pe = logK / n_electrons + log_Q / n_electrons return peCalculate pe (electron activity) for a redox couple.
Parameters
redox_couple:strordict- Redox couple name or reaction dictionary
concentrations:dict- Species concentrations {species: concentration}
pH:float, default7.0- Solution pH
T:float, default298.15- Temperature in Kelvin
P:float, default1.0- Pressure in bar
Returns
float- pe value
def pe_from_eh(self, eh: float, T: float = 298.15) ‑> float-
Expand source code
def pe_from_eh(self, eh: float, T: float = 298.15) -> float: """ Convert Eh (redox potential) to pe. Parameters ---------- eh : float Redox potential (Eh) in Volts T : float, default 298.15 Temperature in Kelvin Returns ------- float Electron activity (pe) """ # pe = F * Eh / (RT * ln(10)) RT_F = 8.314 * T / 96485 pe = eh / (RT_F * 2.302585) return peConvert Eh (redox potential) to pe.
Parameters
eh:float- Redox potential (Eh) in Volts
T:float, default298.15- Temperature in Kelvin
Returns
float- Electron activity (pe)