Setting up a steady-state pwp from rainfall

FLAC3D
1

Hey guys, I want to obtain a steady-state pwp contours from net discharge. I’ve never done this type of simulation before, so I’m totally lost. My code below is error free (syntax) but when I run the model, the flow ratio doesn’t evolve at all, it remains static and I don’t understand why.

I was looking at a similar example (FLAC though) and the discharge was positive, so I guess for infiltration the positive value is also for FLAC3D.

Another question, is it mandatory to have an initial head in the model? I was expecting to get one as a consequence of the infiltration. I had set up a head before, but again, the flow-ratio didn’t experiment any change. Please, any observation will be greatly appreciated.

Also, the ‘Top’ and ‘Bottom’ had been generated before through ‘Face skin’ command, and I used ‘Prop_py.f3sav’ for mechanical calculations before, so I don’t know.

; ------------------------------------------------------------------------------
; 1) Load/restore model,
; ------------------------------------------------------------------------------
model restore 'prop_py.f3sav'
model large-strain off
fish automatic-create off
geom import 'fault_hydro.dxf'
zone group 'fault_hydro' range geometry-distance 'fault_hydro' gap 1
; ------------------------------------------------------------------------------
; 2) Configure fluid model
; ------------------------------------------------------------------------------
model configure fluid

; ------------------------------------------------------------------------------
; 3) Assign fluid constitutive model and orientation
; ------------------------------------------------------------------------------
zone fluid cmodel assign anisotropic
zone fluid property dip 90 dip-direction 90 rotation 0  ; Vertical bedding plane (strike N-S)

; ------------------------------------------------------------------------------
; 4) Basic fluid properties: porosity, Biot (coefficient here), etc.
; ------------------------------------------------------------------------------
; Porosity
zone fluid property porosity 0.05                              ; Default for all
zone fluid property porosity 0.1 range group 'weathered'       ; Weathered zone
zone fluid property porosity 0.15 range group 'fault_hydro'       ; Weathered zone


; CORRECTED Biot modulus 
zone fluid property biot 0.8
zone fluid property biot 0.9 range group 'weathered'
zone fluid property biot 1   range group 'fault_hydro'

; ------------------------------------------------------------------------------
; 5) Global variables for water properties & gravity
; ------------------------------------------------------------------------------
[global mu_water = 1e-3
[global rho_f    = 1000
[global g        = 9.81

; ------------------------------------------------------------------------------
; 6) Base rock permeability (anisotropic)
;    permeability-1 = along-strike (N-S)
;    permeability-2 = vertical (bedding dip=90)
;    permeability-3 = cross-strike (E-W)
; ------------------------------------------------------------------------------
[global K_along_strike  = 1e-8]     ; m/s
[global K_vertical      = 1e-7]     ; m/s (10x higher vertical)
[global K_cross_strike  = 1e-8]     ; m/s

[global kappa_along = (K_along_strike * mu_water) / (rho_f * g)
[global kappa_vert  = (K_vertical     * mu_water) / (rho_f * g)
[global kappa_cross = (K_cross_strike * mu_water) / (rho_f * g)

zone fluid property permeability-1 [kappa_along] permeability-2 [kappa_vert] permeability-3 [kappa_cross]

; ------------------------------------------------------------------------------
; 7) Weathered zones (higher permeability)
; ------------------------------------------------------------------------------
[global K_vert_weathered = 1e-6]  ; m/s
[global kappa_vert_weathered = (K_vert_weathered * mu_water) / (rho_f * g)
zone fluid property permeability-2 [kappa_vert_weathered] range group 'weathered'

; ------------------------------------------------------------------------------


fish define set_fault_zone_permeability
    global K_h_fault = 1e-6    ; Horizontal permeability (m/s)
    global K_v_fault = 1e-5    ; Vertical permeability (m/s)
    local mu_water = 1e-3
    local rho_f    = 1000
    local g        = 9.81

    ; Convert to FLAC3D's intrinsic permeability (m²)
    local kappa_h_fault = (K_h_fault * mu_water) / (rho_f * g)
    local kappa_v_fault = (K_v_fault * mu_water) / (rho_f * g)

    loop foreach local zOne zone.list
        if zone.isgroup(zOne, "fault_hydro") = true then
            ; permeability-1 = along-strike (N-S)
            ; permeability-2 = vertical
            ; permeability-3 = cross-strike (E-W)
            zone.fluid.prop(zOne, "permeability-1") = kappa_h_fault
            zone.fluid.prop(zOne, "permeability-2") = kappa_v_fault
            zone.fluid.prop(zOne, "permeability-3") = kappa_h_fault
        endif
    end_loop
end
@set_fault_zone_permeability

; ------------------------------------------------------------------------------
; 9) Initialize fluid properties in the model
; ------------------------------------------------------------------------------
zone initialize fluid-density 1000
zone gridpoint initialize fluid-tension 0.0

; ------------------------------------------------------------------------------
; 10)Initialize fluid head/pressure 
; ------------------------------------------------------------------------------
; zone face apply head 95 range group 'Bottom' position-z 44.87 95

; ------------------------------------------------------------------------------
; 11) Rain boundary condition: discharge top surface
; ------------------------------------------------------------------------------
fish define dis_vert(zp, side)
    dis_vert = math.abs(zone.face.normal(zp, side)->z)
end

[global q_rain = 3.17e-9]  ; m/s (rain flux)
zone face apply discharge [q_rain] range group 'Top'

; ------------------------------------------------------------------------------
; 12) Draining base boundary
; ------------------------------------------------------------------------------
zone face apply pore-pressure 0 range group 'Bottom'

; ------------------------------------------------------------------------------
; 13) Solve steady-state fluid flow
; ------------------------------------------------------------------------------
model mechanical active off
model fluid active on
model fluid timestep automatic
model solve fluid ratio-flow 1e-3
model save 'insitu0_wet.f3sav'

