Tensional rupture appears when upscaling one_way.dat

Hi everyone!

I am facing some troubles with a simple one-way coupling problem, where hydrostatic pressures systematically trigger tensional rupture in faults.

To understand the problem, I used the 3DEC one_way.dat example from the documentation.
Joint Fluid Flow — 3DEC 7.0 documentation.

The default script highlights shear rupture along the fault plane.

Now, if I scale up the model dimensions from metric to kilometric, tensional rupture occurs at the pressure front.

I cannot figure out the physical problem behind this. Gravity would change something but is not activated here.
Am I missing something in 3DEC parameters?

Thanks a lot for your help and suggestions!

Regards,

Théophile

PS : the one_way.dat I slightly modified to have comparable mesh cell dimensions

model new
; one-way coupled simulation
model config fluid
model large-strain off

def parameters
    ; Default
    ; local edge = 1
    ; Upscaled
    local edge = 1000

    global xmin = -edge
    global xmax = edge
    global ymin = -edge
    global ymax = edge
    global zmin = -edge
    global zmax = edge
    
    global mesh_edge = edge/10.
end
[parameters]


block create brick [xmin], [xmax], [ymin], [ymax], [zmin], [zmax]
block cut joint-set dip 0 dip-direction 90

block zone generate edgelength [mesh_edge]
block zone cmodel assign elastic
block zone property density 2500 shear 2e9 bulk 5e9 

block contact jmodel assign mohr
block contact property stiffness-normal 1e10 stiffness-shear 1e10 

block fluid property bulk 3.8e5
block fluid property density 1000
block fluid property viscosity 1e-3

flowplane vertex property aperture-initial 1e-4 ...
                          aperture-minimum 1e-4 aperture-maximum 1e-4 

;
block gridpoint apply velocity-x 0 range position-x [xmin]
block gridpoint apply velocity-x 0 range position-x [xmax]
block gridpoint apply velocity-y 0 range position-y [ymin]
block gridpoint apply velocity-y 0 range position-y [ymax]
block gridpoint apply velocity-z 0 range position-z [zmin]
block gridpoint apply velocity-z 0 range position-z [zmax]
;
block insitu stress -1e7 -1e7 -1e7 0 0 0

model fluid active off
model mechanical active on
block fluid property bulk 0.0

model history mechanical unbalanced-maximum
model cycle 100

model fluid active on
model mechanical active off
;
flowknot apply pore-pressure 1.1e7 range position-x [xmin]
flowknot apply pore-pressure 0.0 range position-x [xmax]
;
block fluid property bulk 3.8e5

history delete
model his fluid time-total
flowknot history pore-pressure position 0 0 0
model cycle 1000
model fluid active off
model mechanical active on
block fluid property bulk 0.0
model cycle 1000
program return

Hi,
It’s because the pore pressure gradient along the x-direction is not the same.

Hi, thanks a lot for your answer!

The gradient was indeed different, but rectification does not change the issue:

Tensional rupture remains when testing various upscaling factors (e.g., x10 rather than x1000).

Refining the mesh does not solve the issue.

The problem might come from localized effects of pore pressure close to the joint where the closest cells might concentrate most of the deformation, hence leading to excessive normal displacements and eventually to tensional rupture. Increasing the block moduli wipes out tensional rupture, which tends to confirm this assumption but I cannot be 100% sure at this moment.

Regards,

Théophile

PS: a revised version of the script. Use use_default = true for the default one_way.dat model, use_default = false for the upscaled one, and use_default = false + generate new options for the upscaled and refined one.

model new
; one-way coupled simulation
model config fluid
model large-strain off

def parameters
    local use_default = false
    
    ; Default
    local edge_default = 1
    ; Upscaled
    local edge_upscaled = 1000
    
    global pressure_xmin = 1.1e7
    global pressure_xmax = 0
    edge = edge_upscaled
    
    if use_default
        edge = edge_default
        pressure_xmax = pressure_xmin + edge_default/edge_upscaled*(pressure_xmax-pressure_xmin)
    end_if

    global xmin = 0.
    global xmax = 2*edge
    global ymin = -edge
    global ymax = edge
    global zmin = -edge
    global zmax = edge
    
    global mesh_edge = edge/10.
    global mesh_edge_min = edge/20.
    global mesh_edge_max = edge/8.
end
[parameters]


block create brick [xmin], [xmax], [ymin], [ymax], [zmin], [zmax]
block cut joint-set dip 0 dip-direction 90

;block zone generate edgelength [mesh_edge]
block zone size edge [mesh_edge_max] range position-z [zmin] by block-bface
block zone size edge [mesh_edge_max] range position-z [zmax] by block-bface
block zone size edge [mesh_edge_min] range position-z -0. by block-bface
block zone generate-new minimum-edge [mesh_edge_min] maximum-edge [mesh_edge_max] gradation-volume 0.2

block zone cmodel assign elastic
block zone property density 2500 shear 2e9 bulk 5e9 

block contact jmodel assign mohr
block contact property stiffness-normal 1e10 stiffness-shear 1e10 

block fluid property bulk 3.8e5
block fluid property density 1000
block fluid property viscosity 1e-3

flowplane vertex property aperture-initial 1e-4 ...
                          aperture-minimum 1e-4 aperture-maximum 1e-4 

;
block gridpoint apply velocity-x 0 range position-x [xmin]
block gridpoint apply velocity-x 0 range position-x [xmax]
block gridpoint apply velocity-y 0 range position-y [ymin]
block gridpoint apply velocity-y 0 range position-y [ymax]
block gridpoint apply velocity-z 0 range position-z [zmin]
block gridpoint apply velocity-z 0 range position-z [zmax]
;
block insitu stress -1e7 -1e7 -1e7 0 0 0

model fluid active off
model mechanical active on
block fluid property bulk 0.0

model history mechanical unbalanced-maximum
model cycle 100

model fluid active on
model mechanical active off
;
flowknot apply pore-pressure [pressure_xmin] range position-x [xmin]
flowknot apply pore-pressure [pressure_xmax] range position-x [xmax]
;
block fluid property bulk 3.8e5

history delete
model his fluid time-total
flowknot history pore-pressure position 0 0 0
model cycle 5000
model fluid active off
model mechanical active on
block fluid property bulk 0.0
model cycle 5000
program return

Have a look at the displacements and stresses in the two models. And also the normal stresses on the joints. They are completely different in the two models. I’m not exactly sure why - maybe because you have a much larger volume of rock able to accommodate the stress change caused by the pore pressure in the joint. Because the top and bottom boundaries are fixed, it’s harder for the joint to open in the small model. I guess you would expect the same response if the boundaries were really far away from the fault.

Thanks for your advice!
I shall have a look soon.

Regards,

Théophile

I could be wrong, but i would concur with @jhazzard, since the larger model should also have a larger crack aperture as a result of the much longer crack itself (think of it as a large lever vs. a small lever). There should even be analytical solutions for similar situations (I think there was one in a UDEC6 example).

Anyways with larger aperture (normal displacement), but the same normal stiffnesses you will possibly exceed the normal displacement for tensile failure in your large model, but not in the small one.

Hi!

Sorry about my (very) late answer.

Thanks a lot for your suggestions.
Interestingly, I haven’t found an analytical solution so far, it doesn’t seem to be that much of an easy problem.

Anyway, and since I found a workaround for my computations, I did not dig further into this case.
I shall eventually have a deeper look, and your advice would really help.

So, thanks again!

Cheers

Theophile