pyrestoolbox - pypi Package Compare versions

Comparing version 2.2.3 to 2.5.0

pyrestoolbox/brine/_lib_salting_library.py

+"""
+Library of Sechenov (salting-out) correlations for gases in NaCl brine.
+
+All functions return ks on a LOG10 basis with molality (mol/kg H2O) scale:
+
+    log10(S_water / S_brine) = ks * m_NaCl
+
+Sources implemented:
+    1. S&W Eq 8    — Soreide & Whitson 1992, generic (Tb-based), all gases
+    2. Dubessy     — Dubessy et al. 2005, extended Sechenov for CO2 and H2S
+    3. Akinfiev    — Akinfiev et al. 2016, Pitzer model for H2S (-> effective ks)
+    4. Li et al.   — Li, Zhang, Luo & Li 2015, Pitzer for CH4, C2H6, C3H8, nC4H10
+    5. Mao & Duan  — Mao & Duan 2006, Pitzer model for N2
+    6. Duan & Sun  — Duan & Sun 2003, Pitzer model for CO2
+
+Not implemented (too complex for Sechenov extraction):
+    - Springer et al. 2014 (OTC-25295-MS): MSE framework for H2S/CO2.
+      Requires speciation + Debye-Huckel + virial + UNIQUAC simultaneously.
+      Cannot be reduced to an analytical Sechenov form.
+
+Convention:  ks > 0 means salting-out (lower solubility in brine).
+"""
+
+import numpy as np
+
+LN10 = np.log(10.0)  # 2.302585...
+
+
+# ====================================================================
+# 1. Soreide & Whitson 1992, Equation 8  (generic, Tb-based)
+# ====================================================================
+# Boiling points (K) for S&W gases
+TB_K = {
+    'H2': 20.3, 'N2': 77.4, 'CO2': 194.7, 'H2S': 213.6,
+    'CH4': 111.6, 'C2H6': 184.6, 'C3H8': 231.1,
+    'iC4H10': 261.4, 'nC4H10': 272.7, 'iC5H12': 301.0,
+    'nC5H12': 309.2, 'nC6H14': 341.9, 'nC7H16': 371.6,
+    'nC8H18': 398.8, 'nC10H22': 447.3,
+}
+
+
+def ks_sw_eq8(T_K, gas_or_Tb):
+    """Soreide-Whitson 1992 Equation 8 Sechenov coefficient.
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin.
+    gas_or_Tb : str or float
+        Gas name (e.g. 'H2', 'CH4') or boiling point in K.
+
+    Returns
+    -------
+    ks : float or array
+        Sechenov coefficient (log10 basis, molality scale).
+    """
+    if isinstance(gas_or_Tb, str):
+        Tb = TB_K[gas_or_Tb]
+    else:
+        Tb = float(gas_or_Tb)
+    T_F = (np.asarray(T_K, dtype=float) - 273.15) * 9.0 / 5.0 + 32.0
+    return (0.13163 + 4.45e-4 * Tb - 7.692e-4 * T_F
+            + 2.6614e-6 * T_F**2 - 2.612e-9 * T_F**3)
+
+
+# ====================================================================
+# 2. Dubessy et al. (2005) — Extended Sechenov for CO2 and H2S
+#    Oil & Gas Sci. Technol. - Rev. IFP, Vol. 60, No. 2, pp. 339-355
+#    Table 9 coefficients.
+#
+#    log10(K_brine/K_water) = m*b1(T) + m^2*b2(T) + m^3*b3(T)
+#    b_i(T) = sum(B_ij * T^j, j=0..4)
+#
+#    Effective ks = b1(T) + m*b2(T) + m^2*b3(T)
+# ====================================================================
+
+# CO2: column scaling 1e0, 1e-2, 1e-4, 1e-8, 1e-11
+_DUB_CO2_B1 = np.array([3.114712456, -2.7655585e-2, 0.9176713976e-4,
+                         -12.78795941e-8, 6.2704268351e-11])
+_DUB_CO2_B2 = np.array([-2.05637458, 2.081980200e-2, -0.765857702e-4,
+                         12.011325315e-8, -6.790343083e-11])
+_DUB_CO2_B3 = np.array([0.253424331, -0.26047432e-2, 0.0972580216e-4,
+                         -1.551654794e-8, 0.8948557284e-11])
+
+# H2S: column scaling 1e0, 1e-2, 1e-4, 1e-7, 1e-10
+_DUB_H2S_B1 = np.array([-12.4617636, 12.69373100e-2, -4.791540697e-4,
+                         7.9817223650e-7, -4.931093145e-10])
+_DUB_H2S_B2 = np.array([5.327383011, -5.82779828e-2, 2.3650333285e-4,
+                         -4.207913036e-7, 2.7628521914e-10])
+_DUB_H2S_B3 = np.array([-0.75715275, 0.831927411e-2, -0.338668040e-4,
+                         0.6037602785e-7, -0.397049836e-10])
+
+
+def _poly4(c, T):
+    """Evaluate 4th-order polynomial c[0] + c[1]*T + ... + c[4]*T^4."""
+    T = np.asarray(T, dtype=float)
+    return c[0] + c[1]*T + c[2]*T**2 + c[3]*T**3 + c[4]*T**4
+
+
+def ks_dubessy_co2(T_K, m_NaCl=0.0):
+    """Dubessy et al. 2005 effective Sechenov for CO2-NaCl (log10).
+
+    Valid: T <= 543 K (270 C), P <= 300 bar, m <= 6. R^2 = 0.83.
+    """
+    return (_poly4(_DUB_CO2_B1, T_K) + m_NaCl * _poly4(_DUB_CO2_B2, T_K)
+            + m_NaCl**2 * _poly4(_DUB_CO2_B3, T_K))
+
+
+def ks_dubessy_h2s(T_K, m_NaCl=0.0):
+    """Dubessy et al. 2005 effective Sechenov for H2S-NaCl (log10).
+
+    Valid: T <= 523 K (250 C), P <= 150 bar, m <= 6. R^2 = 0.74.
+    """
+    return (_poly4(_DUB_H2S_B1, T_K) + m_NaCl * _poly4(_DUB_H2S_B2, T_K)
+            + m_NaCl**2 * _poly4(_DUB_H2S_B3, T_K))
+
+
+# ====================================================================
+# 3. Akinfiev, Majer & Shvarov (2016) — Pitzer model for H2S-NaCl
+#    Chemical Geology, 424, 1-11. DOI: 10.1016/j.chemgeo.2016.01.006
+#
+#    Pitzer activity coefficient (Eq. 16):
+#      ln(gamma_s) = 2*m_s*lam_ss + 3*m_s^2*tau_sss
+#                  + 2*m_e*B_se + 3*m_e^2*C_see + 6*m_s*m_e*C_sse
+#
+#    For dilute H2S (m_s << m_e), the effective Sechenov is:
+#      ks_eff ≈ (2*B_se + 6*m_s*C_sse) / ln(10)   [on log10/molality]
+#
+#    But the rigorous approach uses the full recommended solubility tables
+#    or the complete Pitzer framework.
+# ====================================================================
+
+# Binary H2S-H2O self-interaction: Eq. (18)
+# lambda_ss = a_lam + b_lam * (100/(T-228)) + c_lam * (T/(T-760))
+_AK_A_LAM = -0.19515
+_AK_B_LAM = 0.102822
+_AK_C_LAM = -0.033726
+_AK_TAU_SSS = 0.004900
+
+# Ternary H2S-H2O-NaCl: Eq. (20)
+# B_se = b_B * (100/(T-228)) + c_B * (T/(T-760))
+_AK_B_B = 0.03568
+_AK_C_B = -0.02354
+_AK_C_SEE = 0.0       # set to zero by authors
+_AK_C_SSE = 0.002558
+
+
+def _akinfiev_lambda_ss(T_K):
+    """H2S-H2S self-interaction parameter (Akinfiev Eq. 18)."""
+    T = np.asarray(T_K, dtype=float)
+    return _AK_A_LAM + _AK_B_LAM * (100.0 / (T - 228.0)) + _AK_C_LAM * (T / (T - 760.0))
+
+
+def _akinfiev_B_se(T_K):
+    """H2S-NaCl interaction parameter B_se (Akinfiev Eq. 20)."""
+    T = np.asarray(T_K, dtype=float)
+    return _AK_B_B * (100.0 / (T - 228.0)) + _AK_C_B * (T / (T - 760.0))
+
+
+def akinfiev_ln_gamma_h2s(T_K, m_h2s, m_NaCl):
+    """Natural log of H2S activity coefficient in H2S-H2O-NaCl (Eq. 16).
+
+    Parameters
+    ----------
+    T_K : float
+        Temperature in Kelvin (283-573 K).
+    m_h2s : float
+        H2S molality (mol/kg H2O).
+    m_NaCl : float
+        NaCl molality (mol/kg H2O), 0-6.
+
+    Returns
+    -------
+    ln_gamma : float
+        Natural log of activity coefficient on Henry scale.
+    """
+    lam = _akinfiev_lambda_ss(T_K)
+    B = _akinfiev_B_se(T_K)
+    return (2.0 * m_h2s * lam + 3.0 * m_h2s**2 * _AK_TAU_SSS
+            + 2.0 * m_NaCl * B + 3.0 * m_NaCl**2 * _AK_C_SEE
+            + 6.0 * m_h2s * m_NaCl * _AK_C_SSE)
+
+
+def ks_akinfiev_h2s(T_K, m_NaCl=1.0, m_h2s_approx=0.1):
+    """Effective Sechenov coefficient for H2S-NaCl from Akinfiev Pitzer (log10).
+
+    Computes the effective ks from the Pitzer activity coefficient difference
+    between saline and fresh solutions, at a specified approximate H2S molality.
+
+    For dilute H2S, ln(gamma_saline/gamma_fresh) ≈ 2*m_e*B_se + 6*m_s*m_e*C_sse
+    so ks ≈ (2*B_se + 6*m_s*C_sse) / ln(10).
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin (283-573 K).
+    m_NaCl : float
+        NaCl molality for computing the effective ks. Default 1 m.
+    m_h2s_approx : float
+        Approximate H2S molality (affects C_sse contribution). Default 0.1 m.
+
+    Returns
+    -------
+    ks : float or array
+        Effective Sechenov coefficient (log10 basis, molality scale).
+    """
+    B = _akinfiev_B_se(T_K)
+    # ln(gamma_brine) - ln(gamma_fresh) at given m_h2s:
+    #   = 2*m_NaCl*B_se + 3*m_NaCl^2*C_see + 6*m_h2s*m_NaCl*C_sse
+    # Dividing by m_NaCl gives the ln-based Sechenov:
+    ks_ln = 2.0 * B + 3.0 * m_NaCl * _AK_C_SEE + 6.0 * m_h2s_approx * _AK_C_SSE
+    return ks_ln / LN10
+
+
+def ks_akinfiev_h2s_from_tables(T_K, m_NaCl):
+    """Effective ks for H2S from Akinfiev 2016 recommended solubility tables.
+
+    Uses linearly interpolated Table 3 (pure water) and Table 4 (brine)
+    values at 5 MPa to compute ks = log10(m_water/m_brine) / m_NaCl.
+
+    This is the most reliable method — uses the model output directly,
+    avoids any approximation in the Pitzer-to-Sechenov conversion.
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin. Supported: 298-523 K.
+    m_NaCl : float
+        Must be one of: 1, 2, 4, or 6 mol/kg.
+
+    Returns
+    -------
+    ks : float or array
+        Effective Sechenov coefficient (log10 basis).
+    """
+    # Tables 3 & 4 at P = 5 MPa (mid-range, well within experimental coverage)
+    # T_K:                298.15  323.15  373.15  423.15  473.15  523.15
+    _T = np.array([298.15, 323.15, 373.15, 423.15, 473.15, 523.15])
+    _m0 = np.array([np.nan, 2.026, 1.651, 1.271, 0.956, 0.300])  # pure water, 5 MPa
+    # NaCl = 1 m
+    _m1 = np.array([np.nan, 1.752, 1.437, 1.110, 0.831, 0.260])
+    # NaCl = 2 m
+    _m2 = np.array([np.nan, 1.523, 1.259, 0.975, 0.726, 0.226])
+    # NaCl = 4 m
+    _m4 = np.array([np.nan, 1.168, 0.983, 0.763, 0.561, 0.171])
+    # NaCl = 6 m
+    _m6 = np.array([np.nan, 0.907, 0.779, 0.606, 0.438, 0.130])
+
+    brine_tables = {1: _m1, 2: _m2, 4: _m4, 6: _m6}
+    if m_NaCl not in brine_tables:
+        raise ValueError(f"m_NaCl must be 1, 2, 4, or 6. Got {m_NaCl}")
+
+    m_brine = brine_tables[m_NaCl]
+    # Compute ks at each table temperature (skip 298.15 where data missing at 5 MPa)
+    valid = ~np.isnan(_m0) & ~np.isnan(m_brine)
+    ks_pts = np.full_like(_T, np.nan)
+    ks_pts[valid] = np.log10(_m0[valid] / m_brine[valid]) / m_NaCl
+
+    # Interpolate to requested temperatures
+    T_valid = _T[valid]
+    ks_valid = ks_pts[valid]
+    T_K = np.asarray(T_K, dtype=float)
+    scalar = T_K.ndim == 0
+    T_K = np.atleast_1d(T_K)
+    result = np.interp(T_K, T_valid, ks_valid)
+    return float(result[0]) if scalar else result
+
+
+# ====================================================================
+# 4. Li, Zhang, Luo & Li (2015) — Pitzer model for CH4-C2H6-C3H8-nC4H10
+#    Applied Geochemistry. DOI: 10.1016/j.apgeochem.2015.02.006
+#
+#    Eq. 15 (for NaCl, with lambda_{i-Cl-} = 0):
+#      ln(gamma_i) = 2 * m * lambda_{i-Na+}(T,P) + m^2 * zeta_{i-Na+-Cl-}
+#
+#    Effective Sechenov (log10):
+#      ks = (2*lambda + zeta*m) / ln(10)
+#
+#    Validity: 273-473 K, 1-1000 bar, 0-6 m NaCl
+#    NOTE: C3H8 and nC4H10 have very limited brine data (1 atm only).
+# ====================================================================
+
+def _li2015_lambda_ch4(T_K, P_bar=100.0):
+    """CH4-Na+ Pitzer lambda (Li et al. 2015 Table 7). T in K, P in bar."""
+    T = np.asarray(T_K, dtype=float)
+    P = float(P_bar)
+    return (-5.7066455e-1 + 7.2997588e-4 * T + 1.5176903e2 / T
+            + 3.1927112e-5 * P - 1.642651e-5 * P / T)
+
+
+def _li2015_lambda_c2h6(T_K, P_bar=100.0):
+    """C2H6-Na+ Pitzer lambda (Li et al. 2015 Table 7). T in K, P in bar."""
+    T = np.asarray(T_K, dtype=float)
+    P = float(P_bar)
+    return (-2.143686 + 2.598765e-3 * T + 4.6942351e2 / T
+            - 4.6849541e-5 * P - 8.4616602e-10 * P**2 * T
+            + 1.095219e-6 * P * T)
+
+
+def _li2015_lambda_c3h8(T_K, P_bar=100.0):
+    """C3H8-Na+ Pitzer lambda (Li et al. 2015 Table 7). T-only (no P dep)."""
+    T = np.asarray(T_K, dtype=float)
+    return 0.513068 - 0.000958 * T
+
+
+def _li2015_lambda_nc4h10(T_K, P_bar=100.0):
+    """nC4H10-Na+ Pitzer lambda (Li et al. 2015 Table 7). T-only (no P dep)."""
+    T = np.asarray(T_K, dtype=float)
+    return 0.52862384 - 1.0298104e-3 * T
+
+
+# Zeta (third virial) — all constants, no T or P dependence
+_LI2015_ZETA = {
+    'CH4':    -2.9990084e-3,
+    'C2H6':   -1.0165947e-2,
+    'C3H8':   -0.007485,
+    'NC4H10':  0.0206946,    # NOTE: positive (unlike the others)
+}
+
+_LI2015_LAMBDA = {
+    'CH4':    _li2015_lambda_ch4,
+    'C2H6':   _li2015_lambda_c2h6,
+    'C3H8':   _li2015_lambda_c3h8,
+    'NC4H10': _li2015_lambda_nc4h10,
+}
+
+
+def ks_li2015(T_K, gas, m_NaCl=0.0, P_bar=100.0):
+    """Effective Sechenov for light HCs from Li et al. 2015 Pitzer (log10).
+
+    ln(gamma) = 2*m*lambda(T,P) + m^2*zeta
+    ks = (2*lambda + zeta*m) / ln(10)
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin (273-473 K).
+    gas : str
+        'CH4', 'C2H6', 'C3H8', or 'nC4H10'.
+    m_NaCl : float
+        NaCl molality. Affects effective ks through zeta term.
+    P_bar : float
+        Pressure in bar (affects CH4 and C2H6 lambda). Default 100.
+
+    Returns
+    -------
+    ks : float or array
+        Effective Sechenov coefficient (log10 basis, molality scale).
+    """
+    g = gas.upper().replace('N', 'N')  # normalize
+    if g == 'NC4H10' or g == 'NC4':
+        g = 'NC4H10'
+    elif g not in _LI2015_LAMBDA:
+        raise ValueError(f"Li 2015 covers CH4, C2H6, C3H8, nC4H10. Got {gas}")
+
+    lam = _LI2015_LAMBDA[g](T_K, P_bar)
+    zeta = _LI2015_ZETA[g]
+    return (2.0 * lam + zeta * m_NaCl) / LN10
+
+
+# ====================================================================
+# 5. Mao & Duan (2006) — Pitzer model for N2-NaCl
+#    Fluid Phase Equilibria, 248, 103-114.
+#    DOI: 10.1016/j.fluid.2006.07.020
+#
+#    Eq. 8 (for NaCl, with lambda_{N2-Cl-} = 0):
+#      ln(gamma_N2) = 2 * m * lambda_{N2-Na+}(T,P) + m^2 * xi_{N2-Na+-Cl-}
+#
+#    Effective Sechenov (log10):
+#      ks = (2*lambda + xi*m) / ln(10)
+#
+#    Validity: 273-400 K, 1-600 bar, 0-6 m NaCl
+# ====================================================================
+
+def _mao2006_lambda_n2(T_K, P_bar=100.0):
+    """N2-Na+ Pitzer lambda (Mao & Duan 2006 Table 3). T in K, P in bar."""
+    T = np.asarray(T_K, dtype=float)
+    P = float(P_bar)
+    return (-2.4434074 + 0.0036351795 * T + 447.47364 / T
+            - 0.000013711527 * P + 0.0000071037217 * P**2 / T)
+
+
+_MAO2006_XI_N2 = -0.0058071053  # constant, no T or P dependence
+
+
+def ks_mao2006_n2(T_K, m_NaCl=0.0, P_bar=100.0):
+    """Effective Sechenov for N2 from Mao & Duan 2006 Pitzer model (log10).
+
+    ln(gamma_N2) = 2*m*lambda(T,P) + m^2*xi
+    ks = (2*lambda + xi*m) / ln(10)
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin (273-400 K).
+    m_NaCl : float
+        NaCl molality. Affects effective ks through xi term.
+    P_bar : float
+        Pressure in bar (1-600 bar). Default 100.
+
+    Returns
+    -------
+    ks : float or array
+        Effective Sechenov coefficient (log10 basis, molality scale).
+    """
+    lam = _mao2006_lambda_n2(T_K, P_bar)
+    return (2.0 * lam + _MAO2006_XI_N2 * m_NaCl) / LN10
+
+
+# ====================================================================
+# 6. Duan & Sun (2003) — Pitzer model for CO2-NaCl
+#    Chemical Geology, 193, 257-271.
+#    DOI: 10.1016/S0009-2541(02)00263-2
+#
+#    Same Pitzer framework as Mao 2006 (N2) — Pitzer et al. (1984) form:
+#      Par(T,P) = c1 + c2*T + c3/T + c4*T^2 + c5/(630-T)
+#               + c6*P + c7*P*ln(T) + c8*P/T + c9*P/(630-T)
+#               + c10*P^2/(630-T)^2 + c11*T*ln(P)
+#
+#    lambda_{CO2-Na+} uses c1,c2,c3,c8,c9 (5 coeffs)
+#    zeta_{CO2-Na+-Cl-} uses c1,c2,c8,c9 (4 coeffs)
+#    lambda_{CO2-Cl-} = 0 (set to zero)
+#
+#    ln(gamma_CO2) = 2*m*lambda + m^2*zeta
+#    Effective Sechenov: ks = (2*lambda + zeta*m) / ln(10)
+#
+#    Validity: 273-533 K, 0-2000 bar, 0-4.3 m NaCl
+# ====================================================================
+
+def _duan2003_lambda_co2(T_K, P_bar=100.0):
+    """CO2-Na+ Pitzer lambda (Duan & Sun 2003 Table 2). T in K, P in bar."""
+    T = np.asarray(T_K, dtype=float)
+    P = float(P_bar)
+    return (-0.411370585 + 6.07632013e-4 * T + 97.5347708 / T
+            - 0.0237622469 * P / T + 0.0170656236 * P / (630.0 - T))
+
+
+def _duan2003_zeta_co2(T_K, P_bar=100.0):
+    """CO2-Na+-Cl- Pitzer zeta (Duan & Sun 2003 Table 2). T in K, P in bar."""
+    T = np.asarray(T_K, dtype=float)
+    P = float(P_bar)
+    return (3.36389723e-4 - 1.98298980e-5 * T
+            + 2.12220830e-3 * P / T - 5.24873303e-3 * P / (630.0 - T))
+
+
+def ks_duan2003_co2(T_K, m_NaCl=0.0, P_bar=100.0):
+    """Effective Sechenov for CO2 from Duan & Sun 2003 Pitzer model (log10).
+
+    ln(gamma_CO2) = 2*m*lambda(T,P) + m^2*zeta(T,P)
+    ks = (2*lambda + zeta*m) / ln(10)
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin (273-533 K).
+    m_NaCl : float
+        NaCl molality (0-4.3 m). Affects effective ks through zeta term.
+    P_bar : float
+        Pressure in bar (0-2000 bar). Default 100.
+
+    Returns
+    -------
+    ks : float or array
+        Effective Sechenov coefficient (log10 basis, molality scale).
+    """
+    lam = _duan2003_lambda_co2(T_K, P_bar)
+    zeta = _duan2003_zeta_co2(T_K, P_bar)
+    return (2.0 * lam + zeta * m_NaCl) / LN10
+
+
+# ====================================================================
+# Convenience: unified interface
+# ====================================================================
+
+def ks_library(T_K, gas, source='sw_eq8', m_NaCl=0.0, **kwargs):
+    """Unified Sechenov coefficient lookup (log10 basis, molality scale).
+
+    Parameters
+    ----------
+    T_K : float or array
+        Temperature in Kelvin.
+    gas : str
+        Gas name: 'H2', 'N2', 'CH4', 'C2H6', 'C3H8', 'CO2', 'H2S', etc.
+    source : str
+        'sw_eq8'        — Soreide & Whitson 1992 Eq 8 (all gases)
+        'dubessy'       — Dubessy et al. 2005 (CO2, H2S only)
+        'akinfiev'      — Akinfiev et al. 2016 Pitzer (H2S only)
+        'akinfiev_table'— Akinfiev et al. 2016 from tables (H2S only)
+        'li2015'        — Li et al. 2015 Pitzer (CH4, C2H6, C3H8, nC4H10)
+        'mao2006'       — Mao & Duan 2006 Pitzer (N2 only)
+        'duan2003'      — Duan & Sun 2003 Pitzer (CO2 only)
+    m_NaCl : float
+        NaCl molality (needed for non-linear Sechenov models).
+
+    Returns
+    -------
+    ks : float or array
+    """
+    gas = gas.upper()
+    if gas == 'NC4H10':
+        pass  # already normalized
+    elif gas == 'NC4':
+        gas = 'NC4H10'
+    src = source.lower()
+
+    if src == 'sw_eq8':
+        return ks_sw_eq8(T_K, gas)
+    elif src == 'dubessy':
+        if gas == 'CO2':
+            return ks_dubessy_co2(T_K, m_NaCl)
+        elif gas == 'H2S':
+            return ks_dubessy_h2s(T_K, m_NaCl)
+        else:
+            raise ValueError(f"Dubessy only covers CO2 and H2S, not {gas}")
+    elif src == 'akinfiev':
+        if gas == 'H2S':
+            return ks_akinfiev_h2s(T_K, m_NaCl, **kwargs)
+        else:
+            raise ValueError(f"Akinfiev 2016 only covers H2S, not {gas}")
+    elif src == 'akinfiev_table':
+        if gas == 'H2S':
+            return ks_akinfiev_h2s_from_tables(T_K, m_NaCl)
+        else:
+            raise ValueError(f"Akinfiev tables only for H2S, not {gas}")
+    elif src == 'li2015':
+        return ks_li2015(T_K, gas, m_NaCl, **kwargs)
+    elif src == 'mao2006':
+        if gas == 'N2':
+            return ks_mao2006_n2(T_K, m_NaCl, **kwargs)
+        else:
+            raise ValueError(f"Mao & Duan 2006 only covers N2, not {gas}")
+    elif src == 'duan2003':
+        if gas == 'CO2':
+            return ks_duan2003_co2(T_K, m_NaCl, **kwargs)
+        else:
+            raise ValueError(f"Duan & Sun 2003 only covers CO2, not {gas}")
+    else:
+        raise ValueError(f"Unknown source: {source}")
+
+
+# ====================================================================
+# Main: comparison table
+# ====================================================================
+if __name__ == '__main__':
+    temps_C = [25, 50, 75, 100, 125, 150, 200, 250]
+    temps_K = [T + 273.15 for T in temps_C]
+
+    print("=" * 80)
+    print("Sechenov Coefficient Library — All values on log10 / molality basis")
+    print("  ks > 0 means salting-out;  log10(S_water/S_brine) = ks * m_NaCl")
+    print("=" * 80)
+
+    # --- CO2 comparison ---
+    print("\n--- CO2: S&W vs Dubessy vs Duan 2003 ---")
+    print(f"{'T(C)':>6}  {'S&W Eq8':>8}  {'Dub m=0':>8}  {'Dub m=1':>8}  "
+          f"{'Dub m=4':>8}  {'Duan m=0':>9}  {'Duan m=1':>9}  {'Duan m=4':>9}")
+    print("-" * 82)
+    for Tc, Tk in zip(temps_C, temps_K):
+        sw = ks_sw_eq8(Tk, 'CO2')
+        d0 = ks_dubessy_co2(Tk, 0)
+        d1 = ks_dubessy_co2(Tk, 1)
+        d4 = ks_dubessy_co2(Tk, 4)
+        if Tk <= 533:
+            dn0 = ks_duan2003_co2(Tk, 0, 100)
+            dn1 = ks_duan2003_co2(Tk, 1, 100)
+            dn4 = ks_duan2003_co2(Tk, 4, 100)
+            print(f"{Tc:>6}  {sw:8.4f}  {d0:8.4f}  {d1:8.4f}  {d4:8.4f}  "
+                  f"{dn0:9.4f}  {dn1:9.4f}  {dn4:9.4f}")
+        else:
+            print(f"{Tc:>6}  {sw:8.4f}  {d0:8.4f}  {d1:8.4f}  {d4:8.4f}  "
+                  f"{'n/a':>9}  {'n/a':>9}  {'n/a':>9}")
+
+    # --- H2S comparison ---
+    print("\n--- H2S ---")
+    print(f"{'T(C)':>6}  {'S&W Eq8':>8}  {'Dub m=0':>8}  {'Dub m=1':>8}  "
+          f"{'Akin m=1':>9}  {'Akin m=4':>9}  {'AkTbl m=1':>10}  {'AkTbl m=4':>10}")
+    print("-" * 90)
+    for Tc, Tk in zip(temps_C, temps_K):
+        sw = ks_sw_eq8(Tk, 'H2S')
+        d0 = ks_dubessy_h2s(Tk, 0)
+        d1 = ks_dubessy_h2s(Tk, 1)
+        ak1 = ks_akinfiev_h2s(Tk, m_NaCl=1.0, m_h2s_approx=0.1)
+        ak4 = ks_akinfiev_h2s(Tk, m_NaCl=4.0, m_h2s_approx=0.1)
+        # Table-based (only for T >= 323.15 K)
+        if Tk >= 323.15:
+            akt1 = ks_akinfiev_h2s_from_tables(Tk, 1)
+            akt4 = ks_akinfiev_h2s_from_tables(Tk, 4)
+            print(f"{Tc:>6}  {sw:8.4f}  {d0:8.4f}  {d1:8.4f}  "
+                  f"{ak1:9.4f}  {ak4:9.4f}  {akt1:10.4f}  {akt4:10.4f}")
+        else:
+            print(f"{Tc:>6}  {sw:8.4f}  {d0:8.4f}  {d1:8.4f}  "
+                  f"{ak1:9.4f}  {ak4:9.4f}  {'n/a':>10}  {'n/a':>10}")
+
+    # --- Light HC comparison: S&W vs Li 2015 ---
+    print("\n--- Light Hydrocarbons: S&W Eq 8 vs Li et al. 2015 Pitzer ---")
+    print("    Li 2015 at P=100 bar, m=0 (infinite dilution) and m=2")
+    hc_gases = ['CH4', 'C2H6', 'C3H8', 'nC4H10']
+    print(f"{'T(C)':>6}", end="")
+    for g in hc_gases:
+        print(f"  {'SW '+g:>11}  {'Li m=0':>7}  {'Li m=2':>7}", end="")
+    print()
+    print("-" * 118)
+    for Tc, Tk in zip(temps_C, temps_K):
+        print(f"{Tc:>6}", end="")
+        for g in hc_gases:
+            sw = ks_sw_eq8(Tk, g)
+            li0 = ks_li2015(Tk, g, m_NaCl=0.0, P_bar=100.0)
+            li2 = ks_li2015(Tk, g, m_NaCl=2.0, P_bar=100.0)
+            print(f"  {sw:11.4f}  {li0:7.4f}  {li2:7.4f}", end="")
+        print()
+
+    # --- All S&W gases ---
+    print("\n--- S&W Eq 8 for all gases (m-independent) ---")
+    gases = ['H2', 'N2', 'CH4', 'C2H6', 'C3H8', 'nC4H10', 'CO2', 'H2S']
+    header = f"{'T(C)':>6}" + "".join(f"  {g:>8}" for g in gases)
+    print(header)
+    print("-" * (8 + 10 * len(gases)))
+    for Tc, Tk in zip(temps_C, temps_K):
+        row = f"{Tc:>6}"
+        for g in gases:
+            row += f"  {ks_sw_eq8(Tk, g):8.4f}"
+        print(row)
+
+    # --- N2 comparison: S&W vs Mao & Duan 2006 ---
+    print("\n--- N2: S&W Eq 8 vs Mao & Duan 2006 Pitzer ---")
+    print("    Mao 2006 at P=100 bar, various m_NaCl")
+    print(f"{'T(C)':>6}  {'S&W Eq8':>8}  {'Mao m=0':>8}  {'Mao m=1':>8}  "
+          f"{'Mao m=2':>8}  {'Mao m=4':>8}")
+    print("-" * 56)
+    for Tc, Tk in zip(temps_C, temps_K):
+        if Tk > 400:
+            # Mao 2006 valid to 400 K only
+            sw = ks_sw_eq8(Tk, 'N2')
+            print(f"{Tc:>6}  {sw:8.4f}  {'n/a':>8}  {'n/a':>8}  {'n/a':>8}  {'n/a':>8}")
+        else:
+            sw = ks_sw_eq8(Tk, 'N2')
+            m0 = ks_mao2006_n2(Tk, m_NaCl=0.0, P_bar=100.0)
+            m1 = ks_mao2006_n2(Tk, m_NaCl=1.0, P_bar=100.0)
+            m2 = ks_mao2006_n2(Tk, m_NaCl=2.0, P_bar=100.0)
+            m4 = ks_mao2006_n2(Tk, m_NaCl=4.0, P_bar=100.0)
+            print(f"{Tc:>6}  {sw:8.4f}  {m0:8.4f}  {m1:8.4f}  {m2:8.4f}  {m4:8.4f}")
+
+    # --- N2: Pressure sensitivity (Mao 2006) ---
+    print("\n--- N2: Mao 2006 ks at different pressures (m=0) ---")
+    print(f"{'T(C)':>6}  {'1 bar':>8}  {'50 bar':>8}  {'100 bar':>8}  "
+          f"{'200 bar':>8}  {'500 bar':>8}  {'S&W':>8}")
+    print("-" * 62)
+    for Tc, Tk in zip(temps_C, temps_K):
+        if Tk > 400:
+            sw = ks_sw_eq8(Tk, 'N2')
+            print(f"{Tc:>6}  {'n/a':>8}  {'n/a':>8}  {'n/a':>8}  {'n/a':>8}  {'n/a':>8}  {sw:8.4f}")
+        else:
+            vals = [ks_mao2006_n2(Tk, 0, P) for P in [1, 50, 100, 200, 500]]
+            sw = ks_sw_eq8(Tk, 'N2')
+            print(f"{Tc:>6}" + "".join(f"  {v:8.4f}" for v in vals) + f"  {sw:8.4f}")
+
+    # --- CH4: Pressure sensitivity (Li 2015) ---
+    print("\n--- CH4: Li 2015 ks at different pressures (m=0) ---")
+    print(f"{'T(C)':>6}  {'1 bar':>8}  {'50 bar':>8}  {'100 bar':>8}  "
+          f"{'200 bar':>8}  {'500 bar':>8}  {'S&W':>8}")
+    print("-" * 62)
+    for Tc, Tk in zip(temps_C, temps_K):
+        vals = [ks_li2015(Tk, 'CH4', 0, P) for P in [1, 50, 100, 200, 500]]
+        sw = ks_sw_eq8(Tk, 'CH4')
+        print(f"{Tc:>6}" + "".join(f"  {v:8.4f}" for v in vals) + f"  {sw:8.4f}")
+
+    # --- H2S model comparison summary ---
+    print("\n--- H2S: Akinfiev table ks at multiple pressures (log10) ---")
+    print("  Showing that ks is largely P-independent (as expected)")
+    # Tables at 1 MPa and 10 MPa for m=2
+    T_tab = [323.15, 373.15, 423.15, 473.15, 523.15]
+    m0_1 = [0.596, 0.297, 0.136, np.nan, np.nan]     # pure water, 1 MPa
+    m2_1 = [0.467, 0.241, 0.111, np.nan, np.nan]      # 2m NaCl, 1 MPa
+    m0_5 = [2.026, 1.651, 1.271, 0.956, 0.300]        # pure water, 5 MPa
+    m2_5 = [1.523, 1.259, 0.975, 0.726, 0.226]        # 2m NaCl, 5 MPa
+    m0_10 = [2.100, 2.703, 2.752, 2.622, 2.045]       # pure water, 10 MPa
+    m2_10 = [1.577, 1.980, 1.969, 1.834, 1.411]       # 2m NaCl, 10 MPa
+    m0_20 = [2.237, 3.065, 4.407, 5.395, 5.458]       # pure water, 20 MPa
+    m2_20 = [1.676, 2.224, 3.039, 3.659, 3.677]       # 2m NaCl, 20 MPa
+
+    print(f"{'T(C)':>6}  {'1 MPa':>8}  {'5 MPa':>8}  {'10 MPa':>8}  {'20 MPa':>8}")
+    print("-" * 46)
+    for i, Tk in enumerate(T_tab):
+        Tc = Tk - 273.15
+        vals = []
+        for mw, mb in [(m0_1[i], m2_1[i]), (m0_5[i], m2_5[i]),
+                        (m0_10[i], m2_10[i]), (m0_20[i], m2_20[i])]:
+            if np.isnan(mw) or np.isnan(mb) or mb <= 0:
+                vals.append("   n/a  ")
+            else:
+                ks = np.log10(mw / mb) / 2.0
+                vals.append(f"{ks:8.4f}")
+        print(f"{Tc:>6.0f}  {'  '.join(vals)}")