image_2025-03-24_114156037
image_2025-03-24_1141560371154×1009 54.7 KB

1 Like
2

I am also struggling with a similar problem.

Especially with flac3d ver 9.1, the problem can be easily solved with a very simple command line “zone fluid steady-state”.

With flac3d ver 9.0 and below, it gets harder with multiple steps to troubleshoot and, if you’re lucky, a solution or results that you can’t understand like yours.

Can’t you update itasca to use the steady-state command with ver 9.0?

Or I believe someone can show me how to solve this problem.

Thank you.

1 Like
3

Hey, thanks for the feedback.

I’m indeed using the version 9.0 of the software, I’ll try with the command ‘steady-state’ but it seems it’s not available in this version since I haven’t seen it in the contextual help.

Hopefully Itascans will share their point of view on this. And if it’s a bug, it should get fixed.

Best,

4

@Felipe,
It’s not a bug. To take advantage to the new developments, such as those in the fluid logic (including the zone fluid steady-state command) you need to purchase (or switch) to version 9.3 (our subscription product). FLAC3D 9.0 is our perpetual product, which does not have these new developments.

As for your model, something isn’t right as I see the timestep being 8e19 or something like that, it’s blocked out in the image. Looking at your data file I am not sure what is wrong at this time.

5

Yes, it is not a bug, but the part Felipe mentioned should have worked properly, but it did not converge, so it was called a bug.

I know that the “steady state” command is not available in the current permanent version 9.0.

Also, looking at the data file above, the rainfall problem was also considered, but the unsaturated behavior is incomplete in 9.0. This makes it difficult to solve the problem in a general way.
Nevertheless, I wish the “steady state” command was included for free in version 9.0 to make steady flow analysis easier without the unsaturated command, but itasca’s policy says that it will be difficult.

In addition to what dblanksma said, there are many new commands introduced in version 9.1 that make the problem easier. There is a command for steady flow and unsaturated behavior (the “unsaturated cutoff” command) and an automatic scaling command for the fluid bulk modulus, which seems to be almost complete in version 9.3.

In this case, if you use fish and show me an example, it will help me to solve the problem in flac2d 9.0, flac3d 9.0. But I don’t know how to do it. So it’s hard to use 9.3 version, and if users who use 9.0 version provide fish or python code and example, it will be solved. What do you think? Is updating to 9.3 version really the only solution?

1 Like
6

Well indeed. From my original post the missing line was the one about the fluid modulus ‘zone gridpoint initialize fluid-modulus 2e9’ and then the timestep got reduced in a significant manner.

But well, the simulation will take AGES and I just want a general way to estimate steady-state pore water pressures, to use it for downstream mechanical simulations. This is a middle-sized model with ~1.2M zones and in mechanical it’s really fast.

Is there any other way? let’s discuss about options.

I don’t want to use the archaic method of Ru coefficients… and sadly I can’t import any pwp contours since I don’t have them available.

7

If you are only trying to get to steady state, then you don’t need to use the true fluid modulus. See equation 12 on this page: Solving Flow-Only and Coupled-Flow Problems — Itasca Software 9.0 documentation

You can set the fluid modulus orders of magnitude smaller, and solve for fluid only to get to the steady state relatively quickly.

8

In addition to what @jhazzard said, since this is only a flow only calculation and as far as I can tell it’s fully saturated (correct me if I am wrong) you can use an implicit solver to get very quick calculation speeds - see the zone fluid implicit command in addition to the documentation @jhazzard suggested.

