proteus.outgas.lavatmos¶

def read_in_element_fracs_normalized(input_path):
    """constructs a dataframe with number element fractions of each element as computed by lavatmos

    input:
    - inputpath: path to the directory where the element_abundance folder which contains the lavatmos output elements is located
        this file contains element abundances in fastchem format where ej=10^(xj-12)
    - time: time of the current timestep at which the vapourisation is computed

    output:
    - dataframe with normalised element fractions ej/etot

    """
    # read file (skip comment line, split on whitespace)

    df = pd.read_csv(
        input_path,
        comment='#',  # ignore lines starting with #
        sep=r'\s+',
        header=None,  # split on any whitespace
    )

    # make first column the headers and second column the data row
    abundance_dict = dict(zip(df[0], df[1]))
    for key in abundance_dict.keys():
        if abundance_dict[key] != 0.0:
            # lavatmos abundances are scaled by 1e20
            abundance_dict[key] = 10 ** (abundance_dict[key] - 12) / 1e20
        else:
            abundance_dict[key] = 0.0

    # add missing elements
    for e in element_list:
        if e not in abundance_dict.keys():
            abundance_dict[e] = 0.0

    total = sum(abundance_dict.values())

    if abs(total) > 1e-30:
        norm_dict = {k: v / total for k, v in abundance_dict.items()}
    else:
        norm_dict = {k: 0 for k in abundance_dict}

    return norm_dict

def run_lavatmos(
    dirs: dict, config: Config, hf_row: dict, volatile_fracs: dict, first_iter: bool
):
    """

    This function imports the LavAtmos package from `LAVA_DIR` in-process and
    runs the melt-vapour equilibrium solve directly (no subprocess or
    container). The result is written to a csv file for diagnostics; the
    caller (`run_vapourisation`) instead reads the element
    abundances and FastChem output that LavAtmos writes as a side effect.

    """

    # Get paths from environment variables
    try:
        paths = paths_importer(dirs)
    except ValueError as e:
        UpdateStatusfile(dirs, 27)
        raise RuntimeError(str(e))

    # Import lavatmos
    lavatmos_dir = str(paths.lavatmos_dir)
    if lavatmos_dir not in sys.path:
        sys.path.append(lavatmos_dir)
    try:
        import lavatmos3
    except ImportError as e:
        UpdateStatusfile(dirs, 27)
        raise RuntimeError(
            f"Failed to import LavAtmos from '{paths.lavatmos_dir}'. "
            'Please check your LAVA_DIR environment variable.'
        ) from e

    Magma = set_magmaproperties(config, hf_row, volatile_fracs, dirs)

    # Import melt composition
    melt_comp_fname = os.path.join(paths.lava_comps, Magma.melt_comp_name + '.csv')
    melt_comp_df = pd.read_csv(melt_comp_fname, names=['spec', 'abund'])
    melt_comp = {}
    for i in melt_comp_df.index:
        melt_comp[melt_comp_df['spec'].loc[i]] = float(melt_comp_df['abund'].loc[i])

    # Pressure used for melt activities [bar]
    P_melt = config.outgas.lavatmos.P_melt

    # Guess for fO2 [bar]
    fO2_initial_guess = 10 ** hf_row['fO2_vapourise_derived']

    # Lavatmos tolerance
    xatol = config.outgas.lavatmos.xatol

    system = lavatmos3.melt_vapor_system(paths)
    lavatmos_output = system.vaporise(
        Magma.T_surf,
        Magma.P_volatile,
        melt_comp,
        volatile_fracs,
        melt_fraction=Magma.melt_fraction,
        P_melt=P_melt,
        fO2_initial_guess=fO2_initial_guess,
        fO2_tries_from_last=bool(not first_iter),
        xatol=xatol,
    )

    # Save results
    output_name = f'{Magma.run_name}.csv'
    lavatmos_output.to_csv(os.path.join(paths.output_dir, output_name))

