%---- Run the great oxidation model for EOS 460 ----
%
% *****NEW COMMENTS****
% Additionally, some parts of example code have changed - these are clearly
% marked.
%
% the reccomended ode solver for this model is ode15s
%
% For a model experiment that we will run, we need the option of passing
% vectors of N and r, and corresponding to time vectors tN and tr, as
% inputs to our model function.
% the first line of that function script will then be:
% function dy = gtoxmodel(t,y,N,r,tN,tr)
%
% We will then need to have our code automatically distinguish as to
% whether N and r were passed as scalars or vectors, and if vectors
% interpolate to the value at time t. The following code would do that:
% if ~isscalar(N)
% N = interp1(tN,N,t,'linear'); % tt and (supplied) NN are vectors. t is time in function.
% end
% if ~isscalar(r)
% r = interp1(tr,r,t,'linear'); % tt and (supplied) rr are vectors. t is time in function.
% end
%
% Note that if N and r are scalar, then tN and tr can simply be ignored,
% but they still need to be passed as inputs. Empty arrays could be used,
% e.g.
% tN = []; tr = [];
%
% N and r must not be global.
%
% given the extra input variables, we need a slightly different calling
% syntax for the ode solver. Use:
% [t,y] = ode23s(fun,tspan,y0,options,N,r,tN,tr);
clear all
%% std setup
% parameters not to change
global d_ar d_mt s p beta w
mia = 1.773e20; % number of moles in atmosphere
solb_o = 0.0013; % solubility of oxygen
d_ar = 100e-6*mia./solb_o;
d_mt = 2e-6*mia./solb_o;
beta = 0.001; % carbon burial fraction *****CHANGED****
w = 6e-9; % weathering parameter
s = 3.7e-5; % hydrogen escape coefficient
p = [0.0030084, -0.1655405, 3.2305351, -25.8343054, 71.5397861]; % methane oxidation paprameters
% standard versions of parameters ******CHANGED**** (not global)
Nstd = 3.75e15;
rstd = 1e11;
% model to use
fun = @gtoxmodel;
% ode options
options = []; %*****NEW*****
%% run model and plot for a few different values of N
% some suggested N values are shown.
% > make plots of the mixing ratios of methane and oxygen, and the size of
% the organic carbon reservoir
% model to use
fun = @gtoxmodel;
% parameters
Nlist = Nstd*[1/3 1 3];
r = rstd;
% N and r are constant, so can use an empty time vector
tr = [];
tN = [];
% time span
% setup figure
for ii = 1:length(Nlist)
N = Nlist(ii);
% initial conditions
M0 = mia*1e-6;
o0 = mia*1e-6;
C0 = ; % start at steady state or it will be a mess.
% model evaluation
[t,y] = ode23s(fun,tspan,y0,options,N,r,tN,tr); %****NEW****
% add to the figure
end
%% Dynamic runs to show the great oxidation
% We will run two versions, one with fixed change in the boundary
% conditions N and r, and another with noise added to them. With noise
% which is added randomly, running the model multiple times will give
% different results.
% time span
tref = 4.5e9;
tspan = [0 tref];
% some stadard parts to use for r
tr = 0:1e8:tref; % time vector
r0 = 10e11; % initial value for r
fr = -0.99; % factor by which r should change
% constant N
tN = tr;
N = Nstd*ones(size(tN));
% initial conditions
M0 = mia*1e-4;
o0 = mia*1e-6;
C0 = beta*(Nstd+r0)/w; % start at steady state or it will be a mess.
y0 = [M0,o0,C0];
%----r decreasing with time, smooth variation
r = r0*(1 + fr*tr/tref);
% model evaluation
% plot a figure
% show [M, O, C, r, N] versus time
% ----r decreasing with time, N constant, with red noise added to r and N
a = 1e-1;
r = r+ a*r.*randn(size(r));
N = N+ a*N.*randn(size(N));
% model evaluation
% plot a figure
% show [M, O, C, r, N] versus time
%% make a version of the bistability diagram
% We're looking to make something like figure 2 in the paper.
% An easy way to do this is to run the model for a sufficient length of time,
% and taking the values of state variables at the end of integration to be
% the steady state.
% suggested initial condition mixing ratios of oxygen are lowO2=1e-8, highO2=1e-2
% and correspnding methane of lowO2=1e-4, highO2=1e-6 *****NEW****
% model to use
fun = @gtoxmodel;
%----varying N----
% parameters
Nlist = logspace(14,16,30); % 30 points evenly spaced in log space between 1e14 and 1e16
r = rstd;
%----varying r----