using equilibrium partitioning to calculate water concentration to use with the SoLiquidus function from the VBRc

vbrc_katz_with_partitioning.m

addpath(getenv('vbrdir'))
vbr_init

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% equilibrium partitioning
%
% solidii are in terms of volatile concentration in the
% fluid phase. Can use equilibrium partitioning to estimate
%
% mass conservation:
%
% c_bulk = c_s * (1 - phi) + c_f * phi
%
% where:
%   c_bulk : bulk concentration
%   c_s : concentration in solid
%   c_f : concentration in fluid
%   phi : volume fraction of fluid (porosity)
%
% define equilibrium partition coefficient
%   k_d = c_s / c_f
% or, in terms of c_s:
%   c_s = k_d * c_f
%
% substitute c_s, do some algebra:
%   c_bulk = k_d * c_f * (1 - phi) + c_f * phi
%   c_bulk = k_d * c_f + c_f * phi * (1 - k_d )
%
% to arrive at:
%
%   c_f = c_bulk / (k_d + phi * (1 - k_d))
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


P_Pa = 2.2 * 1e9;
k_d = 1e-2;

phi_vals = [0, 0.0001, 0.001, 0.005 0.01];

figure()
for iphi = 1:numel(phi_vals)

    phi = phi_vals(iphi);
    H2O_bulk_PPM = linspace(0,3000, 100);
    H2O_F_PPM =  H2O_bulk_PPM ./ (k_d + phi * (1 - k_d));
    H2O_F_wt = H2O_F_PPM * 1e-4; % convert from PPM to wt percent
    CO2_F_wt = zeros(size(H2O_bulk_PPM));

    Solidus = SoLiquidus(P_Pa, H2O_F_wt, CO2_F_wt, 'katz');

    hold all
    labelstr =  ["phi ", num2str(phi)];
    plot(H2O_bulk_PPM, Solidus.Tsol, 'linewidth', 2, 'displayname',labelstr)

end
xlabel('H2O in solid (PPM)')
ylabel('Solidus (C)')
title('2.2 GPa, phi in fraction not percent')
legend()
set(gca,  "fontsize", 12)