Peng_Robinson_Fugacity_and_Global_Warming.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "provenance": [],
      "authorship_tag": "ABX9TyPe4mjtWk3qrJ/OvtelAPcN",
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/walkerjian/Physics/blob/main/Peng_Robinson_Fugacity_and_Global_Warming.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fugacity is a concept used in thermodynamics to describe the \"escaping tendency\" of a substance from its current phase. It is often used in the context of phase equilibria and can be thought of as an effective or \"corrected\" pressure that takes into account non-ideal behavior of real gases. The term \"fugacity\" is derived from the Latin word \"fugere,\" which means \"to flee.\"\n",
        "\n",
        "In the context of global warming, fugacity can be applied to understand the behavior of greenhouse gases in the atmosphere. The Clausius-Clapeyron equation relates the vapor pressure of a substance to its temperature and is used to describe phase transitions, especially between liquid and vapor phases.\n",
        "\n",
        "**Clausius-Clapeyron Equation:**\n",
        "The differential form of the Clausius-Clapeyron equation is given by:\n",
        "\n",
        "$$\n",
        "\\frac{d \\ln P}{dT} = \\frac{\\Delta H_{vap}}{RT^2}\n",
        "$$\n",
        "Where:\n",
        "- $ P $ is the vapor pressure.\n",
        "- $ T $ is the absolute temperature.\n",
        "- $ R $ is the universal gas constant.\n",
        "- $ \\Delta H_{vap} $ is the enthalpy of vaporization.\n",
        "\n",
        "**Example Calculation:**\n",
        "Suppose we want to determine the change in vapor pressure of water with temperature, given that the enthalpy of vaporization for water is approximately $ 40.79 \\, \\text{kJ/mol} $\n",
        "and the vapor pressure at\n",
        "$ 100^\\circ \\text{C} \\, (373.15 \\text{K})$\n",
        "is\n",
        "$ 1 \\text{atm} $.\n",
        "\n",
        "Using the Clausius-Clapeyron equation, the change in vapor pressure with a small change in temperature $ dT $ can be calculated as:\n",
        "$$\n",
        "d \\ln P = \\frac{\\Delta H_{vap}}{RT^2} dT\n",
        "$$\n",
        "\n",
        "Integrating this equation between two temperatures $T_{1}$ and $T_2$ gives:\n",
        "$$\n",
        "\\ln \\left( \\frac{P_2}{P_1} \\right) = \\frac{\\Delta H_{vap}}{R} \\left( \\frac{1}{T_1} - \\frac{1}{T_2} \\right)\n",
        "$$\n",
        "\n",
        "Given $P_1 = 1 \\, \\text{atm} $ at $ T_1 = 373.15 \\, \\text{K} $, we can calculate $P_2$ at a different temperature $T_2$ using the above equation.\n",
        "\n",
        "This calculation can provide insights into how the vapor pressure of water (or other greenhouse gases) changes with temperature, which is crucial for understanding the dynamics of global warming.\n",
        "\n",
        "It's worth noting that while the Clausius-Clapeyron equation provides a fundamental relationship between temperature and vapor pressure, fugacity is a more general concept that can be applied to any phase and any substance, not just vapors. In the context of global warming, fugacity can be used to understand the behavior of gases in the atmosphere, especially when they deviate from ideal behavior."
      ],
      "metadata": {
        "id": "lwQUoWGQQBm4"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Is it not true that the carrying capacity of the atmosphere for water vapour increases exponentially as the temperature rises, and thus in conjunction with fugacity, indicates a scenario of accelerating global warming?\n",
        "\n"
      ],
      "metadata": {
        "id": "QPcXG7vtX9Fb"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Yes, you're correct. The carrying capacity of the atmosphere for water vapor does increase with rising temperatures, and this relationship is often described using the Clausius-Clapeyron equation. The equation shows that the saturation vapor pressure (essentially the maximum amount of water vapor the air can hold) increases roughly exponentially with temperature. This means that for each degree of warming, the air can hold a progressively larger amount of water vapor.\n",
        "\n",
        "Water vapor itself is a potent greenhouse gas. As the temperature rises and the atmosphere can hold more water vapor, this increased water vapor concentration amplifies the greenhouse effect, leading to further warming. This is known as a positive feedback mechanism. Here's a simplified breakdown:\n",
        "\n",
        "1. **Initial Warming**: Due to any reason, such as increased CO₂ concentrations.\n",
        "2. **Increased Water Vapor Capacity**: As the atmosphere warms, it can hold more water vapor.\n",
        "3. **More Water Vapor**: The warmer atmosphere leads to increased evaporation and thus more water vapor in the air.\n",
        "4. **Enhanced Greenhouse Effect**: Water vapor is a greenhouse gas, so its increased concentration enhances the greenhouse effect.\n",
        "5. **Further Warming**: This enhanced greenhouse effect causes further warming, which in turn allows the atmosphere to hold even more water vapor, and the cycle continues.\n",
        "\n",
        "This water vapor feedback is one of the primary reasons why small changes in factors like CO₂ concentrations can lead to significant changes in global temperatures. The feedback mechanism accelerates the rate of global warming.\n",
        "\n",
        "However, it's essential to note that there are other feedback mechanisms, both positive and negative, in the climate system. For instance, increased cloud cover due to higher water vapor concentrations can reflect more sunlight back to space, providing a cooling effect. The net effect on global warming depends on the balance of all these feedbacks. But the water vapor feedback is generally considered one of the most potent positive feedbacks in the climate system."
      ],
      "metadata": {
        "id": "jYaLscpsYECr"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### How does fugacity affect global warming as temperatures increase?\n",
        "\n"
      ],
      "metadata": {
        "id": "Rvz8NsQiYPzS"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fugacity is a measure of the \"escaping tendency\" or \"effective pressure\" of a substance from its current phase. In the context of gases in the atmosphere, fugacity can be thought of as a corrected pressure that accounts for non-ideal behavior of real gases. While fugacity itself doesn't directly cause global warming, it plays a role in understanding the behavior of greenhouse gases in the atmosphere, especially under varying conditions.\n",
        "\n",
        "Here's how fugacity relates to global warming as temperatures increase:\n",
        "\n",
        "1. **Non-Ideal Behavior of Gases**: As the concentration of greenhouse gases (GHGs) in the atmosphere increases, their behavior may deviate from the ideal gas law, especially at higher pressures or in the presence of other gases. Fugacity helps account for these deviations.\n",
        "\n",
        "2. **Partitioning Between Phases**: Many GHGs can exist in multiple phases in the environment (e.g., carbon dioxide in the atmosphere and dissolved in the ocean). Fugacity can be used to understand the partitioning of these gases between different phases. For instance, as temperatures rise, the solubility of gases like CO₂ in the ocean decreases, leading to a higher fugacity in the aqueous phase. This can result in more CO₂ being released from the oceans into the atmosphere, further intensifying the greenhouse effect.\n",
        "\n",
        "3. **Feedback Mechanisms**: The concept of fugacity can be applied to water vapor, a significant greenhouse gas. As temperatures increase, the carrying capacity (or saturation vapor pressure) of the atmosphere for water vapor also increases. This means the atmosphere can hold more water vapor at higher temperatures. Since water vapor has a higher fugacity at higher temperatures, more of it will tend to remain in the gaseous phase rather than condensing. This increased water vapor concentration amplifies the greenhouse effect, leading to further warming.\n",
        "\n",
        "4. **Transport and Diffusion**: Fugacity can also be used in models to understand the transport and diffusion of GHGs in the atmosphere and between other environmental compartments. These models can provide insights into how GHGs spread and how their concentrations change over time and space.\n",
        "\n",
        "In summary, while fugacity itself doesn't directly cause global warming, it's a useful concept in understanding the behavior and distribution of greenhouse gases in the environment, especially under changing temperature conditions. By accounting for the non-ideal behavior of gases and their interactions with other phases, fugacity provides a more accurate picture of the dynamics of GHGs in the context of global warming."
      ],
      "metadata": {
        "id": "5tfHPCerYTu3"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "I take it that there has not been much research on how fugacity affects greenhouse effect?"
      ],
      "metadata": {
        "id": "XvCYeViEYbBA"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "The concept of fugacity is fundamental in thermodynamics and is widely used in chemical engineering and physical chemistry to describe the behavior of substances in various phases, especially when deviations from ideal behavior occur. However, when it comes to the broader topic of the greenhouse effect and global warming, the primary focus has typically been on the concentration, radiative properties, and feedback mechanisms associated with greenhouse gases (GHGs) rather than their fugacities per se.\n",
        "\n",
        "A few reasons for this are:\n",
        "\n",
        "1. **Scale of Interest**: Climate science often deals with large-scale processes and phenomena. While fugacity is crucial for understanding microscale or localized phase behaviors, the macroscopic effects on climate are more influenced by the total concentrations of GHGs, their radiative forcing properties, and the feedback mechanisms in the Earth's climate system.\n",
        "\n",
        "2. **Ideal Behavior Approximation**: For many atmospheric applications, gases can be approximated as behaving ideally, especially at the relatively low pressures found in the atmosphere. In such cases, the distinction between pressure and fugacity becomes less critical.\n",
        "\n",
        "3. **Complexity of Climate System**: The Earth's climate system is influenced by a multitude of factors, including GHG concentrations, solar radiation, cloud cover, albedo, oceanic processes, and more. While fugacity can play a role in understanding the partitioning of gases between phases (e.g., between the atmosphere and oceans), it's just one of many factors that need to be considered.\n",
        "\n",
        "That said, in specific contexts, especially where phase equilibria are important (like CO₂ dissolution in the oceans or methane clathrate stability), the concept of fugacity can be relevant. It might be used in specialized studies or models that look at these processes in detail.\n",
        "\n",
        "In summary, while fugacity is a valuable concept in thermodynamics and can be relevant in specific contexts related to the greenhouse effect, the broader field of climate science tends to focus on other more directly impactful factors when studying global warming and the greenhouse effect."
      ],
      "metadata": {
        "id": "igUnFBNCYfea"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Derive an expression for fugacity\n",
        "\n"
      ],
      "metadata": {
        "id": "6UKpYmkqYmeU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Fugacity is a thermodynamic property used to describe the \"escaping tendency\" of a substance from its current phase. It's especially useful for real gases that do not behave ideally. The concept of fugacity bridges the gap between the ideal gas law and the real behavior of gases.\n",
        "\n",
        "For an ideal gas, the fugacity $(f)$ is equal to its partial pressure $(P)$. However, for real gases, this isn't the case, and we need to introduce a fugacity coefficient $(\\phi)$ to account for the deviation from ideality:\n",
        "\n",
        "$ f = \\phi P $\n",
        "\n",
        "The fugacity coefficient is defined as:\n",
        "\n",
        "$ \\phi = \\frac{f}{P} $\n",
        "\n",
        "Where:\n",
        "- $f$ is the fugacity of the gas.\n",
        "- $P$ is the pressure of the gas.\n",
        "\n",
        "For real gases, the relationship between fugacity, pressure, and the fugacity coefficient can be derived from the equation of state and the Gibbs-Duhem equation. One common approach is to use the van der Waals equation of state or other similar equations to account for the non-ideality of gases.\n",
        "\n",
        "**Derivation using the van der Waals equation:**\n",
        "\n",
        "The van der Waals equation for a real gas is:\n",
        "\n",
        "$\\left( P + \\frac{a}{V^2} \\right) (V - b) = RT $\n",
        "\n",
        "Where:\n",
        "- $P$ is the pressure.\n",
        "- $V$ is the molar volume.\n",
        "- $R$ is the universal gas constant.\n",
        "- $T$ is the absolute temperature.\n",
        "- $a$ and $b$ are van der Waals constants specific to the gas.\n",
        "\n",
        "To derive an expression for fugacity, we would integrate the differential form of the van der Waals equation along an isotherm, considering both the ideal and real behaviors of the gas. This process involves integrating the equation with respect to volume and then exponentiating to get the fugacity in terms of pressure, volume, and temperature.\n",
        "\n",
        "The resulting expression will relate the fugacity of the gas to its pressure, molar volume, temperature, and the van der Waals constants. The exact form of this expression can be quite complex and may require numerical methods for evaluation.\n",
        "\n",
        "In practice, the fugacity coefficient $( \\phi )$ is often determined experimentally or using more advanced equations of state, and the fugacity is then calculated using the relationship $ f = \\phi P $."
      ],
      "metadata": {
        "id": "hCvH3ITXYqce"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Deriving an expression for fugacity using the Peng-Robinson equation of state, which is one of the more advanced and commonly used equations of state for real gases.\n",
        "The Peng-Robinson equation is given by:\n",
        "\n",
        "$ P = \\frac{RT}{V_m - b} - \\frac{a \\alpha}{V_m^2 + 2bV_m - b^2} $\n",
        "\n",
        "Where:\n",
        "- $ P $ is the pressure.\n",
        "- $ R $ is the universal gas constant.\n",
        "- $ T $ is the absolute temperature.\n",
        "- $ V_m $ is the molar volume.\n",
        "- $ a $ and $ b $ are Peng-Robinson constants.\n",
        "- $ \\alpha $ is a temperature-dependent function.\n",
        "\n",
        "To derive an expression for fugacity, we'll start with the differential form of the Gibbs-Duhem equation:\n",
        "\n",
        "$ d\\ln(f) = \\frac{V_m}{RT}dP $\n",
        "\n",
        "Integrating this equation along an isotherm gives:\n",
        "\n",
        "$ \\ln\\left(\\frac{f_2}{f_1}\\right) = \\frac{1}{RT} \\int_{P_1}^{P_2} V_m dP $\n",
        "\n",
        "Where:\n",
        "- $ f_1 $ and $ f_2 $ are the fugacities at pressures $ P_1 $ and $ P_2 $, respectively.\n",
        "\n",
        "Using the Peng-Robinson equation, we can express $ V_m $ as a function of $ P $ and then integrate.\n",
        "\n",
        "However, the integration is non-trivial and typically requires numerical methods. For simplicity, we'll outline the general approach without diving into the detailed integration:\n",
        "\n",
        "1. Solve the Peng-Robinson equation for $ V_m $ as a function of $ P $ and $ T $. This will give three roots, but we're interested in the root corresponding to the gas phase.\n",
        "2. Substitute this expression for $ V_m $ into the integral.\n",
        "3. Integrate with respect to $ P $ from $ P_1 $ to $ P_2 $.\n",
        "4. Use the resulting expression to calculate the ratio $ \\frac{f_2}{f_1} $.\n",
        "\n",
        "The final expression will relate the fugacity at two different pressures based on the Peng-Robinson equation of state.\n",
        "\n",
        "Given the complexity of the Peng-Robinson equation, exact analytical solutions are challenging to obtain, and numerical methods are often employed in practice."
      ],
      "metadata": {
        "id": "UBQLE5JPbcLj"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "A numerical solution using the Peng-Robinson equation of state to compute the molar volume $ V_m $ as a function of $ P $ and $ T $, and then integrate to determine the fugacity ratio $ \\frac{f_2}{f_1} $.\n",
        "\n",
        "```python\n",
        "import numpy as np\n",
        "from scipy.optimize import fsolve\n",
        "from scipy.integrate import quad\n",
        "\n",
        "# Peng-Robinson constants for a hypothetical gas (you can replace these with actual values for a specific gas)\n",
        "R = 8.314  # Universal gas constant [J/(mol*K)]\n",
        "a = 3.6  # [Pa*m^6/mol^2]\n",
        "b = 0.04  # [m^3/mol]\n",
        "Tc = 500  # Critical temperature [K]\n",
        "Pc = 3e6  # Critical pressure [Pa]\n",
        "omega = 0.225  # Acentric factor\n",
        "\n",
        "# Temperature-dependent function alpha\n",
        "def alpha(T):\n",
        "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
        "    return (1 + kappa*(1 - np.sqrt(T/Tc)))**2\n",
        "\n",
        "# Peng-Robinson equation of state\n",
        "def PR_eqn(Vm, P, T):\n",
        "    return P - R*T/(Vm - b) + a*alpha(T)/((Vm**2) + 2*b*Vm - b**2)\n",
        "\n",
        "# Compute molar volume Vm for a given P and T\n",
        "def compute_Vm(P, T):\n",
        "    # Initial guess for Vm (ideal gas law)\n",
        "    Vm_guess = R*T/P\n",
        "    Vm = fsolve(PR_eqn, Vm_guess, args=(P, T))\n",
        "    return Vm[0]\n",
        "\n",
        "# Integrand for the Gibbs-Duhem equation\n",
        "def integrand(P, T):\n",
        "    Vm = compute_Vm(P, T)\n",
        "    return Vm / (R*T)\n",
        "\n",
        "# Compute the fugacity ratio for given pressures P1 and P2 at temperature T\n",
        "def compute_fugacity_ratio(P1, P2, T):\n",
        "    integral, _ = quad(integrand, P1, P2, args=(T))\n",
        "    return np.exp(integral)\n",
        "\n",
        "# Example usage\n",
        "P1 = 1e5  # Initial pressure [Pa]\n",
        "P2 = 2e5  # Final pressure [Pa]\n",
        "T = 300  # Temperature [K]\n",
        "\n",
        "fugacity_ratio = compute_fugacity_ratio(P1, P2, T)\n",
        "print(f\"Fugacity ratio f2/f1 at T = {T} K: {fugacity_ratio}\")\n",
        "\n",
        "```\n",
        "\n",
        "**Documentation**:\n",
        "\n",
        "- We define the Peng-Robinson constants for a hypothetical gas. These values can be replaced with actual values for a specific gas.\n",
        "- The function `alpha(T)` computes the temperature-dependent function $\\alpha$ using the acentric factor $( \\omega )$.\n",
        "- The function `PR_eqn(Vm, P, T)` represents the Peng-Robinson equation of state.\n",
        "- The function `compute_Vm(P, T)` computes the molar volume $( V_m )$ for a given pressure $( P )$ and temperature $( T )$ using the `fsolve` function from `scipy.optimize`.\n",
        "- The function `integrand(P, T)` represents the integrand for the Gibbs-Duhem equation.\n",
        "- The function `compute_fugacity_ratio(P1, P2, T)` computes the fugacity ratio $( \\frac{f_2}{f_1} )$ for given pressures $( P1 )$ and $( P2 )$ at temperature $( T )$ using the `quad` function from `scipy.integrate`.\n",
        "- Finally, we provide an example usage for the code, computing the fugacity ratio for given pressures and temperature.\n",
        "\n",
        "This code provides a numerical approach to compute the fugacity ratio using the Peng-Robinson equation of state. You can adjust the constants and parameters as needed for a specific gas and conditions."
      ],
      "metadata": {
        "id": "6wlMNFmegI5F"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### Set up a test case considering the gases in the atmosphere and how fugacity is affected by, say, a 2 degree C climate change, with concordant increases in rainfall and humidity.\n",
        "\n"
      ],
      "metadata": {
        "id": "HqCCzxTEiSaI"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "To set up a test case considering the gases in the atmosphere and the effects of a 2°C climate change, we'll need to make several assumptions and simplifications. Here's a plan:\n",
        "\n",
        "1. **Gases to Consider**: The primary greenhouse gases in the atmosphere are water vapor (H₂O), carbon dioxide (CO₂), methane (CH₄), nitrous oxide (N₂O), and ozone (O₃). For this test case, we'll focus on water vapor and CO₂, as they are the most significant contributors to the greenhouse effect.\n",
        "\n",
        "2. **Initial Conditions**:\n",
        "   - Current average global temperature: $ T_1 $\n",
        "   - Current average global CO₂ concentration: $ P_{CO2,1} $\n",
        "   - Current average global water vapor concentration (expressed as relative humidity): $ RH_1 $\n",
        "\n",
        "3. **Changed Conditions**:\n",
        "   - Increased temperature: $ T_2 = T_1 + 2 $°C\n",
        "   - Increased CO₂ concentration due to feedback mechanisms (assumption): $ P_{CO2,2} $\n",
        "   - Increased relative humidity due to increased evaporation and rainfall (assumption): $ RH_2 $\n",
        "\n",
        "4. **Calculations**:\n",
        "   - Compute the fugacity of water vapor and CO₂ at $ T_1 $ and $ T_2 $ using the Peng-Robinson equation of state.\n",
        "   - Compare the fugacities to understand the change.\n",
        "\n",
        "Let's set up the test case:\n",
        "\n",
        "1. **Initial Conditions**:\n",
        "   - $ T_1 = 288 $ K (15°C, average global temperature)\n",
        "   - $ P_{CO2,1} = 400 $ ppm (approximate current concentration)\n",
        "   - $ RH_1 = 50\\% $ (assumed average global relative humidity)\n",
        "\n",
        "2. **Changed Conditions**:\n",
        "   - $ T_2 = 290 $ K (17°C)\n",
        "   - $ P_{CO2,2} = 410 $ ppm (assumed increase)\n",
        "   - $ RH_2 = 52\\% $ (assumed increase due to increased evaporation and rainfall)\n",
        "\n",
        "3. **Peng-Robinson Constants**:\n",
        "   - For water vapor and CO₂, we'll need the critical temperature, critical pressure, and acentric factor. These values can be found in literature or thermodynamic databases.\n",
        "\n",
        "4. **Calculations**:\n",
        "   - Use the provided Python code to compute the fugacity of water vapor and CO₂ at $ T_1 $ and $ T_2 $.\n",
        "\n",
        "Given the complexity of the task and the need for specific constants, we'll start by gathering the Peng-Robinson constants for water vapor and CO₂. Let's retrieve these constants.\n",
        "\n",
        "It seems that Wolfram|Alpha couldn't provide the Peng-Robinson constants directly. However, I can provide typical values for these constants from standard literature:\n",
        "\n",
        "For **Water Vapor (H₂O)**:\n",
        "- $ T_c $ (Critical Temperature) = 647.1 K\n",
        "- $ P_c $ (Critical Pressure) = 22.064 MPa\n",
        "- $ \\omega $ (Acentric Factor) = 0.344\n",
        "\n",
        "For **Carbon Dioxide (CO₂)**:\n",
        "- $ T_c $ = 304.2 K\n",
        "- $ P_c $ = 7.38 MPa\n",
        "- $ \\omega $ = 0.225\n",
        "\n",
        "With these constants, we can proceed to compute the fugacity of water vapor and CO₂ at the initial and changed conditions using the Python code provided earlier. However, it's worth noting that the Peng-Robinson equation of state might not be the most accurate for all conditions, especially for water vapor near its critical point. More specialized equations of state or activity coefficient models might be needed for high accuracy."
      ],
      "metadata": {
        "id": "t-lwG12wik6f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "Computing the fugacity of water vapor and CO₂ at the initial conditions (\\( T_1 \\)) and the changed conditions (\\( T_2 \\)).\n",
        "\n",
        "Here's the Python code for the test case:\n",
        "\n",
        "```python\n",
        "import numpy as np\n",
        "from scipy.optimize import fsolve\n",
        "\n",
        "# Universal gas constant [J/(mol*K)]\n",
        "R = 8.314\n",
        "\n",
        "# Peng-Robinson constants for Water Vapor (H2O)\n",
        "a_H2O = 0.45724 * (R**2 * 647.1**2) / 22.064e6\n",
        "b_H2O = 0.07780 * R * 647.1 / 22.064e6\n",
        "omega_H2O = 0.344\n",
        "\n",
        "# Peng-Robinson constants for CO2\n",
        "a_CO2 = 0.45724 * (R**2 * 304.2**2) / 7.38e6\n",
        "b_CO2 = 0.07780 * R * 304.2 / 7.38e6\n",
        "omega_CO2 = 0.225\n",
        "\n",
        "# Temperature-dependent function alpha\n",
        "def alpha(T, Tc, omega):\n",
        "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
        "    return (1 + kappa*(1 - np.sqrt(T/Tc)))**2\n",
        "\n",
        "# Peng-Robinson equation of state\n",
        "def PR_eqn(Vm, P, T, a, b, Tc, omega):\n",
        "    a_alpha = a * alpha(T, Tc, omega)\n",
        "    return P - R*T/(Vm - b) + a_alpha/((Vm**2) + 2*b*Vm - b**2)\n",
        "\n",
        "# Compute molar volume Vm for a given P, T, and gas\n",
        "def compute_Vm(P, T, a, b, Tc, omega):\n",
        "    Vm_guess = R*T/P\n",
        "    Vm = fsolve(PR_eqn, Vm_guess, args=(P, T, a, b, Tc, omega))\n",
        "    return Vm[0]\n",
        "\n",
        "# Compute the fugacity for given P, T, and gas\n",
        "def compute_fugacity(P, T, a, b, Tc, omega):\n",
        "    Vm = compute_Vm(P, T, a, b, Tc, omega)\n",
        "    Z = P*Vm/(R*T)\n",
        "    ln_phi = Z - 1 - np.log(Z) - (a/(2*np.sqrt(2)*b*R*T)) * np.log((Vm + (1+np.sqrt(2))*b)/(Vm + (1-np.sqrt(2))*b))\n",
        "    return np.exp(ln_phi) * P\n",
        "\n",
        "# Initial and changed conditions\n",
        "T1 = 288  # Initial temperature [K]\n",
        "T2 = 290  # Changed temperature [K]\n",
        "P_CO2_1 = 400e-6 * 1e5  # Initial CO2 pressure [Pa]\n",
        "P_CO2_2 = 410e-6 * 1e5  # Changed CO2 pressure [Pa]\n",
        "\n",
        "# Compute fugacities\n",
        "fugacity_H2O_1 = compute_fugacity(0.5*1e5, T1, a_H2O, b_H2O, 647.1, omega_H2O)  # Assuming 50% relative humidity\n",
        "fugacity_H2O_2 = compute_fugacity(0.52*1e5, T2, a_H2O, b_H2O, 647.1, omega_H2O)  # Assuming 52% relative humidity\n",
        "fugacity_CO2_1 = compute_fugacity(P_CO2_1, T1, a_CO2, b_CO2, 304.2, omega_CO2)\n",
        "fugacity_CO2_2 = compute_fugacity(P_CO2_2, T2, a_CO2, b_CO2, 304.2, omega_CO2)\n",
        "\n",
        "print(f\"Fugacity of H2O at T1: {fugacity_H2O_1} Pa\")\n",
        "print(f\"Fugacity of H2O at T2: {fugacity_H2O_2} Pa\")\n",
        "print(f\"Fugacity of CO2 at T1: {fugacity_CO2_1} Pa\")\n",
        "print(f\"Fugacity of CO2 at T2: {fugacity_CO2_2} Pa\")\n",
        "```\n",
        "\n",
        "This code will compute the fugacities of water vapor and CO₂ at the initial and changed conditions using the Peng-Robinson equation of state. The results will give an indication of how the \"escaping tendencies\" of these gases change with the 2°C temperature increase and the assumed changes in their concentrations."
      ],
      "metadata": {
        "id": "vlx5XcdmkDD_"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "from scipy.optimize import fsolve\n",
        "\n",
        "# Universal gas constant [J/(mol*K)]\n",
        "R = 8.314\n",
        "\n",
        "# Peng-Robinson constants for Water Vapor (H2O)\n",
        "a_H2O = 0.45724 * (R**2 * 647.1**2) / 22.064e6\n",
        "b_H2O = 0.07780 * R * 647.1 / 22.064e6\n",
        "omega_H2O = 0.344\n",
        "\n",
        "# Peng-Robinson constants for CO2\n",
        "a_CO2 = 0.45724 * (R**2 * 304.2**2) / 7.38e6\n",
        "b_CO2 = 0.07780 * R * 304.2 / 7.38e6\n",
        "omega_CO2 = 0.225\n",
        "\n",
        "# Temperature-dependent function alpha\n",
        "def alpha(T, Tc, omega):\n",
        "    kappa = 0.37464 + 1.54226*omega - 0.26992*omega**2\n",
        "    return (1 + kappa*(1 - np.sqrt(T/Tc)))**2\n",
        "\n",
        "# Peng-Robinson equation of state\n",
        "def PR_eqn(Vm, P, T, a, b, Tc, omega):\n",
        "    a_alpha = a * alpha(T, Tc, omega)\n",
        "    return P - R*T/(Vm - b) + a_alpha/((Vm**2) + 2*b*Vm - b**2)\n",
        "\n",
        "# Compute molar volume Vm for a given P, T, and gas\n",
        "def compute_Vm(P, T, a, b, Tc, omega):\n",
        "    Vm_guess = R*T/P\n",
        "    Vm = fsolve(PR_eqn, Vm_guess, args=(P, T, a, b, Tc, omega))\n",
        "    return Vm[0]\n",
        "\n",
        "# Compute the fugacity for given P, T, and gas\n",
        "def compute_fugacity(P, T, a, b, Tc, omega):\n",
        "    Vm = compute_Vm(P, T, a, b, Tc, omega)\n",
        "    Z = P*Vm/(R*T)\n",
        "    ln_phi = Z - 1 - np.log(Z) - (a/(2*np.sqrt(2)*b*R*T)) * np.log((Vm + (1+np.sqrt(2))*b)/(Vm + (1-np.sqrt(2))*b))\n",
        "    return np.exp(ln_phi) * P\n",
        "\n",
        "# Initial and changed conditions\n",
        "T1 = 288  # Initial temperature [K]\n",
        "T2 = 290  # Changed temperature [K]\n",
        "P_CO2_1 = 400e-6 * 1e5  # Initial CO2 pressure [Pa]\n",
        "P_CO2_2 = 410e-6 * 1e5  # Changed CO2 pressure [Pa]\n",
        "\n",
        "# Compute fugacities\n",
        "fugacity_H2O_1 = compute_fugacity(0.5*1e5, T1, a_H2O, b_H2O, 647.1, omega_H2O)  # Assuming 50% relative humidity\n",
        "fugacity_H2O_2 = compute_fugacity(0.52*1e5, T2, a_H2O, b_H2O, 647.1, omega_H2O)  # Assuming 52% relative humidity\n",
        "fugacity_CO2_1 = compute_fugacity(P_CO2_1, T1, a_CO2, b_CO2, 304.2, omega_CO2)\n",
        "fugacity_CO2_2 = compute_fugacity(P_CO2_2, T2, a_CO2, b_CO2, 304.2, omega_CO2)\n",
        "\n",
        "print(f\"Fugacity of H2O at T1: {fugacity_H2O_1} Pa\")\n",
        "print(f\"Fugacity of H2O at T2: {fugacity_H2O_2} Pa\")\n",
        "print(f\"Fugacity of CO2 at T1: {fugacity_CO2_1} Pa\")\n",
        "print(f\"Fugacity of CO2 at T2: {fugacity_CO2_2} Pa\")\n"
      ],
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "2e9o_7XdQC0a",
        "outputId": "d4fb7d98-70c7-4ef3-c3f6-a76c57122b12"
      },
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Fugacity of H2O at T1: 49738.79343980095 Pa\n",
            "Fugacity of H2O at T2: 51721.36756280549 Pa\n",
            "Fugacity of CO2 at T1: 39.99988940357763 Pa\n",
            "Fugacity of CO2 at T2: 40.999885401803105 Pa\n"
          ]
        }
      ]
    }
  ]
}