pyrestoolbox/plyasunov/__init__.py

+"""
+Plyasunov apparent molar volume model for dissolved gases in water.
+
+Provides V_phi(gas, T_K, P_MPa) and gas_mw(gas) for:
+    H2, N2, CH4, CO2, C2H6, C3H8, NC4H10, H2S
+
+Based on:
+    Plyasunov, A.V. Fluid Phase Equilibria (2019, 2020, 2021) Parts I-III.
+    IAPWS-IF97 Region 1 for pure water properties.
+"""
+
+from .plyasunov_model import V_phi, gas_mw, V2_inf, B12, A12_inf, MW_GAS

pyrestoolbox/plyasunov/iapws_if97.py

+"""
+IAPWS-IF97 Region 1: Compressed liquid water properties.
+
+Standalone implementation with no external dependencies beyond numpy.
+
+Provides:
+    - rho_if97(T, P): pure water density in kg/m3
+    - kappa_T_if97(T, P): isothermal compressibility in 1/MPa
+    - rho_and_kappa(T, P): both in a single call
+
+Valid range (Region 1):
+    273.15 K <= T <= 623.15 K  (0-350 C)
+    P_sat(T) <= P <= 100 MPa
+
+Reference:
+    Wagner, W. et al. (2000). "The IAPWS Industrial Formulation 1997
+    for the Thermodynamic Properties of Water and Steam."
+    ASME J. Eng. Gas Turbines Power, 122(1), 150-182.
+
+Units: T in K, P in MPa
+"""
+
+# Constants
+R_SPECIFIC = 461.526e-6  # Specific gas constant for water [MPa*m3/(kg*K)]
+P_STAR = 16.53           # Reference pressure [MPa]
+T_STAR = 1386.0          # Reference temperature [K]
+
+# Region 1 coefficients (Table 2 of IAPWS-IF97)
+# Each row: (I_i, J_i, n_i)
+_REGION1_IJN = [
+    (0,  -2,   0.14632971213167e+00),
+    (0,  -1,  -0.84548187389013e+00),
+    (0,   0,  -0.37563603672040e+01),
+    (0,   1,   0.33855169168385e+01),
+    (0,   2,  -0.95791963387872e+00),
+    (0,   3,   0.15772038513228e+00),
+    (0,   4,  -0.16616417199501e-01),
+    (0,   5,   0.81214629983568e-03),
+    (1,  -9,   0.28319080123804e-03),
+    (1,  -7,  -0.60706301565874e-03),
+    (1,  -1,  -0.18990068218419e-01),
+    (1,   0,  -0.32529748770505e-01),
+    (1,   1,  -0.21841717175414e-01),
+    (1,   3,  -0.52838357969930e-04),
+    (2,  -3,  -0.47184321073267e-03),
+    (2,   0,  -0.30001780793026e-03),
+    (2,   1,   0.47661393906987e-04),
+    (2,   3,  -0.44141845330846e-05),
+    (2,  17,  -0.72694996297594e-15),
+    (3,  -4,  -0.31679644845054e-04),
+    (3,   0,  -0.28270797985312e-05),
+    (3,   6,  -0.85205128120103e-09),
+    (4,  -5,  -0.22425281908000e-05),
+    (4,  -2,  -0.65171222895601e-06),
+    (4,  10,  -0.14341729937924e-12),
+    (5,  -8,  -0.40516996860117e-06),
+    (8, -11,  -0.12734301741682e-08),
+    (8,  -6,  -0.17424871230634e-09),
+    (21, -29, -0.68762131295531e-18),
+    (23, -31,  0.14478307828521e-19),
+    (29, -38,  0.26335781662795e-22),
+    (30, -39, -0.11947622640071e-22),
+    (31, -40,  0.18228094581404e-23),
+    (32, -41, -0.93537087292458e-25),
+]
+
+
+def _gamma_derivatives(pi, tau):
+    """
+    Compute gamma_pi and gamma_pipi for Region 1.
+
+    gamma = sum( n_i * (7.1 - pi)^I_i * (tau - 1.222)^J_i )
+
+    Returns:
+        (gamma_pi, gamma_pipi)
+    """
+    a = 7.1 - pi
+    b = tau - 1.222
+
+    gp = 0.0
+    gpp = 0.0
+
+    for I, J, n in _REGION1_IJN:
+        bJ = b ** J
+        if I == 0:
+            continue
+        elif I == 1:
+            gp += -n * bJ
+        else:
+            aI1 = a ** (I - 1)
+            gp += n * (-I) * aI1 * bJ
+            gpp += n * I * (I - 1) * (aI1 / a) * bJ
+
+    return gp, gpp
+
+
+def rho_if97(T, P):
+    """
+    Pure water density from IAPWS-IF97 Region 1.
+
+    Parameters:
+        T: temperature in K (273.15 - 623.15)
+        P: pressure in MPa (up to 100)
+
+    Returns:
+        density in kg/m3
+    """
+    pi = P / P_STAR
+    tau = T_STAR / T
+    gp, _ = _gamma_derivatives(pi, tau)
+    return P_STAR / (R_SPECIFIC * T * gp)
+
+
+def kappa_T_if97(T, P):
+    """
+    Isothermal compressibility of pure water from IAPWS-IF97 Region 1.
+
+    Parameters:
+        T: temperature in K (273.15 - 623.15)
+        P: pressure in MPa (up to 100)
+
+    Returns:
+        kappa_T in 1/MPa
+    """
+    pi = P / P_STAR
+    tau = T_STAR / T
+    gp, gpp = _gamma_derivatives(pi, tau)
+    return -gpp / (gp * P_STAR)
+
+
+def rho_and_kappa(T, P):
+    """
+    Compute both density and compressibility in a single call.
+
+    Returns:
+        (rho [kg/m3], kappa_T [1/MPa])
+    """
+    pi = P / P_STAR
+    tau = T_STAR / T
+    gp, gpp = _gamma_derivatives(pi, tau)
+    rho = P_STAR / (R_SPECIFIC * T * gp)
+    kappa = -gpp / (gp * P_STAR)
+    return rho, kappa

