% Load datasets

% Load all ACC TE data from Evans et al. (2020) [Li, B, Na, Mg, Mn, Sr, Ba, U]
load("Evans_et_al_2020_ACC_TE_data.mat")

% Inorganic aragonite Ds at 25oC (preferred, min, max)
Darag = [2.60E-3 1.86E-03 3.28E-03]; % Li: average & min/max: (Castillo Alvarez et al. 2025)
Darag(2,:) = [5.4 1.5 14.4].*10^-3; % B: average & min/max: (Holcomb et al., 2016)
Darag(3,:) = [2.7028e-04 3.9811e-06 6.87E-04]; % Na: average: Kitano et al. (1975; @ normal SW [Na] (2950/23*40.078/395000)/(10500/23*40.078/380)), min/max (Füger et al., 2019 (min); Okumura & Kitano, 1986)
Darag(4,:) = [0.8E-3 0.39E-3 1.57E-3]; % Mg: (Castillo Alvarez et al. 2025 mean @ Omega = ~10, range from same study)
Darag(5,:) = [0.86 0.82 0.9]; % Mn (all from Raiswell & Brimblecombe, 1977)
Darag(6,:) = [1.015 0.92 1.15]; % Sr: average (average of Zhong & Mucci, 1989), min: (Zhong & MUcci, 1989), max: (Gaetani & Cohen, 2006)
Darag(7,:) = [1.4 1.3 2.25]; % Ba: average (Gaetani & Cohen, 2006 & Dietzel et al. 2004), min: (Dietzel et al., 2004), max: (Gaetani & Cohen, 2006)
Darag(8,:) = nan(1,3); % U ?? 2 ppm in corals, 3.6 ppm in diagenetic aragonite (Gvirtzman et al., 1973)

% Load coral data compilation
% All mmol/mol except Li, B, Ba
CoralData = readtable('Foster_Evans_Tab_S1.csv');

% Two-stage Rayleigh fractionation model:
% i) ACC precipitated from an open reservoir
% ii) ACC separated from that reservoir, some portion is re-dissolved, and some undergoes solid state transformation

% Concentration of the above TE in seawater (umol or mmol)
SWconcs = [39.5 ; 416 ; 485 ; 52.8 ; 6.5e-1 ; 0.091 ; 0.0378 ; 1.5e-2]; % uM uM mM mM uM mM uM uM
SWCa = 10.3;

% ACC trace element distribution coefficients from Evans et al (2020)
% Li, B, Na, Mg, Mn, Sr, Ba, U
% Averaged over all experiments in seawater with no additives
DACC = mean((TEallpHEx(1:26,:)./1000)./(SWall(1:26,:)./CaAll(1:26).*1000),'omitnan')';
%DACC = (1-0.5.^DACC)./(1-0.5);

% Coral distribution coefficients (COCs)
coralDsCOC = [CoralData.Li_COC_umol_mol_ CoralData.B_COC_umol_mol_ CoralData.Na_COC_mmol_mol_ ...
    CoralData.Mg_COC_mmol_mol_ nan(height(CoralData),1) ...
    CoralData.Sr_COC_mmol_mol_ CoralData.Ba_COC_umol_mol_ nan(height(CoralData),1)];