def run_vapourisation(dirs: dict, config: Config, hf_row: dict, first_iter: bool):
    """

    This function runs the ThermoEngine module LavAtmos. Vapourisation of
    refractory species from the melt surface is computed from the magma
    temperature and the volatile surface pressure. This is a phase change out
    of the melt, distinct from the volatile outgassing (solubility) step that
    runs before it.

    Parameters:
        dirs : dict
            Dictionary of directories
        config : Config
            Configuration object
        hf_row : dict
            Dictionary of helpfile variables, at this iteration only

    """
    paths = paths_importer(dirs)

    # set element fractions in atmosphere for lavatmos run.
    input_eles = ['H', 'C', 'N', 'S', 'O'] + noble_gases
    hf_row['M_vaps'] = 0.0

    # lavatmos takes in the abundance fractions of element not mass fractions so divide by atomic number
    molfracs = {}
    nfrac = {'P': 0.0}
    total_mols = 0.0
    for e in input_eles:
        molfracs[e] = hf_row[e + '_kg_atm'] / species_lib[e].weight
        total_mols += molfracs[e]
    for e in input_eles:
        if total_mols > 0:
            if e == 'H':
                nfrac[e] = max(molfracs[e] / total_mols, H_NUMFRAC_FLOOR)
            else:
                nfrac[e] = molfracs[e] / total_mols
        else:
            if e == 'H':
                nfrac[e] = H_NUMFRAC_FLOOR
            else:
                nfrac[e] = 0.0
    log.debug('volatile element fractions going as input to lavatmos : %s', nfrac)
    log.debug('volatile pressure given to lavatmos : %.4e', hf_row['P_surf'])

    # running lavatmos
    run_lavatmos(dirs, config, hf_row, nfrac, first_iter)

    # convert the element abundances from lavatmos file to element fractions, normalized to unity
    element_fracs = read_in_element_fracs_normalized(paths.element_abundance_output)
    log.debug('element fraction after running lavatmos: %s' % element_fracs)

    # read in boa chemistry from last iteration of fastchem and lavatmos
    output_fc = paths.fastchem3_output
    if os.path.exists(output_fc):
        mmr_path = os.path.join(output_fc, 'boa_chem.dat')
    else:
        UpdateStatusfile(dirs, 27)
        raise RuntimeError('cannot find fastchem output from lavatmos loop!')

    # update abundances in output file for next calliope run
    new_atmos_abundances = pd.read_csv(mmr_path, sep=r'\s+')

    # make sure that the column names are consistent with the rest of the code
    if '#p(bar)' in new_atmos_abundances.columns:
        new_atmos_abundances.rename(columns={'#p(bar)': 'Pbar'}, inplace=True)
    if 'm(u)' in new_atmos_abundances.columns:
        new_atmos_abundances.rename(columns={'m(u)': 'mu'}, inplace=True)

    mu_combined = new_atmos_abundances['mu'][0]

    # Mean particle mass of the combined (volatile + rock-vapour) atmosphere.
    kg_pp_new = mu_combined / const_Nav * 1e-3  # kg per particle
    log.debug('new kg per particle: %.4e' % kg_pp_new)

    # New atmospheric mass from the hydrostatic relation
    area = 4.0 * np.pi * hf_row['R_int'] ** 2
    P_surf_new_Pa = new_atmos_abundances['Pbar'][0] * 1e5  # bar -> Pa
    M_atmo_new = P_surf_new_Pa * area / hf_row['gravity']  # kg

    # Split the new total surface pressure into a volatile part (P_vol) and a
    # rock-vapour part (P_vap). P_vap is the excess of the LavAtmos+FastChem total
    # over the volatile-only pressure that went in; it cannot be negative.
    P_surf_new = new_atmos_abundances['Pbar'][0]
    P_vap_new = P_surf_new - hf_row['P_surf']  # rock-vapour pressure contribution
    if P_vap_new < 0.0:
        log.warning(
            'rock-vapour pressure computed negative (%.3e bar); clamping to zero', P_vap_new
        )
        P_vap_new = 0.0
    log.debug('pressure of vapourised rock species: %.4f bar' % P_vap_new)
    log.debug('pressure of volatiles before vapourisation: %.4f bar' % hf_row['P_surf'])

    hf_row['P_vap'] = P_vap_new
    hf_row['P_vol'] = P_surf_new - P_vap_new  # == old P_surf when P_vap_new >= 0

    # Update the tracked atmospheric MMW of the combined volatile + rock-vapour gas
    hf_row['atm_kg_per_mol'] = mu_combined * 1e-3

    mmw_elements = 0
    for e in element_fracs.keys():
        mmw_elements += element_fracs[e] * species_lib[e].weight

    for vol in gas_list:
        if vol in species_lib.keys():
            vol_key = species_lib[vol].fc_name
        else:
            vol_key = vol
        if vol_key not in new_atmos_abundances.columns:
            # FastChem did not emit this species. This happens for noble gases
            # absent from the FastChem element data (e.g. Kr/Xe are not in the
            # bundled logK set). Keep the species' prior value instead of
            # crashing, and note it.
            log.warning(
                'species %s (%s) not present in FastChem output; keeping prior value',
                vol,
                vol_key,
            )
            continue
        new_pp = new_atmos_abundances[vol_key][0] * new_atmos_abundances['Pbar'][0]
        hf_row[vol + '_vmr'] = new_atmos_abundances[vol_key][0]
        hf_row[vol + '_bar'] = new_pp

        hf_row[vol + '_kg_atm'] = (
            new_atmos_abundances[vol_key][0]
            * M_atmo_new
            * species_lib[vol].weight
            / mu_combined
        )  # kg
        hf_row[vol + '_mol_atm'] = hf_row[vol + '_kg_atm'] / hf_row['atm_kg_per_mol']

        if vol in vap_list:
            hf_row[vol + '_kg_total'] = (
                0.0  # ensures that elements Na and K, si, Ti are not added to the planet mass computations in interior model
            )
            # but saved separately later in M_vaps
            hf_row[vol + '_mol_total'] = 0.0
        elif vol in noble_gases:
            # A noble gas shares its string with the matching element symbol, so
            # its `_kg_total` slot is set by escape. Here, set the mol total from
            # the kg total so that it's counted, but don't overwrite it.
            hf_row[vol + '_mol_total'] = float(hf_row.get(vol + '_kg_total', 0.0)) / (
                species_lib[vol].weight * 1e-3
            )
        else:
            hf_row[vol + '_kg_total'] = (
                hf_row[vol + '_kg_atm'] + hf_row[vol + '_kg_solid'] + hf_row[vol + '_kg_liquid']
            )
            hf_row[vol + '_mol_total'] = (
                hf_row[vol + '_mol_atm']
                + hf_row[vol + '_mol_solid']
                + hf_row[vol + '_mol_liquid']
            )

    Omass_after_vapourise = 0.0
    for e in element_list:
        log.debug('element frac:  %s,  %s', e, element_fracs[e])
        if e in input_eles:
            if e == 'O':
                # O_kg_atm must stay representative of the volatile outgassing
                # step alone, because it sets the escape rate when
                # reservoir=outgas
                log.info(
                    'Oxygen in atmosphere before vapourisation: %.5e', hf_row[e + '_kg_atm']
                )
                # The post-vapourisation O mass is tracked separately, since it
                # feeds the AGNI input VMRs computed with FastChem
                Omass_after_vapourise += (
                    element_fracs[e] * M_atmo_new * species_lib[e].weight / mmw_elements
                )
                log.info(
                    'Oxygen in atmosphere after vapourisation: %.5e', Omass_after_vapourise
                )
            else:
                hf_row[e + '_kg_atm'] = (
                    element_fracs[e] * M_atmo_new * species_lib[e].weight / mmw_elements
                )
        else:
            if e not in gas_list:
                hf_row[e + '_kg_atm'] = (
                    element_fracs[e] * M_atmo_new * species_lib[e].weight / mmw_elements
                )
            hf_row['M_vaps'] += (
                element_fracs[e] * M_atmo_new * species_lib[e].weight / mmw_elements
            )  # don't use hf_row[e + '_kg_atm'] because this only accounts for e.g. Si in atomic gas form not total Si

    hf_row['M_vaps'] += (
        Omass_after_vapourise - hf_row['O_kg_atm']
    )  # add vapourised oxygen mass to the rock-vapour element mass
    log.debug(
        'added oxygen from vapourisation to initial O budget in atmosphere [kg]: %.4e ',
        Omass_after_vapourise - hf_row['O_kg_atm'],
    )
    # Oxygen released by the vapourisation step, held separately from the
    # volatile step's O_kg_atm.
    hf_row['O_vapourised_kg'] = Omass_after_vapourise - hf_row['O_kg_atm']

    # Oxygen fugacity of the combined atmosphere. FastChem re-equilibrated the
    # volatile inventory (nfrac above, built from the *_kg_atm the volatile
    # backend wrote) together with the vapourised rock, so its O2 mixing ratio
    # already carries the oxygen of both sources; the partial pressure is that
    # mixing ratio times the combined surface pressure P_vol + P_vap. The
    # volatile-step O2 partial pressure is therefore NOT added on top of this:
    # doing so would count the volatile oxygen twice.
    vmr_O2 = new_atmos_abundances['O2'][0]
    if vmr_O2 < 1e-20:
        log.debug('O2 mixing ratio %.3e below 1e-20; flooring it for the log', vmr_O2)
        vmr_O2 = 1e-20
    psurf = max(hf_row['P_vol'] + hf_row['P_vap'], 1e-10)  # bar; == FastChem Pbar
    pO2 = vmr_O2 * psurf  # bar
    log.debug('combined O2 partial pressure: %.3e bar (VMR %.3e)', pO2, vmr_O2)
    log10_fO2 = np.log10(pO2)

    # OxygenFugacity(T, fO2_shift=0) returns the absolute buffer value
    # IW(T); the shift relative to that buffer is log10_fO2 - IW(T).
    iw_buffer = OxygenFugacity(model=config.outgas.lavatmos.fO2_buffer_model)
    hf_row['fO2_vapourise_derived'] = log10_fO2
    hf_row['fO2_vapourise_shift_IW_derived'] = log10_fO2 - iw_buffer(hf_row['T_magma'])
    hf_row['P_surf'] = new_atmos_abundances['Pbar'][0]
    hf_row['M_atm'] = M_atmo_new
    # Refresh the volatile-only atmospheric mass against the recombined
    # species masses written above. M_atm now includes the rock vapour, so
    # the two columns diverge from here on.
    hf_row['M_vol_atm'] = sum(float(hf_row.get(s + '_kg_atm', 0.0)) for s in vol_gas_list)

    log.debug(
        'log10 fO2 shift compared to IW buffer: %.6f' % hf_row['fO2_vapourise_shift_IW_derived']
    )

    mask = [hf_row[s + '_vmr'] for s in vap_list]
    for i in np.argsort(mask)[::-1]:
        s = vap_list[i]
        _p = hf_row[s + '_bar']
        _x = hf_row[s + '_vmr']
        _s = '    %-6s     = %-9.2f bar (%.2e VMR)' % (s, _p, _x)
        if _p > 0.01:
            log.info(_s)
        else:
            # don't spam log with species of negligible abundance
            log.debug(_s)

    # print total pressure and mmw
    log.info('    total      = %-9.2f bar' % hf_row['P_surf'])
    log.info('    mmw        = %-9.5f g mol-1' % (hf_row['atm_kg_per_mol'] * 1e3))