%% Clear version of Pressure Diffusion in 2D - Italy %% Thanushika Gunatilake & Stephen A. Miller clear all close all %% Load Data Italy_AVN_EQData=csvread('Italy_AVN_EQData.csv'); save('Italy_AVN_EQData.mat','Italy_AVN_EQData'); load('Italy_AVN_EQData.mat') Lat = Italy_AVN_EQData(1:4369,2); Long= Italy_AVN_EQData(1:4369,3); M=Italy_AVN_EQData(1:4369,11); depth=Italy_AVN_EQData(1:4369,4); yr = datetime(Italy_AVN_EQData(1:4369,5),Italy_AVN_EQData(1:4369,6),Italy_AVN_EQData(1:4369,7),Italy_AVN_EQData(1:4369,8),Italy_AVN_EQData(1:4369,9),Italy_AVN_EQData(1:4369,10)); absoluteDays = days(yr-yr(1)); relativeDays = diff(yr); relativeDays=[0;days(relativeDays)]; day = relativeDays; data_days_check_plot=absoluteDays; %% Domain parameters L = 40*1000; % Length of domain [m] H = 15*1000; % Height of domain [m] nx = 300; % Number of gridpoints in x-direction ny = 300; % Number of gridpoints in y-direction dy = H/ny; % Spacing between gridpoints in y-direction dx = L/nx; % Spacing between gridpoints in x-direction [x2d,y2d] = meshgrid(0:dx:(nx-1)*dx,(ny-1)*dy:-dy:0); %Computational domain x2d = x2d'./1000; y2d = y2d'./1000; material = ones(nx,ny); % Model geometry beta_fluid = 1e-10; % Fluid compressibility [1/Pa] beta_phi = 1e-8; % Pore compressibility [1/Pa] porosity_rock = 0.03; % Rock porosity rho_fluid = 1000; % Density of fluid [kg/m3] rho_rock = 2700; % Density of rocks [kg/m3] g = 10; % Gravitational acceleration [m2/s] Geotherm_grad = 20; % Geothermal gradient [K/km] Stress_state = 2; % Stress state: 1-Extension, 2-Compression, 3-Transtension num_faults = 6; % Specify the number of discrete faults in the model (current maximum 5) Dir_name = '2D_output_1';% Directory where the data should be saved Day = 3600*24; % Number of seconds per day dt = Day; % Timestep in seconds tend = 351; % Overall sumulated time in days (%tend=1000,Save_timestep=2000;) tot_time = (tend*Day); % Total nunmber of timesteps Num_timestep = 0; % number of time steps simulated Save_timestep = 5; % 0.01*2000=20Bilder, How many timesteps are there between 2 savings? 0.05*20 = 1Bild pro Tag Save_data = true; % Save the data or not? Save_tiff_pictures = true; % Save pictures as tiff? error_level = 1e-2; error_iterations = realmax; n_iter = 0; %% Faults Faults = repmat(struct('x0', 10000, 'y0', -4000, 'yend', -11000, 'thickness', 150, 'dip', 60), num_faults, 6); Faults(1).x0 = 28*1000; % x-coordinate of faults upper limit Faults(1).y0 = -2200; % y-coordinate of faults upper limit Faults(1).yend = -12500; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(1).thickness = 0.1*1000; % Thickness of the fault [m] Faults(1).dip = 135; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] Faults(2).x0 = 14.78*1000; % x-coordinate of faults upper limit Faults(2).y0 = -4569; % y-coordinate of faults upper limit Faults(2).yend = -7300; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(2).thickness = 0.1*1000; % Thickness of the fault [m] Faults(2).dip = 60; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] Faults(3).x0 = 11.1*1000; % x-coordinate of faults upper limit Faults(3).y0 = -6326; % y-coordinate of faults upper limit Faults(3).yend = -12600; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(3).thickness = 0.1*1000; % Thickness of the fault [m] Faults(3).dip = 45; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] Faults(4).x0 = 22.4*1000; % x-coordinate of faults upper limit Faults(4).y0 = -6150; % y-coordinate of faults upper limit Faults(4).yend = -9850; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(4).thickness = 0.05*1000; % Thickness of the fault [m] Faults(4).dip = 109; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] Faults(5).x0 = 19*1000; % x-coordinate of faults upper limit Faults(5).y0 = -6250; % y-coordinate of faults upper limit Faults(5).yend = -7250; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(5).thickness = 0.08*1000; % Thickness of the fault [m] Faults(5).dip = 80; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] Faults(6).x0 = 18.3*1000; % x-coordinate of faults upper limit Faults(6).y0 = -8552; % y-coordinate of faults upper limit Faults(6).yend = -10350; % y-coordinate of faults lower limit (x-coordinate is calculated) Faults(6).thickness = 0.05*1000; % Thickness of the fault [m] Faults(6).dip = 75; % Dip of the fault. Measured anti-clock wise from the x-axis [deg] %% Boundary Conditions Bound_type_hor = {'Neumann' ,'Neumann' }; Bound_val_hor = [0 , 0]; Bound_type_ver = {'Neumann' ,'Dirichlet' }; Bound_val_ver = [0 0]; %% Material Material_mat(1,:) = [ 60 0 0.01 1e-14 35 ]; Material_mat(2,:) = [ 60 0 0.05 1e-12 100 ]; Material_mat(3,:) = [ 60 0 0.05 1e-12 100 ]; Material_mat(4,:) = [ 60 0 0.05 1e-12 100 ]; Material_mat(5,:) = [ 60 1e-6 0.05 1e-12 100 ]; Material_mat(6,:) = [ 60 1e-6 0.05 1e-12 100 ]; Material_mat(7,:) = [ 60 1e-6 0.05 1e-12 100 ]; mat_matrix=1; mat_fault2=2; mat_fault3=3; mat_fault4=4; mat_fault5=5; mat_fault6=6; mat_fault7=7; %% Create Domain for i=1:nx for j=1:ny material(i,j)=mat_matrix; end end %% Set fault indentifier for k=1:num_faults if Faults(k).y0 <= Faults(k).yend error(['The end of fault %i should be deeper than the start!',k]) end x0 = Faults(k).x0; y0 = Faults(k).y0; xend = x0-abs(Faults(k).yend+Faults(k).y0) / tan(-Faults(k).dip*pi/180.0); yend = Faults(k).yend; for i=1:nx for j=1:ny xpos = i*dx; ypos = -(ny-j)*dy; dist_to_fault = abs((yend-y0)*xpos-(xend-x0)*ypos+xend*y0-yend*x0)/sqrt((yend-y0)^2+(xend-x0)^2); if (dist_to_fault <= Faults(k).thickness && ypos >= Faults(k).yend && ypos <= Faults(k).y0) material(i,j)=1+k; end end end end P_fluid_hydro = rho_fluid*g*(y2d+.1)*1000; % Hydrostatic fluid pressure Sigma_1 = rho_rock*g*(y2d+.1)*1000; % Lithostatic rock pressure if Stress_state==1 Sigma_3 = Sigma_1; Sigma_1 = 2.*Sigma_3; elseif Stress_state==2 Sigma_3 = 0.62*Sigma_1; Sigma_1 = Sigma_1; elseif Stress_state==3 Sigma_3 = Sigma_1; Sigma_1 = 1.5*Sigma_3; else error('Stress state unknown') end %% Variables P_source = zeros(nx,ny); % Fluid source pressure field (impermeable source is handled in initial conditions) P_source2 = zeros(nx,ny); P_source3 = zeros(nx,ny); P_source4 = zeros(nx,ny); parameter P_old = zeros(nx,ny); % Solution field kc0 = ones(nx,ny); % Permeability matrix (m^2) sig_star=zeros(nx,ny); % Stress norm factor matrix Temp = -y2d./1e3*Geotherm_grad + 0; % Temperature, assuming a constant gradient and 0 degrees on top Fi = porosity_rock.*ones(nx,ny); %Porosity matrix Beta = Fi.*(beta_phi + beta_fluid);% beta coefficient used for the diffusion equation. Visc = ones(nx,ny)*10^(-3)-0.9*(y2d./max(max(y2d))).*ones(nx,ny)*10^(-3); %visc in Pa-s (20 dec C=1e-3 Pa-s) theal = zeros(nx,ny); % Healing time. Basically stores the number of timesteps since healing became active. trig_tog = zeros(nx,ny); % Stores if the point has already been triggered. plot_trig = zeros(nx,ny); %Stores the points that were triggered in the current timestep parameter; %% Mechanic (Plot initial failure condition) Sigma_normal = (Sigma_1 + Sigma_3 - 2.*(P_old + P_fluid_hydro))./2+((Sigma_1 - Sigma_3)./2).*(cos(2.*Dip_angle.*pi/180)); Tau = ((Sigma_1-Sigma_3)./2).*sin(2.*Dip_angle.*pi./180); fric_fail=0.48; %% Simulation Begin time_count=1; time = 0; k_mean = NaN(floor(tot_time/dt)+1,1); p_mean = NaN(floor(tot_time/dt)+1,1); Q = NaN(floor(tot_time/dt)+1,1); total_trig = zeros(floor(tot_time/dt)+1,1); time_trig = zeros(floor(tot_time/dt)+1,1); P_source0 = P_source; while time<tot_time count_trig=1; time =time+dt; n_iter= 0; if time_count==int16(66/dt*Day) || time_count==int16(136/dt*Day) || time_count==int16(198/dt*Day) theal = 0.*theal; theal2 = 0.*theal2; end sigma_n=(Sigma_1+Sigma_3-2.*(P_old + P_fluid_hydro))./2+((Sigma_1-Sigma_3)./2).*(cos(2.*Dip_angle.*pi/180)); tau_a=((Sigma_1-Sigma_3)./2).*sin(2.*Dip_angle.*pi./180); for i=1:nx for j=1:ny if (material(i,j)==1 || material(i,j)==2 ||material(i,j)==3 ||material(i,j)==4 ... || material(i,j)==5 || material(i,j)==6 || material(i,j)==7) if ((tau_a(i,j)./sigma_n(i,j)>fric_fail||tau_a(i,j)./sigma_n(i,j)<-fric_fail)&&trig_tog(i,j)==0) trig_tog(i,j)=1; plot_trig(i,j)= plot_trig(i,j)+1; count_trig=count_trig+1; theal(i,j)=0; end theal(i,j)=theal(i,j)+dt; theal2(i,j)=theal(i,j)/24/3600; if (trig_tog(i,j)==1) if time_count<=int16(65/dt*Day) alpha = 0.03; P_source(i,j)=P_source0(i,j).*exp(-alpha.*theal2(i,j)); kc0(i,j)=1e-14+kc0(i,j).*exp(-alpha.*theal2(i,j)); P_source2(i,j) = P_source(i,j); elseif time_count>int16(65/dt*Day) && time_count<=int16(135/dt*Day) alpha2 = 0.0567; P_source(i,j)=P_source2(i,j).*exp(-alpha2.*theal2(i,j)); kc0(i,j)=5e-14+kc0(i,j).*exp(-alpha2.*theal2(i,j)); P_source3(i,j) = P_source(i,j); else alpha4 = 0.01; P_source(i,j)=P_source3(i,j).*exp(-alpha4.*theal2(i,j)); kc0(i,j)=1e-14+kc0(i,j).*exp(-alpha4*theal2(i,j)); end end end end end K = kc0.*exp(-Sigma_normal./1e6./sig_star); k_mu = (1./Visc).*K; [P_new] = p2d(P_old, Beta, k_mu, P_source, nx, ny, dx, dy, dt, Bound_type_hor, Bound_val_hor, Bound_type_ver, Bound_val_ver); error_iterations = max(max(abs(P_old-P_new)))./max(max(abs(P_new))); dP_dt = (P_new-P_old)/(dt/Day); P_old = P_new; n_iter = n_iter + 1; Q = (P_source./Beta); disp(['The error of iterations # ', num2str(n_iter),' = ',num2str(error_iterations)]); time_count=time_count+1; end