coralDsCOC = coralDsCOC./(repmat(SWconcs'./SWCa,size(coralDsCOC,1),1))./[1000 1000 1000 1000 1000 1000 1000 1];

% Coral distribution coefficients (fibres)
coralDsFibre = [CoralData.Li_fibre_umol_mol_ CoralData.B_fibre_umol_mol_ CoralData.Na_fibre_mmol_mol_ ...
    CoralData.Mg_fibre_mmol_mol_ nan(height(CoralData),1) ...
    CoralData.Sr_fibre_mmol_mol_ CoralData.Ba_fibre_umol_mol_ nan(height(CoralData),1)];
coralDsFibre = coralDsFibre./(repmat(SWconcs'./SWCa,size(coralDsFibre,1),1))./[1000 1000 1000 1000 1000 1000 1000 1];

% Distribution coefficients into corals
% Dcoral = [0.00249 ; NaN ; 0.000446 ; 0.00087 ; ...
%     3.3 ; 1.05 ; 1.5 ; NaN];
Dcoral = mean(coralDsFibre,'omitnan');
Dcoral(5) = 0.9; % add Mn D, from Foerster et al. (2025)
% DcoralRange = [0.001 ; NaN ; 0.00005 ; 0.00015 ; ...
%     1 ; 0.1 ; 0.3 ; NaN];
DcoralRange = std(coralDsCOC,'omitnan');
DcoralRange(5) = 0.1; %from Foerster et al. (2025)

F = 0.001:0.001:1-0.001;
mgPml = (-2:0.01:1)';
% conversion between miligrams and milimoles accounting for Mg/Ca and H2O
% of ACC
ACCmgPmmol = (DACC(4)*5300/1030)/(DACC(4)*5300/1030+1)*24.305+...
    (1-(DACC(4)*5300/1030)/(DACC(4)*5300/1030+1)).*40.078+...
    15.9994*3+12.011+(2+15.9994)*0.85;
mmolPml = 10.^(mgPml)./ACCmgPmmol;

% Amount of Ca and TE in ACC
TEmmolPml = repmat(mmolPml,1,size(DACC,1)).*repmat(DACC',size(mmolPml,1),1).*...
    repmat(SWconcs',size(mmolPml,1),1)./SWCa;
TEmmolPml = [zeros(1,size(TEmmolPml,2)) ; TEmmolPml]; % Add reference case that zero Ca is added by ACC
CammolPml = mmolPml.*(1-(DACC(4)*5300/1030)/(DACC(4)*5300/1030+1));
CammolPml = [0 ; CammolPml]; % Add reference case that zero Ca is added by ACC

% percentage calcite derived from solid state transformation
% fraction ranges from 0 to 1 between mg/ml ratios of 0 to 1
pSS = zeros(size(CammolPml));
pSS(201:301) = 0:1/100:1;
pSS(302:size(pSS,1),1) = 1;
pSS = pSS.*1;
pSS = abs(pSS-1);

% Apparent distribution coefficients calculated from modified reservoir in
% which some ACC is dissolved and then removal takes place via calcite
% precipitation
VacCa = CammolPml.*1000+10.3;
VacTE = NaN(size(CammolPml,1),1);
VacTECa = NaN(size(CammolPml,1),1);
DTE = NaN(size(CammolPml,1),size(F,2),8);
DTEnoSS = NaN(size(CammolPml,1),size(F,2),8);
for j = 1:8
    VacTE(:,j) = TEmmolPml(:,j).*1000+SWconcs(j);
    VacTECa(:,j) = (VacTE(:,j)./VacCa)./(SWconcs(j)/10.3);
    DTE(:,:,j) = (1 - repmat(F,size(VacCa,1),1).^Darag(j,1))./(1 - repmat(F,size(VacCa,1),1))...
        .*repmat(VacTECa(:,j),1,size(F,2));
    DTEnoSS(:,:,j) = DTE(:,:,j);
    DTE(:,:,j) = DTE(:,:,j).*repmat(pSS,1,size(DTE,2)) + ...
        abs(1 - repmat(pSS,1,size(DTE,2))).*repmat(DACC(j),size(DTE,1),size(DTE,2));
end


% coral data inferred F values assuming no ACC
coralFCaFibre = NaN(size(coralDsFibre));
for i = 1:height(CoralData)
    for j = 1:size(coralDsFibre,2)
        if ~isnan(coralDsFibre(i,j))
            tempLoc = [abs(DTE(1,:,j)'-coralDsFibre(i,j)) (1:1:size(DTE,2))'];
            tempLoc = sortrows(tempLoc,1);
            % If they disagree by more than 5%, assume Rayleigh model can't match data
            if abs(tempLoc(1,1)/coralDsFibre(i,j))*100>5
                coralFCaFibre(i,j) = NaN;
            else
                coralFCaFibre(i,j) = F(tempLoc(1,2));
            end
        end
    end
end

coralFCaCOC = NaN(size(coralDsCOC));
for i = 1:height(CoralData)
    for j = 1:size(coralDsCOC,2)
        if ~isnan(coralDsCOC(i,j))
            tempLoc = [abs(DTE(1,:,j)'-coralDsCOC(i,j)) (1:1:size(DTE,2))'];
            tempLoc = sortrows(tempLoc,1);
            % If they disagree by more than 5%, assume Rayleigh model can't match data
            if abs(tempLoc(1,1)/coralDsCOC(i,j))*100>5
                coralFCaCOC(i,j) = NaN;
            else
                coralFCaCOC(i,j) = F(tempLoc(1,2));
            end
        end
    end
end

%
%% Manuscript figures
% Figure 1: Aragonite compared to ACC, and an example Rayleigh model

close(figure(1))
H = figure(1);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[19 6.5],...
    'Units','centimeter','color',[1 1 1],'Position',[30 15 19 6.5])

cmap = parula(8);
t = tiledlayout(1,3);

t1 = nexttile;
hold on
for j = 1:7
    fill([Darag(j,3) Darag(j,2) Darag(j,2) Darag(j,3)],...
        [min((TEallpHEx(1:26,j)./1000)./(SWall(1:26,j)./CaAll(1:26).*1000)) ...
        min((TEallpHEx(1:26,j)./1000)./(SWall(1:26,j)./CaAll(1:26).*1000)) ...
        max((TEallpHEx(1:26,j)./1000)./(SWall(1:26,j)./CaAll(1:26).*1000)) ...
        max((TEallpHEx(1:26,j)./1000)./(SWall(1:26,j)./CaAll(1:26).*1000))],...
        cmap(j,:),'facealpha',0.2,'edgecolor',cmap(j,:))
end
plot([1e-5 1000],[1e-5 1000],'--k','linewidth',0.5)

set(t1,'yscale','log','xscale','log','fontsize',8,'box','on')
xticks([1e-4 1e-2 1 100])
ylim([1e-5 500])
xlim([1e-6 500])
xlabel('aragonite D_x range')
ylabel('ACC D_x range')

t2 = nexttile;
hold on
for j = 1:7
    if j~=5
        fill([Darag(j,2) Darag(j,2) Darag(j,3) Darag(j,3)],...
            [min(coralDsFibre(:,j)) max(coralDsFibre(:,j))  max(coralDsFibre(:,j)) min(coralDsFibre(:,j))],...
            cmap(j,:),'facealpha',0.2,'edgecolor',cmap(j,:))
    else
        fill([Darag(j,2) Darag(j,2) Darag(j,3) Darag(j,3)],...
            [Dcoral(5)-DcoralRange(5) Dcoral(5)+DcoralRange(5) Dcoral(5)+DcoralRange(5) Dcoral(5)-DcoralRange(5)],...
            cmap(j,:),'facealpha',0.2,'edgecolor',cmap(j,:))
    end
end
plot([1e-5 1000],[1e-5 1000],'--k','linewidth',0.5)

set(t2,'yscale','log','xscale','log','fontsize',8,'box','on')
xticks([1e-4 1e-2 1 100])
ylim([1e-5 500])
xlim([1e-6 500])
xlabel('aragonite D_x range')
ylabel('coral D_x range')

t3 = nexttile;
plot(1-F,DTE(1,:,4),'color',cmap(4,:))
hold on
plot(1-F,(1 - F.^Darag(4,2))./(1 - F),'--','color',cmap(4,:))
plot(1-F,(1 - F.^Darag(4,3))./(1 - F),'--','color',cmap(4,:))
scatter(1-coralFCaFibre(:,4),coralDsFibre(:,4),10,cmap(4,:),'filled')
xlabel('fraction Ca utilised (f_{Ca})')
ylabel('apparent distribution coefficient')
set(t3,'fontsize',8)

text(t1,Darag(1,2),0.6*min((TEallpHEx(1:26,1)./1000)./(SWall(1:26,1)./CaAll(1:26).*1000)),'Li','color',cmap(1,:),'fontsize',7,'fontweight','bold')
text(t1,Darag(2,2)*1.1,1.5*max((TEallpHEx(1:26,2)./1000)./(SWall(1:26,2)./CaAll(1:26).*1000)),'B','color',cmap(2,:),'fontsize',7,'fontweight','bold')
text(t1,Darag(3,2),1.5*max((TEallpHEx(1:26,3)./1000)./(SWall(1:26,3)./CaAll(1:26).*1000)),'Na','color',cmap(3,:),'fontsize',7,'fontweight','bold')
text(t1,Darag(4,2),1.5*max((TEallpHEx(1:26,4)./1000)./(SWall(1:26,4)./CaAll(1:26).*1000)),'Mg','color',cmap(4,:),'fontsize',7,'fontweight','bold')
text(t1,Darag(5,2)*1.3,0.8*max((TEallpHEx(1:26,5)./1000)./(SWall(1:26,5)./CaAll(1:26).*1000)),'Mn','color',cmap(5,:),'fontsize',7,'fontweight','bold')
text(t1,0.7*Darag(6,2),0.55*min((TEallpHEx(1:26,6)./1000)./(SWall(1:26,6)./CaAll(1:26).*1000)),'Sr','color',cmap(6,:),'fontsize',7,'fontweight','bold')
text(t1,Darag(7,3)*1.2,1.2*min((TEallpHEx(1:26,7)./1000)./(SWall(1:26,7)./CaAll(1:26).*1000)),'Ba','color',cmap(7,:),'fontsize',7,'fontweight','bold')

text(t2,Darag(1,2),0.6*min(coralDsFibre(:,1)),'Li','color',cmap(1,:),'fontsize',7,'fontweight','bold')
text(t2,Darag(2,2)*1.1,0.6*min(coralDsFibre(:,2)),'B','color',cmap(2,:),'fontsize',7,'fontweight','bold')
text(t2,Darag(3,2),0.6*min(coralDsFibre(:,3)),'Na','color',cmap(3,:),'fontsize',7,'fontweight','bold')
text(t2,Darag(4,2)*0.9,1.7*max(coralDsFibre(:,4)),'Mg','color',cmap(4,:),'fontsize',7,'fontweight','bold')
text(t2,Darag(5,2)*0.2,1*(Dcoral(5)+DcoralRange(5)),'Mn','color',cmap(5,:),'fontsize',7,'fontweight','bold')
text(t2,0.65*Darag(6,2),0.3*max(coralDsFibre(:,6)),'Sr','color',cmap(6,:),'fontsize',7,'fontweight','bold')
text(t2,Darag(7,3)*1.2,0.7*max(coralDsFibre(:,7)),'Ba','color',cmap(7,:),'fontsize',7,'fontweight','bold')

text(t1,1.5e-6,250,'A','fontsize',9,'fontweight','bold')
text(t2,1.5e-6,250,'B','fontsize',9,'fontweight','bold')
text(t3,0.02,0.012*0.96,'C','fontsize',9,'fontweight','bold')

set(t,"TileSpacing","compact")
set(t,"Padding","tight")

print(H,'-depsc','-r600','Figure1_Dcomp_Rayleigh_ex_R1')

%
% Figure 2: Example expansion of the model to include ACC

close(figure(2))
H = figure(2);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[9 7],...
    'Units','centimeter','color',[1 1 1],'Position',[30 15 9 7])

Hp = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,6)) - log10(DTE(1,:,6)));
set(Hp,'linestyle','none')
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,6).*SWconcs(6)/SWCa.*1000,[7 8 8.9 9 9.1],'--r','linewidth',0.5);
clabel(c,h,'fontsize',5,'color','r')

[c1 h1] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,6).*SWconcs(6)/SWCa.*1000,...
    round([Dcoral(6) Dcoral(6)].*SWconcs(6)/SWCa.*1000,2),'-k',...
    'linewidth',1.5);
