%Cellular automata / slider-block model of earthquake rupture with %asperities. Translated from John Hooker's Excel code and based on the %model of Huang et al. (1992). Coded into Matlab by David Oakley with %extensions (v4 onward) by JH. %This version reverses the dip direction and so has normal-sense %displacement %LANF v02 created July 2023 to incorporate flow law %LANF v03 created march 2024 to separate "asperitization" from %strengthening and pore/perm reduction. Without this, if it takes a longer %time for an asperity to strengthen, then it also takes a longer time for %it to plug up pore-perm %v03 incorporates flow law appropriate for anhydrite (Hooker et al, 2024, %submitted) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %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. Row 1 at the bottom for % normal fault version % 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 %% %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 for P = 1 %fluid_prod_mag = 10^-(O+3); %defperm = 10^-(P+5); %fluid_prod_mag_stats(O,P) = fluid_prod_mag; %defperm_stats(O,P) = defperm; %% %Model parameters. nx = 100; %Number of cells in the x (strike) direction. ny = 30; %Number of cells in the y (dip) direction. n_incs = 2000000;%10100000; %Number of increments to run through. burn_in = 1000000; %Number of increments to ignore at the beginning. Useful for 'rand' or 'flat' initial conditions. init_conditions = 'rand'; %Options for this are 'file','rand',or 'flat'. init_loc_max = 1; %If init_condistions is 'rand', this is the maximum initial loc value. max_slip = 0; max_slip_loc = 0; %% %%%Simulation options mapruptures = 1; frame_int = 2000; %Number of increments between movie frames. rupture_snap =6; % Magnitude above which rupture is mapped aftershocks = 0; aftershock_count = aftershocks; %%%Spring material properties phi = 1.05; %Ratio of static to dynamic friction. default 1.5 D_phi_asp = 0; %Change in ratio of static to dynamic friction for asperities. For V6, now applied to all cells, so consider that if we start using it. kc = 5; %Coupling spring constant (non-asperity). kl = 1; %Leading spring constant. alpha = kc/kl; %Ratio of coupling and leading spring constants. %%%Temperature (and asperity) parameters non_asp_str = 1; %Failure stress of non-asperity cells. megapascals max_dstrength = 0.01; %Max strength added to asperities. megapascals tensile_strength_factor = 0.1; %factor to multiple max_dstrength by to get tensile strength R = 8.314e-3; %Gas constant. Ea = 54; %Activation energy for nucleation. C = ones(ny,nx)*1e-2; %Mineral mucleation rate factor G = 48; %Energy term for mineral strengthening. asp_temp = 100;%temper (C) for full asperitization asp_time = 1e8;%time for full asperitization, s (3e7 s/year) Tmin = 25; %Temperature at the top of the fault. Tmax = 400; %Temperature at the bottom of the fault. T = repmat((Tmax:-(Tmax-Tmin)/(ny-1):Tmin)',1,nx); %Array of temperatures, constant along strike of the slab but decreasing down dip. %%%Model geometry parameters L = 60; %Length of interface in km. Default from Syracuse et al., 2010, = 225C/130km dip = 20 * pi/180; %interface dip in radians (enter degrees as first term) (used to determine overburden) block_height = 1; %block height in meters block_area = (L/ny)^2*1e6; %Area of one block, in square meters. plate_rate = 0.001; %Plate motion rate in m/yr. plate_rate_persec = plate_rate/365/24/3600;%Plate motion rate in m/s rupture_snap = 10^((rupture_snap-2/3*log10(32e16)+10.7)*1.5); dt = 100000;% time (s) per model increment. 10^6 s is about 11.5 days time_increment = dt; %%%Fluid parameters defperm = 1e-6;%in m/s defK = defperm/9.8;%convert perm into seconds for calc faultmap = zeros(ny,nx);%% implement heterogeneity in the hanging wall faultmap(15,15:35) = 1; faultmap(15,65:85) = 1; faultKfactor = 0; %factor of 10 by which fault permeability is increased (decreased for negative number) defK = ones(ny,nx).*defK.*10.^[faultmap.*faultKfactor];% UN/COMMENT TO MAKE DEFK ACCOUNT FOR FAULT K = defK; minKfactor = 0.99; %factor by which a sealed cell's permeability is reduced (1 means reduced to zero) defpore = 0.3; %default porosity for unhealed cells porosity = defpore.*ones(ny,nx); minporefactor = 0; % 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 = 1;%1 = sea-floor constant pressure boundary condition at top; 0 otherwise sea_floor_pressure = 0; %sea floor pressure, MPa, for sea floor boundary condition overburden = ones(ny,nx)*9.8*(2400-1100)/1000000*sin(dip)*L*1000/ny.*(ny:-1:1)';% effective overburden fault_ceiling = ones(ny,nx); fault_ceiling(faultmap == 1) = 0.01; fluidmean = 1; %ny*0.1; %Depth of max fluid production fluidsd = 1; %Controls height of fluid production zone fluid_prod_mag = 1e-2; %Controls amount of fluid production. %units are kg/s/mwidthx)/m(height); scales with block volume but not porosity 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; fluid_cycle_scaler = 0;%0 for constant fluid production %%%Flow law parameters F_thickness = 1; %%%User%%% F_H = -27; F_M = 5e-7.*exp(0.0094*(overburden-fluid_pressure)); F_S = F_M.*exp(-F_H./(R.*(T+273.15))); F_E = 30; F_A = ones(ny,nx)*1e-17; F_d = 1e-5; %%%User%%% Grain size in m %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%Initialize arrays that need to be initialized. 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); %fluid_mass = importdata('init_mass.mat'); end %% %placeholder variables %asp_area = zeros(1,nruptures); %Asperity area in each step. rupture_slip = zeros(1,n_incs); %Total slip in each rupture. %max_slip = zeros(1,n_incs); %maximum slip on a cell for each rupture. %plate_displ = zeros(1,n_incs); %Total plate displacement between ruptures. 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. %next_fail_disp = 0; %Initialize it to zero. 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); %fluid_loss_frame = zeros(ny,nx,1); fail = zeros(ny,nx); F_strainrate_frame = zeros(ny,nx,1); rupture_init_X = zeros(1,n_incs); rupture_init_Y = zeros(1,n_incs); %rupture_init_Xn = zeros(1,n_incs);%location of each rupture surface normalized to the bounds of the rupture %rupture_init_Yn = zeros(1,n_incs); rupture_time = zeros(1,n_incs); %% %Begin simulation for n = 1:n_incs %Run through the total number of increments. time_increment = dt*(10^(-log10(defperm)-6))^(n<(burn_in*0.9));%defined burn_in_dt_factor, see notes on may 28 2020 %time_increment = dt*10000^(n<2000000); if min(min(fluid_mass)) == 0 disp('Error: A cell has run out of fluid.'); break end % if n ~= 1 slip_surplus(n) = sum(sum(loc)); %Slip surplus (+) or deficit (-) at the beginning of the next rupture. %slip_surplus(n) = loc(fluidmean,nx/2);%FEB 2020: tracking a single cell at L:R center and height of max production %rupture slip is always negative, so we use -rupture_slip to make a surplus of fault slip positive. 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 %minstress = overburden - 2*stress/sin(2*dip); 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 %next_fail_disp = min(min(fail)); %Plate displacement needed for the next failure. 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(ny,:) = sea_floor_pressure;%keep top row at default fluid_mass(ny,:) = porosity(ny,:).*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); %adds to the cumulative slip of all cells that failed. %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 if rupture_area(n) > 0 %max_slip(n) = min(min(cum_slip)); %Used to track max slipping cell per rupture; uncomment with cum_slip as desired max_slip(end+1) = -min(min(cum_slip)); max_slip_loc(end+1) = find(cum_slip == min(min(cum_slip))); end if -rupture_slip(n)*block_area>=rupture_snap rupturemap_frame(:,:,size(rupturemap_frame,3)+1) = -cum_slip;%rupture_map.*-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; %V33 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(:,:); F_strainrate_frame(:,:,n/frame_int+1-burn_in/frame_int) = F_strainrate(:,:); %fluid_loss_frame(:,:,n/frame_int+1-burn_in/frame_int) = fluid_loss(:,:); % particles_frame(:,1,n/frame_int+1-burn_in/frame_int) = particles(:,1); % particles_frame(:,2,n/frame_int+1-burn_in/frame_int) = particles(:,2); 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_M = 5e-7.*exp(0.0094*(overburden-fluid_pressure)); F_S = F_M.*exp(-F_H./(R.*(T+273.15))); F_strainrate = F_S.*F_A.*stress./F_d.^3./(T+273.15).*exp(-F_E./(R.*(T+273.15))); 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;%updated failure criterion, V35 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_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-burn_in*dt/3600/24/365; 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);%*mperunit;%m per block; REMOVE save('results.mat') end end %end %% %Make plots. if O*P == 1 %Plot the asperity area and rupture sizes over time. figure(1) subplot(2,1,1) %plot(rupture_time(rupture_time>0),Mw_finite,'Marker','.','MarkerFaceColor',[0 0 0],'LineStyle','none') 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]) %xlabel('Cumulative Plate Displacement (m)') ylabel('Magnitude') subplot(2,1,2) %plot(rupture_time(rupture_time>0),slip_surplus(rupture_time>0)) 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')