%% CrustMaker Instructions
% ========================================================================
% 0 - Build the base data files
% A)METAMORPHIC/PERPLE_X DATA CODE: run Perple_X under a given set of
% conditions to calculate expected Vp, Vs, etc. for a given rock
% composition at a range of temperatures and pressures.
%
% B)'DataS' BUILDING CODE: explains structure and formatting of
% seismic inputs
%
% C)SUBSETS AND FILTERING: Select a "subset" of X rock samples (default is 500) because it
% becomes very time consuming to run CrustMaker on a geochemical set
% much larger than that. This code will create a vector Matlab file
% that contains X indices (numbers) which are randomly selected. The
% randomly selected numbers correspond to randomly selected samples
% from your desired sample set. The data is then "filtered" based on
% desired criteria to get rid of samples that would be considered
% identical within uncertainty, so as not to bias the data towards a
% formation that was simply oversampled. Criteria for filtering is
% based on differences in composition (major oxide content) and
% Perple_X-calculated Vp, Vs, and Vp/Vs ratio at each T/P interval. The
% Matlab file created will be a 1x200 structure in the form of
% 'todelete(1:100).depth(1:100).samp(1:length(full sample set))' logical.
% Samples given a "1" are "duplicates" and will be deleted in the final
% CrustMaker considerations.
%
% D)OVERLAP CALCULATIONS: Generates the "S1 X regime X" files used by CrustMaker.
%
%
% ========================================================================
%
% == Data files used: ==
% dataStecton21.mat - seismic profiles for each location/tectonic regime
% MajorOxidesGran20a.mat - oxide data for geochemical samples of interest
% GranulitesAllTemps20a.mat - Perple_X calculated seismic velocities for
% said oxide data
% == Data files created: ==
% silicasubsetgran_500.mat - indexing for a random subset of samples, should
% you wish to downsize your geochemical sample set
% Sdelgran3_500.mat - an indexing file which tells CrustMaker which samples
% are indistinguishable from each other and shouldn't be counted twice
% == Non-Matlab functions used: ==
% pressureconvert.m
% tconvert.m
% pxtpinvert.m
%
%% (A) METAMORPHIC/PERPLE_X DATA CODE
% Run Perple_X. This program will calculate
% various physical properties and mineralogies expected under various
% pressure and temperature conditions given a bulk rock composition.
%% (B) 'DataS' BUILDING CODE
%{
Your "dataS" file should contain a structure of the format:
dataS(n).file.(physicalproperty_pdf), where "n" is the number of stations or tectonic
provinces you are considering, and "physical property_pdf" is an ensemble
of 1-D seismic profiles to produce an average seismic distribtion for a
given station or tectonic province. Examples are "Vp_pdf", "Vs_pdf", and
"ratio_pdf" for the Vp, Vs, and Vp/Vs distributions. Dimensions should be
depth x iterations, so if your profiles start at 1 km depth and reach
60 km depth, and you have 4000 iterations, then your
"dataS(n).file.(physicalproperty_pdf)" variable should be a 60x4000 matrix.
For the global CrustMaker model, we have produced the "dataSbuild" script,
which draws on the Mooney 2015 geophysical database to produce average Vp, Vs, and
Vp/Vs distributions for different tectonic provinces.
%}
load('dataStecton21.mat')
%% (C) SUBSETS AND FILTERING
% Make a subset of geochemical sample data by choosing the metamorphic grad
% and a) randomly selecting 500 samples or b) selecting up to 500 samples
% based on location. 500 seems to be a sufficient number of samples to
% capture data distribution (mean and median stable and not significantly
% different from full distribution mean and median) without overloading the
% system and having the calculations take too long.
%=========================================================================%
%%% -----------------------------Subsets------------------------------ %%%
% Make the subset defining vector
today = date;
%--------- Choosing the sample type ---------%
%---- Granulite ----%
load('MajorOxidesGran20a.mat')
Granulites = MajorOxidesGran20a;
k = 500;
a = randsample(length(Granulites),k);
save(['silicasubsetgran_500_',today],'a','-v7.3')
%---- Amphibolite ----%
load('MajorOxidesAmph20a.mat')
k = 500;
a = randsample(length(MajorOxidesAmph20a),k);
save(['silicasubsetamph_500_',today],'a','-v7.3')
%%% --------------------- Filtering for deletion --------------------- %%%
%---- Granulite ----%
load(['silicasubsetgran_500_',today])
load('GranulitesAllTemps20a.mat')
oxides = Granulites; % change this based on which rocks you're using
data = oxides;
n = length(data(:,1));
distancew = NaN(n,n);
tic
% Filter based on distance between samples in oxide space
for i = 1:length(data(1:n,1))
for j = 1:length(data(1:n,1))
a1 = data(i,1);
b1 = data(i,2);
c1 = data(i,3);
d1 = data(i,4);
e1 = data(i,5);
f1 = data(i,6);
g1 = data(i,7);
h1 = data(i,9);
a2 = data(j,1);
b2 = data(j,2);
c2 = data(j,3);
d2 = data(j,4);
e2 = data(j,5);
f2 = data(j,6);
g2 = data(j,7);
h2 = data(j,9);
distancew(i,j) = sqrt((a1-a2)^2 + (b1-b2)^2 + (c1-c2)^2 + (d1-d2)^2 + (e1-e2)^2 + (f1-f2)^2 + (g1-g2)^2 + (h1-h2)^2);
end
end
toc
distancew(distancew==0)=10;
distancew(distancew<10)=NaN; % mark as identical if total distance between any two samples in oxide space is < 3 wt.%
oxdel = sum(distancew,2);
todelete = [];
counter = 0;
% Filter based on if velocities or ratios are indistinguishable within
% uncertainty at each temparture (1:200 standard) and pressure (1:100
% standard) interval
velocity = GranulitesAllTemps20a;
for R = 1:200
for depth = 1:100
A1 = round(velocity(R).file.Vs_pdf(depth,:),2); % 0.01 km/s uncertainty
A2 = round(velocity(R).file.Vp_pdf(depth,:),2); % 0.01 km/s uncertainty
A3 = round(velocity(R).file.ratio_pdf(depth,:),3); % 0.001 km/s uncertainty
A4 = ~isnan(oxdel);
[C1,IA1,IC1] = unique(A1,'stable');
[C2,IA2,IC2] = unique(A2,'stable');
[C3,IA3,IC3] = unique(A3,'stable');
todelete(R).depth(depth).samp = ~ismember(1:length(A1),[IA1; IA2; IA3; A4]);
end
end
save(['Sdelgran3_500_',today],'todelete','-v7.3')
%---- Amphibolite ----%
load(['silicasubsetamph_500_',today])
load('AmphibolitesAllTemps20a.mat')
oxides = MajorOxidesAmph20a; % change this based on which rocks you're using
data = oxides;
n = length(data(:,1));
distancew = NaN(n,n);
tic
% Filter based on distance between samples in oxide space
for i = 1:length(data(1:n,1))
for j = 1:length(data(1:n,1))
a1 = data(i,1);
b1 = data(i,2);
c1 = data(i,3);
d1 = data(i,4);
e1 = data(i,5);
f1 = data(i,6);
g1 = data(i,7);
h1 = data(i,9);
a2 = data(j,1);
b2 = data(j,2);
c2 = data(j,3);
d2 = data(j,4);
e2 = data(j,5);
f2 = data(j,6);
g2 = data(j,7);
h2 = data(j,9);
distancew(i,j) = sqrt((a1-a2)^2 + (b1-b2)^2 + (c1-c2)^2 + (d1-d2)^2 + (e1-e2)^2 + (f1-f2)^2 + (g1-g2)^2 + (h1-h2)^2);
end
end
toc
distancew(distancew==0)=10;
distancew(distancew<3)=NaN; % mark as identical if total distance between any two samples in oxide space is < 3 wt.%
oxdel = sum(distancew,2);
todelete = [];
counter = 0;
% Filter based on if velocities or ratios are indistinguishable within
% uncertainty at each temparture (1:200 standard) and pressure (1:100
% standard) interval
velocity = AmphibolitesAllTemps20a;
for R = 1:200
for depth = 1:100
A1 = round(velocity(R).file.Vs_pdf(depth,:),2); % 0.01 km/s uncertainty
A2 = round(velocity(R).file.Vp_pdf(depth,:),2); % 0.01 km/s uncertainty
A3 = round(velocity(R).file.ratio_pdf(depth,:),3); % 0.001 km/s uncertainty
A4 = ~isnan(oxdel);
[C1,IA1,IC1] = unique(A1,'stable');
[C2,IA2,IC2] = unique(A2,'stable');
[C3,IA3,IC3] = unique(A3,'stable');
todelete(R).depth(depth).samp = ~ismember(1:length(A1),[IA1; IA2; IA3; A4]);
end
end
save(['Sdelamph3_500_',today],'todelete','-v7.3')
%% (D) OVERLAP CALCULATIONS FOR VARIOUS REGIONS
% This only needs to be run if you haven't run in for a certain
% location/crust type before.
%
% Overlap calculations calculate the overlap at each P and T of the
% distribution of a velocity or Vp/Vs (or density, hypothetically) at each
% station. This only needs to be done once for each metamorphic x seismic
% model.
%
% Global probability overlap calcs have 8 "stations" each, corresponding to
% 8 different crustal types:
%
%
% For regions with high resolution data, you must do these overlap
% calculations for each station.
%=========================================================================%
% This is the baseline code that is used in the overlap calculations.
%-------------------------------------------------------------------------%
load('dataStecton21.mat')
load('AmphibolitesAllTemps20a.mat')
seis_set = dataS;
met_set = AmphibolitesAllTemps20a;
% Calculating Vs
param = 1;
S1vs = pxtpinvert(met_set,seis_set,param);
% Calculating Vp
param = 2;
S1vp = pxtpinvert(met_set,seis_set,param);
% Calculating Vp/Vs
param = 3;
S1ratio = pxtpinvert(met_set,seis_set,param);
for i = 1:length(seis_set)
S1.station(i).vp = S1vp.station(i).vp;
S1.station(i).vs = S1vs.station(i).vs;
S1.station(i).ratio = S1ratio.station(i).ratio;
end
save('S12021amph','S1','-v7.3')