clabel(c1,h1,'fontsize',6,'color','k')

xlabel('fraction Ca utilised (f_{Ca})')
ylabel('ACC delivery to ECF log_{10}(mg/ml)')

tMid = 0;
clim(gca,[tMid-2 tMid+2])

set(gca,'box','on','layer','top','fontsize',8)
xlim([0 1])

cmap = [bone ; flipud(pink)];
colormap(cmap)

set(gca,'position',[0.13,0.11,0.6,0.815])

H1C = colorbar;
set(H1C,'position',[0.86 0.11 0.02 0.815],...
     'Ticks',-2:0.5:2,'TickLabels',["-2" "-1.5" "-1.0" "-0.5" "0" "0.5" "1.0" "1.5" "2.0"])
ylabel(H1C,'log_{10}(change in D)','rotation',270)

ax2 = axes;
ax2.Visible = "Off";
cmap2 = [zeros(512*2,3) ; gray];
colormap(ax2,cmap2)
H2C = colorbar(ax2);
set(H2C,'position',[0.735 0.11 0.02 0.815],'fontsize',7,'Ticks',[0 4/6 1],'TickLabels',[0 0 1])
ylabel(H2C,'ACC fraction crystallised via solid state transformation','rotation',270)

print(H,'-depsc','-r600','Figure2_example_Rayleigh_ACC_model_v2_R1')

