% OAE2Os_repository.m
% Code author: Matthew M. Jones, PhD (Northwestern University/Smithsonian National Museum of Natural History)
% Jones et al. (in press) "Mid-Cretaceous ocean acidification linked to massive volcanism" Nature Geoscience

%%% This MATLAB script plots Os timeseries data through Cretaceous Ocean Anoxic Event 2 from the Angus Core in Colorado and International Ocean Discovery Program (IODP) Site U1516 in the Mentelle Basin. It then plots output from a simple box model that can be iteratively adjusted to match the measured Os chemostratigraphy. This script follows Tejada et al. (2009).

% establish Angus Core Os time series (from Jones et al 2021 GSA Bulletin) and define time zero in floating time series ("ALIGN" = time adjustment to align with Mentelle Basin time series).
AngusOsTime = [475 320 186 113 61 -12 -25 -35 -51 -59 -71 -81 -91 -108 -235]';
AngusOsDat = [0.683 0.652 0.388 0.167 0.193 0.184 0.168 0.180 0.175 0.247 0.772 0.855 0.863 0.910 0.753]';
AngusOs192Dat = [25.3 49.0 60.9 184.8 567.0 285.9 931.4 1183.9 902.2 632.6 89.7 87.9 87.7 82.0 59.5]';
ALIGN = AngusOsTime(find(AngusOsDat<0.3,1,'last'));
AngusOs = [AngusOsTime-ALIGN,AngusOsDat];
AngusOs192 = [AngusOsTime-ALIGN,AngusOs192Dat];

% IODP Site U1516 Os time series
U1516OsTime = [1214 1057 880 763 676 641 586 507 463 368 303 247 213 127 39 0 -28 -96 -169 -195 -225 -318 -417 -438 -561 -777 -1003 -1166]';
U1516OsDat = [0.55 0.57 0.56 0.13 0.27 0.52 0.55 0.47 0.14 0.18 0.42 0.17 0.31 0.18 0.17 0.43 0.43 0.38 0.62 0.63 0.56 0.66 0.76 0.76 0.73 0.74 0.73 0.75]';
U1516Os192Dat = [10.2    8.1    6.9  169.5    6.7   12.3   11.6   12.9   65.6   17.2   58.8   21.7   31.1 1847.6  273.8   30.5   31.4   21.7   12.7   11.5 17.1   14.0   22.2   35.5   20.5   15.3   12.6   25.3]';
U1516Os = [U1516OsTime,U1516OsDat];
U1516Os192 = [U1516OsTime,U1516Os192Dat];


% Initialize vectors for model
N = 3000; % duration of model run in kiloyears
M = zeros(1,N); % Vector to store the marine Os mass through time
Ro = zeros(1,N); % Vector to store the Os isotopic ratio through time

% Create time vector
TIME = (1:N)-1500; % time vector in Kyr, can be shifted to align with Angus Core datum
AMP = 35;

% Initialize model parameters
M(1) = 13000; % tonnes of Os in Ocean (from Tejada et al 2009; Levasseur et al 1998)
Ro(1) = 0.75; % pre-OAE2 188/187Os (Du Vivier et al 2014)

% Constant Os fluxes in tonnes/kyr
Friv = 1.15*295; % Riverine osmium flux
Fmantle = 0.85*414; % Unradiogenic continental weathering osmium flux pre-OAE2
Fcosmic = 17.6; % Cosmic osmium flux

% End member isotopic compositions
Rriv = 1.4; % isotopic ratio of continental weathered Os; Tejada et al 2009; Levasseur et al. 1999
Rmantle = 0.126; % Levasseur et al. 1999
Rcosmic = 0.126; % Levasseur et al. 1999
Rclip = 0.126;


%%%% SCENARIO - LINEARLY INTERPOLATED USER-DEFINED PULSE OF LIP VOLCANISM %%%%
Famp = AMP*Fmantle;
% ENTER CONTROL POINTS HERE [TIME_1 Famp_1; TIME_2 Famp_2; TIME_3 Famp_3; TIME_N Famp_N]
ControlPts = [min(TIME) 0; -440 0; -200 0.75; -20 2.5; -10 AMP; 50 AMP; 170 3; 280 1.25; 550 0.25; 800 0; max(TIME) 0];
Mag = ControlPts(:,2)*Fmantle;
Control_INDs = find(ismember(TIME,ControlPts(:,1)));
Fclip = interp1(Control_INDs,Mag,1:N);

IND_end = find(TIME==-10,1,'first');
Fclip(IND_end:end) = Fclip(IND_end:end)+0.75*Fmantle;

% run a simple box model of marine osmium cycle
for t = 2:N
    % determine changes in mass and isotopic composition of marine osmium
    % reservoir per time step
    Fsed = 0.05*M(t-1);
    dM = Friv + Fmantle + Fcosmic + Fclip(t) - Fsed;
    
    M(t) = M(t-1)+dM;
    dRo = (Friv*(Rriv-Ro(t-1)) + Fmantle*(Rmantle-Ro(t-1)) + Fclip(t)*(Rclip-Ro(t-1)) + Fcosmic*(Rcosmic-Ro(t-1)))/M(t-1);
    Ro(t) = Ro(t-1)+dRo;
end


% plot geochemical data alongside model output 

% plot panel 1 (right) with Os isotope data and model output. Angus Core Osi data in red, IODP Site U1516 Osi data in blue, model output in black.
figure
set(gcf,'Position',[0,0,500,500])
subplot(1,5,3:5)
p1 = plot(Ro,TIME,'b','LineWidth',3);
set(gca,'FontWeight','bold','YTickLabel',[],'TickDir','out','XColor','r')
axis([0 1.0 min(TIME)-1 max(TIME)+1])
hold on
x = [0.126 0.126];
y = [min(TIME) max(TIME)];
line(x,y,'Color','k','LineStyle','--')
p2 = plot(AngusOs(:,2),AngusOs(:,1)-20,'^r-.','MarkerFaceColor','w');
p3 = plot(U1516Os(:,2),U1516Os(:,1)-20,'or:','MarkerFaceColor','w');
xlabel('^{187/188}Os_{ocean}','FontWeight','bold')
ax1 = gca;
hold on

ax2 = axes('Position',ax1.Position,'XAxisLocation','top','Color','none','YTickLabel',[],'XColor','k','FontWeight','bold','TickDir','out');
hold(ax2,'on')
p5 = plot(ax2,U1516Os192(:,2),U1516Os192(:,1),'ok','MarkerSize',3,'MarkerFaceColor','k');
p6 = plot(ax2,AngusOs192(:,2),AngusOs192(:,1),'^k','MarkerSize',3,'MarkerFaceColor','k');

axis([0 2000 min(TIME)-1 max(TIME)+1])
xlabel('^{192}Os (ppt)')

% left panel plots user-input LIP volcanic Os flux
subplot(1,5,1:2)
p4=plot(Fclip/Fmantle,TIME,'k-','LineWidth',2);
set(gca,'FontWeight','bold','TickDir','out')
hold on
plot(Fclip(ismember(TIME,ControlPts(:,1)))/Fmantle,ControlPts(:,1),'ko','MarkerSize',5,'MarkerFaceColor','k')
ylabel('Time (kyr)')
xlabel({'LIP volcanic activity';'(LIP Os flux / background mantle Os flux)'})

lgd = legend([p4 p1 p3 p2],{'box-model input','box-model output','U1516','Angus'});
lgd.Location = 'southwest';

set(findall(gcf,'-property','FontWeight'),'FontWeight','bold')
set(findall(gcf,'-property','Fontsize'),'Fontsize',10)