9
  • I noticed that you multiplied mu_water when calculating permeability, which seems incorrect and unnecessary.
  • Just as we can use a lower bulk modulus to speed up calculations, we can also assign a higher ratio-flow. Depending on the model’s zone size, values like 0.01, 0.02, or even 0.1 may be sufficient. This can be verified by plotting a gridpoint pore-pressure history—if the curve flattens out, the value can be considered acceptable.
  • Increasing the smallest zone size can significantly improve computation speed, as the allowable time step is proportional to the square of the zone size.
10

Thanks a lot for the great suggestions, however the problem persists. I set up a very low fluid modulus (1e3) which seems to be the lowest possible value, because lower than that the pore-pressure contours sky-rocket to implausible values (~1e18).

@cheng the calculation was correct, but I do agree that it wasn’t necessary since you can work directly with the conductivities. I realized that later.

Thing is, despite of all your observations, the flow-ratio remains static, with no signs whatsoever of evolution or convergence. I also set up some pore-pressures based on zones, because I’d read that maybe my model wasn’t giving any results because there were no gradients at all. But it didn’t work.

Here in this picture you can see that the model updates pwp, however, there’s no sign of convergence:

image
image1199×948 124 KB

I don’t know what to fix. This model as is, should be able to converge.

In this example: Rainfall on a Slope (FLAC2D) — Itasca Software 9.3 documentation

A fish function is implemented to correct the discharge to account for the inclination of the slope, however when I do that in FLAC3D, I get this error and I don’t really understand how to fix it either:

flac3d>[global q_rain = 3.17e-8] ; m/s (rain flux) 3.17e-9
flac3d>zone face apply discharge [q_rain] fish-local dis_vert range group ‘Top’
*** Integer value cannot be less than 1
While executing line 109 columns 25 to 41 of source C:/

My updated, new condensed code is here:

; ------------------------------------------------------------------------------
; 1) Load/restore model,
; ------------------------------------------------------------------------------
model restore 'prop_py.f3sav'
model large-strain off
fish automatic-create off
geom import 'fault_hydro.dxf'
zone group 'fault_hydro' range geometry-distance 'fault_hydro' gap 1
zone fluid zone-based-pp true ;For setting up pore pressure based on zones
; ------------------------------------------------------------------------------
; 2) Configure fluid model
; ------------------------------------------------------------------------------
model configure fluid
zone fluid biot off ; This is the standard anyways
; ------------------------------------------------------------------------------
; 3) Assign fluid constitutive model and orientation
; ------------------------------------------------------------------------------
zone fluid cmodel assign anisotropic
zone fluid property dip 90 dip-direction 90 rotation 0  ; Vertical bedding plane (strike N-S)
; ------------------------------------------------------------------------------
; 4) Basic fluid properties: porosity, Biot (coefficient here), etc.
; ------------------------------------------------------------------------------
; Porosity
zone fluid property porosity 0.05 range group 'slates'                     
zone fluid property porosity 0.1 range group 'weathered'       ; Weathered zone
zone fluid property porosity 0.15 range group 'fault_hydro'       ; Weathered zone
; ------------------------------------------------------------------------------
; 5) Global variables for water properties & gravity
; ------------------------------------------------------------------------------
[global rho_f    = 1000
[global g        = 9.81
; ------------------------------------------------------------------------------
; 6) Base rock conductivity (anisotropic)
;    -1 = along-strike (N-S)
;    -2 = vertical (bedding dip=90)
;    -3 = cross-strike (E-W)
; ------------------------------------------------------------------------------
[global K_along_strike  = 1e-7]     ; m/s
[global K_vertical      = 1e-6]     ; m/s (10x higher vertical)
[global K_cross_strike  = 1e-7]     ; m/s

zone fluid property hydraulic-conductivity-1 [K_along_strike]...
hydraulic-conductivity-2 [K_vertical]...
hydraulic-conductivity-3 [K_cross_strike] range group 'slates'
; ------------------------------------------------------------------------------
; 7) Weathered zones (higher permeability)
; ------------------------------------------------------------------------------
[global K_vert_weathered = 1e-5]  ; m/s
[global K_1_weathered = 1e-6
[global K_3_weathered = 1e-6
zone fluid property hydraulic-conductivity-1 [K_along_strike]...
hydraulic-conductivity-2 [K_vert_weathered]...
hydraulic-conductivity-3 [K_3_weathered] range group 'weathered'
; ------------------------------------------------------------------------------
fish define set_fault_zone_conductivity
    global K_h_fault = 1e-5    ; m/s
    global K_v_fault = 1e-4    ; m/s
    
    loop foreach local zOne zone.list
        if zone.isgroup(zOne, "fault_hydro") = true then
            ; permeability-1 = along-strike (N-S)
            ; permeability-2 = vertical
            ; permeability-3 = cross-strike (E-W)
            zone.fluid.prop(zOne, "hydraulic-conductivity-1") = K_h_fault
            zone.fluid.prop(zOne, "hydraulic-conductivity-2") = K_v_fault
            zone.fluid.prop(zOne, "hydraulic-conductivity-3") = K_h_fault
        endif
    end_loop
end
@set_fault_zone_conductivity
; ------------------------------------------------------------------------------
; 9) Initialize fluid properties in the model
; ------------------------------------------------------------------------------
zone initialize fluid-density 1e3
zone gridpoint initialize fluid-tension 0.0
zone gridpoint initialize fluid-modulus 1e3 ;Low to speed up convergence.

; ------------------------------------------------------------------------------
; 10)Initialize fluid head/pressure 
; ------------------------------------------------------------------------------
; zone face apply head 95 range group 'Bottom' position-z 44.87 95
[global water_level = 95.0  ; Water table elevation in meters

; Loop over all zones and assign pore pressure where below water table
fish define set_water_pore_pressure
    ; Loop over all zones in the model
    loop foreach local zoneID zone.list
        ; retrieve the z-coordinate of the zone’s centroid
        local z_centre = zone.pos.z(zoneID)
        
        if (z_centre < water_level) then
            ; calculate pore water pressure: p = density * gravity * (water_level - z_centre)
            zone.pp(zoneID) = g * rho_f * (water_level - z_centre)
        else
            ; zones above the water table get zero pore pressure
            zone.pp(zoneID) = 0.0
        endif
    end_loop
end
@set_water_pore_pressure

; ------------------------------------------------------------------------------
; 11) Rain boundary condition: discharge top surface
; ------------------------------------------------------------------------------
fish define dis_vert(zp, side)
    dis_vert = math.abs(zone.face.normal(zp, side)->z)