pyrestoolbox/plyasunov/plyasunov_model.py

+"""
+Plyasunov A12-infinity model for computing V2-infinity (infinite-dilution
+partial molar volume) of dissolved gases in water.
+
+Covers 8 gases using unified equation form from:
+  - Part II (Plyasunov 2020): CO2, C2H6, C3H8, n-C4H10
+  - Part III (Plyasunov 2021): H2S
+  - Part IV (Plyasunov 2021): H2, N2, CH4
+
+Master equation:
+    A12_inf = 1 + rho1* * [a0 + a1*rho1* + a2*(rho1*)^2 + ... + a5*(rho1*)^5]
+
+Recovery of V2_inf:
+    V2_inf(T,P) = A12_inf(T, rho1*) * kappa_T(T,P) * R * T
+
+Units:
+    T in K, P in MPa, rho1* in kg/m3, V2_inf in cm3/mol, kappa_T in 1/MPa
+
+References:
+    Plyasunov, A.V. (2019). "Values of the apparent molar volumes V_phi and
+    the osmotic coefficients phi of NaCl(aq) at infinite dilution...Part I:
+    Aqueous solutions of H2, N2, and CH4." Fluid Phase Equilibria, 496, 43-51.
+
+    Plyasunov, A.V. (2020). "...Part II: CO2, C2H6, C3H8, n-C4H10."
+    Fluid Phase Equilibria, 523, 112757.
+
+    Plyasunov, A.V. (2021). "...Part III: H2S."
+    Fluid Phase Equilibria, 530, 112883.
+"""
+
+import numpy as np
+from .water_properties import rho_w, kappa_T, TC_WATER, MW_WATER
+
+# Constants
+R_CM3_MPA = 8.314462   # cm3*MPa/(mol*K)
+NA = 6.02214076e23     # Avogadro's number
+OMEGA = 1e-3 / MW_WATER  # converts B12 from cm3/mol to m3/kg
+
+
+# ============================================================================
+# B12 CROSS SECOND VIRIAL COEFFICIENTS
+# ============================================================================
+
+def _B12_polynomial(T, coeffs):
+    """
+    B12 for simple fluids (H2, N2, CH4) from Part IV Eq. 5:
+    B12 = sum( ai * (T/100)^bi )  [cm3/mol]
+    coeffs: list of (a, b) tuples, 4 terms
+    """
+    T100 = T / 100.0
+    result = 0.0
+    for a, b in coeffs:
+        result += a * T100**b
+    return result
+
+
+def _B12_CO2(T):
+    """
+    B12 for CO2 from Part II Eq. 8 (Hellmann 2019):
+    B12 = b1 + b2/(T*)^0.5 + b3/T* + b4/(T*)^3 + b5/(T*)^6 + b6/(T*)^(21/2)
+    where T* = T/100, result in cm3/mol
+    """
+    Tstar = T / 100.0
+    b = [15.210, 149.72, -534.54, -2234.6, -13017.0, -39482.0]
+    return (b[0]
+            + b[1] / Tstar**0.5
+            + b[2] / Tstar
+            + b[3] / Tstar**3
+            + b[4] / Tstar**6
+            + b[5] / Tstar**10.5)
+
+
+def _B12_square_well(T, sigma, eps_kB, lam):
+    """
+    B12 from square-well potential (Part II Eq. 9):
+    B12 = (2/3)*pi*NA*sigma12^3 * {1 - (lambda12^3-1)*[exp(epsilon12/(kB*T)) - 1]}
+
+    Parameters:
+        sigma: collision diameter in Angstrom
+        eps_kB: well depth / kB in K
+        lam: well width in collision diameters
+
+    Returns: B12 in cm3/mol
+    """
+    sigma_cm = sigma * 1e-8
+    b_hs = (2.0 / 3.0) * np.pi * NA * sigma_cm**3
+    exp_term = np.exp(eps_kB / T) - 1.0
+    return b_hs * (1.0 - (lam**3 - 1.0) * exp_term)
+
+
+def _B12_group_contribution(T, groups):
+    """
+    B12 for hydrocarbons via group contributions:
+    B12(compound) = sum( n_i * B12(group_i) )
+    groups: list of (n_copies, sigma, eps_kB, lambda) tuples
+    """
+    total = 0.0
+    for n, sigma, eps_kB, lam in groups:
+        total += n * _B12_square_well(T, sigma, eps_kB, lam)
+    return total
+
+
+# Square-well parameters for functional groups (Part II, Table 2)
+_SW_CH3 = (2.788, 283.8, 1.430)
+_SW_CH2 = (2.226, 271.4, 1.430)
+_SW_H2S = (2.85, 650.0, 1.324)  # Part III
+
+# B12 polynomial coefficients (Part IV, Table 1): list of (a, b) tuples
+_B12_POLY_H2 = [(33.047, -0.21), (-250.41, -1.50), (285.42, -2.26), (-186.78, -3.21)]
+_B12_POLY_N2 = [(156.679, -0.33), (-183.541, -0.57), (-194.330, -1.47), (-154.815, -3.66)]
+_B12_POLY_CH4 = [(109.22, -0.2), (-202.52, -0.6), (-235.86, -2.0), (-297.63, -3.0)]
+
+
+def B12(gas, T):
+    """
+    Cross second virial coefficient B12 for gas-water interaction.
+
+    Parameters:
+        gas: gas name string (case-insensitive)
+        T: temperature in K
+
+    Returns:
+        B12 in cm3/mol
+    """
+    gas = gas.upper()
+    if gas == 'H2':
+        return _B12_polynomial(T, _B12_POLY_H2)
+    elif gas == 'N2':
+        return _B12_polynomial(T, _B12_POLY_N2)
+    elif gas == 'CH4':
+        return _B12_polynomial(T, _B12_POLY_CH4)
+    elif gas == 'CO2':
+        return _B12_CO2(T)
+    elif gas == 'C2H6':
+        return _B12_group_contribution(T, [(2, *_SW_CH3)])
+    elif gas == 'C3H8':
+        return _B12_group_contribution(T, [(2, *_SW_CH3), (1, *_SW_CH2)])
+    elif gas in ('NC4H10', 'N-C4H10', 'NC4'):
+        return _B12_group_contribution(T, [(2, *_SW_CH3), (2, *_SW_CH2)])
+    elif gas == 'H2S':
+        return _B12_square_well(T, *_SW_H2S)
+    else:
+        raise ValueError(f"Unknown gas: {gas}")
+
+
+# ============================================================================
+# p_in COEFFICIENT TABLES
+# ============================================================================
+
+# Format: p_in[i][n] where i=1..5, n=0..6
+# Each gas has a 5x7 array of coefficients
+
+# Paper 8 (Part IV, Table 2) - full precision - supersedes Paper 5
+P_IN_H2 = [
+    [2.1202801e-3, -6.037483e-3, 1.0195766e-2, -1.0693934e-2, 5.4991895e-3, -1.0878264e-3, 1.1746219e-6],
+    [-1.1046217e-5, 3.1667503e-5, -5.3803531e-5, 5.6835389e-5, -2.9567299e-5, 5.9819654e-6, -2.019032e-8],
+    [1.9307298e-8, -5.4979556e-8, 9.2491093e-8, -9.7577108e-8, 5.1231461e-8, -1.0724676e-8, 1.2613281e-10],
+    [-1.3252458e-11, 3.7335356e-11, -6.1484868e-11, 6.403194e-11, -3.3547995e-11, 7.4031616e-12, -3.4170589e-13],
+    [3.011170e-15, -8.3139044e-15, 1.3052306e-14, -1.2814023e-14, 6.0470404e-15, -9.8662924e-16, -4.1237636e-17],
+]
+
+P_IN_N2 = [
+    [-2.2070753e-3, 6.6067485e-3, -1.2027859e-2, 1.3728692e-2, -7.7929165e-3, 1.7076684e-3, -2.8424026e-6],
+    [1.1431298e-5, -3.3833752e-5, 6.1522458e-5, -7.1045896e-5, 4.1165829e-5, -9.3369807e-6, 4.7959089e-8],
+    [-1.8029853e-8, 5.3340827e-8, -9.7754731e-8, 1.1514188e-7, -6.8971366e-8, 1.667536e-8, -3.0103987e-10],
+    [1.0436147e-11, -3.0813479e-11, 5.7057178e-11, -6.8981444e-11, 4.3377118e-11, -1.2090314e-11, 9.5828102e-13],
+    [-1.6676726e-15, 4.8355967e-15, -8.7551827e-15, 1.0376018e-14, -6.124985e-15, 1.3543054e-15, -1.6688106e-18],
+]
+
+P_IN_CH4 = [
+    [-1.3023759e-3, 3.8481069e-3, -6.8502936e-3, 7.7002273e-3, -4.3454809e-3, 9.5403855e-4, -1.8741835e-6],
+    [7.3809403e-6, -2.1303297e-5, 3.7100486e-5, -4.1343101e-5, 2.3418494e-5, -5.2627676e-6, 3.0296735e-8],
+    [-1.2366446e-8, 3.5518855e-8, -6.1394616e-8, 6.8386294e-8, -3.9321302e-8, 9.3024856e-9, -1.7821599e-10],
+    [7.5953065e-12, -2.1583027e-11, 3.6743322e-11, -4.050275e-11, 2.3523613e-11, -6.1571572e-12, 4.6261515e-13],
+    [-1.3716224e-15, 3.7669953e-15, -5.924387e-15, 5.8112619e-15, -2.7435547e-15, 3.8223834e-16, 5.3380661e-17],
+]
+
+# Paper 6 (Part II, Table 5) - full precision
+P_IN_CO2 = [
+    [9.1583134e-4, -2.7160016e-3, 5.0839423e-3, -5.8354599e-3, 3.2543133e-3, -7.0687955e-4, 1.4494891e-6],
+    [-3.6106209e-6, 1.1315777e-5, -2.2758954e-5, 2.7797001e-5, -1.6473215e-5, 3.8439408e-6, -2.557928e-8],
+    [5.0157576e-9, -1.5637417e-8, 3.200829e-8, -4.1225984e-8, 2.6278766e-8, -6.8414637e-9, 1.6610365e-10],
+    [-2.0446824e-12, 6.2444294e-12, -1.3130671e-11, 1.9003321e-11, -1.3986032e-11, 4.6916137e-12, -5.1705069e-13],
+    [-1.1881416e-16, 4.4574039e-16, -7.9103055e-16, 2.5194274e-16, 2.4022155e-16, -9.7895994e-17, -1.4094429e-17],
+]
+
+P_IN_C2H6 = [
+    [-5.1475443e-4, 1.729277e-3, -3.6536764e-3, 4.6194269e-3, -2.7977608e-3, 6.2657963e-4, -7.9712465e-8],
+    [3.5428879e-6, -1.0888811e-5, 2.1555401e-5, -2.6562891e-5, 1.5872915e-5, -3.5330389e-6, 1.9105807e-9],
+    [-7.076262e-9, 2.150819e-8, -4.1847563e-8, 5.0828275e-8, -3.0091634e-8, 6.7020044e-9, -1.7913067e-11],
+    [5.7671272e-12, -1.7389573e-11, 3.3600113e-11, -4.0403444e-11, 2.3837864e-11, -5.457811e-12, 1.0241543e-13],
+    [-1.6279128e-15, 4.8764202e-15, -9.3154249e-15, 1.1023937e-14, -6.3736955e-15, 1.4100215e-15, -1.4735524e-17],
+]
+
+P_IN_C3H8 = [
+    [2.2892124e-3, -6.4740679e-3, 1.1030983e-2, -1.1944689e-2, 6.4760838e-3, -1.3894223e-3, 3.3011188e-6],
+    [-1.029341e-5, 3.0090181e-5, -5.3316193e-5, 5.9603698e-5, -3.35325e-5, 7.5922124e-6, -5.4814576e-8],
+    [1.5037316e-8, -4.4050725e-8, 7.8975079e-8, -9.0810405e-8, 5.3609476e-8, -1.3301616e-8, 3.3278243e-10],
+    [-7.3401228e-12, 2.1306209e-11, -3.7787336e-11, 4.4705147e-11, -2.8278765e-11, 8.5764552e-12, -9.1713991e-13],
+    [5.5117215e-16, -1.4000747e-15, 1.7856523e-15, -1.6585346e-15, 6.3567082e-16, 7.7519099e-17, -7.8311547e-17],
+]
+
+P_IN_NC4H10 = [
+    [2.5603794e-3, -6.9646157e-3, 1.0987245e-2, -1.1110696e-2, 5.7384581e-3, -1.2233119e-3, 4.9218547e-6],
+    [-8.9460753e-6, 2.5173896e-5, -4.1659657e-5, 4.4713209e-5, -2.5204428e-5, 6.1085451e-6, -8.1017199e-8],
+    [7.0204532e-9, -1.9465263e-8, 3.2883394e-8, -4.0132208e-8, 2.8071928e-8, -9.1449293e-9, 4.8274281e-10],
+    [3.6497249e-12, -1.1330901e-11, 2.0340955e-11, -1.7139410e-11, 2.6365251e-12, 3.4703436e-12, -1.2475357e-12],
+    [-3.7951579e-15, 1.1392044e-14, -2.075511e-14, 2.2507396e-14, -1.20847e-14, 2.7449667e-15, -1.3968834e-16],
+]
+
+# Paper 7 (Part III) - full precision from Table 5
+P_IN_H2S = [
+    [6.7014672e-4, -1.9403321e-3, 3.5360019e-3, -4.0122786e-3, 2.2524974e-3, -4.8843038e-4, 1.6044152e-6],
+    [-3.1805655e-6, 9.5800671e-6, -1.7994973e-5, 2.1087746e-5, -1.2352096e-5, 2.8686985e-6, -2.7620761e-8],
+    [5.0089935e-9, -1.4885595e-8, 2.7733075e-8, -3.3416171e-8, 2.0773681e-8, -5.4162153e-9, 1.7656758e-10],
+    [-2.5739894e-12, 7.422722e-12, -1.3339226e-11, 1.677997e-11, -1.1669169e-11, 4.0370617e-12, -5.6862157e-13],
+    [1.6189502e-16, -3.2073034e-16, 1.4144345e-16, -7.3004629e-17, 5.7935525e-17, -3.11107e-18, 1.3146726e-18],
+]
+
+# Map gas names to p_in tables
+P_IN_MAP = {
+    'H2': P_IN_H2,
+    'N2': P_IN_N2,
+    'CH4': P_IN_CH4,
+    'CO2': P_IN_CO2,
+    'C2H6': P_IN_C2H6,
+    'C3H8': P_IN_C3H8,
+    'NC4H10': P_IN_NC4H10,
+    'N-C4H10': P_IN_NC4H10,
+    'NC4': P_IN_NC4H10,
+    'H2S': P_IN_H2S,
+}
+
+# Molecular weights of dissolved gases (g/mol)
+MW_GAS = {
+    'H2': 2.01588,
+    'N2': 28.0134,
+    'CH4': 16.0425,
+    'CO2': 44.0095,
+    'C2H6': 30.069,
+    'C3H8': 44.096,
+    'NC4H10': 58.122,
+    'N-C4H10': 58.122,
+    'NC4': 58.122,
+    'H2S': 34.081,
+}
+
+
+# ============================================================================
+# CORE MODEL
+# ============================================================================
+
+def _compute_ai(p_in_row, T, i_index):
+    """
+    Compute temperature-dependent coefficient a_i from p_in coefficients.
+
+    Eq. 34: ai = pi0*theta^(-0.5) + pi1*theta^(-1) + pi2*theta^(-2) + ...
+                 + pi5*theta^(-5) + pi6*theta^(-n)
+    where theta = T/Tc, n = 9(i=1), 8(i=2), 7(i=3), 6(i=4,5)
+    i_index: 1-based index (1-5)
+    """
+    theta = T / TC_WATER
+    exponents = [-0.5, -1.0, -2.0, -3.0, -4.0, -5.0]
+    n_map = {1: 9, 2: 8, 3: 7, 4: 6, 5: 6}
+    last_exp = -n_map[i_index]
+
+    result = 0.0
+    for k in range(6):
+        result += p_in_row[k] * theta**exponents[k]
+    result += p_in_row[6] * theta**last_exp
+
+    return result
+
+
+def A12_inf(gas, T, rho1_star):
+    """
+    Compute A12_inf = V2_inf/(kappa_T*R*T) at given T and water density.
+
+    Eq. 32: A12_inf = 1 + rho1*[a0 + a1*rho1* + ... + a5*(rho1*)^5]
+    where a0 = 2*Omega*B12
+
+    Parameters:
+        gas: gas name string
+        T: temperature in K
+        rho1_star: pure water density in kg/m3
+
+    Returns:
+        A12_inf (dimensionless)
+    """
+    gas = gas.upper()
+    if gas not in P_IN_MAP:
+        raise ValueError(f"Unknown gas: {gas}. Available: {list(P_IN_MAP.keys())}")
+
+    p_in = P_IN_MAP[gas]
+
+    a0 = 2.0 * OMEGA * B12(gas, T)
+    a_coeffs = [a0]
+    for i in range(5):
+        a_coeffs.append(_compute_ai(p_in[i], T, i + 1))
+
+    rho = rho1_star
+    poly_sum = 0.0
+    rho_power = 1.0
+    for k in range(6):
+        poly_sum += a_coeffs[k] * rho_power
+        rho_power *= rho
+
+    return 1.0 + rho * poly_sum
+
+
+def V2_inf(gas, T, P):
+    """
+    Compute V2_inf (infinite-dilution partial molar volume) at given T, P.
+
+    V2_inf = A12_inf * kappa_T * R * T
+
+    Parameters:
+        gas: gas name string (case-insensitive)
+        T: temperature in K
+        P: pressure in MPa
+
+    Returns:
+        V2_inf in cm3/mol
+    """
+    rho = rho_w(T, P)
+    kt = kappa_T(T, P)
+    a12 = A12_inf(gas, T, rho)
+    return a12 * kt * R_CM3_MPA * T
+
+
+def V_phi(gas, T, P):
+    """
+    Apparent molar volume V_phi at infinite dilution.
+    Alias for V2_inf.
+
+    Parameters:
+        gas: gas name string (case-insensitive).
+             Supported: H2, N2, CH4, CO2, C2H6, C3H8, NC4H10, H2S
+        T: temperature in K
+        P: pressure in MPa
+
+    Returns:
+        V_phi in cm3/mol
+    """
+    return V2_inf(gas, T, P)
+
+
+def gas_mw(gas):
+    """Return molecular weight of gas in g/mol."""
+    gas = gas.upper()
+    if gas not in MW_GAS:
+        raise ValueError(f"Unknown gas: {gas}")
+    return MW_GAS[gas]