%
% Figure 3: Impact of ACC on Rayleigh distillation models across TE systems

close(figure(3))
H = figure(3);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[18.5 11],...
    'Units','centimeter','color',[1 1 1],'Position',[30 15 18.5 11])

H1 = axes(); %Li
H1p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,1)) - log10(DTE(1,:,1)));
set(H1p,'linestyle','none')
set(H1,'box','on','position',[0.05 0.57 0.22 0.38],'fontsize',8)
ylabel('log ACC/SW ratio (mg/ml)')
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,1).*SWconcs(1)/SWCa.*1000,...
    round([Dcoral(1) Dcoral(1)].*SWconcs(1)/SWCa.*1000,0),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')

H2 = axes(); %Na
H2p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,3)) - log10(DTE(1,:,3)));
set(H2p,'linestyle','none')
set(H2,'box','on','position',[0.31 0.57 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,3).*SWconcs(3)/SWCa.*1000,...
    round([Dcoral(3) Dcoral(3)].*SWconcs(3)/SWCa.*1000,0),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')

H31 = axes();
set(H31,'position',[0.57 0.57 0.22 0.38],'fontsize',8)
H31.XColor = 'w';
H31.YAxisLocation = 'right';
ylim([-2 1])
yticks(-2:1:1)
set(H31,'yticklabels',[round(((CammolPml(2)*1000)/(CammolPml(2)*1000 + 10.6))*100,1) ...
    round(((CammolPml(101)*1000)/(CammolPml(101)*1000 + 10.6))*100,0) round(((CammolPml(201)*1000)/(CammolPml(201)*1000 + 10.6))*100,0) ...
    round(((CammolPml(301)*1000)/(CammolPml(301)*1000 + 10.6))*100,0)])
H31l = ylabel('skeletal Ca dervied from ACC (%)','rotation',-90);
pos = get(H31l,'position');
set(H31l,'position',pos+[0.1 0 0])
H3 = axes(); %Mg
H3p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,4)) - log10(DTE(1,:,4)));
set(H3p,'linestyle','none')
set(H3,'box','on','position',[0.57 0.57 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,4).*SWconcs(4)/SWCa.*1000,...
    round([Dcoral(4) Dcoral(4)].*SWconcs(4)/SWCa.*1000,1),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')

H4 = axes(); %Mn
H4p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,5)) - log10(DTE(1,:,5)));
set(H4p,'linestyle','none')
set(H4,'box','on','position',[0.05 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,5).*SWconcs(5)/SWCa.*1000,...
    round([Dcoral(5) Dcoral(5)].*SWconcs(5)/SWCa.*1000,2),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')
ylabel('log ACC/SW ratio (mg/ml)')
xlabel('fraction Ca utilised (f_{Ca})')

H5 = axes(); %Sr
H5p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,6)) - log10(DTE(1,:,6)));
set(H5p,'linestyle','none')
set(H5,'box','on','position',[0.31 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,6).*SWconcs(6)/SWCa.*1000,...
    round([Dcoral(6) Dcoral(6)].*SWconcs(6)/SWCa.*1000,1),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')
xlabel('fraction Ca utilised (f_{Ca})')

H61 = axes();
set(H61,'position',[0.57 0.09 0.22 0.38],'fontsize',8)
H61.XColor = 'w';
H61.YAxisLocation = 'right';
ylim([-2 1])
yticks(-2:1:1)
set(H61,'yticklabels',[round(((CammolPml(2)*1000)/(CammolPml(2)*1000 + 10.6))*100,1) ...
    round(((CammolPml(101)*1000)/(CammolPml(101)*1000 + 10.6))*100,0) round(((CammolPml(201)*1000)/(CammolPml(201)*1000 + 10.6))*100,0) ...
    round(((CammolPml(301)*1000)/(CammolPml(301)*1000 + 10.6))*100,0)])
H61l = ylabel('skeletal Ca dervied from ACC (%)','rotation',-90);
pos = get(H61l,'position');
set(H61l,'position',pos+[0.1 0 0])
H6 = axes(); %Ba
H6p = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,7)) - log10(DTE(1,:,7)));
set(H6p,'linestyle','none')
set(H6,'box','on','position',[0.57 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,7).*SWconcs(7)/SWCa.*1000,...
    round([Dcoral(7) Dcoral(7)].*SWconcs(7)/SWCa.*1000,1),'-k',...
    'linewidth',1.5);
clabel(c,h,'fontsize',6,'color','k')
xlabel(H6,'fraction Ca utilised (f_{Ca})')

