Added Pearson Iceland Tectonics inputs

Pearson_and_Loveless_Tectonics_2021_Inputs/Supporting Functions/GPSCoordinates.m

@@ -0,0 +1,89 @@
+function [anglesbetween, misfit, results, rakes] = GPSCoordinates(faults, faultsHotspot, d, bc, azimuth, hotspotContribution, GPSData)
+%GPScomparison evaluates model results in comparison with strain axes
+%derived from GPS velocity data 
+%Inputs: 
+% faults: meshed fault map       faultsHotspot: meshed fault map and spherical mesh hotspot
+% d: nelx3 array of boundary conditions & bc: NEL-by-3 array specifying type
+% of boundary conditions, see tribemx for specifications 
+% azimuth: angle of uniaxial remote stress, defined clockwise from due east
+% hotspot Contribution: amount of element-normal dislocation on the hotspot
+% in mm/yr
+% GPSData: comparison data, including: stationsTotal (UTM coordinates of GPS stations, column 1: 
+% Eastings, column 2: Northings), GPS (trends of axes of maximum and minimum extensional strain calculated
+% directly from GPS velocities, 1:max, 2:min), and stationsList(UTM coordinates, broken down by region in order:
+% NW, CW, SW, NC, CC, SC, NE, CE, SE)
+
+
+%Outputs: 
+%anglesbetween: angular misfits between axis of maximum extension for each
+%comparison data point
+%misfit: overall average misfit for model 
+%results:  axes of maximum extension for each observation coordinate defined for
+%the model
+%rakes: overall rake calculated for each fault from the model 
+
+%set up grid of observation points at each GPS station location
+stationsTotal=GPSData{1};
+GPS=GPSData{2};
+stationsList=GPSData{3};
+xyutm=GPSData{4};
+obs.x=stationsTotal(:, 1); obs.y=stationsTotal(:, 2); obs.z=0*obs.x;
+obs.v="d";
+obs.rc=[xyutm(1,1), xyutm(1,2), 0]; %reference coordinates, close to the plate boundary
+
+%run tribemx to calculate displacements at each station
+d(sum(faults.nEl)+1:end, 3)=hotspotContribution;
+[REMSrotate, ~]=calculateREMS(azimuth);
+[slip, trac, out] = tribemx(faultsHotspot, d, bc, REMSrotate, obs);
+disp2=[out.u(:,1), out.u(:,2)];
+faultsHotspot=rake(faultsHotspot, slip);
+rakes=faultsHotspot.sliprake;
+
+
+start=0; type=1; plotpar=0;
+array=zeros(9, 1);  arrayAlt=zeros(9, 1); misfitAlt=zeros(9, 1);
+misfit=zeros(9, 1); anglesbetween=zeros(9,1); results=zeros(9,1);
+for k=1:9  %run for each of the 9 regions
+    %calculate strains based on displacement
+    stationsNew=stationsList{1, k};
+    nStations=length(stationsNew);
+    par = [500e3, nStations, 250e3];
+    [~,~,~,pstrains,~] = GridStrain(stationsNew,disp2(start+1:start+length(stationsNew), :),type,par,plotpar);
+    start=start+length(stationsNew);
+    %rename
+    e1mag = pstrains(1, 1, :); e1mag = e1mag(:);
+    e1tre = pstrains(1, 2, :); e1tre = e1tre(:);
+    e1plu = pstrains(1, 3, :); e1plu = e1plu(:);
+    e2mag = pstrains(2, 1, :); e2mag = e2mag(:);
+    e2tre = pstrains(2, 2, :); e2tre = e2tre(:);
+    e2plu = pstrains(2, 3, :); e2plu = e2plu(:);
+    e3mag = pstrains(3, 1, :); e3mag = e3mag(:);
+    e3tre = pstrains(3, 2, :); e3tre = e3tre(:);
+    e3plu = pstrains(3, 3, :); e3plu = e3plu(:);
+
+    %permutations
+    e1horizmag = e1mag; % Default assumption is that e1 horizontal is e1
+    e1horiztre = e1tre;
+    e1horizplu = e1plu;
+    vert_e1 = e1plu == pi/2; % Logical index of any e1 plunges that are vertical
+    e1horizmag(vert_e1) = e2mag(vert_e1); % Reassign those values to be e2 values
+    e1horiztre(vert_e1) = e2tre(vert_e1);
+    e1horizplu(vert_e1) = e2plu(vert_e1);
+    e3horizmag = e3mag; % Default assumption is that e3 horizontal is e3
+    e3horiztre = e3tre;
+    e3horizplu = e3plu;
+    vert_e3 = e3plu == pi/2; % Logical index of any e3 plunges that are vertical
+    e3horizmag(vert_e3) = e2mag(vert_e3); % Reassign those values to be e2 values
+    e3horiztre(vert_e3) = e2tre(vert_e3);
+    e3horizplu(vert_e3) = e2plu(vert_e3);
+    array(k, 1)=rad2deg(e1horiztre);
+    arrayAlt(k, 1)=mod(rad2deg(e1horiztre)+180, 360);
+    misfit(k, 1)=mod(abs(array(k, 1)-GPS(k, 1)), 180);
+    misfitAlt(k, 1)=mod(abs(arrayAlt(k, 1)-GPS(k, 1)), 180);
+    anglesbetween(k, 1)=min(misfit(k,1), misfitAlt(k,1));
+    results(k,1)=array(k,1);
+end
+results(:,2)=zeros(9,1);
+misfit=mean(anglesbetween);
+end
+

