Miscible CO₂ Injection Convergence Problem and Capillary Pressure Issues

[QUIUQ]

Problem Description:
I encountered convergence issues while simulating CO₂ injection using the MOOSE framework. The problem setup is as follows:
Simulation Conditions:
Axisymmetric cylindrical space injection with isothermal conditions and two-phase miscible gas (CO₂) injection.

When using a constant capillary pressure (pc = 0), the injection simulation runs successfully for 1 year.
However, when I use the capillary pressure model [pc], with the PorousFlowCapillaryPressureVG type, and the following parameters:

[pc]  
    type = PorousFlowCapillaryPressureVG  
    alpha = 0.5  
    m = 0.45  
    sat_lr = 0.2  
[]  

I found that during the simulation, the pressure falls below the range of the tabulated data. The formation pressure is around 3e7, but the simulated pressure goes out of range.

Additionally, when using the following parameters:

[pc]  
    type = PorousFlowCapillaryPressureVG  
    alpha = 1e-5  
    m = 0.5  
    sat_lr = 0.2  
[]  

The program quickly fails during the advection process from the first row of the grid on the left to the right side. I'm not sure why this happens. Even with alpha = 1e-5, which should be a reasonable value in practice, the simulation still fails.

I suspect the issue is with the capillary pressure settings or that the settings for the miscible CO₂ model are incorrect.
Moreover, I noticed that during injection, the overpressure of water does not increase significantly, while in TOUGH, the overpressure increases much more compared to the MOOSE simulation.
And here is the link to my code. https://github.com/QUIUQ/CO2-injection-MOOSE

[Mesh]
    type = GeneratedMesh
    dim = 2
    ny = 25
    nx = 50
    ymax = 100
    xmin = 0.1
    xmax = 1200
    bias_x = 1.05
    coord_type = RZ
    rz_coord_axis = Y
[]

[GlobalParams]
    PorousFlowDictator = dictator
    gravity = '0 -9.81 0'
    temperature_unit = Celsius
[]

[Variables]
    [pgas]
        #scaling = 1e-7
    []
    [zi]
        initial_condition = 0
    []
[]

[ICs]
    [pp_initial]
        type = FunctionIC
        variable = pgas
        function = pp_ic
    []
[]

[Functions]
    [pp_ic]
        type = SolutionFunction
        solution = soln
    []
    [injection_rate]
        type = ParsedFunction
        symbol_values = injection_area
        symbol_names = area
        expression = '-5/area'
    []
[]


[UserObjects]
    [soln]
        type = SolutionUserObject
        mesh = initial_condition.e
        system_variables = porepressure
    []
    [dictator]
        type = PorousFlowDictator
        porous_flow_vars = 'pgas zi'
        number_fluid_phases = 2
        number_fluid_components = 2
    []
    [pc]
        type = PorousFlowCapillaryPressureVG
        alpha = 1e-5
        m = 0.5
        sat_lr = 0.3
    []
    [pc0]
        type = PorousFlowCapillaryPressureConst
        pc = 0
    []
    [fs]
        type = PorousFlowBrineCO2
        brine_fp = brine
        co2_fp = co2
        capillary_pressure = pc
    []
[]

[AuxVariables]
    [temperature]
        initial_condition = 50
    []
    [xnacl]
        initial_condition = 0.1
    []
    [pressure_gas]
        order = CONSTANT
        family = MONOMIAL
    []
    [pressure_water]
        order = CONSTANT
        family = MONOMIAL
    []
    [saturation_gas]
        order = CONSTANT
        family = MONOMIAL
    []
    [saturation_water]
        order = CONSTANT
        family = MONOMIAL
    []
    [density_water]
        order = CONSTANT
        family = MONOMIAL
    []
    [density_gas]
        order = CONSTANT
        family = MONOMIAL
    []
    [x0_water]
        order = CONSTANT
        family = MONOMIAL
    []
    [x0_gas]
        order = CONSTANT
        family = MONOMIAL
    []
    [x1_water]
        order = CONSTANT
        family = MONOMIAL
    []
    [x1_gas]
        order = CONSTANT
        family = MONOMIAL
    []
    [cap_pre]
        order = CONSTANT
        family = MONOMIAL
    []
