;Model to test the effects of the contact between an impermeable and a permeable layers. This is to test the ideas produced by the fractal study
; compl3 has an increased bulk modulus of the perm. layer compared to compl2 which has all properties of the two layers equal except for permeab.
;in comparisson with compl4 this one has a sloping boundary between the two layers, 
;plus I decreased the cohesion and the tensile strength of the imperm. layer.

; this is an intermediate velocity and the results here are very similar to those at the 
;slowest velocity in Compl5 plot 050 which corresponds roughlyh but not exactly to the same
;amount of deformation (here it is 430m*2) there it is ca. 300m*2

; compl8 equal compl5 intermediate velocity but shear and bulk have been further increased in the imperm. layer. and
;decreased p at top to 1 kbar
; compl9  kept the same moduli and Mohr-Coulomb parameters, reduced the veloc by half and added a permability function from
; Nerys's pump, 
; compl15b is similar to the very successful compl10 but improved resolution to 300*150.
; compl17 is like 15b but de3creased perm to 1e-18 since in 15b with 1e-16 the results were
;significantly different than compl10 which had 1e-18 


tit
FLAC Modelling Input File
;
config extra 4 gw
grid  200   100 
gen 0.00e0,0.00e0 0.00e0,1.00e4 2.8e4,1.00e4 2.8e4,0.0e0 i=1,201 j=1,101
mark i=1
mark j=1
mark i = 201 
mark j = 101


***************************************
model mohr
water bulk 2e9 dens 1000
prop poro=0.3 ; to be elaborated below
;; 1 mech step = 1 fluid step; modify if necessary!!!
set nmech 1 ngw 10 flow on
***************************************

tab 1 0, 0, 0, 10000, 12500, 10000, 15500, 5000, 28000, 5000, 28000, 0

gen tab 1


;
; Properties for each sub-regions: 
pro dens= 2720  bulk=4.0e9 shear=2.45e9 reg=20, 20
pro coh=3.00e7 tens=1.5e7 reg=20, 20
pro fric=40 dil= 10  reg=20, 20                     ; changed fric and dil compared to all previour "compl" runs
pro perm=1.00e-14 porosity=0.2 reg=20, 20

; Properties for each sub-regions: 
pro dens= 2720  bulk=8.e10 shear=2.45e10  reg=180, 80 ; bulk and shear are further increased compared to compl5 
pro coh=3.00e7 tens=1.5e7  reg=180, 80
pro fric=40 dil= 10   reg=180, 80
pro perm=1.00e-18 porosity=0.2 reg=180, 80


;
apply syy = -1.e8 j= 101 
*****
ini sxx=-3.66800e8 var=0, 2.668e8 i=1, 201 j=1,101
ini syy=-3.66800e8 var=0, 2.668e8 i=1, 201 j=1,101
ini szz=-3.66800e8 var=0, 2.668e8 i=1, 201 j=1,101
ini sxy=0e6


fix x i=1
fix x i= 201 
fix y j=1

;
def hydro_pp
loop i (1,igp)
loop j (1,jgp)
if j = jgp Then
sss2 = 3.67e7
;sss2 = 7.35e7
else
aaa = abs(y(i, j) - y(i, jgp))
sss2 = 1000 * 9.8 * aaa+3.67e7
;sss2 = 1000 * 9.8 * aaa+7.35 e7 ;for 2kb
end_if
gpp(i, j) = sss2
end_loop
end_loop
end
hydro_pp
*
;from Nerys's pump
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 perm_yield
while_stepping
    incrval = 50.0
    incrval2 = 20.0
    incrval3 = 10.0
    loop i (1,200)
        loop j (1,100)
               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
              end_if
              end_if
              end_if
              end_loop
    end_loop
end
perm_yield
;
ini pp=3.67e7 j=101
fix pp j=101

ini sat 1
fix sat
;
set gravity=9.8
;
hist 1 unbal nstep 1

set small
set mech on fl off
step 200
ini   xd 0    yd 0
ini xvel 0  yvel 0
;plot hold hist 1 
set mech on fl on
step 400  ; this can be increased if you think the unbalanced forces haven't bottomed out
ini   xd 0    yd 0 
ini xvel 0  yvel 0

set large
;plot hold hist 1


;initial velocity command should be here ie ini xvel 0.004 etc;


;stress loading time
set nmech 1 ngw 1 flow on
ini xvel  0.001 var=-0.002 ,0 i=1, 201  j=1, 101 
step 5000
save compl17_0001.sav
step 5000
save compl17_0002.sav



set nmech 1 ngw 10 flow on


;ini xvel  0.001 var=-0.002 ,0 i=1, 301  j=1, 101 ; halved the velocity compared to compl8

;set mech off
apply dis= 1e-10 i= 1 201 j= 1


def files
    if nn<10 then
     savefile = 'compl17_00' + string(nn) + '.sav'
    else
    if nn>9 then
     if nn<17 then
     savefile = 'compl17_0' + string(nn) + '.sav'
    else
    if nn>17 then
     return
    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 10000
	end_command
      ;vsi_perm
   ;end_loop
files
end_loop
end
step_loop




return