; This is the input file for model 10, but can be seen as a general input file. The elastic moduli, shear and bulk, are 
; varied, as well as the contact angle between the low and high permeability fields. Another important variable is
; the increase in permeability as a result of yielding.

; 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 comparison 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 roughly but not exactly to the same
; amount of deformation (here it is 430m*2) there it is ca. 300m*2


tit
FLAC Modelling Input File
;
config extra 4 gw
grid  100   100 
gen 0.00e0,0.00e0 0.00e0,1.00e4 3.0e4,1.00e4 3.0e4,0.0e0 i=1,101 j=1,101
mark i=1
mark j=1
mark i = 101 
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, 13500, 10000, 16500, 5000, 30000, 5000, 30000, 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=75, 75 ; bulk and shear are further increased compared to compl5 
pro coh=3.00e7 tens=1.5e7 reg=75, 75
pro fric=40 dil= 10  reg=75, 75
pro perm=1.00e-18 porosity=0.2 reg=75, 75

plo hold perm blo

; the following five lines impose a pressure at the top of the model, simulating burial of the model at 1 kbar in this case,
; and initialize the model with a lithostatic gradient from 1 kbar at the top to 3.668 at the base of the model, 
; which corresponds to a depth of 10 km for the rock density used in these models of 2720 kg/m3.

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

fix x i=1
fix x i= 101 
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
*

; the fish function below controls permeability changes as a function of yielding.

def perm_yield

while_stepping
    incrval = 50.0
    incrval2 = 20.0
    incrval3 = 10.0
    loop i (1,100)
        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;


;
set nmech 1 ngw 1 flow on
ini xvel  0.002 var=-0.004 ,0 i=1, 101  j=1, 101 
step 5000
save compl10_0001.sav
step 5000
save compl10_0002.sav



set nmech 1 ngw 10 flow on


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

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


def files
    if nn<10 then
     savefile = 'compl10_00' + string(nn) + '.sav'
    else
    if nn>9 then
     if nn<100 then
     savefile = 'compl10_0' + string(nn) + '.sav'
    else
    if nn>99 then
     savefile = 'compl10_' + 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 10000
	end_command
      ;vsi_perm
   ;end_loop
files
end_loop
end
step_loop




return