pyrestoolbox/plyasunov/water_properties.py

+"""
+Pure water properties for the Plyasunov model.
+
+Provides:
+    - rho_w(T, P): pure water density in kg/m3
+    - kappa_T(T, P): isothermal compressibility in 1/MPa
+
+Backend: IAPWS-IF97 Region 1 (no external dependencies).
+
+Parameters:
+    T: temperature in Kelvin
+    P: pressure in MPa
+"""
+
+from .iapws_if97 import rho_if97, kappa_T_if97
+
+# Constants
+MW_WATER = 18.015268   # g/mol, molar mass of water (IAPWS)
+TC_WATER = 647.096     # K, critical temperature of water
+
+
+def rho_w(T, P):
+    """
+    Pure water density.
+
+    Parameters:
+        T: temperature in K
+        P: pressure in MPa
+
+    Returns:
+        density in kg/m3
+    """
+    return rho_if97(float(T), float(P))
+
+
+def kappa_T(T, P):
+    """
+    Isothermal compressibility of pure water.
+
+    Parameters:
+        T: temperature in K
+        P: pressure in MPa
+
+    Returns:
+        kappa_T in 1/MPa
+    """
+    return kappa_T_if97(float(T), float(P))

pyrestoolbox/tests/test_unified_brine_design.py

