/
fstrat_solver.m
343 lines (235 loc) · 14.5 KB
/
fstrat_solver.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
%% This code uses sulfur isotope data from ice cores to solve for the fraction of stratospheric sulfate in a sulfate peak
%% Written by Andrea Burke
%% When using cite Burke et al. (2023) "High sensitivity of summer temperatures to stratospheric sulfur
%% loading from volcanoes in the Northern Hemisphere." Proceedings of the National Academy of Sciences (PNAS).
%% Read in data table
numVars = 14; % number of columns of data to read in
varNames = {'Core','Eruption','Type','BotDepth','TopDepth','Age','Volume', 'Conc','d34S', 'd34S_err','d33S', 'd33S_err','D33S', 'D33S_err'} ;
varTypes = {'char', 'char', 'char', 'double','double','double','double','double','double','double','double','double','double','double'};
data_range = 'A4:N122';
opts = spreadsheetImportOptions('NumVariables',numVars,...
'VariableNames',varNames,...
'VariableTypes', varTypes,...
'DataRange', data_range);
imported_data = readtable('Burke_2023_PNAS.xlsx', opts);
%% Choose which event(s) to analyze and set up Monte Carlo simulation parameters
cores = {'Tunu'}; % choose from 'Tunu', 'B40', 'NGRIP'
eruptions = {'UE 1453'};% list eruption(s) to consider
niters = 10000;
thresholdMIF = 0.1; % This is the number above which a sample is considered to have stratospheric sulfate
d34trange = 15; % range of d34S of tropospheric volcanic sulfate
d34tmin = -5; % minimum value of the range of d34S of volcanic sulfate
stratmin = 0; %starting minimum value of the d34S stratospheric
stratmax = 30; %starting maximum value of the d34S stratospheric
d34stratrestricted = true; % true if the strictest restrictions on d34S stratospheric are kept such that there is a monotonic decrease for d34Sstrat. This should be set to false for Huaynaputina
%% Run analyses core by core, and first determine what the background and uncertainty for the core should be
for kk = 1:length(cores)
% find all data from a given core, and make sure it is sorted from
% deepest to shallowest
core_ind = find(strcmp(imported_data.Core(:,1),cores(kk)));
core_data = sortrows(imported_data(core_ind,:), {'TopDepth'}, 'descend');
%Determine the how variable the background is for the core by
% taking standard deviation of concentration and d34S of all background
% samples in the core
if not(strcmp(cores(kk),'NGRIP'))
all_bkgd_ind = find(strcmp(core_data.Type(:,1), 'bkgd')); %indices of all background samples
bkgd_err = std(core_data.Conc(all_bkgd_ind)); % standard deviation of concentration of all background samples in core
d34bkgd_err = std(core_data.d34S(all_bkgd_ind)); %standard deviation of d34S of all background samples in core
else
bkgd_err = 17; %standard deviation of NGRIP background samples run in STAiG lab
d34bkgd_err = 1.3; %standard deviation of NGRIP background samples run in STAiG lab
end
for jj = 1:length(eruptions) % treat each eruption event separately
% find all samples from a given event and save them in working matrix D
eruption_ind = find(strcmp(core_data.Eruption(:,1),eruptions(jj)));
D = core_data(eruption_ind, :);
%find background sample(s)
indbkgd = find(strcmpi(D.Type(:), 'bkgd'));
if indbkgd % if there are background sample(s), calculate the mean of concentration and d34S
bkgd = mean(D.Conc(indbkgd)); %average background concentration for this eruption
d34bkgd = mean(D.d34S(indbkgd)); % average d34S of background for this eruption
elseif strcmp(cores(kk),'Tunu') && strcmp(eruptions(jj), 'UE 540') % for the case of 540 in Tunu there is no appropriate background so use the background from 536
bkgd = 20.2; %concentration
d34bkgd = 14.65; %isotope value
disp('using Tunu 536 background for Tunu 540 event')
else
% if there is not a background
disp('No background for:')
disp(cores(kk))
disp(eruptions(jj))
end
indevent = find(strcmpi(D.Type(:), 'event')); % Find all of the samples classified as having volcanic sulfate (i.e.'event')
j = 1; %start at first event sample
while j <= length(indevent)
if D.D33S(indevent(j)) < thresholdMIF % if there is no measurable MIF signal
%Just calculate the fvolc and its isotopic composition
% load in the isotope data for this sample, and create Monte
% Carlo matrices
d34M = randn(niters,1)*D.d34S_err(indevent(j))/2+D.d34S(indevent(j));
d33M = randn(niters,1)*D.d33S_err(indevent(j))/2+D.d33S(indevent(j));
% Calculate the Monte Carlo matrices for the isotopic
% composition of background sulfate
d34b = randn(niters,1)*d34bkgd_err+d34bkgd;
d33b = 1000 * ((1 + d34b./1000).^0.515 - 1);
% calculate the fraction background and fraction volcanic and create Monte Carlo
% matrices
bkgdMC = randn(niters,1)*bkgd_err+bkgd;
fb = bkgdMC./D.Conc(indevent(j)); %MC matrix for background fraction
indexneg = find(fb(:,1)<0); %any negative fb limited to 0
fb(indexneg) = 0;
index = find(fb(:,1)>1); % find any impossible background fractions (greater than 1)
fb(index) = 1; % must be 1
fvolc = 1-fb; %calculate fraction volcanic
% calculate isotopic composition of volcanic sulfate
for i = 1:niters
if fvolc(i)~=0
d34volc(i) = (d34M(i) - fb(i).*d34b(i))./fvolc(i);
d33volc(i) = (d33M(i) - fb(i).*d33b(i))./fvolc(i);
D33volc(i) = (d33volc(i)/1000 - ((d34volc(i)./1000+1).^0.515 -1))*1000;
else
d34volc(i) = NaN;
d33volc(i) = NaN;
D33volc(i) = NaN;
end
end
% store median and standard deviation
D.fvolc{indevent(j)} = median(fvolc,'omitnan');
D.fvolc_err{indevent(j)} = std(fvolc,'omitnan');
D.d34volc{indevent(j)} = median(d34volc,'omitnan');
D.d34volc_err{indevent(j)} = std(d34volc,'omitnan');
D.d33volc{indevent(j)} = median(d33volc,'omitnan');
D.d33volc_err{indevent(j)} = std(d33volc,'omitnan');
D.D33volc{indevent(j)} = median(D33volc,'omitnan');
D.D33volc_err{indevent(j)} = std(D33volc,'omitnan');
j = j+1; %go to next sample
else %if there is measureable MIF signal
firststrat_ind = indevent(j); %store this index as the index for the first sample with stratospheric sulfur
% load in the isotope data for this sample, and create Monte
% Carlo matrices
d34M = randn(niters,1)*D.d34S_err(indevent(j))/2+D.d34S(indevent(j));
d33M = randn(niters,1)*D.d33S_err(indevent(j))/2+D.d33S(indevent(j));
% Calculate the Monte Carlo matrices for the isotopic
% composition of background sulfate
d34b = randn(niters,1)*d34bkgd_err+d34bkgd;
d33b = 1000 * ((1 + d34b./1000).^0.515 - 1);
%Calculate the Monte Carlo matrices for the isotopic
%composition of tropospheric volcanic sulfate
d34t = rand(niters,1)*d34trange+d34tmin;
d33t = 1000 * ((1 + d34t./1000).^0.515 - 1);
% calculate the fraction background and fraction volcanic and create Monte Carlo
% matrices
bkgdMC = randn(niters,1)*bkgd_err+bkgd;
fb = bkgdMC./D.Conc(indevent(j)); %MC matrix for background fraction
indexneg = find(fb(:,1)<0); %any negative fb limited to 0
fb(indexneg) = 0;
index = find(fb(:,1)>1); % find any impossible background fractions (greater than 1)
fb(index) = 1; % must be 1
fvolc = 1-fb; %calculate fraction volcanic
% calculate isotopic composition of volcanic sulfate
for i = 1:niters
if fvolc(i)~=0
d34volc(i) = (d34M(i) - fb(i).*d34b(i))./fvolc(i);
d33volc(i) = (d33M(i) - fb(i).*d33b(i))./fvolc(i);
D33volc(i) = (d33volc(i)/1000 - ((d34volc(i)./1000+1).^0.515 -1))*1000;
else
d34volc(i) = NaN;
d33volc(i) = NaN;
D33volc(i) = NaN;
end
end
% store median and standard deviation
D.fvolc{indevent(j)} = median(fvolc,'omitnan');
D.fvolc_err{indevent(j)} = std(fvolc,'omitnan');
D.d34volc{indevent(j)} = median(d34volc,'omitnan');
D.d34volc_err{indevent(j)} = std(d34volc,'omitnan');
D.d33volc{indevent(j)} = median(d33volc,'omitnan');
D.d33volc_err{indevent(j)} = std(d33volc,'omitnan');
D.D33volc{indevent(j)} = median(D33volc,'omitnan');
D.D33volc_err{indevent(j)} = std(D33volc,'omitnan');
l = randn(niters,1)*0.006+0.608; %allow there to be uncertainty on lambda
[solutions] = fstrat_MC(d34M, d33M,fb, d34b,d33b, d34t, d33t, l,stratmin, stratmax);
D.solutions{indevent(j)} = solutions;
D.fstrat{indevent(j)} = median(solutions(:,1));
D.fstrat_err{indevent(j)} = std(solutions(:,1));
D.d34Sstrat{indevent(j)} = median(solutions(:,2));
D.d34Sstrat_err{indevent(j)} = std(solutions(:,2));
j = length(indevent) +1; % get out of the while loop
end
end
%Now go through the rest of the event samples
for ii = firststrat_ind+1:indevent(end)
MCsolutions = datasample(D.solutions{ii-1},niters);
% load in the isotope data for this sample, and create Monte
% Carlo matrices
d34M = randn(niters,1)*D.d34S_err(ii)/2+D.d34S(ii);
d33M = randn(niters,1)*D.d33S_err(ii)/2+D.d33S(ii);
D33M = (d33M./1000 - ((d34M./1000+1).^0.515-1))*1000; %calculate the D33 of each MC iteration
% Calculate the Monte Carlo matrices for the isotopic
% composition of background sulfate
d34b = randn(niters,1)*d34bkgd_err+d34bkgd;
d33b = 1000 * ((1 + d34b./1000).^0.515 - 1);
%Calculate the Monte Carlo matrices for the isotopic
%composition of tropospheric volcanic sulfate
%d34t = MCsolutions(:,3); %pass through d34S tropospheric options consistent with previous sample solutions
d34t = rand(niters,1)*d34trange+d34tmin;
d33t = 1000 * ((1 + d34t./1000).^0.515 - 1);
% calculate the fraction background and fraction volcanic and create Monte Carlo
% matrices
bkgdMC = randn(niters,1)*bkgd_err+bkgd;
fb = bkgdMC./D.Conc(ii); %MC matrix for background fraction
indexneg = find(fb(:,1)<0); %any negative fb limited to 0
fb(indexneg) = 0;
index = find(fb(:,1)>1); % find any impossible background fractions (greater than 1)
fb(index) = 1; % must be 1
fvolc = 1-fb; %calculate fraction volcanic
% calculate isotopic composition of volcanic sulfate
for i = 1:niters
if fvolc(i)~=0
d34volc(i) = (d34M(i) - fb(i).*d34b(i))./fvolc(i);
d33volc(i) = (d33M(i) - fb(i).*d33b(i))./fvolc(i);
D33volc(i) = (d33volc(i)/1000 - ((d34volc(i)./1000+1).^0.515 -1))*1000;
else
d34volc(i) = NaN;
d33volc(i) = NaN;
D33volc(i) = NaN;
end
end
% store median and standard deviation
D.fvolc{ii} = median(fvolc,'omitnan');
D.fvolc_err{ii} = std(fvolc,'omitnan');
D.d34volc{ii} = median(d34volc,'omitnan');
D.d34volc_err{ii} = std(d34volc,'omitnan');
D.d33volc{ii} = median(d33volc,'omitnan');
D.d33volc_err{ii} = std(d33volc,'omitnan');
D.D33volc{ii} = median(D33volc,'omitnan');
D.D33volc_err{ii} = std(D33volc,'omitnan');
l = MCsolutions(:,4); %pass through lamda options consistent with previous sample solutions
if ~strcmp(eruptions(jj), 'Huaynaputina') %for non-Huaynaputina samples
stratmax_input = MCsolutions(:,2);%each subsequent sample must have a d34Sstrat lower than the previous sample
stratmin_input = -60*ones(niters,1); % the minimum basic constraint on d34S strat is -60 per mil
if d34stratrestricted %if want to include strictest restrictions on d34S strat
negindices = find(D33M <= 0); % indices if the measured D33S of the Monte Carlo simulation is negative
posindices = find(D33M >0); % % indices if the measured D33S of the Monte Carlo simulation is positive
stratmin_input(posindices,1) = MCsolutions(posindices,3); %if D33S is positive the minimum d34S can be is the starting value of d34S tropospheric
stratmax_input(negindices,1) = min(MCsolutions(negindices,2),MCsolutions(negindices,3)); %if D33S is negative the maximum d34S must be lower than both the starting tropospheric sulfate and the previous d34S stratospheric sample
end
else
stratmax_input = stratmax; %for Huaynaputina, because two eruptions involved, we relax this constraint and keep the max to 30 per mil
stratmin_input = -60;
end
[solutions] = fstrat_MC(d34M, d33M,fb, d34b,d33b, d34t, d33t, l,stratmin_input, stratmax_input);
D.solutions{ii} = solutions;
D.fstrat{ii} = median(solutions(:,1));
D.fstrat_err{ii} = std(solutions(:,1));
D.d34Sstrat{ii} = median(solutions(:,2));
D.d34Sstrat_err{ii} = std(solutions(:,2));
end
eruptionsplit = split(eruptions{jj});
if d34stratrestricted
filename = strcat(cores(kk),'_',eruptionsplit(end),'_restricted');
else
filename = strcat(cores(kk),'_',eruptionsplit(end));
end
save(filename{1},'D')
end
end