text(H1,0.03,0.84,'Li','fontweight','bold','fontsize',9,'color','k')
text(H2,0.03,0.84,'Na','fontweight','bold','fontsize',9,'color','k')
text(H3,0.03,0.84,'Mg','fontweight','bold','fontsize',9,'color','w')
text(H4,0.03,0.84,'Mn','fontweight','bold','fontsize',9,'color','w')
text(H5,0.03,0.84,'Sr','fontweight','bold','fontsize',9,'color','k')
text(H6,0.03,0.84,'Ba','fontweight','bold','fontsize',9,'color','k')

tMid = 0;
clim(H1,[tMid-2 tMid+2])
clim(H2,[tMid-2 tMid+2])
clim(H3,[tMid-2 tMid+2])
clim(H4,[tMid-2 tMid+2])
clim(H5,[tMid-2 tMid+2])
clim(H6,[tMid-2 tMid+2])

set(H1,'box','on','layer','top')
set(H2,'box','on','layer','top')
set(H3,'box','on','layer','top')
set(H4,'box','on','layer','top')
set(H5,'box','on','layer','top')
set(H6,'box','on','layer','top')

cmap = [bone ; flipud(pink)];
colormap(cmap)

H1C = colorbar;
set(H1C,'position',[0.925 0.09 0.015 0.85],...
     'Ticks',-2:0.5:2,'TickLabels',["-2" "-1.5" "-1.0" "-0.5" "0" "0.5" "1.0" "1.5" "2.0"])
ylabel(H1C,'log_{10}(change in D)','rotation',270)

ax2 = axes;
ax2.Visible = "Off";
cmap2 = [zeros(512*2,3) ; gray];
colormap(ax2,cmap2)
H2C = colorbar(ax2);
set(H2C,'position',[0.85 0.57 0.01 0.38],'fontsize',7,'Ticks',[0 4/6 9/12 10/12 11/12 1],'TickLabels',[0 0 0.25 0.5 0.75 1])
ylabel(H2C,'fraction solid state ACC transformation','rotation',270)

ax3 = axes;
ax3.Visible = "Off";
cmap2 = [zeros(512*2,3) ; gray];
colormap(ax3,cmap2)
H3C = colorbar(ax3);
set(H3C,'position',[0.85 0.09 0.01 0.38],'fontsize',7,'Ticks',[0 4/6 9/12 10/12 11/12 1],'TickLabels',[0 0 0.25 0.5 0.75 1])
ylabel(H3C,'fraction solid state ACC transformation','rotation',270)

print(H,'-dpdf','-r600','Figure3_full_ACC_Rayleigh_model_R1')
%print(H,'-depsc','-r600','FigureS1_full_ACC_Rayleigh_model_ACC_D05_R1')

%
% Figure 4: COC/fibre ratios in context of this model

close(figure(4))
H = figure(4);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[16 10],...
    'Units','centimeter','color',[1 1 1],'Position',[0 0 16 10])
cmap = parula(12);
t = tiledlayout(2,3);

% regression
tempX = CoralData.Mg_fibre_mmol_mol_;
tempY = CoralData.Mg_COC_mmol_mol_;
tempY(tempY>15) = NaN;
tempex = isnan(tempX) | isnan(tempY);
tempX = tempX(~tempex);
tempY = tempY(~tempex);
modelfun = @(b,x) b(1)+b(2).*x; 
beta0 = [0.5 0.5];
mdl1 = fitnlm(tempX,tempY,modelfun,beta0);
% plot
t1 = nexttile;
plot([0 15],[0 15],':','color',[0.5 0.5 0.5],'linewidth',0.5)
hold on
scatter(CoralData.Mg_fibre_mmol_mol_,CoralData.Mg_COC_mmol_mol_,10,[0.9 0.9 0.9],'filled','markeredgecolor','k','linewidth',0.25)
plot((round(min(tempX),1):0.1:round(max(tempX),1)),...
    mdl1.Coefficients{1,1}+mdl1.Coefficients{2,1}.*...
    (round(min(tempX),1):0.1:round(max(tempX),1)),'-','color','k',...
    'linewidth',1.)
[beta,R,J,COVB,MSE] = nlinfit(tempX,tempY,modelfun,beta0);
[ypred,delta] = nlpredci(modelfun,...
    (round(min(tempX),1):0.1:round(max(tempX),1)),beta,R,'Covar',COVB,...
    'MSE',MSE,'SimOpt','off');
lower = ypred - delta;
upper = ypred + delta;
plot((round(min(tempX),1):0.1:round(max(tempX),1)),...
    [lower;upper],'--','color','k','linewidth',0.5)

xlabel('fibre Mg/Ca (mmol/mol)')
ylabel('COC Mg/Ca (mmol/mol)')

% regression
tempX = CoralData.Sr_fibre_mmol_mol_;
tempY = CoralData.Sr_COC_mmol_mol_;
tempY(tempY>15) = NaN;
tempex = isnan(tempX) | isnan(tempY);
tempX = tempX(~tempex);
tempY = tempY(~tempex);
modelfun = @(b,x) b(1)+b(2).*x; 
beta0 = [0.5 0.5];
mdl2 = fitnlm(tempX,tempY,modelfun,beta0);
% plot
t2 = nexttile(4);
plot([0 15],[0 15],':','color',[0.5 0.5 0.5],'linewidth',0.5)
hold on
scatter(CoralData.Sr_fibre_mmol_mol_,CoralData.Sr_COC_mmol_mol_,10,[0.9 0.9 0.9],'filled','markeredgecolor','k','linewidth',0.25)
plot((round(min(tempX),1):0.1:round(max(tempX),1)),...
    mdl2.Coefficients{1,1}+mdl2.Coefficients{2,1}.*...
    (round(min(tempX),1):0.1:round(max(tempX),1)),'-','color','k',...
    'linewidth',1.)
