proteus.outgas.binodal¶

def apply_binodal_h2(hf_row: dict, config: Config) -> None:
    """Apply Rogers+2025 binodal to partition H2 between atm and mantle.

    Modifies hf_row in place: redistributes H2_kg_atm, H2_kg_liquid,
    and updates H2_bar and H2_vmr accordingly. Reads ``T_magma`` [K],
    ``P_surf`` [bar], ``M_mantle`` [kg], ``R_int`` [m] and ``gravity``
    [m s-2] from the row, and skips the step when ``T_magma``, ``M_mantle``,
    ``R_int`` or ``gravity`` is not positive.

    Does nothing if H2 is not included or has zero total mass. It does not
    read ``outgas.h2_binodal``: the outgassing wrapper calls it only when that
    flag is true.

    Parameters
    ----------
    hf_row : dict
        Current helpfile row with volatile mass and state keys.
    config : Config
        PROTEUS configuration object.
    """
    from zalmoxis.binodal import rogers2025_suppression_weight

    # Check if H2 is included and has meaningful mass
    if not config.outgas.calliope.is_included('H2'):
        return

    H2_kg_total = float(hf_row.get('H2_kg_total', 0.0))
    if H2_kg_total <= 0:
        return

    # Get state variables
    T_magma = float(hf_row.get('T_magma', 0.0))
    P_surf_Pa = float(hf_row.get('P_surf', 0.0)) * 1e5  # bar -> Pa
    M_mantle = float(hf_row.get('M_mantle', 0.0))
    R_int = float(hf_row.get('R_int', 0.0))
    gravity = float(hf_row.get('gravity', 0.0))

    if T_magma <= 0 or M_mantle <= 0 or R_int <= 0 or gravity <= 0:
        log.debug('Skipping binodal: missing state variables')
        return

    # Compute mass fractions for the binary H2-MgSiO3 system.
    # w_H2 is the total H2 mass relative to the mantle + H2 total.
    w_H2 = H2_kg_total / (M_mantle + H2_kg_total)
    w_sil = 1.0 - w_H2

    # Compute suppression weight (sigma)
    sigma = rogers2025_suppression_weight(
        P_Pa=P_surf_Pa,
        T_K=T_magma,
        w_H2=w_H2,
        w_sil=w_sil,
    )

    # Partition H2 between liquid mantle (dissolved) and atmosphere
    H2_kg_liquid_new = sigma * H2_kg_total
    H2_kg_atm_new = (1.0 - sigma) * H2_kg_total

    # Store previous values for logging
    H2_kg_atm_old = float(hf_row.get('H2_kg_atm', 0.0))
    H2_kg_liquid_old = float(hf_row.get('H2_kg_liquid', 0.0))

    # Update hf_row
    hf_row['H2_kg_liquid'] = H2_kg_liquid_new
    hf_row['H2_kg_atm'] = H2_kg_atm_new
    # Solid H2 is zero (H2 does not partition into solid silicate)
    hf_row['H2_kg_solid'] = 0.0

    # The element-level hydrogen totals that the outgassing module wrote
    # before binodal partitioning must follow the H2 just moved. H2 is pure
    # hydrogen by mass (M_H2 = 2 M_H), so the hydrogen mass shifted equals the
    # H2 mass shifted. Without this the *_kg_atm element ratios would still
    # count the relocated hydrogen in the atmosphere.
    dH_atm = H2_kg_atm_new - H2_kg_atm_old
    dH_liquid = H2_kg_liquid_new - H2_kg_liquid_old
    # A round-off negative becomes 0; NaN or a larger negative is kept for the
    # impact check to refuse, since it marks a ledger defect.
    H_atm = float(hf_row.get('H_kg_atm', 0.0)) + dH_atm
    round_off = _H_ATM_ROUND_OFF_REL * max(abs(H2_kg_atm_old), abs(H2_kg_atm_new))
    floor = math.isfinite(round_off) and -round_off <= H_atm < 0.0
    hf_row['H_kg_atm'] = 0.0 if floor else H_atm
    hf_row['H_kg_liquid'] = float(hf_row.get('H_kg_liquid', 0.0)) + dH_liquid

    # Recompute H2 partial pressure from atmospheric mass
    # P_H2 = m_H2 * g / (4 * pi * R^2)
    area = 4.0 * math.pi * R_int**2
    if area > 0 and gravity > 0:
        hf_row['H2_bar'] = H2_kg_atm_new * gravity / area / 1e5  # Pa -> bar
    else:
        hf_row['H2_bar'] = 0.0

    # Recompute total surface pressure (sum of all partial pressures)
    from proteus.utils.constants import gas_list

    P_total = sum(float(hf_row.get(s + '_bar', 0.0)) for s in gas_list)
    hf_row['P_surf'] = P_total

    # Recompute VMRs
    if P_total > 0:
        for s in gas_list:
            hf_row[s + '_vmr'] = float(hf_row.get(s + '_bar', 0.0)) / P_total
    else:
        for s in gas_list:
            hf_row[s + '_vmr'] = 0.0

    # Update molar quantities
    import scipy.constants as const

    N_av = const.Avogadro
    hf_row['H2_mol_liquid'] = H2_kg_liquid_new / _MU_H2 / N_av
    hf_row['H2_mol_atm'] = H2_kg_atm_new / _MU_H2 / N_av
    hf_row['H2_mol_solid'] = 0.0
    hf_row['H2_mol_total'] = H2_kg_total / _MU_H2 / N_av

    # Recompute atmospheric MMW
    _recompute_atm_mmw(hf_row)

    # Note: sigma was evaluated at the pre-redistribution P_surf. A fixed-point
    # iteration (recompute sigma at updated P_surf, re-partition, repeat) would
    # be more self-consistent but the correction is small: the H2 partial
    # pressure shift is typically < 1% of total P_surf, so one pass suffices.
    log.info(
        'Binodal H2: sigma=%.4f, H2_atm %.2e->%.2e kg, '
        'H2_liquid %.2e->%.2e kg (T=%.0f K, P=%.1f bar)',
        sigma,
        H2_kg_atm_old,
        H2_kg_atm_new,
        H2_kg_liquid_old,
        H2_kg_liquid_new,
        T_magma,
        float(hf_row['P_surf']),
    )