+#!/usr/bin/env python3
+"""
+Test scripts for the unified brine density & viscosity framework design.
+
+These are exploratory/validation tests, NOT part of the main test suite yet.
+They help understand the current implementations and validate design decisions.
+
+Run with: PYTHONPATH=/home/mark/projects python3 tests/test_unified_brine_design.py
+"""
+
+import sys
+import numpy as np
+
+sys.path.insert(0, '/home/mark/projects')
+
+from pyrestoolbox import brine
+
+# ============================================================================
+# R12: Compare Spivey Pitzer CH4 density vs Garcia-form CH4 density
+# ============================================================================
+# The Spivey/McCain approach (brine_props) computes CH4-saturated brine density
+# using Pitzer interaction parameter derivatives (Eq 4.22-4.24).
+#
+# The Garcia approach uses V_phi(T) polynomial + Eq 18 mixing rule.
+# For CO2: V_phi = 37.51 - 0.09585*T + 8.74e-4*T^2 - 5.044e-7*T^3
+#
+# We need to understand: what effective V_phi does the Spivey method imply
+# for CH4, and how does it compare to literature values?
+# ============================================================================
+
+def extract_spivey_ch4_vphi():
+    """
+    Extract the effective apparent molar volume of dissolved CH4
+    implied by the Spivey/McCain approach at various temperatures.
+
+    The Spivey method computes:
+        Vmch4b = R*T * (du/dp + 2*m*dlambda/dp + m^2*deta/dp)  [Eq 4.22]
+    where Vmch4b is partial molar volume in cm3/mol at the given T,P,salinity.
+
+    We'll extract Vmch4b at various conditions and compare to published
+    V_phi^inf values for CH4 in water.
+    """
+    from pyrestoolbox.brine.brine import _Eq41, _RHOW_T70_ARR, _EWT_ARR, _FWT_ARR
+    from pyrestoolbox.brine.brine import _DM2T_ARR, _DM32T_ARR, _DM1T_ARR, _DM12T_ARR
+    from pyrestoolbox.brine.brine import _EMT_ARR, _FM32T_ARR, _FM1T_ARR, _FM12T_ARR
+
+    print("=" * 80)
+    print("R12: Spivey Pitzer CH4 Partial Molar Volume vs Temperature")
+    print("=" * 80)
+    print()
+
+    # Spivey Pitzer interaction parameter coefficients for CH4
+    u_arr = [0, 8.3143711, -7.2772168e-4, 2.1489858e3, -1.4019672e-5,
+             -6.6743449e5, 7.698589e-2, -5.0253331e-5, -30.092013, 4.8468502e3, 0]
+    lambda_arr = [0, -0.80898, 1.0827e-3, 183.85, 0, 0, 3.924e-4, 0, 0, 0, -1.97e-6]
+
+    print(f"{'T(°F)':>8} {'T(°C)':>8} {'P(psia)':>8} {'Wt%':>6} {'Vmch4b':>10} {'rho_brine':>10} {'rho_ch4sat':>12} {'delta_rho%':>10}")
+    print("-" * 90)
+
+    results = []
+
+    for degf in [100, 150, 200, 250, 300]:
+        for p in [2000, 5000, 8000]:
+            for wt in [0, 3, 10]:
+                degc = (degf - 32) / 1.8
+                degk = degc + 273
+                Mpa = p * 0.00689476
+                m = 1000 * (wt / 100) / (58.4428 * (1 - wt / 100)) if wt > 0 else 0
+
+                # Compute du/dp and dlambda/dp (pressure derivatives of Pitzer params)
+                dudptm = (u_arr[6] + u_arr[7] * degk + u_arr[8] / degk + u_arr[9] / (degk * degk))
+                dlambdadptm = lambda_arr[6] + 2 * lambda_arr[10] * Mpa
+
+                # Vmch4b = R*T*(du/dp + 2*m*dlambda/dp)  [cm3/mol]
+                Vmch4b = 8.314467 * degk * (dudptm + 2 * m * dlambdadptm)
+
+                # Also get full brine_props results for density comparison
+                bw, lden, visw, cw, rsw = brine.brine_props(p=p, degf=degf, wt=wt, ch4_sat=1.0)
+                bw0, lden0, visw0, cw0, rsw0 = brine.brine_props(p=p, degf=degf, wt=wt, ch4_sat=0.0)
+
+                delta_pct = (lden - lden0) / lden0 * 100
+
+                results.append({
+                    'degf': degf, 'degc': degc, 'p': p, 'wt': wt,
+                    'Vmch4b': Vmch4b, 'rho_nogas': lden0, 'rho_ch4sat': lden,
+                    'delta_pct': delta_pct
+                })
+
+                if wt == 0 and p == 5000:  # Print subset
+                    print(f"{degf:>8.0f} {degc:>8.1f} {p:>8.0f} {wt:>6.0f} {Vmch4b:>10.2f} {lden0:>10.5f} {lden:>12.5f} {delta_pct:>10.3f}")
+
+    print()
+    print("Key observations:")
+    print("  - Vmch4b is the Spivey/Pitzer partial molar volume of CH4 in brine (cm3/mol)")
+    print("  - Garcia CO2 V_phi at 25°C = 37.51 cm3/mol; literature CH4 V_phi ~37 cm3/mol")
+    print("  - delta_rho% shows density change from dissolved CH4 (should be NEGATIVE for CH4)")
+    print()
+
+    # Summary statistics
+    vmch4_values = [r['Vmch4b'] for r in results if r['wt'] == 0]
+    print(f"  Vmch4b range (fresh water): {min(vmch4_values):.2f} to {max(vmch4_values):.2f} cm3/mol")
+
+    delta_values = [r['delta_pct'] for r in results]
+    print(f"  Density change range: {min(delta_values):.3f}% to {max(delta_values):.3f}%")
+
+    return results
+
+
+def compare_garcia_co2_implementations():
+    """
+    Verify that the Garcia CO2 density correction in CO2_Brine_Mixture
+    produces consistent results. This is validation for TASK D6.
+    """
+    print()
+    print("=" * 80)
+    print("D6: Garcia CO2 density implementation verification")
+    print("=" * 80)
+    print()
+
+    print(f"{'T(°C)':>8} {'P(bar)':>8} {'ppm':>8} {'xCO2':>10} {'rho_co2sat':>12} {'rho_brine':>10} {'delta%':>8}")
+    print("-" * 75)
+
+    test_cases = [
+        (100, 85, 0),      # Moderate conditions, fresh water
+        (200, 85, 0),      # Higher pressure
+        (300, 85, 0),      # High pressure
+        (100, 85, 30000),  # 3% NaCl
+        (200, 85, 30000),
+        (300, 150, 0),     # High T high P
+        (200, 150, 50000), # High salinity
+    ]
+
+    for pbar, degc, ppm in test_cases:
+        try:
+            mix = brine.CO2_Brine_Mixture(pres=pbar, temp=degc, ppm=ppm, metric=True)
+            xCO2 = mix.x[0]
+            rho_sat = mix.bDen[0]
+            rho_brine = mix.bDen[1]
+            delta = (rho_sat - rho_brine) / rho_brine * 100
+
+            print(f"{degc:>8.0f} {pbar:>8.0f} {ppm:>8.0f} {xCO2:>10.5f} {rho_sat:>12.5f} {rho_brine:>10.5f} {delta:>8.3f}")
+        except Exception as e:
+            print(f"{degc:>8.0f} {pbar:>8.0f} {ppm:>8.0f} {'ERROR':>10} {str(e)[:40]}")
+
+    print()
+    print("Key observations:")
+    print("  - CO2 INCREASES brine density (positive delta)")
+    print("  - Effect is ~1-3% at typical reservoir conditions")
+    print("  - Effect decreases with temperature (V_phi increases)")
+
+
+def compare_viscosity_corrections():
+    """
+    Compare viscosity with and without dissolved gas corrections.
+    Shows magnitude of effect and validates the Islam-Carlson CO2 correction.
+    """
+    print()
+    print("=" * 80)
+    print("Viscosity correction magnitude comparison")
+    print("=" * 80)
+    print()
+
+    print("--- CO2 effect on brine viscosity (Islam-Carlson) ---")
+    print(f"{'T(°C)':>8} {'P(bar)':>8} {'ppm':>8} {'xCO2':>10} {'vis_co2':>10} {'vis_brine':>10} {'increase%':>10}")
+    print("-" * 75)
+
+    for pbar, degc, ppm in [(200, 85, 0), (200, 85, 30000), (300, 85, 0), (300, 150, 0), (100, 50, 0)]:
+        try:
+            mix = brine.CO2_Brine_Mixture(pres=pbar, temp=degc, ppm=ppm, metric=True)
+            xCO2 = mix.x[0]
+            vis_sat = mix.bVis[0]
+            vis_brine = mix.bVis[1]
+            increase = (vis_sat - vis_brine) / vis_brine * 100
+
+            print(f"{degc:>8.0f} {pbar:>8.0f} {ppm:>8.0f} {xCO2:>10.5f} {vis_sat:>10.4f} {vis_brine:>10.4f} {increase:>10.2f}")
+        except Exception as e:
+            print(f"{degc:>8.0f} {pbar:>8.0f} {ppm:>8.0f} {'ERROR':>10} {str(e)[:40]}")
+
+    print()
+    print("--- CH4 effect on brine viscosity (currently: NONE in Spivey) ---")
+    print(f"{'T(°F)':>8} {'P(psi)':>8} {'wt%':>6} {'vis_ch4sat':>12} {'vis_nogas':>10} {'change%':>10}")
+    print("-" * 60)
+
+    for degf, p, wt in [(200, 5000, 0), (200, 5000, 3), (250, 5000, 0), (150, 3000, 0)]:
+        bw1, den1, vis1, cw1, rsw1 = brine.brine_props(p=p, degf=degf, wt=wt, ch4_sat=1.0)
+        bw0, den0, vis0, cw0, rsw0 = brine.brine_props(p=p, degf=degf, wt=wt, ch4_sat=0.0)
+        change = (vis1 - vis0) / vis0 * 100
+
+        print(f"{degf:>8.0f} {p:>8.0f} {wt:>6.0f} {vis1:>12.4f} {vis0:>10.4f} {change:>10.2f}")
+
+    print()
+    print("Key observations:")
+    print("  - CO2: Islam-Carlson gives ~5-15% viscosity increase at saturation")
+    print("  - CH4: Currently 0% change (Spivey doesn't correct for dissolved CH4)")
+    print("  - Literature suggests CH4 effect is ~6% max (SPE-14211)")
+    print("  - For N2, H2: effect is negligible due to low solubility")
+
+
+def demonstrate_garcia_generalization():
+    """
+    Demonstrate how the Garcia mixing rule generalizes to any gas.
+    Uses the existing CO2 implementation as the template and shows
+    what the unified function signature would look like.
+    """
+    print()
+    print("=" * 80)
+    print("Demonstration: Garcia mixing rule with different gases")
+    print("=" * 80)
+    print()
+
+    # Garcia Eq 18: rho = (1 + MW_gas/MW_brine * x/(1-x)) / (V_phi * x/((1-x)*MW_brine) + 1/rho_brine)
+    # This is IDENTICAL for any gas — only V_phi and MW change
+
+    def garcia_density(rho_brine, T_celsius, x_gas, MW_brine, MW_gas, vphi_fn):
+        """Generic Garcia Eq 18"""
+        x_not_gas = 1.0 - x_gas
+        x_ratio = x_gas / x_not_gas
+        mw_ratio = MW_gas / MW_brine
+        vphi = vphi_fn(T_celsius)
+        return (1.0 + mw_ratio * x_ratio) / (vphi * x_ratio / MW_brine + 1.0 / rho_brine)
+
+    # V_phi polynomials: V_phi(T) = a0 + a1*T + a2*T^2 + a3*T^3 where T in °C
+    # CO2: from Garcia (2001) — these are the KNOWN coefficients
+    vphi_co2 = lambda T: 37.51 + T * (-0.09585 + T * (8.74e-4 + T * (-5.044e-7)))
+
+    # CH4: PLACEHOLDER — approximate from literature ~37 cm3/mol at 25°C
+    # Real coefficients need to be fitted from Hnedkovsky (1996) or Plyasunov (2019)
+    # For now, use a rough estimate based on similarity to CO2 V_phi shape
+    vphi_ch4_placeholder = lambda T: 37.0 + T * (0.10) + T**2 * (1.0e-4)  # PLACEHOLDER ONLY
+
+    # H2S: PLACEHOLDER — V_phi ~35 cm3/mol at 25°C (from extending_garcia doc)
+    vphi_h2s_placeholder = lambda T: 35.0 + T * (-0.05)  # PLACEHOLDER ONLY
+
+    MW = {'CO2': 44.01, 'CH4': 16.043, 'H2S': 34.082, 'N2': 28.014, 'H2': 2.016}
+
+    # Test conditions
+    T_celsius = 85.0
+    rho_brine = 1.005  # typical brine density g/cm3
+    MW_brine = 18.15   # slightly saline brine MW
+
+    print(f"Conditions: T = {T_celsius}°C, rho_brine = {rho_brine} g/cm3")
+    print()
+    print(f"{'Gas':>6} {'MW':>8} {'V_phi':>8} {'MW/V_phi':>8} {'x_gas':>8} {'rho_sat':>10} {'delta%':>8} {'Effect':>12}")
+    print("-" * 80)
+
+    for gas_name, vphi_fn, mw in [
+        ('CO2', vphi_co2, MW['CO2']),
+        ('CH4', vphi_ch4_placeholder, MW['CH4']),
+        ('H2S', vphi_h2s_placeholder, MW['H2S']),
+    ]:
+        x_gas = 0.02  # typical dissolved mole fraction
+        vphi_val = vphi_fn(T_celsius)
+        mw_over_vphi = mw / vphi_val
+        rho_sat = garcia_density(rho_brine, T_celsius, x_gas, MW_brine, mw, vphi_fn)
+        delta = (rho_sat - rho_brine) / rho_brine * 100
+        effect = "increases" if delta > 0 else "decreases"
+
+        print(f"{gas_name:>6} {mw:>8.3f} {vphi_val:>8.2f} {mw_over_vphi:>8.3f} {x_gas:>8.4f} {rho_sat:>10.5f} {delta:>8.3f} {effect:>12}")
+
+    print()
+    print("Note: CH4 and H2S V_phi values are PLACEHOLDERS. Real coefficients needed.")
+    print("Note: MW/V_phi > ~1.0 → increases density; MW/V_phi < ~1.0 → decreases density")
+    print()
+
+    # Mixed gas demonstration
+    print("--- Mixed gas example ---")
+    print("If multiple gases dissolved, use mole-fraction-weighted V_phi and MW:")
+
+    x_co2 = 0.015
+    x_ch4 = 0.005
+    x_total = x_co2 + x_ch4
+    y_co2 = x_co2 / x_total  # fraction of CO2 among dissolved gases
+    y_ch4 = x_ch4 / x_total
+
+    vphi_eff = y_co2 * vphi_co2(T_celsius) + y_ch4 * vphi_ch4_placeholder(T_celsius)
+    mw_eff = y_co2 * MW['CO2'] + y_ch4 * MW['CH4']
+
+    vphi_eff_fn = lambda T: vphi_eff  # constant for this demo
+    rho_mixed = garcia_density(rho_brine, T_celsius, x_total, MW_brine, mw_eff, vphi_eff_fn)
+    delta_mixed = (rho_mixed - rho_brine) / rho_brine * 100
+
+    print(f"  CO2: x={x_co2}, CH4: x={x_ch4}, total: x={x_total}")
+    print(f"  V_phi_eff = {vphi_eff:.2f} cm3/mol, MW_eff = {mw_eff:.3f} g/mol")
+    print(f"  rho_mixed = {rho_mixed:.5f} g/cm3 (delta = {delta_mixed:.3f}%)")
+
+
+def print_unified_api_design():
+    """Print the proposed unified API design."""
+    print()
+    print("=" * 80)
+    print("PROPOSED UNIFIED API DESIGN")
+    print("=" * 80)
+    print("""
+Option A: Single function approach
+-----------------------------------
+brine.brine_density(
+    p,              # Pressure (psia)
+    degf,           # Temperature (°F)
+    wt=0,           # Salt wt% (0-100)
+    x_ch4=0,        # Dissolved CH4 mole fraction in brine
+    x_co2=0,        # Dissolved CO2 mole fraction in brine
+    x_h2s=0,        # Dissolved H2S mole fraction in brine
+    x_n2=0,         # Dissolved N2 mole fraction in brine
+    x_h2=0,         # Dissolved H2 mole fraction in brine
+    x_c2h6=0,       # Dissolved C2H6 mole fraction in brine
+    x_c3h8=0,       # Dissolved C3H8 mole fraction in brine
+) -> float          # Density in g/cm3
+
+brine.brine_viscosity(
+    p, degf, wt=0,
+    x_ch4=0, x_co2=0, x_h2s=0, x_n2=0, x_h2=0, x_c2h6=0, x_c3h8=0,
+) -> float          # Viscosity in cP
+
+brine.brine_properties(
+    p, degf, wt=0,
+    x_ch4=0, x_co2=0, x_h2s=0, x_n2=0, x_h2=0, x_c2h6=0, x_c3h8=0,
+) -> dict           # All properties: density, viscosity, Bw, compressibility
+
+Notes:
+- Gas-free brine density from Spivey (shared engine, already exists)
+- Dissolved gas density correction from Garcia Eq 18 with per-gas V_phi(T)
+- Base viscosity from Mao-Duan (already exists)
+- Viscosity correction: multiplicative per-gas (Islam-Carlson form)
+- Dissolved gas mole fractions come from the CALLER (could be from
+  Spycher-Pruess for CO2, Spivey for CH4, future Soreide-Whitson for
+  multi-gas, or any external source)
+- This DECOUPLES density/viscosity from solubility models
+
+Option B: Keep existing APIs, add unified function alongside
+------------------------------------------------------------
+- Keep brine_props() as-is (backward compatible, CH4-focused)
+- Keep CO2_Brine_Mixture as-is (backward compatible, CO2-focused)
+- Add new unified functions that accept arbitrary compositions
+- Existing functions could internally call the unified functions
+  once validated
+
+RECOMMENDED: Option B for initial implementation
+  - No breaking changes
+  - Allows side-by-side validation
+  - Migrate users gradually
+  - Eventually deprecate old functions
+""")
+
+
+if __name__ == '__main__':
+    # Run all explorations
+    extract_spivey_ch4_vphi()
+    compare_garcia_co2_implementations()
+    compare_viscosity_corrections()
+    demonstrate_garcia_generalization()
+    print_unified_api_design()

pyrestoolbox/tests/test_viscosity_scaling.py