[]

[Kernels]
    [mass0]
        type = PorousFlowMassTimeDerivative
        variable = pgas
        fluid_component = 0
    []
    [flux0]
        type = PorousFlowAdvectiveFlux
        variable = pgas
        fluid_component = 0
    []
    [mass1]
        type = PorousFlowMassTimeDerivative
        variable = zi
        fluid_component = 1
    []
    [flux1]
        type = PorousFlowAdvectiveFlux
        variable = zi
        fluid_component = 1
    []
[]

[AuxKernels]
    [pressure_water]
        type = PorousFlowPropertyAux
        variable = pressure_water
        property = pressure
        phase = 0
        execute_on = timestep_end
    []
    [pressure_gas]
        type = PorousFlowPropertyAux
        variable = pressure_gas
        property = pressure
        phase = 1
        execute_on = timestep_end
    []
    [saturation_water]
        type = PorousFlowPropertyAux
        variable = saturation_water
        property = saturation
        phase = 0
        execute_on = timestep_end
    []
    [saturation_gas]
        type = PorousFlowPropertyAux
        variable = saturation_gas
        property = saturation
        phase = 1
        execute_on = timestep_end
    []
    [density_water]
        type = PorousFlowPropertyAux
        variable = density_water
        property = density
        phase = 0
        execute_on = timestep_end
    []
    [density_gas]
        type = PorousFlowPropertyAux
        variable = density_gas
        property = density
        phase = 1
        execute_on = timestep_end
    []
    [x1_water]
        type = PorousFlowPropertyAux
        variable = x1_water
        property = mass_fraction
        phase = 0
        fluid_component = 1
        execute_on = timestep_end
    []
    [x1_gas]
        type = PorousFlowPropertyAux
        variable = x1_gas
        property = mass_fraction
        phase = 1
        fluid_component = 1
        execute_on = timestep_end
    []
    [x0_water]
        type = PorousFlowPropertyAux
        variable = x0_water
        property = mass_fraction
        phase = 0
        fluid_component = 0
        execute_on = timestep_end
    []
    [x0_gas]
        type = PorousFlowPropertyAux
        variable = x0_gas
        property = mass_fraction
        phase = 1
        fluid_component = 0
        execute_on = timestep_end
    []
    [capillary_pressure]
        type = PorousFlowPropertyAux
        variable = cap_pre
        property = capillary_pressure
        execute_on = timestep_end
    []
[]

[BCs]
    [methane_injection]
        type = PorousFlowSink
        boundary = left
        variable = zi
        flux_function = injection_rate
        fluid_phase = 1
    []
    [outer_pressure_fixed]
        type = FunctionDirichletBC
        variable = pgas
        boundary = right
        function = pp_ic
    []
[]

[FluidProperties]
    [water]
        type = Water97FluidProperties
    []
    [watertab]
        type = TabulatedBicubicFluidProperties
        fp = water
    []
    [brine]
        type = BrineFluidProperties
        water_fp = watertab
    []
    [co2sw]
        type = CO2FluidProperties
    []
    [co2]
        type = TabulatedBicubicFluidProperties
        fp = co2sw
    []

[]

[Materials]
    [temperature]
        type = PorousFlowTemperature
        temperature = temperature
    []
    [brineco2]
        type = PorousFlowFluidState
        gas_porepressure = pgas
        z = zi
        temperature = temperature
        xnacl = xnacl
        capillary_pressure = pc
        fluid_state = fs
    []
    [porosity]
        type = PorousFlowPorosityConst
        porosity = 0.1
    []
    [permeability]
        type = PorousFlowPermeabilityConst
        permeability = '1.5833e-13 0 0  0 1.5833e-13 0  0 0 1.5833e-13' #160.43 mD
    []
    [relperm_water]
        type = PorousFlowRelativePermeabilityVG
        m = 0.45946
        phase = 0
        s_res = 0.300
        sum_s_res = 0.3
      []
      [relperm_gas]
        type = PorousFlowRelativePermeabilityBC
        phase = 1
        s_res = 0.0
        sum_s_res = 0.3
        lambda = 2
        nw_phase = true
      []