[beta,R,J,COVB,MSE] = nlinfit(tempX,tempY,modelfun,beta0);
[ypred,delta] = nlpredci(modelfun,...
    (round(min(tempX),1):0.1:round(max(tempX),1)),beta,R,'Covar',COVB,...
    'MSE',MSE,'SimOpt','off');
lower = ypred - delta;
upper = ypred + delta;
plot((round(min(tempX),1):0.1:round(max(tempX),1)),...
    [lower;upper],'--','color','k','linewidth',0.5)

xlabel('fibre Sr/Ca (mmol/mol)')
ylabel('COC Sr/Ca (mmol/mol)')

% panel C calculations
yRayleigh = [DTE(101,50,4) DTE(101,100,4) DTE(101,200,4) DTE(101,400,4) DTE(101,999,4)].*SWconcs(4)/SWCa*1000;
xRayleigh = [DTE(101,50,6) DTE(101,100,6) DTE(101,200,6) DTE(101,400,6) DTE(101,999,6)].*SWconcs(6)/SWCa*1000;
txtRayleigh = {'0.95','0.9','0.8','0.6','0'};

yACC = [DTE(1,550,4) DTE(151,550,4) DTE(201,550,4) DTE(202,550,4) DTE(251,550,4)].*SWconcs(4)/SWCa*1000;
xACC = [DTE(1,550,6) DTE(151,550,6) DTE(201,550,6) DTE(202,550,6) DTE(251,550,6)].*SWconcs(6)/SWCa*1000;
txtACC = {'-2','-0.5','0','0.01','0.5'};

t3 = nexttile(2,[2 2]);

plot(DTE(101,:,4).*SWconcs(4)/SWCa*1000,DTE(101,:,6).*SWconcs(6)/SWCa*1000,'color',cmap(3,:))
hold on
scatter(yRayleigh,xRayleigh,10,cmap(3,:),'filled')
text(yRayleigh+0.08,xRayleigh-0.03,txtRayleigh,'color',cmap(3,:),'fontsize',6)

plot(DTEnoSS(:,550,4).*SWconcs(4)/SWCa*1000,DTEnoSS(:,550,6).*SWconcs(6)/SWCa*1000,'--','color',cmap(6,:))
plot(DTE(:,550,4).*SWconcs(4)/SWCa*1000,DTE(:,550,6).*SWconcs(6)/SWCa*1000,'color',cmap(6,:))
scatter(yACC,xACC,10,cmap(6,:),'filled')
text(yACC-0.7,xACC-0.02,txtACC,'color',cmap(5,:),'fontsize',6)

coralMgMean = [mean(CoralData.Mg_fibre_mmol_mol_(CoralData.deep_sea~=1),'omitnan') mean(CoralData.Mg_COC_mmol_mol_(CoralData.deep_sea~=1),'omitnan')];
coralMgSTD = [std(CoralData.Mg_fibre_mmol_mol_(CoralData.deep_sea~=1),'omitnan') std(CoralData.Mg_COC_mmol_mol_(CoralData.deep_sea~=1),'omitnan')];
coralSrMean = [mean(CoralData.Sr_fibre_mmol_mol_(CoralData.deep_sea~=1),'omitnan') mean(CoralData.Sr_COC_mmol_mol_(CoralData.deep_sea~=1),'omitnan')];
coralSrSTD = [std(CoralData.Sr_fibre_mmol_mol_(CoralData.deep_sea~=1),'omitnan') std(CoralData.Sr_COC_mmol_mol_(CoralData.deep_sea~=1),'omitnan')];
%errorbar(coralMgMean',coralSrMean'.*1.05,coralSrSTD',coralSrSTD',coralMgSTD',coralMgSTD','o','markersize',10,'color','k',...
%    'markerfacecolor',[0.7 0.7 0.7],'linewidth',0.5,'capsize',0)
 scatter(coralMgMean,coralSrMean*1.05,...
     40,[0.7 0.7 0.7],'filled','markeredgecolor','k','linewidth',0.5)
%plot([1 3],[1 3.*0.0826+1],'--','color',[0.3 0.3 0.3])

xlabel('Mg/Ca (mmol/mol)')
ylabel('Sr/Ca (mmol/mol)')
axis([0 20 8.6 9.7])

text(t1,1,8,['y = ',num2str(round(mdl1.Coefficients{2,1},2)),'x + ',num2str(round(mdl1.Coefficients{1,1},2))],'fontsize',6)
text(t1,1,6.5,{['R^2 = ',num2str(round(mdl1.Rsquared.Ordinary,2))],'p <0.01'},'fontsize',6)
text(t2,7.5,11.25,['y = ',num2str(round(mdl2.Coefficients{2,1},2)),'x + ',num2str(round(mdl2.Coefficients{1,1},2))],'fontsize',6)
text(t2,7.5,10.7,{['R^2 = ',num2str(round(mdl2.Rsquared.Ordinary,2))],'p <0.01'},'fontsize',6)

