%% Model description
% This model conceptualizes seafloor hydrothermal fluids as the end
% product of a sequential two-box model comprising low-temperature and 
% high-temperature reactions.
%
% Within each box, fluids and rock fully equilibrate as fresh rock is
% iteratively added to a fluid-rock system until a final total rock mass
% fraction (i.e. fluid:rock ratio) is achieved. K concentrations are 
% controlled by a constant partition coefficient between fluids and rock. 
% Isotopic equilibrium is controlled by an equilibrium isotope 
% fractionation factor. 
%
% These reactions are parameterized by the equations: 
% 1) 
% K_final = (K_rk / D) - (K_rk / D - K_initial).* exp(-D ./ WR);
% where:
%    K_rk is the K concentration of fresh rock
%    D is the equilibrium K partition coefficient
%    K_initial is the K concentration of the initial fluid
%    WR is the integrated water:rock mass ratio
% and 2)
% del41K_final = (del41K_rk + DEL41K) - (del41K_rk + DEL41K - del41K_initial) .* exp(-D ./ WR);
% where:
%    del41K_rk is the del41K value of fresk rock
%    DEL41K is the equilibrium isotope fractionation factor between rock 
%    and aqueous fluid.
%    del41K_initial is the del41K value of the initial fluid
%
% Analyzing these equations, fluid compositions evolve along linear
% trajectories in [K , del41K] space, from the inital fluid composition
% [K_initial , del41K_initial] to a rock-buffered composition
% [K_rk / D , del41K_rk + DEL41K]
% Progress along this path is governed by the term exp(-D ./ WR), which is
% common to both equations.
% Distinct partition coefficients and isotope fractionation are assigned 
% for low- and high-temperature boxes, resulting in an inflection 
% separating low-temperature and high-temperature fluid evolution paths.
 
%% Contraints Imposed by Measured Values
% Fluids begin as seawater and evolve toward a low temperature
% rock-buffered fluid before evolving toward a high-temperature
% rock-buffered fluid. Constraints are imposed by the starting seawater
% composition and the final composition of Lost City vent fluids: 
% [10.46 mmol/kg , 0.05 per mille].  
%
% Concentration of K in seawater (mmol/kg)
K_sw =  10.1; 
% delta 41 K of seawater (per mille) 
del41K_sw = 0.13; 
% Concentration of K in Lost City vent fluid (mmol/kg)
K_LC =  10.46; 
% delta 41 K of Lost City vent fluid (per mille) 
del41K_LC = 0.0462; 
%
%% Calculation of high-temperature fluids passing through Lost City vent fluids
% The composition of the high-temperature rock-buffered fluid is a function of:
%    1) The concentration of K in fresh rock (mmol/kg)
K_rk =  20.0; % 20 mmol/kg gabbro (D-MORB, Gale et al., 2013).
% K_rk =  1.0; % Maximum for 1 mmol/kg peridotite. (Wang and Ionov, EPSL 2023). 
%    2) The delta 41 K value of fresh rock (per mille) 
del41K_rk = -0.45;
%    3) The equilibrium K partition coefficient: D = K_rock/K_fluid
D_highT = 0.1; % This value is from 350 C basalt alteration experiments (Evans et al., 2023).
%    4) The high-temperature equilibrium K isotope fractionation factor:
DEL41K_highT = 0.0; % This value is unknown, but set at 0.0 per mille.

% The composition of the high-temperature rock-buffered fluid.
K_RB_highT = K_rk/D_highT;
del41K_RB_highT = del41K_rk + DEL41K_highT;

% Slopes of high-temperature solutions lines.
m_highT = ((del41K_RB_highT - del41K_LC) ./ (K_RB_highT - K_LC));
% Y intercepts of high-temperature solutions lines.
b_highT = del41K_LC - m_highT .* K_LC;

%% Calculation of low-temperature fluids passing through seawater
% The low-temperature equilibrium K partition coefficient: 
% D = K_rock/K_fluid
D_lowT = 2.7; % This value is from low-temperature experiments (Seyfried et al., 1979).
% The low-temperature equilibrium K isotope fractionation factor:
DEL41K_lowT = [0.0 : 0.1 : 0.5]; % Values ranging from 0 to 0.5 per mille of Liu et al. (2021)

% The composition of the low-temperature rock-buffered fluid.
K_RB_lowT = K_rk/D_lowT;
del41K_RB_lowT = del41K_rk + DEL41K_lowT;

% Slopes of low-temperature solutions lines, constrained by seawater
% composition.
m_lowT = ((del41K_RB_lowT - del41K_sw) ./ (K_RB_lowT - K_sw));
% Y intercepts of high temperature solutions lines.
b_lowT = del41K_sw - m_lowT .* K_sw;

%% Calculation of altered seawater composition
% Compositions of low-temperature altered seawater are intersections
% between low-temperature and high-temperature solutions.
K_asw = (b_lowT - b_highT) ./ (m_highT - m_lowT);
del41K_asw = m_lowT .* K_asw + b_lowT;
del41K_asw_check = m_highT .* K_asw + b_highT;
scatter(K_asw , del41K_asw, 'k')

%% Calculation of low-temperature water:rock ratio
% The expression exp(-D_lowT./WR) scales between 0 and 1. The 
% expression = 1 at the initial fluid composition and 0 at the final 
% rock-buffered fluid composition.