[]

[Preconditioning]
    [andy]
      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
    end_time = 31536000
    nl_abs_tol = 1e-7
    nl_rel_tol = 1e-5
    dtmax = 1e5
    nl_max_its = 15
    automatic_scaling = true
    #compute_scaling_once = true
    [TimeStepper]
      type = IterationAdaptiveDT
      dt = 1e1
      growth_factor = 1.25
    []
  []

  [Postprocessors]
    [mass_ph0]
      type = PorousFlowFluidMass
      fluid_component = 0
      execute_on = 'initial timestep_end'
    []
    [mass_ph1]
      type = PorousFlowFluidMass
      fluid_component = 1
      execute_on = 'initial timestep_end'
    []
    [injection_area]
      type = AreaPostprocessor
      boundary = left
      execute_on = initial
    []
  []

  [Outputs]
    execute_on = 'initial timestep_end'
    exodus = true
    perf_graph = true
    [csv]
        type = CSV
        sync_times = '2.592e6 8.64e6 31536000'
        sync_only = true
    []
  []

[GiudGiud]

hello

the tabulatedFluidProperties have other options than throwing for the out of bounds cases
https://mooseframework.inl.gov/source/fluidproperties/TabulatedBicubicFluidProperties.html

can you try those for your case?

[QUIUQ]

Thank you, @GiudGiud

However, I believe the main issue lies elsewhere. I have used the capillary pressure parameters from the restart examples, so I believe the parameters are correct. https://mooseframework.inl.gov/modules/porous_flow/restart.html

The change I made was to the fluid properties, replacing immiscible methane with miscible CO₂ gas. However, the simulation fails quickly （the dt become extremely small）. Upon examining the results, I found that the simulation stalls when the CO₂ gas fills the left mesh (the injection zone) and stops when attempting to advect to the right mesh.

Is there something wrong with my input file? This issue is really frustrating, and I hope I can get some help. and my input file https://github.com/QUIUQ/CO2-injection-MOOSE/blob/main/co2_injection_vg.i

[GiudGiud]

Have you tried some of the solutions from this page? what were the results?
https://mooseframework.inl.gov/moose/application_usage/failed_solves.html

[GiudGiud]

I dont see anything obviously wrong in the input.

does this simulation solve with other fluid properties? Maybe more simple ones? ( replacing immiscible methane/ did those work?)

[QUIUQ]

Thank you, @GiudGiud
I have successfully run simulations using immiscible methane with the capillarypressureVG model and an alpha value of 1e-5.

Additionally, I tried using an alpha value of 1 for miscible CO₂, which also worked.

However, when I set alpha to 1e-5 for miscible CO₂, the simulation fails to advect to the next mesh.

I suspect that the alpha value might be the issue.

An alpha of 1e-5 is considered a reasonable engineering value, isn't it?

Could it be that a smaller alpha increases the capillary pressure, preventing the gas from entering adjacent meshes?

[GiudGiud]

An alpha of 1e-5 is considered a reasonable engineering value, isn't it?

I would not know

from alpha = 1 to alpha = 1e-5, what changes in the numerical solve?
If we can see an impact on the scaling of the equations, or the conditioning of the matrix, or the distance between the solution and the initial condition, then this could be a clue to how to fix this

[cpgr]

I think it is still the issue that is noted in #23197 (reply in thread)

The workaround in that case was to make a second (constant) capillary pressure user object and use that in the material. It doesn't affect the results, only the derivatives and hence has helped convergence.

[QUIUQ]

Thank you, @GiudGiud and @cpgr , for your responses and suggestions.
I appreciate your help! The workaround with the constant capillary pressure user object has helped with convergence. I also realized that I should set the pc_max, because 𝑝_water = 𝑝_gas − 𝑝𝑐, and if 𝑝𝑐 is too high, 𝑝_water might fall below the range of the tabulated flow. I'm not sure if my understanding is correct.
I'll continue learning the MOOSE Framework. Thank you again for your guidance and support!