Potential error in PorousFlowUnsaturated Action: neglect residual_saturation for CapillaryPressureVG

[GiudGiud]

Description

See below

How to reproduce

Run input file attached below

Impact

Some simulations cannot be performed with expected accuracy

Discussed in https://github.com/idaholab/moose/discussions/27693

^{Originally posted by amoure11 May 23, 2024} Hi,

I might have found an error in the PorousFlowUnsaturated Action: I ran a simple example with constant fixed pressure P=-1E6, and with the following PorousFlowUnsaturated parameters:

  relative_permeability_exponent = 3
  relative_permeability_type = Corey
  residual_saturation = 0.1
  van_genuchten_alpha = 1E-6
  van_genuchten_m = 0.6

According to the van Genuchten expression S_eff = (1+(-alpha*P)^(1/(1-m))^(-m), I would get S_eff = 0.6597 According to the effective saturation expression: S_eff = (S-S_res)/(1-S_res), I would get S = 0.6937

I open the output (exodus) file, and saturation is ~0.66, which corresponds to S_eff, and not S. Thus, I am assuming that the residual saturation (S_res) is not taken into account when computing the 'CapillaryPressureVG' Could you clarify this?

FYI, I ran the same example without using the Action 'PorousFlowUnsaturated', and the saturation value in the output file is ~0.69.

---------------------------- Example ----------------
[Mesh]
  [generate_mesh]
    type = GeneratedMeshGenerator
    dim = 3
    nx = 20
    xmin = -800
    xmax = 800
    ny = 20
    ymin = -800
    ymax = 800
    nz = 15
    zmin = -1200
    zmax = 0
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '-800 -800 -800'
    top_right = '800 800 -400'
    input = generate_mesh
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'y<-799.9 & x>-200 & x<200'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
    initial_condition = -1E6
  []
[]

[PorousFlowUnsaturated]
  porepressure = porepressure
  coupling_type = Hydro
  gravity = '0 0 0'
  fp = the_simple_fluid
  relative_permeability_exponent = 3
  relative_permeability_type = Corey
  residual_saturation = 0.1
  van_genuchten_alpha = 1E-6
  van_genuchten_m = 0.6
  use_displaced_mesh = false
[]

[AuxVariables]
  [saturation_aux]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [saturation_aux]
    type = PorousFlowPropertyAux
    variable = saturation_aux
    property = saturation
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = SimpleFluidProperties
    bulk_modulus = 2E9
    viscosity = 1.0E-3
    density0 = 1000
    thermal_expansion = 0.0002
    cp = 4194
    cv = 4186
    porepressure_coefficient = 0
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []
[]

[Preconditioning]
  active = preferred_but_might_not_be_installed
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  automatic_scaling = true
  end_time = 1E6
  dt = 1E6
  nl_abs_tol = 1E-10
  nl_rel_tol = 1E-5
[]

[Outputs]
  exodus = true
[]