end

[global q_rain = 3.17e-8]  ; m/s (rain flux) 3.17e-9
zone face apply discharge [q_rain] range group 'Top'
; ------------------------------------------------------------------------------
; 12) Draining base boundary
; ------------------------------------------------------------------------------
zone face apply pore-pressure 0 range group 'Bottom'
; ------------------------------------------------------------------------------
; 13) Solve steady-state fluid flow
; ------------------------------------------------------------------------------
model mechanical active off
model fluid active on
model fluid timestep fix 1e4 ;decreased on purpose to have meaningful pwp contours
model solve fluid ratio-flow 1e-2
model save 'insitu0_wet.f3sav'

Please any advice is highly appreciated.

11

When reviewing your initial data file, I realized I missed mentioning a few important points.

  1. Hydraulic Conductivity Over Permeability?
    I overlooked the manual’s recommendation to use hydraulic conductivity rather than permeability when approaching fluid flow. Fortunately, FLAC3D 9.0 supports the necessary commands to apply this correctly. You’ve implemented this well in your updated data file.

  2. Biot Theory
    Initially, your data file had Biot theory enabled, but it was later disabled. Based on my experience with basic examples, keeping Biot theory on tends to yield more intuitive results. This is because, in a fully coupled analysis, simply increasing permeability without considering other factors can lead to incorrect results or even divergence. Using Biot theory alongside the bulk modulus provides more consistent outcomes. However, your initial data file was missing a fluid bulk modulus equivalent to the Biot bulk modulus in the “zone gridpoint initialise biot” command.

  3. Enabling the Implicit Method
    Instead of enabling the implicit method, you only modified the time step. However, enabling the implicit method generally improves stability and convergence, making it a better approach.

  4. Utilizing FLAC3D Geometry for Initial Conditions
    FLAC3D’s powerful geometry features can be used to set the initial groundwater table as a boundary condition. Since FLAC3D 9.0 lacks a steady-state command, you may need to use something like “water set geometry~~” to define initial pore pressure conditions. Then, using “zone water clear” allows you to simulate pore pressure changes over time for mechanical analysis or groundwater level variations.

  5. Using the dis_vert FISH Function
    It looks like you intended to use the dis_vert FISH function but didn’t actually apply it. According to the manual, implementing it as
    zone face apply drain [q_rain] fish-local dis_vert range group ‘Top’
    would be a more effective approach.

finally, I recently tested FLAC2D 9.0 while working on a similar example to “Rainfall on a Slope” from FLAC2D, as documented in Itasca Software 9.3. However, I encountered a phenomenon where I couldn’t replicate the same results due to the differences in unsaturated soil characteristics and varying water heads on both sides. I don’t know

1 Like
12

I suggest that you send the complete data file (not from a save file because we do not know which you have done) to flac3dsupport@itascacg.com so we can run it, duplicate the issue, and determine the cause. Also, please specify the version and subversion you are using.

1 Like
13

Mail sent. Hopefully we can solve this problem. Thanks.

I’ll update the results

Best,