text(t3,9,9.4,'rate','fontsize',6)
text(t3,13,8.82,'f_{Ca} @ 10^{-1} mg/ml ACC','fontsize',6,'color',cmap(3,:))
text(t3,12,9.02,{'log_{10}(ACC/SW) @ f_{Ca}=0.45','(solid state crystallisation','begins at 10^0)'},'fontsize',6,'color',cmap(6,:))
text(t3,2,9.5,{'ACC/SW','no solid state','crystallisation'},'fontsize',6,'color',cmap(6,:))
text(t3,coralMgMean(1)-2.9,coralSrMean(1)*1.05-0.03,'mean fibre','fontsize',6,'color',[0.3 0.3 0.3],'fontweight','bold')
text(t3,coralMgMean(2)+0.3,coralSrMean(2)*1.05+0.035,'mean COC','fontsize',6,'color',[0.3 0.3 0.3],'fontweight','bold')

text(t1,8*0.025,15*0.94,'A','fontsize',9,'fontweight','bold')
text(t2,5*0.025+7,5.5*0.94+7,'B','fontsize',9,'fontweight','bold')
text(t3,20*0.0125,1.1*0.975+8.6,'C','fontsize',9,'fontweight','bold')

set(t1,'fontsize',8,'box','on','layer','top','xlim',[0 8],'ylim',[0 15])
set(t2,'fontsize',8,'box','on','layer','top','xlim',[7 12],'ylim',[7 12.5])
set(t3,'fontsize',8,'box','on','layer','top')

set(t,"TileSpacing","compact")
set(t,"Padding","compact")


% starting from ACC/SW = 0
arrowX = [DTE(201,750,4).*SWconcs(4)/SWCa*1000+0.3 DTE(201,750,4).*SWconcs(4)/SWCa*1000+0.3*DTE(201,750,4).*SWconcs(4)/SWCa*1000*10^-4/10^-4.7+0.3];
arrowY = [DTE(201,750,6).*SWconcs(6)/SWCa*1000+0.02 DTE(201,750,6).*SWconcs(6)/SWCa*1000+0.3*DTE(201,750,6).*SWconcs(6)/SWCa*1000*((1.13/0.97 - 1)/2)+0.02];

% starting from mean fibre
%arrowX = [coralMgMean(1)+0.5 coralMgMean(1)+0.45*coralMgMean(1)*10^-4/10^-4.7+0.5];
%arrowY = [coralSrMean(1)*1.05+0.02 coralSrMean(1)*1.05+0.45*coralSrMean(1)*1.05*((1.13/0.97 - 1)/2)+0.02];
ax = t3;
[normx, normy] = coord2norm(ax, arrowX, arrowY);
annotation('arrow', normx, normy);

print(H,'-dpdf','-r600','Figure4_COC_fibre_comparison_v2_R1')

%
% Figure S2: Impact of ACC on Rayleigh distillation models across TE systems, no solid state

close(figure(3))
H = figure(3);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[18.5 11],...
    'Units','centimeter','color',[1 1 1],'Position',[30 15 18.5 11])

H1 = axes(); %Li
H1p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,1)) - log10(DTEnoSS(1,:,1)));
set(H1p,'linestyle','none')
set(H1,'box','on','position',[0.05 0.57 0.22 0.38],'fontsize',8)
ylabel('log ACC/SW ratio (mg/ml)')
xlim([0 1])
hold on
contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,1)),[log10(Dcoral(1)) log10(Dcoral(1))],'-k',...
    'linewidth',1.5)

H2 = axes(); %Na
H2p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,3)) - log10(DTEnoSS(1,:,3)));
set(H2p,'linestyle','none')
set(H2,'box','on','position',[0.31 0.57 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
[c h] = contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,3)),[log10(Dcoral(3)) log10(Dcoral(3))],'-k',...
    'linewidth',1.5);

H31 = axes();
set(H31,'position',[0.57 0.57 0.22 0.38],'fontsize',8)
H31.XColor = 'w';
H31.YAxisLocation = 'right';
ylim([-2 1])
yticks(-2:1:1)
set(H31,'yticklabels',[round(((CammolPml(2)*1000)/(CammolPml(2)*1000 + 10.6))*100,1) ...
    round(((CammolPml(101)*1000)/(CammolPml(101)*1000 + 10.6))*100,0) round(((CammolPml(201)*1000)/(CammolPml(201)*1000 + 10.6))*100,0) ...
    round(((CammolPml(301)*1000)/(CammolPml(301)*1000 + 10.6))*100,0)])
H31l = ylabel('skeletal Ca dervied from ACC (%)','rotation',-90);
pos = get(H31l,'position');
set(H31l,'position',pos+[0.1 0 0])
H3 = axes(); %Mg
H3p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,4)) - log10(DTEnoSS(1,:,4)));
set(H3p,'linestyle','none')
set(H3,'box','on','position',[0.57 0.57 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,4)),[log10(Dcoral(4)) log10(Dcoral(4))],'-k',...
    'linewidth',1.5)

H4 = axes(); %Mn
H4p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,5)) - log10(DTEnoSS(1,:,5)));
set(H4p,'linestyle','none')
set(H4,'box','on','position',[0.05 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,5)),[log10(Dcoral(5)) log10(Dcoral(5))],'-k',...
    'linewidth',1.5)
ylabel('log ACC/SW ratio (mg/ml)')
xlabel('F Ca utilised')

H5 = axes(); %Sr
H5p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,6)) - log10(DTEnoSS(1,:,6)));
set(H5p,'linestyle','none')
set(H5,'box','on','position',[0.31 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,6)),[log10(Dcoral(6)) log10(Dcoral(6))],'-k',...
    'linewidth',1.5)
xlabel('F Ca utilised')