Pearson_and_Loveless_Tectonics_2021_Inputs/Supporting Functions/calculateREMS.m

@@ -0,0 +1,11 @@
+function[REMS, remsinitial]=calculateREMS(alpha)
+%calculate uniaxial remote stress tensor at a given azimuth, given as an
+%angle counterclockwise from due East
+rotation=zeros(2,2);
+rotation(1,1)=cosd(alpha); rotation(2,2)=cosd(alpha);
+rotation(2,1)=sind(alpha); rotation(1,2)=-sind(alpha);
+a=[6e6 0; 0 0];
+remsinitial=rotation'*a*rotation;
+REMS=zeros(1,6);
+REMS(1)=remsinitial(1,1); REMS(2)=remsinitial(2,2); REMS(4)=remsinitial(1,2);
+end

Pearson_and_Loveless_Tectonics_2021_Inputs/Supporting Functions/gmshfaults_Windows.m

@@ -0,0 +1,53 @@
+function allstring = gmshfaults_Windows(topdir, botdir, sz, path)
+% gmshfaults  Meshes faults in 3D using Gmsh.
+%   gmshfaults(BOTDIR, ELSIZE) meshes all faults whose bottom coordinates
+%   are contained in text files within the folder BOTDIR, using the 
+%   element size specified using ELSIZE.  Geometry (.geo) and mesh (.msh)
+%   files are saved to the bottom coordinate directory.
+%
+%   ALL = gmshfaults(...) returns a string to the command line that can be
+%   used to read all faults in using ReadPatches.m, i.e.
+%   >> p = ReadPatches(ALL);
+% 
+%
+
+% Hard code full path to top coordinates directory here
+%topdir = ['~' filesep, 'Desktop', filesep,'chile', filesep,'Manual_Faultloc_All', filesep,'Manual_Faultloc_Originals', filesep,'Manual_Faultloc_16km'];
+
+% List all bottom coordinate files
+files = dir([botdir filesep '*.txt']);
+
+for i = 1:size(files, 1)
+   % Read bottom coordinates
+   botcoords = load([botdir filesep files(i).name]);
+   % Get number of coordinates
+   nc = size(botcoords, 1);
+   % Read top coordinates
+   topcoords = load([topdir filesep files(i).name]);
+   % Write the Gmsh geometry file
+   geofile = [botdir filesep files(i).name(1:end-4) '.geo'];
+   fid = fopen(geofile, 'w');
+   fprintf(fid, 'charl = %g;\n', sz);
+   fprintf(fid, 'Point(%g) = {%g, %g, %g, charl};\n', [1:2*nc; [topcoords(:, 1); flipud(botcoords(:, 1))]'; [topcoords(:, 2); flipud(botcoords(:, 2))]'; [topcoords(:, 3); flipud(botcoords(:, 3))]']);
+   fprintf(fid, 'CatmullRom(%g) = {%g:%g};\n', [1:2; [1 nc+1; nc 2*nc]]);
+   fprintf(fid, 'CatmullRom(%g) = {%g,%g};\n', [3:4; [nc 2*nc; nc+1 1]]);
+   fprintf(fid, 'Line Loop(1) = {1, 3, 2, 4};\nRuled Surface(1) = {1};\n');
+   fclose(fid);
+   fid = fopen('infile.txt','rt') ;
+
+   % Mesh using Gmsh (default)
+   % system(sprintf('/Applications/Gmsh/Gmsh.app/Contents/MacOS/gmsh -2 %s -o %s.msh -v 0 > junk', geofile, geofile(1:end-4)));
+   % Mesh using Gmsh (Linux)
+%    system(sprintf('gmsh -2 %s -o %s.msh -v 0 > junk', geofile, geofile(1:end-4)));
+    % Mesh using Gmsh (Windows)
+       system(sprintf(strcat(path,' -2 %s -o %s.msh -v 0 > junk'), geofile, geofile(1:end-4)));
+
+end
+
+% Make a space-separated string of all .msh files
+alldir = dir([botdir '/*.msh']); 
+xlsfiles={alldir.name};
+[~,idx]=natsort(xlsfiles);
+alldir=alldir(idx);
+allstring = sprintf('%s ', alldir.name); 
+allstring = allstring(1:end-1); % Remove trailing space

Pearson_and_Loveless_Tectonics_2021_Inputs/Supporting Functions/makeStereonet.m

@@ -0,0 +1,30 @@
+function [] = makeStereonet(results, type)
+%makeStereonet plots an equal area stereonet, plotting axes on top 
+figure
+Stereonet(0,90*pi/180,10*pi/180,1);
+hold on
+results=deg2rad(results);
+for i=1:length(results)
+    %Plot sigma1/P axis/maximum extension (black, filled circle)
+    [xp,yp] = StCoordLine(results(i,1),results(i,2),1);
+    h2=plot(xp,yp,'ko','MarkerFaceColor','k');
+    %Plot T axis (red circle)
+    if size(results, 2)==4 
+        [xp,yp] = StCoordLine(results(i,3),results(i,4),1);
+        h3=plot(xp,yp,'ro','MarkerFaceColor','r');
+    elseif size(results, 2)==6
+        [xp,yp] = StCoordLine(results(i,3),results(i,4),1);
+        h3=plot(xp,yp,'bo','MarkerFaceColor','b');
+        [xp,yp] = StCoordLine(results(i,5),results(i,6),1);
+        h4=plot(xp,yp,'ro','MarkerFaceColor','r');
+    end
+end
+legend('Location','northwest');
+if type==0 || type==2
+    legend([h2 h3],{'P', 'T'});
+elseif type==1 || type==3
+    legend([h2 h3 h4],{'σ1','σ2', 'σ3'});
+else
+    legend(h2,{'Extension'}); 
+end
+end

Pearson_and_Loveless_Tectonics_2021_Inputs/Supporting Functions/meshFaults.m

@@ -0,0 +1,15 @@
+function [faults, faultsHotspot] = meshFaults(filename, hotspot, dips, path)
+%meshFaults takes in fault mesh and spherical hotspot mesh and outputs
+%patch for running in tribemx
+faults=load(filename);
+coords=utmconv(faults(:,[1 2]), 0);
+faults=[coords, faults(:,3)];
+save ('faultsutm.txt', 'faults', '-ascii')
+maketopbotcoords('faultsutm.txt', 10e3);
+rotatefaultcoords_dir('faultsutm_botcoords', dips);
+addpath(strcat(erase(filename, '.txt'), 'utm_botcoords'), strcat(erase(filename, '.txt'), 'utm_botcoords_truedip'), strcat(erase(filename, '.txt'), 'utm_topcoords'));
+faults=gmshfaults_Windows('faultsutm_topcoords', 'faultsutm_botcoords_truedip', 7500, path);
+faults = ReadPatches(faults);
+faults = PatchCoordsx(faults);
+faultsHotspot=patchcat(faults, hotspot);
+end