{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Benchmarking models for equilibrium constants\n",
    "\n",
    "This notebook benchmarks the models for equilibrium constants in VolFe where possible.\n",
    "\n",
    "The following equilibrium constants and models do not have material available in the original papers for benchmarking:\n",
    "- KHOg: 'Ohmoto97' Reaction (d) in Table 1 from Ohmoto & Kerrick (1977)\n",
    "- KHOSg: 'Ohmoto97' Reaction (h) in Table 1 from Ohmoto & Kerrick (1977)\n",
    "- KOSg: 'Ohmoto97' Reaction (f) in Table 1 from Ohmoto & Kerrick (1977)\n",
    "- KOCg: 'Ohmoto97' Reaction (c) in Table 1 from Ohmoto & Kerrick (1977)\n",
    "- KCOHg: 'Ohmoto97' Reaction (e) in Table 1 from Ohmoto & Kerrick (1977)\n",
    "- KCOs: 'Holloway92' Eq (3) KI in Holloway et al. (1992)\n",
    "\n",
    "## Python set-up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import VolFe as vf\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for equilibrium constant of SO3 \n",
    "\n",
    "option = KOSg2, function = KOSg2\n",
    "\n",
    "### 'ONeill22' Eq. (6b) from O'Neill & Mavrogenes (2022)\n",
    "\n",
    "Appendix A. Supplementary material - Supplementary data 2. Tab = Table S6 S redox calculator (sulfate,KSOg2=ONeill22.xlsx).\n",
    "\n",
    "Cell AI12: -17.61\n",
    "\n",
    "Matches to 2 decimal places. Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 'C\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-17.606473217676754"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PT = {\"P\":1000.}\n",
    "PT[\"T\"]=1200.\n",
    "\n",
    "my_models = [[\"KOSg2\", \"ONeill22\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "math.log(1./vf.KOSg2(PT,models=my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Models for equilibrium constant of OCS \n",
    "\n",
    "option = KOCSg, function = KOCSg\n",
    "\n",
    "### 'Moussallam19' Eq. (8) for 2CO2 + OCS ⇄ 3CO + SO2 in Moussallam et al. (2019)\n",
    "\n",
    "From Section 2. Method - Detailed examples of gas oxygen fugacity and equilibrium temperature calculationd - The CO2-CO-OCS-SO2 method\n",
    "\n",
    "\"Given a XSO2 of 0.01, a CO/CO2 molar ratio of 0.043, a CO/OCS molar ratio of 192 and assuming equilibration at atmospheric pressure (∼0.6 bar at Erebus) yields an equilibrium temperature of 1292 K (1019 °C).\"\n",
    "\n",
    "Equation (8):\n",
    "\n",
    "log P = (-15386.45/T) + 9.24403 - log[((xCO/xCO2)^2)(xCO/xOCS)xSO2]\n",
    "\n",
    "K = (-15386.45/T) + 9.24403 = log P + log[((xCO/xCO2)^2)(xCO/xOCS)xSO2]\n",
    "\n",
    "The next cell calculates T from their eq. (17), which they state is 1292 K"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1291.2818121927335"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P = 0.6 # bar\n",
    "xSO2 = 0.01\n",
    "xCO_xCO2 = 0.043\n",
    "xCO_xOCS = 192.\n",
    "\n",
    "-15386.45/(math.log10(P) - 9.24403 + math.log10((xCO_xCO2**2.)*xCO_xOCS*xSO2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below we calculate K for their example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.671610609753634"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "math.log10(P) + math.log10(xCO_xCO2**2.*xCO_xOCS*xSO2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which agrees to the first decimal place with VolFe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-2.6636045625507876"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "PT = {\"P\":0.6}\n",
    "PT[\"T\"]=1019.\n",
    "\n",
    "my_models = [[\"KOCSg\", \"Moussallam19\"]]\n",
    "my_models = vf.make_df_and_add_model_defaults(my_models)\n",
    "math.log10(vf.KOCSg(PT,models=my_models))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "volfe-dev",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}