; this results are very satisfactory and there is no pumping. ; compared to Nerys' pump2 I have halved the veloc, increased ngw from 1 to 10, and increased discharge to 1e-6. * Model of fault pumping with lithostatic pressure + movie + Save files + fluid historys + gwflow histories+pumping fish file. * fluid pressure evolution model * gyfld2istry added to arc2526.dat * FLAC 3.22 DEC alpha * new ti pump2.dat * config gwflow extra 15 grid 100 100 * ge -3500,0 -3500,40 950,40 950,0 rat=0.85,1 i=1,25 j=1,2 ge -3500,40 -3500,1830 950,1830 950,40 rat=0.85,.982 i=1,25 j=2,51 ge -3500,1830 -3500,4960 950,4960 950,1830 rat=0.85,1.038 i=1,25 j=51,100 ge -3500,4960 -3500,5000 950,5000 950,4960 rat=0.85,1 i=1,25 j=100,101 ge 950,0 950,40 2050,40 2050,0 i=25,77 j=1,2 ge 950,40 950,1830 2050,3170 2050,40 rat=1,1 i=25,77 j=2,51 ge 950,1830 950,4960 2050,4960 2050,3170 rat=1,1 i=25,77 j=51,100 ge 950,4960 950,5000 2050,5000 2050,4960 rat=1,1 i=25,77 j=100,101 ge 2050,0 2050,40 6500,40 6500,0 rat=1.18,1 i=77,101 j=1,2 ge 2050,40 2050,3170 6500,3170 6500,40 rat=1.18,0.96339 i=77,101 j=2,51 ge 2050,3170 2050,4960 6500,4960 6500,3170 rat=1.18,1.0183 i=77,101 j=51,100 ge 2050,4960 2050,5000 6500,5000 6500,4960 rat=1.18,1 i=77,101 j=100,101 model mohr-coulomb model ubi i=35 66 j=50,51 def mesh_fis loop i (25,76) rat1 = (i-25.0)/51.0 yrat = 0.982 - (rat1*0.01861) rrat1 = (76 - i)/51.0 yyrat = 1.0183 + (rrat1*0.0197) xl1=x(i,2) yl1=y(i,2) xr1=x(i+1,2) yr1=y(i+1,2) xl2=x(i,51) yl2=y(i,51) xr2=x(i+1,51) yr2=y(i+1,51) xxl1=x(i,51) yyl1=y(i,51) xxr1=x(i+1,51) yyr1=y(i+1,51) xxl2=x(i,100) yyl2=y(i,100) xxr2=x(i+1,100) yyr2=y(i+1,100) ii1 = i ii2 = i+1 command gen xl1,yl1 xl2,yl2 xr2,yr2 xr1,yr1 rat=1,yrat i=ii1,ii2 j=2,51 gen xxl1,yyl1 xxl2,yyl2 xxr2,yyr2 xxr1,yyr1 rat=1,yyrat i=ii1,ii2 j=51,100 end_command end_loop end mesh_fis ini y add -20e3 ; properties ; limestone (E = 7e10 v=0.25) pro bul 4.6667e10 shear 2.8e10 coh 2.0e7 ten 2e6 prop dens 2500 ;prop shear 1e9 bulk 1.6666e9 dens 2700 prop porosity 0.3 perm 1e-10 prop dil 2 fric 30 ; impermeable layer def depth depth1 = -17.2e3 depth2 = -17.8e3 loop i (1,100) loop j (1,100) if y(i,j) < depth1 then if y(i,j) > depth2 then k11(i,j) = 1.0e-15 k22(i,j) = 1.0e-15 friction(i,j) = 35 end_if end_if end_loop end_loop end depth * ;table 1 -0.5,1e-10 -0.1,1e-10 ;table 1 0.0,1e-9 0.2,1e-8 0.5,1e-8 ; fault as a zone with ubiquitous joints pro jfric 5 jcoh 5e6 ja 50.63 jtens 2e6 i=35,66 j=50,51 prop perm 1e-15 i=35,66 j=50,51 ;prop per_table=1 i=35,66 j=50,51 water bulk 2e9 dens 1000 * set large * fix y j=1 fix x i=1 fix x i=101 fix y j=101 * set gravity 9.81 * ;free pp i=1,101 j=101 ;free pp i=1,101 j=1 ;free pp i=1 j=1,101 ;free pp i=101 j=1,101 * apply syy = -3.678e8 j= 101 ***** ini sxx=-4.90400e8 var=0, 1.226e8 i=1, 101 j=1,101 ini syy=-4.90400e8 var=0, 1.226e8 i=1, 101 j=1,101 ini szz=-4.90400e8 var=0, 1.226e8 i=1, 101 j=1,101 ini sxy=0e6 ; def hydro_pp loop i (1,igp) loop j (1,jgp) if j = jgp Then sss2 = 1.4715e8 else aaa = abs(y(i, j) - y(i, jgp)) sss2 = 1000 * 9.8 * aaa+1.4715e8 ;for 15km end_if gpp(i, j) = sss2 end_loop end_loop end hydro_pp * ; ini pp=1.4715e8 j=101 fix pp j=101 ;plo hold pp fill his nstep 1 his unbal set small set mech on flow off step 200 set mech off flow on step 200 set mech on flow on step 200 set large free y j=101 * ini sat 1 fix sat * ;pore pressure histories are for zones ;hist 2 top of fault his pp i 66 j 51 ;hist 3 right side of top of fault his pp i 73 j 37 ;hist 4center of fault his pp i 50 j 51 ;hist 5 center of imperm layer right side of fault his pp i 72 j 34 ;hist 6 bottom right of fault his pp i 73 j 18 ;hist 7 bottom of fault his pp i 31 j 51 ;hist 8 bottom left of fault his pp i 14 j 38 ;hist 9 bottom left hand corner his pp i 4 j 10 ;hist 10 bottom right hand corner his pp i 99 j 7 ;hist 11 under imperm layer rhside his pp i 1 j 58 ;hist 12 under imperm layer lhside his pp i 100 j 21 * ;perm hists ;hist 13 top of fault his ex_3 i 66 j 51 ;hist 14 right side of top of fault his ex_3 i 73 j 37 ;hist 15 center of fault his ex_3 i 50 j 51 ;hist 16 center of imperm layer right side of fault his ex_3 i 72 j 34 ;hist 17 bottom right of fault his ex_3 i 73 j 18 ;hist 18 bottom of fault his ex_3 i 31 j 51 ;hist 19 bottom left of fault his ex_3 i 14 j 38 ;hist 20 bottom left hand corner his ex_3 i 4 j 10 ;hist 21 bottom right hand corner his ex_3 i 99 j 7 ;hist 22 under imperm layer rhside his ex_3 i 1 j 58 ;hist 23 under imperm layer lhside his ex_3 i 100 j 21 ;vsi hists ;hist 2 top of fault his vsi i 66 j 51 ;hist 3 right side of top of fault his vsi i 73 j 37 ;hist 4center of fault his vsi i 50 j 51 ;hist 5 center of imperm layer right side of fault his vsi i 72 j 34 ;hist 6 bottom right of fault his vsi i 73 j 18 ;hist 7 bottom of fault his vsi i 31 j 51 ;hist 8 bottom left of fault his vsi i 14 j 38 ;hist 9 bottom left hand corner his vsi i 4 j 10 ;hist 10 bottom right hand corner his vsi i 99 j 7 ;hist 11 under imperm layer rhside his vsi i 1 j 58 ;hist 12 under imperm layer lhside his vsi i 100 j 21 ;ex_2 hists, for flow through zones ;hist 2 top of fault his ex_2 i 66 j 51 ;hist 3 right side of top of fault his ex_2 i 73 j 37 ;hist 4center of fault his ex_2 i 50 j 51 ;hist 5 center of imperm layer right side of fault his ex_2 i 72 j 34 ;hist 6 bottom right of fault his ex_2 i 73 j 18 ;hist 7 bottom of fault his ex_2 i 31 j 51 ;hist 8 bottom left of fault his ex_2 i 14 j 38 ;hist 9 bottom left hand corner his ex_2 i 4 j 10 ;hist 10 bottom right hand corner his ex_2 i 99 j 7 ;hist 11 under imperm layer rhside his ex_2 i 1 j 58 ;hist 12 under imperm layer lhside his ex_2 i 100 j 21 ; set of fish functions which aim to simulate fault valve behaviour ; then when failure occurs permeability is increased def store_perm loop i (1,izones) loop j (1,jzones) ex_1(i,j) = k11(i,j) end_loop end_loop end store_perm * ;def vsi_perm ; change perm with respect to volume strain ;loop i (35,66) ; loop j (50,51) ; if vsi(i,j)>0.2 then ; k11(i,j)=k11(i,j)*10 ; k12(i,j)=0 ; k22(i,j)=k11(i,j) ; else ; if vsi(i,j)<-0.2 then ; k11(i,j)=k11(i,j)/10.0 ; k12(i,j)=0 ; k22(i,j)=k11(i,j) ; else ; if vsi(i,j)<=0.2 then ; if vsi(i,j)>=-0.2 then ; k11(i,j)=k11(i,j)*(exp(11.513*vsi(i,j))) ; k12(i,j)=0 ; k22(i,j)=k11(i,j) ; end_if ; end_if ; end_if ; end_if ; end_loop ;end_loop ;end * def perm_yield ;take this out for one test while_stepping incrval = 1.e6 incrval2 = 1.e6 incrval3 = 1.e-2 loop i (35,66) loop j (50,51) if state(i,j)=1 then ; at yield in shear or volume k11(i,j) = ex_1(i,j) * incrval2 k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval2 else if state(i,j)=2 then ; yielded in the past, not now k11(i,j) = ex_1(i,j) * incrval3 k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval3 else if state(i,j)=3 then ; yield in tension k11(i,j) = ex_1(i,j) * incrval k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval else if state(i,j)=6 then ; ubiquitous joints at yield in shear or volume k11(i,j) = ex_1(i,j) * incrval2 k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval2 else if state(i,j)=7 then ; ubiquitous joints yielded in the past, not now k11(i,j) = ex_1(i,j) * incrval3 k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval3 else if state(i,j)=8 then ; ubiquitous joints yield in tension k11(i,j) = ex_1(i,j) * incrval k12(i,j) = 0.0 k22(i,j) = ex_1(i,j) * incrval end_if end_if end_if end_if end_if end_if end_loop end_loop end * ; so can now use hist ex_3 to look at perm. def presentperm while_stepping loop i (1,izones) loop j (1,jzones) ex_3(i,j)=(k11(i,j)) end_loop end_loop end * ; FISH function to choose zones for plotting fluid flow def ff_plot2 k=0 kk=0 ll=0 loop i (1,izones) k=k+1 if k=2 then k=1 end_if loop j (1,jzones) ex_5(i,j)=0 ex_6(i,j)=0 kk=kk+1 if kk=2 then kk=1 end_if if k=1 then if kk=1 then ll=ll+1 ex_5(i,j)=xflow(i,j) ex_6(i,j)=yflow(i,j) end_if end_if end_loop end_loop end * ini ex_2=0.0 def flow_nick while_stepping loop i(1,igp) loop j(1,jgp) ex_2(i,j)=ex_2(i,j)+gflow(i,j) end_loop end_loop end flow_nick * save pump2a.sav * set nmech 1 ngw 10 flow on * initial xvel 0.0002 var -0.0004 0 i=1,101 j=1,101 ;this is half the veloc of Nerys pump2 apply dis= 1e-6 i= 1 101 j= 1 ; much faster discharge than Nerys' pump2 ;apply discharge -1e-9 j=101 * set nmech 1 ngw 10 flow on ; Nerys's pump2 had ngw=1 * def files if nn<10 then savefile = 'pump3_00' + string(nn) + '.sav' else if nn>9 then if nn<100 then savefile = 'pump3_0' + string(nn) + '.sav' else if nn>99 then savefile = 'pump3_' + string(nn) + '.sav' end_if end_if end_if end_if COMMAND save @savefile END_COMMAND end define step_loop loop nn (1,60) ;loop n (1,1000) command step 2000 end_command ;vsi_perm ;end_loop files end_loop end step_loop return