%Welcome to MEFISTO, the Mineralization, Earthquake, and Fluid-flow
%Integrated SimulaTOr. For more information, visit www.mefisto.org.
% FROM V2:fluid_prod_mag Controls amount of fluid production. units are kg/s/mwidthx)/m(height); scales with block volume but not porosity.
% NEW FOR V3, April 2023: fluid_prod_mean: mean fluid production, as before
% NEW FOR V3, April 2023: fluid_cycles: number of (sine-function) cycles of
% fluid pressure to go through over the entire duration, including burn-in
% NEW FOR V3, April 2023: fluid_cycle_scaler: controls amplitude of fluid
% production magnitude cycles. 1 takes fluid prod all the way to zero at
% troughs and doubles it at crests. Set between 0 and 1.
% NEW FOR V4, July 2023: Updated flow law, includes pressure-solution
% (Fisher and Hirth, 2024, Science Advances) and dislocation creep (Tokle et al., 2019,
% EPSL)
%NEW for V05: March 2024. Separate "asperitization" from
%strengthening and pore/perm reduction. This approach allows each *effect*
%of asperity growth to be modeled separately from one another and from the
%asperity growth itself. Previous A replaced by asp_norm, calculated from
%user-variables asp_time and asp_temp, which are reference values
%representing conditions that would allow a cell to reach full
%asperitization.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1. DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1a. Basic geometry variables.
% block_height: Cell height (m). Normal to interface.
% dip: Interface dip (radians).
% dt: time (s) per model increment.
% L: Dip-parallel length of interface (km).
% nx: Number of cells in the x (strike) direction.
% ny: Number of cells in the y (dip) direction.
% plate_rate: Plate convergence rate (m/yr).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1b. Model execution variables.
% aftershocks: Number of ruptures that are mapped after a large event (see
%'rupture_snap') occurs
% burn_in: Number of increments to ignore at the beginning. Intended for
%model to reach steady-state conditions before generating output.
% frame_int: Number of increments between movie frames.
% init_conditions: ['file', 'rand', or 'flat'] Controls the initial distributions of cell locations and asperitized cells.
% init_loc_max = 1; If init_conditions is 'rand', this is the maximum initial loc value.
% mapruptures: [1 or 0] Option of whether to save a map of aftershock
%ruptures.
% n_incs = Number of time increments to run through.
% rupture_snap: Earthquake magnitude above which subsequent ruptures [n =
%aftershocks] are mapped.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1c. Model physical variables.
%1c1. Elasticity
% D_phi_asp: Change in ratio of static to dynamic friction for asperities.
% kc = 3; %Spring constant of coupling springs, which connect neighboring cells.
% kl = 1; Spring constant of leading springs, which connect each cell to the subducting slab.
% max_dstrength: Max shear strength added to asperities (MPa).
% non_asp_str: Failure shear stress of non-asperity cells (MPa).
% phi [>1]: Ratio of static to dynamic friction.
% tensile_strength_factor: ratio of max_dstrength to tensile strength.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1c2. Plasticity
%F_thickness = xx %Thickness of the shear zone, for crystal-plasticity
%purposes (meters)
% Solution parameters (Fisher and Hirth, 2023)
% F_H_sol = 28000+4500ln(fluid_pressure/(200*1e6));
% F_A_sol = 0.0037*exp(F_H_sol/(R*453.*(T+273.15)));
% Pressure-solution parameters (Fisher and Hirth, 2023)
% F_C_ps = F_A_sol*exp(-F_H_sol/(R.*(T+273.15)));
% F_A_ps = 1.5e-24;
% F_Q_ps = 50;
% F_d = 10-5; %Grain size (meters)
% F_strainrate_ps = F_A_ps.*F_C_ps.*stress./F_d^3.*exp(-F_Q_ps/(R.*(T+273.15)));
%Quartz dislocation creep law parameters (Tokle, 2019)
%F_A_dc = 1.1e-12;
%F_Q_dc = 115;
%F_m = 1.2;
%F_fug_water = 5521*exp(-(31280 - 2e-5.*fluid_pressure)/(R.*(T+273.15)));
%F_strainrate_dc = F_A_dc.*stress^3.*F_fug_water^F_m.*exp(-F_Q_dc/(R.*(T+273.15)));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1c3. Hydraulics
% defperm: Default hydraulic conductivity of cells (m/s)
% minKfactor [0 1): factor by which a sealed cell's permeability is reduced (1 means reduced to zero)
% defpore = 0.01; default porosity for unhealed cells (unitless).
% minporefactor [0 1): factor by which a sealed cell's porosity is reduced
% water_comp = 4.6e-4;%per megapascal. use 4.6e-4 for water
% sea_floor [0 or 1]: option for sea-floor constant pressure boundary
% condition at top; otherwise top is no-flow boundary.
% sea_floor_pressure: sea floor pressure, MPa, for sea floor constant-pressure boundary condition
% fluid_prod_mag: Controls amount of fluid production. units are kg/s/mwidthx)/m(height); scales with block volume but not porosity.
% fluidmean: Depth of max fluid production (units of cell rows).
% fluidsd: Standard deviation of fluid production zone. Controls height of fluid production zone, which is normally distributed parallel to dip and constant parallel to strike.
% faultmap: vary hanging-wall default permeability or maximum fluid
% pressure, simulating fault barriers or valves
% fluid_cycles, fluid_cycle_scaler: vary the wavelength and amplitude of
% fluid generation over the course of the simulation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%1c4. Temperature and mineral controls.
% asp_time (seconds): Time it would take to reach complete effects of asperity
% growth for a cell at asp_temp.
% asp_temp (celsius): Temperature it would take to reach full asperity growth
% in asp_time.
% C: Asperity nucleation factor.
% Ea: Activation energy for asperity nucleation.
% G: Activation energy for strengthening.
% R: Gas constant.
% Tmin: Temperature at the top of the slab (celsius).
% Tmax [>Tmin]: Temperature at the bottom of the slab (celsius).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
cc = clock;
tic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Use this section if you wan to run MEFISTO multiple times and make a
%matrix of ouctomes
%Enter parameters to vary. Be sure to comment them out below if you want
%them to vary; otherwise they will be overwritten below.
%If you don't want anything to vary, then just enter loop sizes = 1 here
for O = 1:1
for P = 1:1
fluid_prod_mag = 10^-(O-1);
defperm = 10^(P-11);
fluid_prod_mag_stats(O,P) = fluid_prod_mag;
defperm_stats(O,P) = defperm;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Model parameters.Change variables followed by '%%%User%%%' as desired.
nx = 30; %%%User%%%
ny = 30; %%%User%%%
n_incs = 2200000;%%%User%%%
burn_in = 2000000; %%%User%%%
init_conditions = 'rand'; %%%User%%%
init_loc_max = 1; %%%User%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Simulation options
mapruptures = 1;%%%User%%%
frame_int = 50;%%%User%%%
rupture_snap =1; %%%User%%%
aftershocks = 0;%%%User%%%
aftershock_count = aftershocks;%%%User%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Spring material properties
phi = 1.5; %%%User%%%
D_phi_asp = 0; %%%User%%%
kc = 3; %%%User%%%
kl = 1; %%%User%%%
alpha = kc/kl;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Temperature and mineral parameters
non_asp_str = 1; %%%User%%%
max_dstrength = 0.01; %%%User%%%
tensile_strength_factor = 0.1; %%%User%%%
Ea = 54; %%%User%%%
C = 1e5; %%%User%%%
asp_temp = 200; %%%User%%%
asp_time = 1e7; %%%User%%% 1e7
G = 54; %%%User%%%
Tmin = 20; %%%User%%%
Tmax = 350; %%%User%%%
R = 8.314e-3;
T = repmat((Tmin:(Tmax-Tmin)/(ny-1):Tmax)',1,nx);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Model geometry parameters
dip = 10 * pi/180; %%%User%%%
L = (Tmax-Tmin)/10/sin(dip);%%%User%%%
block_height = 1; %%%User%%%
plate_rate = 0.05; %%%User%%%
dt = 100000; %%%User%%%
plate_rate_persec = plate_rate/365/24/3600;
rupture_snap = 10^((rupture_snap-2/3*log10(32e16)+10.7)*1.5);
block_area = (L/ny)^2*1e6;
time_increment = dt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Fluid parameters
defperm = 1e-9; %%%User%%%
defK = defperm/9.8; %%%User%%%
faultmap = zeros(ny,nx);
%faultmap(22,12:18) = 1;%row [ny] is at the bottom
faultKfactor = 0; %%User%%
defK = ones(ny,nx).*defK.*10.^(faultmap.*faultKfactor);
K = ones(ny,nx).*defK;
minKfactor = 0.99; %%%User%%%
defpore = 0.1; %%%User%%%
porosity = defpore.*ones(ny,nx);
minporefactor = 0.5; %%%User%%%
water_comp = 4.6e-4;
sea_floor = 1; %%%User%%%
sea_floor_pressure = 0; %%%User%%%
overburden = ones(ny,nx)*9.8*(2400-1100)/1000000*sin(dip)*L*1000/ny.*(1:ny)';
fault_ceiling = ones(ny,nx);
fault_ceiling(faultmap==1) = 0.01;
fluidmean = ny*0.8; %%%User%%%
fluidsd = 1; %%%User%%%
fluid_prod_mag = 1e-11; %%%User%%%
defmass = defpore*block_area*block_height*1100;
fluid_mass = porosity.*block_area.*block_height.*1100;%mass filling porosity initially, density 1.1 g/cc
fluid_pressure = (defpore./porosity.*fluid_mass./defmass-1)./water_comp;% + overburden*(1.1/2.4);%nonhydrostatic component; hydrostatic component omitted and added to failure crit
fluid_prod_mean = repmat(1/(sqrt(2*pi*fluidsd^2)).*exp(-((1:ny)-fluidmean).^2./(2*fluidsd^2))',1,nx)*sqrt(block_area)*block_height*fluid_prod_mag;%gaussian fluid production distribution with depth.
fluid_cycles = 8; %%User%% Number of cycles of fluid production
fluid_cycle_scaler = 0;%%User%% 0 for constant fluid production
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Flow law parameters
F_thickness = 5; %%%User%%%
F_H_sol = 28000+4500.*log(max((fluid_pressure+overburden*1100/2400)*1e6,0.001)/(200*1e6));
F_A_sol = 0.0037.*exp(-F_H_sol./(R*1000*453.*(T+273.15)));
F_C_ps = F_A_sol.*exp(-F_H_sol./(R*1000.*(T+273.15)));
F_A_ps = 1.5e-24;
F_Q_ps = 50;
F_d = 1e-5; %%%User%%%
F_A_dc = 1.1e-12;
F_Q_dc = 115;
F_m = 1.2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
switch init_conditions
case 'file'
loc = load('loc_init.csv'); %Location of each slider block relative to rigid block.
asp = logical(load('asp_init.csv')); %Logical array of whether or not a cell is an asperity.
case 'rand'
loc = init_loc_max*rand(ny,nx);
asp = logical(round(rand(ny,nx)));
case 'flat'
loc = zeros(ny,nx);
asp = false(ny,nx);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%placeholder variables
rupture_slip = zeros(1,n_incs); %Total slip in each rupture.
slip_surplus = zeros(1,n_incs); %Tracks the slip surplus (+) or deficit (-) at the start of each rupture.
nloops = 0; %Number of times through the while loop involved in each rupture.
rupture_area = zeros(1,n_incs); %Number of cells that slip in each rupture, not counting cells that slip more than once as extra cells.
rupture_map = zeros(ny,nx); %map of cells that slip during any single rupture.
fail_cells = false(ny,nx); %Initially no cells are failing; even if they are, the first run through the while loop will determine that and go from there.
cum_slip = zeros(ny,nx); %Tracks the cumulative slip on each cell.
dstrength = zeros(ny,nx); %added strength of each cell, will be based on thermal exposure since last rupture.
aspage = zeros(ny,nx);%age of each asperity
asp_level = zeros(ny,nx);%degree to which each cell is asperitized (0 to 1)
asp_norm = log(asp_time+1).*(exp(-G./(R*(asp_temp+273.15))));
loc_frame = zeros(ny,nx,n_incs/frame_int+1-burn_in); %Stores a frame of the cumulative slip every frame_int ruptures.
dstrength_frame = zeros(ny,nx,n_incs/frame_int+1-burn_in); %Make a movie of cell strength
fluid_mass_frame = zeros(ny,nx,n_incs/frame_int+1-burn_in);
fluid_pressure_frame = zeros(ny,nx,n_incs/frame_int+1-burn_in);
rupturemap_frame = zeros(ny,nx,1);
aftershock_frame = zeros(ny,nx,1);
asp_level_frame = zeros(ny,nx,1);
%fluid_loss_frame = zeros(ny,nx,1);
F_strainrate_frame = zeros(ny,nx,1);
fail = zeros(ny,nx);
rupture_init_X = zeros(1,n_incs);
rupture_init_Y = zeros(1,n_incs);
rupture_time = zeros(1,n_incs);
max_slip = 0;
max_slip_loc = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Begin simulation
for n = 1:n_incs %Run through the total number of increments.
time_increment = dt*(10^(-log10(defperm)-8))^(n<(burn_in*0.9));% This step forces the model to increase time-steps during the burn-in, so that equilibrium fluid conditions are reached.
if min(min(fluid_mass)) == 0
disp('Error: A cell has run out of fluid.');
break
end
slip_surplus(n) = sum(sum(loc)); %Slip surplus (+) or deficit (-) at the beginning of the next rupture.
if min(min(fail)) <= 0
[a,b] = min(min(fail));
[c,d] = min(min(fail,[],2));
rupture_init_X(n) = b;
rupture_init_Y(n) = d;
rupture_time(n) = time_increment/365/24/3600*(n-1);%yr units
end
while min(min(fail)) <= 0 %1e-15) %Keep running through the rupture until no cells are displaced past their point of failure.
nloops = nloops+1;
stress = loc+alpha*([loc(:,1:nx-1)-loc(:,2:nx),zeros(ny,1)]...
+[loc(1:ny-1,:)-loc(2:ny,:);zeros(1,nx)]...
+[zeros(1,nx);loc(2:ny,:)-loc(1:ny-1,:)]...
+[zeros(ny,1),loc(:,2:nx)-loc(:,1:nx-1)]);%Shear stress on each slider block
fail = dstrength+non_asp_str+0.1*(overburden-fluid_pressure+stress/sin(2*dip)*(cos(2*dip)-1))-stress;%updated V35%(0.6*(minstress-fluid_pressure)+dstrength+non_asp_str)/(1-0.6/sin(2*dip)-0.6*-cos(2*dip)/sin(2*dip))-stress;%updated failure criterion treating stress as shear stress, july 6 2019, V25; V27 changes overburden to minstress
fail_cells = fail<=1e-15; %Logical array telling which cells fail. Changed to 1e-15 to deal with rounding errors.
slip = -2*(stress-1./(phi+D_phi_asp.*asp))/(1+4*alpha); %Distance each block would slip if it were to fail.
rupture_map = or(rupture_map>0,fail_cells>0); % rupture map is a logical of each cell that has slipped in this rupture.
loc(fail_cells) = loc(fail_cells)+slip(fail_cells); %Move failing cells by their (negative) displacement to bring them to failure.
asp(fail_cells) = false; %Reset the cells that failed so they are no longer asperities.
dstrength(fail_cells) = 0; %Reset the added strength to failed cells back to zero.
porosity(fail_cells) = defpore; %Set the porosity of failed cells to the default, unhealed value.
fluid_pressure(fail_cells) = min((overburden(fail_cells)+dstrength(fail_cells)*tensile_strength_factor).*fault_ceiling(fail_cells),(defpore./porosity(fail_cells).*fluid_mass(fail_cells)./defmass-1)./water_comp); % drop pressure in failed cells according to porosity; no flow during propagation
if sea_floor == 1
fluid_pressure(1,:) = sea_floor_pressure;%keep top row at constant pressure
fluid_mass(1,:) = porosity(1,:).*block_area.*block_height.*1100;
end
aspage(fail_cells) = 0;%reset the time since asperitization of each failed cell to zero.
rupture_slip(n) = rupture_slip(n)+sum(slip(fail_cells)); %Add to the total amount of displacement during this rupture.
cum_slip(fail_cells) = cum_slip(fail_cells)+slip(fail_cells);
%cum_slip = cum_slip+fail_cells;
rupture_area(n) = sum(sum(rupture_map)); %Rupture area is the total of all cells that slipped in that rupture.
if nloops>2*nx*ny
disp('Error: Rupture not stopping')
break
end
end
if mapruptures == 1
if n > burn_in
%max_slip(n) = min(min(cum_slip)); %Used to track max slipping cell per rupture; uncomment with cum_slip as desired
if rupture_area(n) > 0
max_slip(end+1) = -min(min(cum_slip));
max_slip_loc(end+1) = find(cum_slip == min(min(cum_slip)));
if aftershock_count < aftershocks
aftershock_frame(:,:,size(rupturemap_frame,3)*aftershocks+aftershock_count) = rupture_map.*aftershock_count;
aftershock_count = aftershock_count + 1;
end
end
if -rupture_slip(n)*block_area>=rupture_snap
rupturemap_frame(:,:,size(rupturemap_frame,3)+1) = -cum_slip;
init_Xf(size(rupturemap_frame,3)) = b;
init_Yf(size(rupturemap_frame,3)) = d;
aftershock_count = 1;
end
end
end
rupture_map = zeros(ny,nx); %Clear the rupture map.
cum_slip = zeros(ny,nx);
nloops = 0;
rn = rand(size(asp)); %
asp(~asp & rn=(overburden+dstrength*tensile_strength_factor).*fault_ceiling;
fluid_loss = zeros(ny,nx);
fluid_loss(crack_cells) = fluid_mass(crack_cells) - (1+fluid_pressure(crack_cells).*water_comp).*defmass.*porosity(crack_cells)./defpore;
if n > burn_in
if mod(n,frame_int) == 0
loc_frame(:,:,n/frame_int+1-burn_in/frame_int) = loc(:,:); %Add to the cumulative slip on each cell that failed.
fluid_mass_frame(:,:,n/frame_int+1-burn_in/frame_int) = fluid_mass(:,:);
fluid_pressure_frame(:,:,n/frame_int+1-burn_in/frame_int) = fluid_pressure(:,:);
dstrength_frame(:,:,n/frame_int+1-burn_in/frame_int) = dstrength(:,:);
%fluid_loss_frame(:,:,n/frame_int+1-burn_in/frame_int) = fluid_loss(:,:);
F_strainrate_frame(:,:,n/frame_int+1-burn_in/frame_int) = F_strainrate(:,:);
asp_level_frame(:,:,n/frame_int+1-burn_in/frame_int) = asp_level(:,:);
end
end
fluid_mass(crack_cells) = fluid_mass(crack_cells) - fluid_loss(crack_cells);%leak excess fluid mass out the roof to limit pressure
loc = loc+time_increment*plate_rate_persec; %Move the driving plate to the point of the next rupture. Should always be +ive outside rupture loop.
F_H_sol = 28000+4500.*log(max(1e6,(fluid_pressure+overburden*1100/1300)*1e6)/(200*1e6));
F_A_sol = 0.0037.*exp(-F_H_sol./(R*453.*(T+273.15)));
F_C_ps = F_A_sol.*exp(-F_H_sol./(R*1000.*(T+273.15)));
F_strainrate_ps = F_A_ps.*F_C_ps.*stress./F_d.^3.*exp(-F_Q_ps./(R*1000.*(T+273.15)));
F_fug_water = 5521.*exp(-(31280 - 2e-5.*(fluid_pressure+overburden*1100/1300).*1e6)./(R*1000.*(T+273.15)));
F_strainrate_dc = F_A_dc.*stress.^3.*F_fug_water.^F_m.*exp(-F_Q_dc./(R*1000.*(T+273.15)));
F_strainrate = F_strainrate_dc + F_strainrate_ps;
loc = loc-F_strainrate.* F_thickness.* time_increment; %Flow law
stress = loc+alpha*([loc(:,1:nx-1)-loc(:,2:nx),zeros(ny,1)]...
+[loc(1:ny-1,:)-loc(2:ny,:);zeros(1,nx)]...
+[zeros(1,nx);loc(2:ny,:)-loc(1:ny-1,:)]...
+[zeros(ny,1),loc(:,2:nx)-loc(:,1:nx-1)]);%Stress on each slider block
fail = dstrength+non_asp_str+0.1*(overburden-fluid_pressure+stress/sin(2*dip)*(cos(2*dip)-1))-stress;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Wrap up
t = toc;
disp(['Finished calculations in ',num2str(t),' seconds'])
%If there is a burn-in period to remove, take out those results from the
%result vectors.
if burn_in>0
rupture_time = rupture_time(burn_in+1:end);
rupture_time = (rupture_time>0) .* (1:(size(rupture_time,2))) * dt/365/24/3600;
rupture_slip = rupture_slip(burn_in+1:end);
rupture_slip = rupture_slip(rupture_time>0);
slip_surplus = slip_surplus(burn_in+1:end);
slip_surplus = slip_surplus(rupture_time>0);
rupture_area = rupture_area(burn_in+1:end);
rupture_area = rupture_area(rupture_time>0);
rupture_init_X = rupture_init_X(burn_in-1:find(rupture_init_X,1,'last')); %did burn_in -1 to synch up with movie
rupture_init_Y = rupture_init_Y(burn_in-1:find(rupture_init_Y,1,'last'));
rupture_init_X = rupture_init_X(rupture_init_X>0);
rupture_init_Y = rupture_init_Y(rupture_init_Y>0);
rupture_time = rupture_time(rupture_time>0);
%rupture_time = rupture_time-rupture_time(1)+1;
end
Mw = 2/3*log10(-rupture_slip*block_area) + 2/3*log10(32 * 1e16) - 10.7;
Mw_finite = Mw(isfinite(Mw));
stress_drop = 32000 * rupture_slip./rupture_area./sqrt(rupture_area*block_area);%megapascales
slip_surplus = slip_surplus./(nx*ny);
Mw_stats(O,P,1:size(Mw_finite,2)) = Mw_finite;
rupture_time_stats(O,P,1:size(rupture_time,2)) = rupture_time;
rupture_slip_stats(O,P,1:size(rupture_slip,2)) = rupture_slip;
slip_surplus_stats(O,P,1:size(slip_surplus,2)) = slip_surplus;
fluid_prod_stats(O,P) = sum(sum(fluid_prod));
%fluid_loss_stats(O,P,:) = permute((sum(sum(fluid_loss_frame./time_increment))),[3 2 1]);
totmass_stats(O,P,:) = permute(sum(sum(fluid_mass_frame(2:end,:,:))),[3 2 1]);
pgrad_stats(O,P,:) = permute(sum(fluid_pressure_frame(round(fluidmean),:,:)-fluid_pressure_frame(2,:,:))*1000000/nx/(L*1000*fluidmean/ny),[3 2 1]);
porosity_frame = defpore - defpore*minporefactor*(dstrength_frame./max_dstrength);
porosity_hist_stats(O,P,:) = permute(mean(mean(porosity_frame)),[3 2 1]);
save('results.mat')
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Make plots.
if O*P == 1
%Plot the asperity area and rupture sizes over time.
figure(1)
subplot(2,1,1)
errorbar(rupture_time(rupture_time>0),Mw_finite,Mw_finite,zeros(size(Mw_finite)),'.k','CapSize',0)
ylim([min(Mw_finite)-0.5 max(Mw_finite)+0.5])
ylabel('Magnitude')
subplot(2,1,2)
plot(rupture_time,slip_surplus)
xlabel('Time (yr)')
ylabel('Average slip deficit (m)')
figure(2)
X = min(Mw_finite):0.1:max(Mw_finite);
for i = 1:size(X,2)
Y(i) = sum(Mw_finite>X(1,i));
end
plot(X,log10(Y),'Marker','o','MarkerEdgeColor',[1 0 0],'LineStyle','none')
xlabel('Magnitude')
ylabel('log_{10} Cumulative number')
axis equal
end
save(num2str(cc(1:5)),'defperm','defpore','dip','Ea','G','C','fluid_prod_mag','minKfactor','minporefactor','alpha','phi','sea_floor')