Solve_uncoupled_Fin3 - 41561_2018_200_MOESM5

%   This script uses 'Definitions_uncoupled.m’ and 

'Equations_uncoupled.m’ to

%   then solve for the time dependent solution of sulfur and carbon 

in the whole ocean

%   based on equations that have:

%   1. CaCO3 output based on a first order rate law with DIC

%   2. The Org C output value based on DIC-Corg fractionation with 

modernish

%   values of 27 permil

%   3. The reduced sulfur output is tied to the organic carbon 

output via

%   the classic Berner ratio of 1-3 wt %

%   4. The oxidized sulfur output (evaporites) is a constant (steady 

state solution based %   on a first order

%   rate law with [SO4].)

% 5. using foram carbonate +1permil as the overall ocean d13C.

%   Set the initial conditions

Definitions_uncoupled;

%   Set the initial condition of the ecf to be identical to seawater

c0 = zeros(4,1);    

%%%%%%%%%%%%%%%  Sulfur Initial Conditions   %%%%%%%%%%%%%%%%%%%%%

                   % Initial isotope ratio of ocean in permil

RatioS_0 = ((delS_0./1000)+1)*stdratio_S;  % Initial isotope ratio 

in mol_34/mol_32

          % Initial total [SO4] in ocean in mol/m^3

% Stot_0 = 2.8;

% c0(S_32) = Stot_0./(1+RatioS_0);         % Initial S_32 [SO4] of 

ocean in mol/m^3

% c0(S_34) = RatioS_0.*c0(S_32);   % Initial S-34 SO4 conc of ocean

%%%%%%%%%%%%%%%  Carbon Initial Conditions   %%%%%%%%%%%%%%%%%%%%%

delC_0 = del_C_foram_start + d13C_offset;                     % 

Initial isotope ratio of ocean in permil

RatioC_0 = ((delC_0./1000)+1)*stdratio_C;  % Initial isotope ratio 

in mol_34/mol_32

Ctot_0 = dic.*1029;          % Initial total [DIC] in ocean in mol/

m^3

c0(C_12) = Ctot_0./(1+RatioC_0);         % Initial C-12 [DIC] of 

ocean in mol/m^3

c0(C_13) = RatioC_0.*c0(C_12);   % Initial C-13 DIC conc of ocean

%%%%%%%%%%%%%%  Now start the looping for concentration and sulfur 

flux in   %%%%%%%%%%%%%%

ExpNum = 0;

% Variable for counting the number of total experiments run

count1 = 0;

global Stot_0 F_32_in del_in_S;

% for k = 2;

%     CS_red_ratio = (32./12).*k;

%     count1 = count1 +1;

%

%       burial_ratio(count1) = CS_red_ratio;

    count2 = 0;

for j = 4;

    Stot_0 = j*2.8;

    count2 = count2 +1;

sulfur_original(count2) = Stot_0;

    c0(S_32) = Stot_0./(1+RatioS_0);         % Initial S_32 [SO4] of 

ocean in mol/m^3

    c0(S_34) = RatioS_0.*c0(S_32);   % Initial S-34 SO4 conc of 

ocean

count3 = 0;

        for g=4;

            del_in_S = g;

            count3 = count3 +1;

            delS_in(count3) = del_in_S;

count4 = 0;

            for i =3.5; 

        F_32_in = 1e12.*i;

% Loop for the SO Mix value

        count4 = count4 +1;

        ExpNum = ExpNum +1;

        sulfur_in(count4) = F_32_in;

            F_34_in = ((del_in_S./1000)+1)*stdratio_S.*F_32_in;

        % eps_S_red = ((del_in_S - delS_0).*(F_32_in+ F_34_in))./

((F_C12_red_start./CS_red_ratio)+F_C12_red_start./

CS_red_ratio.*(((delS_0-30)./1000)+1).*stdratio_S);

       % k_evapnew = (F_32_in + F_34_in - F_C12_red_start./

CS_red_ratio - F_C12_red_start./CS_red_ratio.*(((delS_0+eps_0)./

1000)+1).*stdratio_S)./(Stot_0.*VolOcean);