H61 = axes();
set(H61,'position',[0.57 0.09 0.22 0.38],'fontsize',8)
H61.XColor = 'w';
H61.YAxisLocation = 'right';
ylim([-2 1])
yticks(-2:1:1)
set(H61,'yticklabels',[round(((CammolPml(2)*1000)/(CammolPml(2)*1000 + 10.6))*100,1) ...
    round(((CammolPml(101)*1000)/(CammolPml(101)*1000 + 10.6))*100,0) round(((CammolPml(201)*1000)/(CammolPml(201)*1000 + 10.6))*100,0) ...
    round(((CammolPml(301)*1000)/(CammolPml(301)*1000 + 10.6))*100,0)])
H61l = ylabel('skeletal Ca dervied from ACC (%)','rotation',-90);
pos = get(H61l,'position');
set(H61l,'position',pos+[0.1 0 0])
H6 = axes(); %Ba
H6p = pcolor(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,7)) - log10(DTEnoSS(1,:,7)));
set(H6p,'linestyle','none')
set(H6,'box','on','position',[0.57 0.09 0.22 0.38],'fontsize',8)
xlim([0 1])
hold on
contour(fliplr(F),mgPml,log10(DTEnoSS(2:size(DTEnoSS,1),:,7)),[log10(Dcoral(7)) log10(Dcoral(7))],'-k',...
    'linewidth',1.5)
xlabel(H6,'F Ca utilised')

text(H1,0.03,0.82,'Li','fontweight','bold','fontsize',9,'color','k')
text(H2,0.03,0.82,'Na','fontweight','bold','fontsize',9,'color','k')
text(H3,0.03,0.82,'Mg','fontweight','bold','fontsize',9,'color','w')
text(H4,0.03,0.82,'Mn','fontweight','bold','fontsize',9,'color','w')
text(H5,0.03,0.82,'Sr','fontweight','bold','fontsize',9,'color','k')
text(H6,0.03,0.82,'Ba','fontweight','bold','fontsize',9,'color','k')

tMid = 0;
clim(H1,[tMid-2 tMid+2])
clim(H2,[tMid-2 tMid+2])
clim(H3,[tMid-2 tMid+2])
clim(H4,[tMid-2 tMid+2])
clim(H5,[tMid-2 tMid+2])
clim(H6,[tMid-2 tMid+2])

set(H1,'box','on','layer','top')
set(H2,'box','on','layer','top')
set(H3,'box','on','layer','top')
set(H4,'box','on','layer','top')
set(H5,'box','on','layer','top')
set(H6,'box','on','layer','top')

cmap = [bone ; flipud(pink)];
colormap(cmap)

H1C = colorbar;
set(H1C,'position',[0.925 0.09 0.015 0.85],...
     'Ticks',-2:0.5:2,'TickLabels',["-2" "-1.5" "-1.0" "-0.5" "0" "0.5" "1.0" "1.5" "2.0"])
ylabel(H1C,'log_{10}(change in D)','rotation',270)

ax2 = axes;
ax2.Visible = "Off";
cmap2 = [zeros(512*2,3)];
colormap(ax2,cmap2)
H2C = colorbar(ax2);
set(H2C,'position',[0.85 0.57 0.01 0.38],'fontsize',7,'Ticks',[0 1],'TickLabels',[0 0])
ylabel(H2C,'fraction solid state ACC transformation','rotation',270)

ax3 = axes;
ax3.Visible = "Off";
cmap2 = [zeros(512*2,3)];
colormap(ax3,cmap2)
H3C = colorbar(ax3);
set(H3C,'position',[0.85 0.09 0.01 0.38],'fontsize',7,'Ticks',[0 1],'TickLabels',[0 0])
ylabel(H3C,'fraction solid state ACC transformation','rotation',270)

print(H,'-dpng','-r600','FigureS2_full_ACC_Rayleigh_model_no_solid_state_R1')

%
% Figure S3: Example contour plot showing the degree to which ACC may be more clearly visible in spatially-resolved measurements

close(figure(2))
H = figure(2);
set(H,'PaperUnits','centimeter','PaperPositionMode','auto','papersize',[8.5 7],...
    'Units','centimeter','color',[1 1 1],'Position',[30 15 8.5 7])

Hp = pcolor(fliplr(F),mgPml,log10(DTE(2:size(DTE,1),:,5)) - log10(DTE(1,:,5)));
set(Hp,'linestyle','none')
hold on
[c h] = contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,5),[1 2 5 10 50 100 500],'--r','linewidth',0.5);
clabel(c,h,'fontsize',5,'color','r')

contour(fliplr(F),mgPml,DTE(2:size(DTE,1),:,5).*SWconcs(5)/SWCa.*1000,...
    [Dcoral(5) Dcoral(5)].*SWconcs(5)/SWCa.*1000,'-k',...
    'linewidth',1.5)

xlabel('F Ca utilised')
ylabel('ACC delivery to ECF log_{10}(mg/ml)')

tMid = 0;
clim(gca,[tMid-2 tMid+2])

set(gca,'box','on','layer','top','fontsize',8)
xlim([0 1])
ylim([-2 1])

cmap = [bone ; flipud(pink)];
colormap(cmap)

set(gca,'position',[0.13,0.11,0.7,0.815])

H1C = colorbar;
set(H1C,'position',[0.85 0.11 0.02 0.815],...
     'Ticks',-2:0.5:2,'TickLabels',["-2" "-1.5" "-1.0" "-0.5" "0" "0.5" "1.0" "1.5" "2.0"])
ylabel(H1C,'log_{10}(change in D)','rotation',270)