+"""
+Exploratory analysis: Can dissolved gas viscosity effects be scaled from CO2?
+
+Investigates whether the ratio of free-gas viscosities, molecular weights,
+apparent molar volumes, or "dissolved densities" (MW/V_phi) can predict
+the viscosity effect of dissolving different gases in water/brine.
+
+Known data:
+  - CO2: Islam-Carlson (2012): mu = mu_brine * (1 + 4.65 * x_CO2^1.0134)
+  - CH4: Kumagai & Yokoyama (1998): ~1-6% decrease at saturation conditions
+  - H2S: no published correction
+  - N2, H2: no published correction
+"""
+
+import sys
+import os
+import numpy as np
+
+# Add paths
+sys.path.insert(0, os.path.join(os.path.dirname(__file__), '..'))
+sys.path.insert(0, '/home/mark/projects/pyResToolbox/Garcia_code')
+
+# ============================================================================
+# 1. Islam-Carlson CO2 viscosity effect
+# ============================================================================
+def islam_carlson_factor(x_co2):
+    """Viscosity multiplier from Islam-Carlson (2012)."""
+    if x_co2 <= 0:
+        return 1.0
+    return 1.0 + 4.65 * x_co2 ** 1.0134
+
+print("=" * 70)
+print("1. ISLAM-CARLSON CO2 VISCOSITY EFFECT")
+print("=" * 70)
+print(f"\n  {'x_CO2':>8} {'Factor':>8} {'%change':>8}")
+print("  " + "-" * 28)
+for x in [0.005, 0.01, 0.02, 0.03, 0.04, 0.05]:
+    f = islam_carlson_factor(x)
+    print(f"  {x:8.3f} {f:8.4f} {(f-1)*100:+7.2f}%")
+
+# Note: at x_CO2 = 0.01, effect is ~4.4% INCREASE
+# The exponent ~1.01 means it's nearly linear: ~4.65 * x_CO2
+
+# ============================================================================
+# 2. Free gas viscosities at reservoir conditions
+# ============================================================================
+print("\n" + "=" * 70)
+print("2. FREE GAS VISCOSITIES AT RESERVOIR CONDITIONS")
+print("=" * 70)
+
+# Use pyResToolbox gas_ug for hydrocarbon gases
+# For non-HC gases, use Lee-Gonzalez-Eakin with appropriate MW and Tc/Pc
+from pyrestoolbox import gas as rtb_gas
+
+def lee_gonzalez_eakin(T_R, rho_gcc, MW):
+    """Lee-Gonzalez-Eakin gas viscosity correlation.
+    T_R: temperature in Rankine
+    rho_gcc: gas density in g/cc
+    MW: molecular weight
+    Returns: viscosity in cP
+    """
+    K = ((9.4 + 0.02 * MW) * T_R**1.5) / (209 + 19 * MW + T_R)
+    X = 3.5 + 986.0 / T_R + 0.01 * MW
+    Y = 2.4 - 0.2 * X
+    return 1e-4 * K * np.exp(X * rho_gcc**Y)
+
+
+# Conditions: 150°F (610 R, 338.7 K), 3000 psia (20.68 MPa)
+T_F = 150.0
+P_psia = 3000.0
+T_K = (T_F - 32) * 5/9 + 273.15
+T_R = T_F + 459.67
+P_MPa = P_psia / 145.038
+
+print(f"\n  Conditions: T = {T_F}°F ({T_K:.1f} K), P = {P_psia} psia ({P_MPa:.1f} MPa)")
+
+# For pure CH4 - use pyResToolbox
+try:
+    ug_ch4 = rtb_gas.gas_ug(p=P_psia, sg=0.5537, degf=T_F)  # CH4 sg ~0.5537
+    print(f"\n  CH4 (pyResToolbox, sg=0.554): {float(ug_ch4):.4f} cP")
+except Exception as e:
+    print(f"\n  CH4 via pyResToolbox failed: {e}")
+    ug_ch4 = 0.018  # fallback
+
+# For various pure gases, compute approximate viscosity via LGE
+# Need density - use ideal gas as rough estimate, then correct
+R_gcc = 8.314  # J/(mol*K)
+
+gas_props = {
+    'CH4':    {'MW': 16.043, 'Tc_K': 190.56, 'Pc_MPa': 4.599},
+    'CO2':    {'MW': 44.010, 'Tc_K': 304.13, 'Pc_MPa': 7.377},
+    'H2S':    {'MW': 34.081, 'Tc_K': 373.21, 'Pc_MPa': 8.937},
+    'N2':     {'MW': 28.013, 'Tc_K': 126.19, 'Pc_MPa': 3.396},
+    'H2':     {'MW': 2.016,  'Tc_K': 33.19,  'Pc_MPa': 1.313},
+    'C2H6':   {'MW': 30.069, 'Tc_K': 305.32, 'Pc_MPa': 4.872},
+    'C3H8':   {'MW': 44.096, 'Tc_K': 369.83, 'Pc_MPa': 4.248},
+}
+
+# Simple gas density estimate (ideal gas with compressibility guess)
+def gas_density_approx(MW, T_K, P_MPa, Tc_K, Pc_MPa):
+    """Approximate gas density using Pitzer correlation for Z."""
+    Tr = T_K / Tc_K
+    Pr = P_MPa / Pc_MPa
+    # Simple Z estimate (van der Waals like)
+    if Tr > 1.0:
+        Z = 1.0 + (0.083 - 0.422/Tr**1.6) * Pr / Tr
+    else:
+        Z = 0.3  # liquid-like
+    Z = max(Z, 0.1)
+    rho = P_MPa * MW / (Z * 8.314e-3 * T_K)  # kg/m3
+    return rho / 1000.0, Z  # g/cc, Z
+
+print(f"\n  {'Gas':>8} {'MW':>6} {'Tr':>6} {'Pr':>6} {'Z':>6} {'rho(g/cc)':>10} {'ug_LGE(cP)':>12}")
+print("  " + "-" * 65)
+
+gas_viscosities = {}
+for gas, props in gas_props.items():
+    MW = props['MW']
+    Tc = props['Tc_K']
+    Pc = props['Pc_MPa']
+    Tr = T_K / Tc
+    Pr = P_MPa / Pc
+    rho_gcc, Z = gas_density_approx(MW, T_K, P_MPa, Tc, Pc)
+
+    # LGE viscosity
+    ug = lee_gonzalez_eakin(T_R, rho_gcc, MW)
+    gas_viscosities[gas] = ug
+
+    note = ""
+    if Tr < 1.0:
+        note = " ** SUBCRITICAL - liquid phase **"
+
+    print(f"  {gas:>8} {MW:6.2f} {Tr:6.2f} {Pr:6.2f} {Z:6.3f} {rho_gcc:10.4f} {ug:12.5f}{note}")
+
+print("\n  Note: LGE is calibrated for HC gases. Values for CO2, H2S, N2, H2 are approximate.")
+print("  CO2 and H2S are near/below Tc at 150°F — phase behavior complicates things.")
+
+# ============================================================================
+# 3. V_phi (apparent molar volume) and "dissolved density"
+# ============================================================================
+print("\n" + "=" * 70)
+print("3. APPARENT MOLAR VOLUMES AND 'DISSOLVED DENSITY' AT RESERVOIR CONDITIONS")
+print("=" * 70)
+
+from plyasunov_model import V_phi, gas_mw
+from garcia_mixing import density_change_pct
+
+# At our reservoir conditions
+print(f"\n  Conditions: T = {T_K:.1f} K, P = {P_MPa:.1f} MPa")
+print(f"\n  {'Gas':>8} {'MW':>6} {'V_phi':>8} {'MW/Vphi':>8} {'drho%':>8} {'Direction':>10}")
+print("  " + "-" * 56)
+
+gases = ['H2', 'N2', 'CH4', 'CO2', 'C2H6', 'C3H8', 'NC4H10', 'H2S']
+vphi_data = {}
+
+for gas_name in gases:
+    MW = gas_mw(gas_name)
+    vphi = V_phi(gas_name, T_K, P_MPa)
+    dissolved_density = MW / vphi  # g/cm3 "effective density" of dissolved gas
+
+    _, _, drho_pct = density_change_pct(gas_name, 0.02, T_K, P_MPa)
+
+    direction = "INCREASES" if drho_pct > 0 else "decreases"
+
+    vphi_data[gas_name] = {'MW': MW, 'vphi': vphi, 'dissolved_rho': dissolved_density, 'drho_pct': drho_pct}
+    print(f"  {gas_name:>8} {MW:6.2f} {vphi:8.2f} {dissolved_density:8.4f} {drho_pct:+7.3f}% {direction:>10}")
+
+# Water density for reference
+rho_water = 1000 * 18.015 / (18.015 * 1000 / 1000)  # ~1 g/cc at ambient, but let's get actual
+from water_properties import rho_w, MW_WATER
+rho_w_val = rho_w(T_K, P_MPa) / 1000  # g/cc
+V1_star = MW_WATER / (rho_w(T_K, P_MPa) / 1000)  # cm3/mol... wait
+V1_star = MW_WATER * 1000 / rho_w(T_K, P_MPa)  # cm3/mol
+water_dissolved_density = MW_WATER / V1_star
+print(f"\n  Water: MW={MW_WATER:.3f}, V*1={V1_star:.2f} cm3/mol, MW/V*1={water_dissolved_density:.4f} g/cm3")
+print(f"  Water density at conditions: {rho_w(T_K, P_MPa):.2f} kg/m3 = {rho_w_val:.4f} g/cm3")
+
+# ============================================================================
+# 4. KEY ANALYSIS: Can we predict viscosity effect direction and magnitude?
+# ============================================================================
+print("\n" + "=" * 70)
+print("4. ANALYSIS: SCALING PARAMETERS FOR VISCOSITY EFFECTS")
+print("=" * 70)
+
+print("""
+KNOWN VISCOSITY EFFECTS:
+  CO2: INCREASES water viscosity by ~4.65% per mol% (Islam-Carlson)
+  CH4: DECREASES water viscosity by ~1-6% at saturation (Kumagai-Yokoyama)
+  H2S: Unknown (but chemically similar to H2O - hydrogen bonding)
+  N2, H2, C2H6, C3H8: Unknown
+
+HYPOTHESIS A: Scale by free-gas viscosity ratio (ug_gas / ug_CO2)
+""")
+
+ug_co2 = gas_viscosities.get('CO2', 0.03)
+print(f"  CO2 Islam-Carlson coefficient: a = 4.65")
+print(f"  CO2 free-gas viscosity (LGE approx): {ug_co2:.5f} cP")
+print(f"\n  {'Gas':>8} {'ug(cP)':>10} {'ug/ug_CO2':>10} {'Scaled a':>10} {'%eff @ x=0.02':>14}")
+print("  " + "-" * 58)
+
+for gas in ['CH4', 'CO2', 'H2S', 'N2', 'H2', 'C2H6', 'C3H8']:
+    ug = gas_viscosities.get(gas, 0)
+    ratio = ug / ug_co2 if ug_co2 > 0 else 0
+    scaled_a = 4.65 * ratio
+    eff = scaled_a * 0.02  # approximate effect at x=0.02
+    print(f"  {gas:>8} {ug:10.5f} {ratio:10.4f} {scaled_a:10.3f} {eff*100:+13.2f}%")
+
+print("""
+  PROBLEM: This predicts ALL gases INCREASE viscosity (all a > 0).
+  But CH4 is known to DECREASE water viscosity.
+  Free-gas viscosity ratios cannot predict the SIGN of the effect.
+""")
+
+print("HYPOTHESIS B: Scale by 'dissolved density' ratio (MW/V_phi)")
+print("  Idea: Heavy molecules in small volumes increase viscosity;")
+print("  light molecules in large volumes decrease it.")
+print(f"\n  Water 'density': MW/V1* = {water_dissolved_density:.4f} g/cm3")
+print(f"  If MW_gas/V_phi > MW_water/V1*, dissolved gas is 'heavier' → increases viscosity")
+print(f"  If MW_gas/V_phi < MW_water/V1*, dissolved gas is 'lighter' → decreases viscosity")
+
+print(f"\n  {'Gas':>8} {'MW/Vphi':>8} {'vs H2O':>8} {'Density':>10} {'Visc':>10}")
+print("  " + "-" * 50)
+
+for gas_name in gases:
+    d = vphi_data[gas_name]
+    ratio = d['dissolved_rho'] / water_dissolved_density
+    density_dir = "INCREASES" if d['drho_pct'] > 0 else "decreases"
+    visc_dir = "INCREASES" if d['dissolved_rho'] > water_dissolved_density else "decreases"
+    print(f"  {gas_name:>8} {d['dissolved_rho']:8.4f} {ratio:8.4f} {density_dir:>10} {visc_dir:>10}")
+
+print(f"""
+  The MW/V_phi ratio correctly predicts:
+    - CO2 increases (MW/Vphi = {vphi_data['CO2']['dissolved_rho']:.4f} > water {water_dissolved_density:.4f}) ✓
+    - CH4 decreases (MW/Vphi = {vphi_data['CH4']['dissolved_rho']:.4f} < water {water_dissolved_density:.4f}) ✓
+    - H2S: {'increases' if vphi_data['H2S']['dissolved_rho'] > water_dissolved_density else 'decreases'}
+
+  But the sign correlation is IDENTICAL to the density effect direction.
+  This is expected: both density and viscosity effects stem from the same
+  physics (how the dissolved molecule fits in the water structure).
+""")
+
+# ============================================================================
+# 5. PROPOSED APPROACH: Viscosity-density analogy
+# ============================================================================
+print("=" * 70)
+print("5. PROPOSED APPROACHES")
+print("=" * 70)
+
+print("""
+APPROACH 1: Density-viscosity proportionality
+
+  Observation: For CO2 at x=0.02:
+    - Density increase: ~1.1%
+    - Viscosity increase: ~9.3% (Islam-Carlson)
+    - Ratio: viscosity_effect / density_effect ≈ 8.5
+
+  This suggests viscosity is ~8-9x more sensitive than density.
+
+  For other gases, scale:
+    a_gas = a_CO2 * (drho%_gas / drho%_CO2)
+
+  This gives the right SIGN (negative for CH4, positive for CO2)
+  and uses the density model we already have.
+""")
+
+# Compute this
+drho_co2 = vphi_data['CO2']['drho_pct']
+a_co2 = 4.65
+ic_factor_002 = islam_carlson_factor(0.02) - 1  # fractional viscosity change at x=0.02
+visc_to_density_ratio = ic_factor_002 / (drho_co2 / 100)
+
+print(f"  CO2 at x=0.02: drho = {drho_co2:+.4f}%, dmu = {ic_factor_002*100:+.4f}%")
+print(f"  Viscosity/density sensitivity ratio: {visc_to_density_ratio:.2f}")
+
+print(f"\n  {'Gas':>8} {'drho%':>8} {'Predicted dmu%':>15} {'Predicted a':>12}")
+print("  " + "-" * 48)
+
+for gas_name in gases:
+    d = vphi_data[gas_name]
+    # Scale: dmu/mu = visc_to_density_ratio * drho/rho
+    predicted_dmu_pct = visc_to_density_ratio * d['drho_pct']
+    # Back out 'a' coefficient: dmu/mu ≈ a * x at small x
+    predicted_a = predicted_dmu_pct / (0.02 * 100) if abs(d['drho_pct']) > 1e-6 else 0
+
+    print(f"  {gas_name:>8} {d['drho_pct']:+7.3f}% {predicted_dmu_pct:+14.3f}% {predicted_a:+11.3f}")
+
+print("""
+APPROACH 2: Direct V_phi-based correction (more physical)
+
+  Islam-Carlson is: mu = mu_brine * (1 + a * x^b)
+  where a=4.65, b=1.0134 for CO2.
+
+  Generalize: mu = mu_brine * Product_over_gases(1 + a_i * x_i^b)
+  where a_i scales with how "disruptive" the dissolved molecule is.
+
+  Option A: a_i = a_CO2 * (drho_i% / drho_CO2%)
+    Uses density ratio as proxy. Simple, uses existing model.
+
+  Option B: a_i = a_CO2 * (MW_i/Vphi_i - MW_H2O/V1*) / (MW_CO2/Vphi_CO2 - MW_H2O/V1*)
+    Uses "excess dissolved density" relative to water.
+    More physical basis.
+""")
+
+# Compute Option B
+delta_CO2 = vphi_data['CO2']['dissolved_rho'] - water_dissolved_density
+
+print(f"\n  Option B scaling:")
+print(f"  CO2 'excess dissolved density': {delta_CO2:+.4f} g/cm3")
+print(f"\n  {'Gas':>8} {'MW/Vphi':>8} {'delta':>8} {'delta/dCO2':>10} {'a_i':>8}")
+print("  " + "-" * 48)
+
+for gas_name in gases:
+    d = vphi_data[gas_name]
+    delta = d['dissolved_rho'] - water_dissolved_density
+    ratio = delta / delta_CO2 if abs(delta_CO2) > 1e-10 else 0
+    a_i = a_co2 * ratio
+    print(f"  {gas_name:>8} {d['dissolved_rho']:8.4f} {delta:+8.4f} {ratio:+9.4f} {a_i:+8.3f}")
+
+# ============================================================================
+# 6. CROSS-CHECK: CH4 predicted vs literature
+# ============================================================================
+print("\n" + "=" * 70)
+print("6. CROSS-CHECK: CH4 PREDICTED EFFECT vs LITERATURE")
+print("=" * 70)
+
+# Kumagai & Yokoyama (1998): CH4 decreases viscosity by ~1-6% at saturation
+# Typical CH4 saturation mole fraction in water at reservoir conditions: 0.001-0.003
+# At 150°C, 30 MPa: x_CH4 ~ 0.003 (very rough)
+
+print("""
+  Literature: Kumagai & Yokoyama (1998) report CH4 DECREASES water
+  viscosity by 1-6% depending on T, P, and CH4 concentration.
+
+  Typical CH4 saturation in water: x_CH4 ~ 0.001-0.003 at reservoir conditions
+""")
+
+# Using our density-based scaling:
+for x_ch4 in [0.001, 0.002, 0.003, 0.005, 0.01]:
+    drho_ch4_x = vphi_data['CH4']['drho_pct'] * (x_ch4 / 0.02)  # rough linear scale
+    predicted_dmu = visc_to_density_ratio * drho_ch4_x / 100
+    print(f"  x_CH4={x_ch4:.3f}: predicted dmu = {predicted_dmu*100:+.2f}%")
+
+print("""
+  At x_CH4=0.002-0.003, we predict ~-1.5% to -2.3% viscosity decrease.
+  Kumagai-Yokoyama report 1-6% decrease over a range of conditions.
+  Order of magnitude is reasonable, though our prediction is conservative.
+
+  The larger effects in K&Y may reflect higher CH4 concentrations at
+  their specific conditions and/or additional physical mechanisms.
+""")
+
+# ============================================================================
+# 7. TEMPERATURE DEPENDENCE
+# ============================================================================
+print("=" * 70)
+print("7. TEMPERATURE DEPENDENCE OF SCALING PARAMETERS")
+print("=" * 70)
+
+print(f"\n  V_phi and MW/V_phi at P=30 MPa, various temperatures:")
+print(f"\n  {'T(°C)':>6}  ", end="")
+for gas_name in ['CO2', 'CH4', 'H2S', 'N2', 'H2']:
+    print(f"  {gas_name:>12}", end="")
+print()
+print("  " + "-" * 70)
+
+for label in ['V_phi', 'MW/Vphi']:
+    for T_C in [25, 50, 100, 150, 200]:
+        T_K_local = T_C + 273.15
+        print(f"  {T_C:>5}°C ", end="")
+        for gas_name in ['CO2', 'CH4', 'H2S', 'N2', 'H2']:
+            MW = gas_mw(gas_name)
+            vphi = V_phi(gas_name, T_K_local, 30.0)
+            if label == 'V_phi':
+                print(f"  {vphi:12.2f}", end="")
+            else:
+                print(f"  {MW/vphi:12.4f}", end="")
+        print(f"  ({label})")
+    print()
+
+# Check: does the scaling ratio stay roughly constant with T?
+print("  Scaling ratio (MW/Vphi_gas - MW/V1*) / (MW/Vphi_CO2 - MW/V1*) at P=30 MPa:")
+print(f"\n  {'T(°C)':>6}", end="")
+for gas_name in ['CH4', 'H2S', 'N2', 'H2']:
+    print(f"  {gas_name:>8}", end="")
+print()
+print("  " + "-" * 42)
+
+for T_C in [25, 50, 100, 150, 200]:
+    T_K_local = T_C + 273.15
+    V1 = MW_WATER * 1000 / rho_w(T_K_local, 30.0)
+    wd = MW_WATER / V1
+
+    vphi_co2 = V_phi('CO2', T_K_local, 30.0)
+    delta_co2_local = gas_mw('CO2') / vphi_co2 - wd
+
+    print(f"  {T_C:>5}°C", end="")
+    for gas_name in ['CH4', 'H2S', 'N2', 'H2']:
+        vphi_g = V_phi(gas_name, T_K_local, 30.0)
+        delta_g = gas_mw(gas_name) / vphi_g - wd
+        ratio = delta_g / delta_co2_local if abs(delta_co2_local) > 1e-10 else 0
+        print(f"  {ratio:+8.3f}", end="")
+    print()
+
+print("""
+  If the ratios are roughly constant with temperature, then a single
+  set of scaling coefficients works at all T. If they vary significantly,
+  we need T-dependent coefficients.
+""")
+
+# ============================================================================
+# 8. SUMMARY AND RECOMMENDATION
+# ============================================================================
+print("=" * 70)
+print("8. SUMMARY AND RECOMMENDATION")
+print("=" * 70)
+
+print("""
+FINDINGS:
+  1. Free-gas viscosity ratio CANNOT work: all ratios are positive,
+     but CH4 (and likely N2, H2) DECREASE water viscosity. The free-gas
+     viscosity has no physical connection to how dissolved molecules
+     affect liquid water structure.
+
+  2. The density effect (which we already compute) is a reliable proxy
+     for the viscosity effect SIGN: gases that increase density also
+     increase viscosity, and vice versa. This is physically expected -
+     both are controlled by how the solute molecule occupies space in
+     the hydrogen bond network.
+
+  3. Two scaling approaches that preserve the correct sign:
+
+     (A) Density-ratio scaling: a_i = a_CO2 * (drho%_i / drho%_CO2)
+         Simple, uses existing density model. Predicts CH4 effect in
+         the right ballpark vs Kumagai-Yokoyama.
+
+     (B) Dissolved-density scaling: uses (MW/Vphi - MW_H2O/V1*) ratios.
+         More physical but equivalent to (A) for small x.
+
+     Both give the same qualitative results. (A) is simpler to implement.
+
+  4. The scaling ratios are reasonably stable across temperature (25-200°C),
+     meaning a single set of coefficients should work.
+
+  5. Recommended viscosity model:
+
+     mu = mu_brine * Product_i(1 + a_i * x_i^b)
+
+     where b ≈ 1.01 (from Islam-Carlson), and:
+       a_CO2 = +4.65 (Islam-Carlson, calibrated)
+       a_i   = a_CO2 * (drho%_gas_i / drho%_CO2) for other gases
+
+     Or equivalently: a_i = 4.65 * (density_effect_i / density_effect_CO2)
+     computed at a reference condition from the Plyasunov/Garcia model.
+
+CAVEATS:
+  - Only calibrated against CO2 (Islam-Carlson) with order-of-magnitude
+    check against CH4 (Kumagai-Yokoyama).
+  - Assumes viscosity scales linearly with density effect, which is an
+    approximation. The true relationship may be nonlinear.
+  - H2S is a special case (hydrogen bonding with water). Its viscosity
+    effect may not follow the density-based scaling.
+  - No validation data exists for C2H6, C3H8, n-C4H10, H2, or N2
+    dissolved in water.
+""")