% Relative distance from SW to ASW, scaled from 1 to 0
progress_parameter_lowT = (K_RB_lowT - K_asw) ./ (K_RB_lowT - K_sw);
progress_parameter_lowT_check = (del41K_RB_lowT - del41K_asw) ./ (del41K_RB_lowT - del41K_sw);

% Conversion of relative distance to low-temperature water:rock ratio
WR_lowT = -D_lowT ./ log(progress_parameter_lowT);
WR_lowT_check = -D_lowT ./ log(progress_parameter_lowT_check);
WR_lowT
%% Calculation of high-temperature water:rock ratio
% The expression exp(-D_highT./WR) scales between 0 and 1. The 
% expression = 1 at the initial fluid composition and 0 at the final 
% rock-buffered fluid composition.

% Relative distance from ASW to LC, scaled from 1 to 0
progress_parameter_highT = (K_RB_highT - K_LC) ./ (K_RB_highT - K_asw);
progress_parameter_highT_check = (del41K_RB_highT - del41K_LC) ./ (del41K_RB_highT - del41K_asw);

% Conversion of relative distance to high-temperature water:rock ratio
WR_highT = -D_highT ./ log(progress_parameter_highT);
WR_highT_check = -D_highT ./ log(progress_parameter_highT_check);
WR_highT

%% Plot Solutions
hold on
%% Plot low temperature fluid evolution.
% Water:Rock ranges from Infinite to WR_lowT
plot_intervals_n = 10;
Intervals = [0 : 1/plot_intervals_n : 1];
RW_lowT = 1 ./ WR_lowT ; %Convert to Rock:Water ratios
RW_lowT_plot = Intervals' * RW_lowT;

K_lowT_plot = (K_rk / D_lowT) - (K_rk / D_lowT - K_sw) * exp(-D_lowT .* RW_lowT_plot);
del41K_lowT_plot = ones(1,length(Intervals))' * (del41K_rk + DEL41K_lowT) - ones(1,length(Intervals))' * (del41K_rk + DEL41K_lowT - del41K_sw) .* exp(-D_lowT .* RW_lowT_plot);
plot(cat(1,(ones(1,length(K_asw))' *K_sw)',K_asw),cat(1,(ones(1,length(del41K_asw))' *del41K_sw)',del41K_asw),'Color',[0.6 0.6 0.6],'LineWidth',1)

%% Plot high temperature fluid evolution. 
% Water:Rock ranges from Infinite to 1
X_R = [0:5:50];
X_F = 100 - X_R;
WR = [1,2,3,4,5,10,100];%power(10,[-1:1:2]); %X_F ./ X_R;

K_final = (K_rk / D_highT) - ((K_rk / D_highT - K_asw))'.* exp(-D_highT ./ WR);
del41K_final = (del41K_rk + DEL41K_highT) - ((del41K_rk + DEL41K_highT - del41K_asw))'.* exp(-D_highT ./ WR); 
plot(K_final',del41K_final','k', 'LineWidth',1);
scatter(K_final(end,:)',del41K_final(end,:)',200,'k.')
text(K_final(end,:)' - 1 , del41K_final(end,:)',cellstr(num2str(WR'))+":1", 'HorizontalAlignment', 'right')



%% Plot Lost City Vent Fluid and Seawater
LostCityData = [10.5, 0.03, 0.03; 
                10.5, 0.05, 0.03;
                10.5, 0.06, 0.03;
                10.4, 0.04, 0.02;
                10.5, 0.07, 0.03;
                10.3, 0.03, 0.01;
                10.5, 0.05, 0.02;
                10.5, 0.04, 0.02];

errorbar(K_sw,del41K_sw,0.02,0.02,0.2,0.2,'k','LineStyle','none');
scatter(K_sw,del41K_sw,50,'d','MarkerEdgeColor','k','MarkerFaceColor',[0.6 0.9 1]);
errorbar(LostCityData(:,1),LostCityData(:,2),LostCityData(:,3),LostCityData(:,3),'k','LineStyle','none')
scatter(LostCityData(:,1),LostCityData(:,2),50,'d','MarkerEdgeColor','k','MarkerFaceColor',[1. 0.40 0.40]);
%% Plot low temperature gabbro endpoint.
rectangle('position',[K_rk/D_lowT, del41K_rk, K_rk / 2.5 - K_rk/D_lowT, max(DEL41K_lowT)]);
%text(6.9 , -0.045, "Low-T Gabbro",'Rotation',90);

%% Plot low-temperature peridotite endpoint.
K_peridotite = 1
rectangle('position',[0, del41K_rk, K_peridotite/D_lowT, max(DEL41K_lowT)],'FaceColor',[0.4 0.4 0.4]);
%text(1 , -0.045, "Low-T Peridotite",'Rotation',90);

%% Plot high-temperature peridotite endpoint.
endpoint_minK = K_peridotite/0.2;%5
endpoint_maxK = K_peridotite / D_highT;%10
endpoint_center = 0.5*(endpoint_maxK + endpoint_minK);
endpoint_halflength = endpoint_maxK - endpoint_center;
errorbar(endpoint_center, -0.04, 0 , 0 , endpoint_halflength,endpoint_halflength,'k');
%rectangle('position',[K_peridotite/0.2, del41K_rk + min(DEL41K_highT), K_peridotite / D_highT - K_peridotite/0.2, abs(DEL41K_highT)]);
%text(endpoint_minK , -0.045, ["High-T","Peridotite"]);

%% Adjust axes
xlim([0 , 40])
ylim([-0.05 , 0.15])