c=ode15s(@(t,c) Equations_uncoupled(t,c),[0:3.85e5: 76e6],c0);

ConcMatrix = c.y;

        Result_uncoupled_Fin3(ExpNum).time = c.x;

        Result_uncoupled_Fin3(ExpNum).Stot_0 = Stot_0;

        Result_uncoupled_Fin3(ExpNum).F_32_in = F_32_in;

        %Result_CS_uncoupled_smo(ExpNum).CS_red_ratio = 

CS_red_ratio;

        Result_uncoupled_Fin3(ExpNum).del_in_S = del_in_S;

        Result_uncoupled_Fin3(ExpNum).S32 = c.y(1,:);

        Result_uncoupled_Fin3(ExpNum).S34 = c.y(2,:);

        Result_uncoupled_Fin3(ExpNum).C12 = c.y(3,:);

        Result_uncoupled_Fin3(ExpNum).C13 = c.y(4,:);

        Result_uncoupled_Fin3(ExpNum).Sconc = c.y(1,:) + c.y(2,:);

        Result_uncoupled_Fin3(ExpNum).S_3432 = c.y(2,:)./c.y(1,:);

        Result_uncoupled_Fin3(ExpNum).del_34 = ((c.y(2,:)./

c.y(1,:))./stdratio_S-1).*1000;

        %Result_CS_eps(ExpNum).del_pyr = ((((c.y(2,:)./c.y(1,:))./

stdratio_S-1).*1000 + eps_S_red));

        Result_uncoupled_Fin3(ExpNum).DIC_conc = c.y(3,:) + 

c.y(4,:);

        Result_uncoupled_Fin3(ExpNum).C_1312 = c.y(4,:)./c.y(3,:);

        Result_uncoupled_Fin3(ExpNum).del_13 = ((c.y(4,:)./

c.y(3,:))./stdratio_C-1).*1000;

        Result_uncoupled_Fin3(ExpNum).del_13foram = (((c.y(4,:)./

c.y(3,:))./stdratio_C-1).*1000)-d13C_offset;

        %Result_CS_uncoupled_smo(ExpNum).eps_S_red = eps_S_red;

end

end

end

%end

save Result_uncoupled_Fin3;

%

% time = c.x;

%

%

% DIC_conc = c.y(3,:) + c.y(4,:);

% C_1312 = c.y(4,:)./c.y(3,:);

% del_13 = (C_1312./stdratio_C-1).*1000;

%

%

%

% SO4_conc = c.y(1,:) + c.y(2,:);

% S_3432 = c.y(2,:)./c.y(1,:);

% del_34 = (S_3432./stdratio_S-1).*1000;

% %del_pyr = (del_34 + eps_S_red);

% %

% %

% %

% figure(1)

% subplot(1,3,1)

% plot(log10(time),SO4_conc ,'o')

% ylabel('[SO4] (mol/m^3)')

% xlabel('log years')

% subplot(1,3,2)

% plot(log10(time), del_34, '.')

% ylabel('del 34S')

% xlabel('log years')

% subplot (1,3,3)

% % plot (log10(time), del_pyr, 'x')

% % ylabel('del 34S pyrite')

% % xlabel ('log years')

%

% DIC_conc = c.y(3,:) + c.y(4,:);

% C_1312 = c.y(4,:)./c.y(3,:);

% del_13 = (C_1312./stdratio_C-1).*1000;

% %

% figure(2)

% subplot(1,2,1)

% plot(log10(time),DIC_conc ,'o')

% ylabel('[DIC] (mol/m^3)')

% xlabel('years')

% subplot(1,2,2)

% plot(log10(time), del_13, '.')

% ylabel('del 13C')

% xlabel('years')

%

%

%

%

%

% ResTime = c.y(1,end).*VolOcean./F_32_in