PKG-INFO

 Metadata-Version: 2.4
 Name: pyrestoolbox
-Version: 2.2.3
+Version: 2.5.0
 Summary: pyResToolbox - A collection of Reservoir Engineering Utilities
@@ -5,0 +5,0 @@ Home-page: https://github.com/mwburgoyne/pyResToolbox

pyrestoolbox.egg-info/PKG-INFO

 Metadata-Version: 2.4
 Name: pyrestoolbox
-Version: 2.2.3
+Version: 2.5.0
 Summary: pyResToolbox - A collection of Reservoir Engineering Utilities
@@ -5,0 +5,0 @@ Home-page: https://github.com/mwburgoyne/pyResToolbox

pyrestoolbox.egg-info/SOURCES.txt

@@ -15,2 +15,4 @@ LICENSE
 pyrestoolbox/brine/__init__.py
+pyrestoolbox/brine/_lib_salting_library.py
+pyrestoolbox/brine/_lib_vle_engine.py
 pyrestoolbox/brine/brine.py
@@ -48,2 +50,6 @@ pyrestoolbox/classes/__init__.py
 pyrestoolbox/oil/oil.py
+pyrestoolbox/plyasunov/__init__.py
+pyrestoolbox/plyasunov/iapws_if97.py
+pyrestoolbox/plyasunov/plyasunov_model.py
+pyrestoolbox/plyasunov/water_properties.py
 pyrestoolbox/shared_fns/__init__.py
@@ -61,3 +67,5 @@ pyrestoolbox/shared_fns/shared_fns.py
 pyrestoolbox/tests/test_simtools.py
+pyrestoolbox/tests/test_unified_brine_design.py
+pyrestoolbox/tests/test_viscosity_scaling.py
 pyrestoolbox/validate/__init__.py
 pyrestoolbox/validate/validate.py

pyrestoolbox/docs/brine.rst

 ===================================
-Brine PVT with differing degrees of methane or CO2 saturation
+Brine PVT with differing degrees of methane, CO2, or multicomponent gas saturation
 ===================================
 
-For Brine:Methane (brine_props) - Calculates Brine properties from modified Spivey Correlation per McCain Petroleum Reservoir Fluid Properties pg 160. Includes effect of user specified salt concentration and degree of methane saturation.
+Three brine property models are available:
+
+**brine_props** — Methane-saturated brine using IAPWS-IF97 freshwater density with Spivey salt correction per McCain Petroleum Reservoir Fluid Properties pg 160. Includes effect of user specified salt concentration and degree of methane saturation.
 Returns tuple of (Bw (rb/stb), Density (sg), viscosity (cP), Compressibility (1/psi), Rw GOR (scf/stb))
 
-For Brine:CO2 - Returns a class object with calculated CO2 saturated brine property attributes with the following methods
+**CO2_Brine_Mixture** — CO2-saturated brine via Spycher-Pruess mutual solubility model. Returns a class object with calculated CO2 saturated brine property attributes.
 
-.. list-table:: Method employed for different calculations
+**SoreideWhitson** — Multicomponent gas-saturated brine via Soreide-Whitson (1992) VLE model. Supports mixtures of C1, C2, C3, nC4, CO2, H2S, N2 and H2 in fresh or saline water.
+
+.. list-table:: Method employed for different calculations (CO2_Brine_Mixture)
    :widths: 30 40
@@ -18,10 +22,12 @@ :header-rows: 1
      - Spycher & Pruess (2010), modified SRK Cubic EOS method
-   * - Pure Brine Density
+   * - Pure Water Density
+     - IAPWS-IF97 Region 1 (international reference standard)
+   * - Brine Salinity Correction
      - Spivey et al. (modified), per "Petroleum Reservoir Fluid Property Correlations", (McCain, Spivey & Lenn: Chapter 4)
    * - CO2 Corrected Brine Density
-     - Molar volume of dissolved CO2 estimated with Garcia (2001) equation, used with xCO2 calculated from Spycher & Pruess, and CO2-free brine density from Spivey et al to calculate insitu density
+     - Molar volume of dissolved CO2 estimated with Garcia (2001) equation, used with xCO2 calculated from Spycher & Pruess, and CO2-free brine density to calculate insitu density
    * - Pure Brine viscosity
      - Mao-Duan (2009).
    * - CO2 Corrected Brine Viscosity
-     - "Viscosity Models and Effects of Dissolved CO2", Islam-Carlson (2012) to adjust the pure brine viscosity for xCO2 calculated from Spycher & Pruess.  
+     - "Viscosity Models and Effects of Dissolved CO2", Islam-Carlson (2012) to adjust the pure brine viscosity for xCO2 calculated from Spycher & Pruess.
      
@@ -67,7 +73,7 @@
     >>> print('Rsw:', rsw)
-    Bw: 1.0151710978322923
-    SGw: 0.9950036658248123
-    Visw: 0.4993957925685796
-    Cw: 0.00015465633691558178
-    Rsw: 1.2549011339427625
+    Bw: 1.0152007040056148
+    SGw: 0.9950108179684295
+    Visw: 0.4994004662758671
+    Cw: 0.0001539690974662865
+    Rsw: 1.2540982731813703
 
@@ -156,5 +162,5 @@ pyrestoolbox.brine.CO2_Brine_Mixture
     >>> mix.bw  # Returns [CO2 Saturated, Pure Brine, Freshwater]
-    [1.1085795290443725, 1.0543051245909865, 1.0542061001251017]
+    [1.108578337107381, 1.054302417027164, 1.0542033928155845]
     >>> mix.x  # Returns molar fractions in aqueous phase [xCO2, xH2O]
-    array([0.02431245, 0.95743175])
+    array([0.02431225, 0.95743209])
     
@@ -167,3 +173,3 @@ Usage example for 175 Bara x 85 degC and 0% NaCl brine:
     >>> mix.Rs  # Returns sm3 dissolved CO2 / sm3 Brine
-    24.742923469934272
+    24.743651168969475
 
@@ -177,3 +183,3 @@ pyrestoolbox.brine.make_pvtw_table
 
-Generates a PVTW (water PVT) table over a pressure range using the Spivey correlation (brine_props).
+Generates a PVTW (water PVT) table over a pressure range using brine_props (IAPWS-IF97 freshwater with Spivey salt correction).
 Optionally exports ECLIPSE PVTW keyword file and Excel spreadsheet.
@@ -256,17 +262,226 @@
 
-    SoreideWhitson(**kwargs) -> raises NotImplementedError
+    SoreideWhitson(pres, temp, ppm=0, y_CO2=0, y_H2S=0, y_N2=0, y_H2=0, sg=0.65, metric=True, cw_sat=False) -> class
 
-Placeholder for the Soreide-Whitson (1992) PR-EOS model for gas solubility in water/brine.
+Soreide-Whitson (1992) VLE model for multicomponent gas solubility in water/brine, with Garcia/Plyasunov
+density corrections and calibrated viscosity corrections. Supports gas mixtures containing any combination
+of C1, C2, C3, nC4, CO2, H2S, N2 and H2 in fresh or saline water.
 
-Planned support includes multicomponent gas mixtures (C1, C2, C3, nC4, CO2, H2S, N2, H2)
-solubility in fresh and saline water.
+The hydrocarbon portion of the gas (1 - y_CO2 - y_H2S - y_N2 - y_H2) is automatically split among C1-C4
+based on the gas specific gravity using constrained exponential decay to match the implied HC molecular weight.
 
-Planned return values:
-- Mole fraction of dissolved gas components in aqueous phase
-- Mole fraction of vaporised water in gas phase
-- Water content of gas (stb/mmscf)
-- Gas solubility in water (scf/stb)
+.. list-table:: Method employed for different calculations (SoreideWhitson)
+   :widths: 30 40
+   :header-rows: 1
 
-Reference: Soreide, I. and Whitson, C.H., "Peng-Robinson Predictions for Hydrocarbons,
-CO2, N2, and H2S with Pure Water and NaCl Brine", Fluid Phase Equilibria, 77, 217-240, 1992.
+   * - Property
+     - Calculation method
+   * - Gas-Brine Equilibrium
+     - Soreide-Whitson (1992), Peng-Robinson EOS VLE flash
+   * - Pure Water Density
+     - IAPWS-IF97 Region 1 (international reference standard)
+   * - Brine Salinity Correction
+     - Spivey et al. (modified), per "Petroleum Reservoir Fluid Property Correlations", (McCain, Spivey & Lenn: Chapter 4)
+   * - Gas-Corrected Brine Density
+     - Garcia (2001) Eq. 18 mixing rule with Plyasunov (2019-2021) apparent molar volumes for each dissolved gas
+   * - Pure Brine Viscosity
+     - Mao-Duan (2009)
+   * - Gas-Corrected Brine Viscosity
+     - Per-gas multiplicative corrections: Islam-Carlson (2012) for CO2, Murphy-Gaines (1974) for H2S, Ostermann (SPE-14211, 1985) for CH4
 
+.. list-table:: Inputs
+   :widths: 10 15 40
+   :header-rows: 1
+
+   * - Parameter
+     - Type
+     - Description
+   * - pres
+     - float
+     - Pressure (Bar or psia)
+   * - temp
+     - float
+     - Temperature (deg C or deg F)
+   * - ppm
+     - float
+     - Parts per million (Wt) NaCl equivalent in brine. Default 0
+   * - y_CO2
+     - float
+     - Mole fraction of CO2 in dry gas. Default 0
+   * - y_H2S
+     - float
+     - Mole fraction of H2S in dry gas. Default 0
+   * - y_N2
+     - float
+     - Mole fraction of N2 in dry gas. Default 0
+   * - y_H2
+     - float
+     - Mole fraction of H2 in dry gas. Default 0
+   * - sg
+     - float
+     - Gas specific gravity, used to estimate HC split among C1-C4. Default 0.65
+   * - metric
+     - Boolean
+     - If True, treats input pressure & temperature as metric, otherwise treats as Field units. Default True
+   * - cw_sat
+     - Boolean
+     - If True, will also calculate saturated brine compressibility. Default False
+
+.. list-table:: Results
+   :widths: 10 15 40
+   :header-rows: 1
+
+   * - Class Attribute
+     - Unit
+     - Description
+   * - .x
+     - dict
+     - Dissolved gas mole fractions in aqueous phase, e.g. {'CO2': 0.024, 'CH4': 0.0015}
+   * - .x_total
+     - float
+     - Total dissolved gas mole fraction (sum of all x_i)
+   * - .y
+     - dict
+     - Gas phase compositions (dry basis, normalized)
+   * - .y_H2O
+     - float
+     - Water mole fraction in gas phase
+   * - .water_content
+     - dict
+     - Water content of gas: {'y_H2O': ..., 'stb_mmscf': ..., 'lb_mmscf': ...}
+   * - .bDen
+     - (gm/cm3)
+     - Brine density [Gas Saturated, Gas-Free Brine, Freshwater]
+   * - .bVis
+     - (cP)
+     - Brine viscosity [Gas Saturated, Gas-Free Brine, Freshwater]
+   * - .bVisblty
+     - (1/Bar or 1/psi)
+     - Gas-saturated brine viscosibility
+   * - .bw
+     - (rm3/sm3 or rb/stb)
+     - Brine formation volume factor [Gas Saturated, Gas-Free Brine, Freshwater]
+   * - .Rs
+     - dict (sm3/sm3 or scf/stb)
+     - Per-gas solution ratios, e.g. {'CO2': 15.2, 'CH4': 2.1}
+   * - .Rs_total
+     - (sm3/sm3 or scf/stb)
+     - Total solution gas ratio (sum of all per-gas Rs)
+   * - .Cf_usat
+     - (1/Bar or 1/psi)
+     - Brine undersaturated compressibility
+   * - .Cf_sat
+     - (1/Bar or 1/psi)
+     - Brine saturated compressibility. Requires cw_sat input to be set True to calculate
+   * - .gas_comp
+     - dict
+     - Normalized gas composition used (including estimated HC split from SG)
+   * - .MwBrine
+     - (g/mol)
+     - Molecular weight of gas-free brine
+
+Examples:
+
+Pure CO2 case at 5000 psia x 275 deg F and 3% NaCl brine:
+
+.. code-block:: python
+
+    >>> from pyrestoolbox import brine
+    >>> mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=1.0, metric=False)
+    >>> mix.bDen  # Returns [Gas Saturated, Gas-Free Brine, Freshwater]
+    [0.9733769457162755, 0.968164592979362, 0.9476497407774847]
+    >>> mix.Rs  # Returns per-gas Rs dict (scf/stb)
+    {'CO2': 140.90858294709142}
+    >>> mix.bw  # Returns [Gas Saturated, Gas-Free, Freshwater]
+    [1.0968991160573776, 1.0543023174291248, 1.0542033190822462]
+
+Pure CH4 case (SG=0.554) at 5000 psia x 275 deg F and 3% NaCl brine:
+
+.. code-block:: python
+
+    >>> mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=0, sg=0.554, metric=False)
+    >>> mix.Rs
+    {'CH4': 21.414423540331008}
+    >>> mix.bDen
+    [0.9641137202631425, 0.968164592979362, 0.9476497407774847]
+
+Mixed gas (10% CO2, 5% H2S, SG=0.7) at 200 Bar x 80 degC and 10,000 ppm NaCl:
+
+.. code-block:: python
+
+    >>> mix = brine.SoreideWhitson(pres=200, temp=80, ppm=10000, y_CO2=0.1, y_H2S=0.05, sg=0.7, metric=True)
+    >>> mix.gas_comp  # Estimated gas composition including HC split
+    {'CO2': 0.1, 'H2S': 0.05, 'CH4': 0.8133, 'C2H6': 0.0351, 'C3H8': 0.0015, 'nC4H10': 0.0001}
+    >>> mix.Rs_total  # Total dissolved gas (sm3/sm3)
+    6.32875877743837
+    >>> mix.bDen
+    [0.9854845215724698, 0.9871360082710434, 0.9804911502375318]
+
+Pure CO2 fresh water at 175 Bar x 85 degC with saturated compressibility:
+
+.. code-block:: python
+
+    >>> mix = brine.SoreideWhitson(pres=175, temp=85, ppm=0, y_CO2=1.0, metric=True, cw_sat=True)
+    >>> mix.Rs_total  # sm3 dissolved CO2 / sm3 Brine
+    24.188037633302223
+    >>> mix.Cf_sat
+    0.00016012590421810821
+    >>> mix.water_content
+    {'y_H2O': 0.013982317814779299, 'stb_mmscf': 1.923030543083137, 'lb_mmscf': 673.2379475216799}
+
+References:
+
+- Soreide, I. and Whitson, C.H., "Peng-Robinson Predictions for Hydrocarbons, CO2, N2, and H2S with Pure Water and NaCl Brine", Fluid Phase Equilibria, 77, 217-240, 1992.
+- Garcia, J.E., "Density of Aqueous Solutions of CO2", LBNL Report 49023, 2001.
+- Plyasunov, A.V., "Values of the Apparent Molar Volumes...," Fluid Phase Equilibria, Parts I-IV, 2019-2021.
+- Islam, A.W. and Carlson, E.S. (2012), "Viscosity Models and Effects of Dissolved CO2", Energy & Fuels, 26(8), 5330-5336.
+- Ostermann, R.D. et al. (1985), "The Effect of Dissolved Gas on Reservoir Brine Viscosity", SPE 14211.
+
+Density and Viscosity Calculation Details
+----------------------
+
+**Gas-corrected brine density** uses the Garcia (2001) Eq. 18 mixing rule, which is thermodynamically generic
+and applicable to any dissolved gas or mixture of dissolved gases. For each dissolved gas species, the apparent
+molar volume at infinite dilution V_phi(T,P) is computed using the Plyasunov (2019-2021) A12-infinity model
+with IAPWS-IF97 Region 1 pure water properties. This provides rigorous T- and P-dependent V_phi for all 8
+supported gas species (H2, N2, CH4, CO2, C2H6, C3H8, nC4H10, H2S).
+
+For mixed dissolved gases, mole-fraction-weighted effective V_phi and MW are used:
+
+.. code-block:: text
+
+    V_phi_eff = sum(yi * V_phi_i) / sum(yi)
+    MW_eff    = sum(yi * MW_i) / sum(yi)
+
+    rho = (1 + x_gas * MW_eff / (MW_brine * x_brine)) / (x_gas * V_phi_eff / (MW_brine * x_brine) + 1/rho_brine)
+
+where x_gas is the total dissolved gas mole fraction, x_brine = 1 - x_gas, and rho_brine is the gas-free
+brine density (IAPWS-IF97 freshwater with Spivey salt correction).
+
+Gases with MW/V_phi > 1 (CO2, H2S) increase brine density; gases with MW/V_phi < 1 (CH4, N2, H2, C2H6, C3H8)
+decrease it.
+
+**Gas-corrected brine viscosity** applies multiplicative per-gas corrections to the gas-free Mao-Duan (2009)
+brine viscosity:
+
+.. list-table:: Viscosity correction by gas
+   :widths: 15 40 20
+   :header-rows: 1
+
+   * - Gas
+     - Correction
+     - Reference
+   * - CO2
+     - mu x (1 + 4.65 x xCO2^1.0134)
+     - Islam-Carlson (2012)
+   * - H2S
+     - mu x (1 + 1.5 x xH2S^1.0134)
+     - Murphy-Gaines (1974)
+   * - CH4
+     - mu x max(1.109 - 5.98e-4*T + 1.09e-6*T^2, 1.0), T in degF
+     - Ostermann (SPE-14211, 1985)
+   * - Other
+     - No correction (factor = 1.0)
+     - Negligible effect at typical dissolved concentrations
+
+For gas mixtures, the corrections are applied multiplicatively across all dissolved species.
+

pyrestoolbox/tests/test_brine.py

@@ -190,9 +190,50 @@ #!/usr/bin/env python3
 
-def test_soreide_whitson_not_implemented():
-    """SoreideWhitson should raise NotImplementedError"""
+def test_soreide_whitson_pure_co2():
+    """SoreideWhitson with pure CO2 should produce positive density, viscosity, and Rs"""
+    mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=1.0, metric=False)
+    # Density: gas-saturated should be greater than gas-free (CO2 increases brine density)
+    assert mix.bDen[0] > mix.bDen[1], f"CO2-saturated density {mix.bDen[0]} should exceed gas-free {mix.bDen[1]}"
+    # Viscosity: gas-saturated should be greater than gas-free (CO2 increases viscosity)
+    assert mix.bVis[0] > mix.bVis[1], f"CO2-saturated viscosity {mix.bVis[0]} should exceed gas-free {mix.bVis[1]}"
+    # Rs should be positive
+    assert mix.Rs_total > 0, f"Rs_total should be positive, got {mix.Rs_total}"
+    # Bw should be > 1 (dissolved gas expands brine)
+    assert mix.bw[0] > 1.0, f"Gas-saturated Bw should exceed 1.0, got {mix.bw[0]}"
+    # x should have CO2 key with positive value
+    assert 'CO2' in mix.x, f"x should contain CO2, got {mix.x}"
+    assert mix.x['CO2'] > 0, f"xCO2 should be positive, got {mix.x['CO2']}"
+
+
+def test_soreide_whitson_pure_ch4():
+    """SoreideWhitson with pure CH4 should produce physically reasonable results"""
+    mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=0, sg=0.554, metric=False)
+    # CH4 decreases brine density
+    assert mix.bDen[0] < mix.bDen[1], f"CH4-saturated density {mix.bDen[0]} should be less than gas-free {mix.bDen[1]}"
+    # Rs should be positive
+    assert mix.Rs_total > 0, f"Rs_total should be positive, got {mix.Rs_total}"
+    # Should contain CH4 in dissolved gases
+    assert 'CH4' in mix.x, f"x should contain CH4, got {mix.x}"
+
+
+def test_soreide_whitson_mixed_gas():
+    """SoreideWhitson with mixed gas should return per-gas details"""
+    mix = brine.SoreideWhitson(pres=200, temp=80, ppm=10000, y_CO2=0.1, y_H2S=0.05, sg=0.7, metric=True)
+    # Should have multiple gases in results
+    assert len(mix.x) >= 2, f"Should have at least 2 dissolved gases, got {len(mix.x)}"
+    assert len(mix.Rs) >= 2, f"Should have at least 2 Rs entries, got {len(mix.Rs)}"
+    # Gas composition should be normalized
+    total_y = sum(mix.gas_comp.values())
+    assert abs(total_y - 1.0) < 1e-10, f"Gas composition should sum to 1.0, got {total_y}"
+    # All densities should be positive
+    for d in mix.bDen:
+        assert d > 0, f"All densities should be positive, got {mix.bDen}"
+
+
+def test_soreide_whitson_bad_gas_fractions():
+    """SoreideWhitson should raise ValueError if non-HC fractions exceed 1.0"""
     try:
-        brine.SoreideWhitson()
-        assert False, "Expected NotImplementedError"
-    except NotImplementedError as e:
-        assert "not yet implemented" in str(e).lower(), f"Unexpected error message: {e}"
+        brine.SoreideWhitson(pres=200, temp=80, y_CO2=0.6, y_H2S=0.5, metric=True)
+        assert False, "Expected ValueError for fractions > 1.0"
+    except ValueError:
+        pass
 
@@ -199,0 +240,0 @@

pyrestoolbox/tests/test_doc_examples.py

@@ -380,7 +380,7 @@ #!/usr/bin/env python3
     bw, lsg, visw, cw, rsw = brine.brine_props(p=160, degf=135, wt=1.5, ch4_sat=1.0)
-    assert abs(bw - 1.0151710978322923) / 1.0151710978322923 < RTOL
-    assert abs(lsg - 0.9950036658248123) / 0.9950036658248123 < RTOL
-    assert abs(visw - 0.4993957925685823) / 0.4993957925685823 < RTOL
-    assert abs(cw - 0.00015465242842085548) / 0.00015465242842085548 < RTOL
-    assert abs(rsw - 1.2548933067431638) / 1.2548933067431638 < RTOL
+    assert abs(bw - 1.0152007040056148) / 1.0152007040056148 < RTOL
+    assert abs(lsg - 0.9950108179684295) / 0.9950108179684295 < RTOL
+    assert abs(visw - 0.4994004662758671) / 0.4994004662758671 < RTOL
+    assert abs(cw - 0.0001539690974662865) / 0.0001539690974662865 < RTOL
+    assert abs(rsw - 1.2540982731813703) / 1.2540982731813703 < RTOL
 
@@ -392,3 +392,3 @@ def test_doc_co2_brine_field():
     assert len(mix.bw) == 3
-    assert abs(float(mix.bw[0]) - 1.108579075130017) / 1.108579075130017 < RTOL
+    assert abs(float(mix.bw[0]) - 1.108578337107381) / 1.108578337107381 < RTOL
     assert isinstance(mix.x, np.ndarray)
@@ -400,3 +400,3 @@ assert abs(mix.x[0] - 0.02431225) / 0.02431225 < RTOL
     mix = brine.CO2_Brine_Mixture(pres=175, temp=85)
-    assert abs(mix.Rs - 24.742717860296057) / 24.742717860296057 < RTOL
+    assert abs(mix.Rs - 24.743651168969475) / 24.743651168969475 < RTOL
 
@@ -409,4 +409,33 @@ def test_doc_make_pvtw_table():
         assert key in result, f"Missing key: {key}"
-    assert abs(result['bw_ref'] - 1.0275744009507692) / 1.0275744009507692 < RTOL
+    assert abs(result['bw_ref'] - 1.027589195773527) / 1.027589195773527 < RTOL
 
+def test_doc_sw_pure_co2_field():
+    """brine.rst: SoreideWhitson pure CO2 field units"""
+    mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=1.0, metric=False)
+    assert isinstance(mix.bDen, list) and len(mix.bDen) == 3
+    assert abs(mix.bDen[0] - 0.9733769457162755) / 0.9733769457162755 < RTOL
+    assert abs(mix.Rs['CO2'] - 140.90858294709142) / 140.90858294709142 < RTOL
+    assert abs(mix.bw[0] - 1.0968991160573776) / 1.0968991160573776 < RTOL
+
+def test_doc_sw_pure_ch4_field():
+    """brine.rst: SoreideWhitson pure CH4 field units"""
+    mix = brine.SoreideWhitson(pres=5000, temp=275, ppm=30000, y_CO2=0, sg=0.554, metric=False)
+    assert abs(mix.Rs['CH4'] - 21.414423540331008) / 21.414423540331008 < RTOL
+    assert abs(mix.bDen[0] - 0.9641137202631425) / 0.9641137202631425 < RTOL
+
+def test_doc_sw_mixed_gas_metric():
+    """brine.rst: SoreideWhitson mixed gas metric"""
+    mix = brine.SoreideWhitson(pres=200, temp=80, ppm=10000, y_CO2=0.1, y_H2S=0.05, sg=0.7, metric=True)
+    assert abs(mix.Rs_total - 6.32875877743837) / 6.32875877743837 < RTOL
+    assert abs(mix.bDen[0] - 0.9854845215724698) / 0.9854845215724698 < RTOL
+    assert 'CO2' in mix.gas_comp and 'H2S' in mix.gas_comp and 'CH4' in mix.gas_comp
+
+def test_doc_sw_co2_freshwater_cw_sat():
+    """brine.rst: SoreideWhitson pure CO2 freshwater with Cf_sat"""
+    mix = brine.SoreideWhitson(pres=175, temp=85, ppm=0, y_CO2=1.0, metric=True, cw_sat=True)
+    assert abs(mix.Rs_total - 24.188037633302223) / 24.188037633302223 < RTOL
+    assert abs(mix.Cf_sat - 0.00016012590421810821) / 0.00016012590421810821 < RTOL
+    assert isinstance(mix.water_content, dict)
+    assert abs(mix.water_content['stb_mmscf'] - 1.923030543083137) / 1.923030543083137 < RTOL
+
 # =============================================================================
@@ -413,0 +442,0 @@ # Layer Module Documentation Examples (docs/layer.rst)

pyrestoolbox/tests/test_gas.py

@@ -350,2 +350,29 @@ #!/usr/bin/env python3
 
+def test_bns_z_subcritical_fugacity_selection():
+    """BNS Z-factor must use fugacity-based root selection in the 3-root regime.
+    Pure H2S at 80 degF is sub-critical (Tc_H2S = 672 R = 212 degF).
+    The cubic has 3 real roots between ~330-605 psia; the thermodynamically
+    stable root (lowest fugacity) switches from vapor to liquid at ~340 psia.
+    Without fugacity selection the solver returns the metastable vapor root,
+    producing a discontinuous Z-vs-P curve."""
+    # Reference values from bns.py (SPE-229932-MS reference implementation)
+    # Pressures chosen to span: vapor-only, fugacity crossover, liquid-only
+    ref = {
+        100:  0.9496046861917848,   # Vapor (single root)
+        300:  0.83365711413685,     # Vapor (3 roots, vapor has min fugacity)
+        350:  0.042605680488353545, # Liquid (3 roots, liquid has min fugacity)
+        400:  0.04862660238021889,  # Liquid (3 roots, liquid has min fugacity)
+        600:  0.07256104178906098,  # Liquid (3 roots merging)
+        1000: 0.11977671955327113,  # Liquid (single root)
+    }
+    pressures = np.array(sorted(ref.keys()))
+    z_arr = gas.gas_z(p=pressures, sg=0.8, degf=80, h2s=1.0, zmethod='BNS', cmethod='BNS')
+    for i, p in enumerate(pressures):
+        assert abs(z_arr[i] - ref[p]) < 1e-8, (
+            f"BNS Z at p={p} psia: got {z_arr[i]:.10f}, expected {ref[p]:.10f}"
+        )
+    # Verify monotonic in the liquid region (no discontinuous jump)
+    z_liquid = z_arr[pressures >= 350]
+    assert np.all(np.diff(z_liquid) > 0), "Z should increase monotonically in liquid region"
+
 def test_gas_water_content_salinity():
@@ -352,0 +379,0 @@ """Water content should decrease with salinity"""

setup.cfg

 [metadata]
 name = pyrestoolbox
-version = 2.2.3
+version = 2.5.0
 author = Mark W. Burgoyne
@@ -5,0 +5,0 @@ author_email = [e-mail]

setup.py

@@ -15,3 +15,3 @@ #!/usr/bin/python3
     include_package_data=True,
-    version='2.2.3',  # Ideally should be same as your GitHub release tag version
+    version='2.5.0',  # Ideally should be same as your GitHub release tag version
     packages=find_packages(),
@@ -18,0 +18,0 @@ description='pyResToolbox - A collection of Reservoir Engineering Utilities',

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

Improved metrics

Number of package files

69

13.11%

Lines of code

11289

62.62%

Worsened metrics

Total package byte prevSize

1732526

-6.74%