CO2 multiphase flow coupled with Geochemistry

[Femi-Oyenowo]

Hi,

Please I am trying to simulate sCO2 storage in an aquifer with a complex brine composition (i.e. has NaCl, HCO3, K, Ca etc). I want to simulate the CO2 dissolution into the aqueous phase as well as the subsequent geochemical reactions in the aqueous phase with rock formations. Could you give me some pointers on how to do this?
I tried using the GeoTes example as a template, though the Geochemistry file setup was straightforward, the porous flow file setup has not been.

[GiudGiud]

Hello

This should be a relevant example I think
https://mooseframework.inl.gov/modules/porous_flow/co2_intercomparison.html

Guillaume

[Femi-Oyenowo]

Not exactly. Though it does have some elements necessary for the setup.
The fluid systems contains 7 components and 2 phases. Phase 0 (Aqueous phase) should contain all 7 components, while phase 1 (gaseous phase) should contain only components 0 and 1 (H2O and CO2). There is only one source term and that is for CO2.
The goal is to use the multiapp transfer to send the change in components at each node to a geochemistry app, to simulate the necessary reactions.
This is what my input file looks like:

[Mesh]
  type = GeneratedMesh
  dim = 2
  xmin = 0
  xmax = 100
  ymin = 0
  ymax = 100
  nx = 10
  ny = 10
[]


[GlobalParams]
    PorousFlowDictator = dictator
    gravity = '0 0 -9.8'
    point = '1 1 0'
    # block = 0
  []

[AuxVariables]
  [xnacl]
    initial_condition = 0.0
  []
    [pgasphase]
        order = CONSTANT
        family = MONOMIAL
    []
    [pwater]
        order = CONSTANT
        family = MONOMIAL
        []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [CO2massfrac_aq]
    order = CONSTANT
    family = MONOMIAL
  []
  [H2Omassfrac_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pc]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
    [pwater]
        type = PorousFlowPropertyAux
        variable = pwater
        property = pressure
        phase = 0
        execute_on = timestep_end
    []
    [pgasphase]
        type = PorousFlowPropertyAux
        variable = pgasphase
        property = pressure
        phase = 1
        execute_on = timestep_end
    []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [CO2massfrac_aq]
    type = PorousFlowPropertyAux
    variable = CO2massfrac_aq
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [H2Omassfrac_gas]
    type = PorousFlowPropertyAux
    variable = H2Omassfrac_gas
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
  [pc]
      type = PorousFlowPropertyAux
      variable = pc
      property = capillary_pressure
      execute_on = 'initial timestep_end'
      gas_phase = 1
      liquid_phase = 0
    []
[]

[Variables]
    [pgas]
    []
    [zi]
    initial_condition = 0.2 #the brine contains some carbonate ions to start with hence this value for CO2 fraction
    []
    [f_H]
        initial_condition= -3.25058E-06
    []
    [f_Cl]
        initial_condition= 0.001869559
    []
    [f_SiO2]
        initial_condition= 1.65565E-05
    []
    [f_Al]
        initial_condition=7.e-07
    []
    [f_Ca]
        initial_condition=1.169839
    []
  []
[ICs]
    [pgas]
      type = FunctionIC
      variable = pgas
      function = '(1000*9.81*(100-y))'
    []
  []

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

  [mass3]
    type = PorousFlowMassTimeDerivative
    fluid_component = 3
    variable = f_Cl
  []
  [flux3]
    type = PorousFlowAdvectiveFlux
    fluid_component = 3
    variable = f_Cl
  []

  [mass4]
    type = PorousFlowMassTimeDerivative
    fluid_component = 4
    variable = f_SiO2
  []
  [flux4]
    type = PorousFlowAdvectiveFlux
    fluid_component = 4
    variable = f_SiO2
  []

  [mass5]
    type = PorousFlowMassTimeDerivative
    fluid_component = 5
    variable = f_Al
  []
  [flux5]
    type = PorousFlowAdvectiveFlux
    fluid_component = 5
    variable = f_Al
  []

  [mass6]
    type = PorousFlowMassTimeDerivative
    fluid_component = 6
    variable = f_Ca
  []
  [flux6]
    type = PorousFlowAdvectiveFlux
    fluid_component = 6
    variable = f_Ca
  []
[]

[UserObjects]
    [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi f_H f_Cl f_SiO2 f_Al f_Ca'
    number_fluid_phases = 2
    number_fluid_components = 7
    []
    [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 0.491
    m = 0.45946
    sat_lr = 0.3
    []
    [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
    []
    [injected_mass]
    type = PorousFlowSumQuantity
    []
[]

[FluidProperties]
    [co2sw]
    type = CO2FluidProperties
    []
    [co2]
    type = TabulatedBicubicFluidProperties
    fp = co2sw
    fluid_property_file = Sco2_fluid_properties.csv
    interpolated_properties = 'density enthalpy viscosity'
    temperature_min = 300
    temperature_max = 400
    pressure_min = 10e6
    pressure_max = 80e6
    num_T = 50
    num_p = 200
    []
    [water]
    type = Water97FluidProperties
    []
    [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    fluid_property_file = water_fluid_properties.csv
    interpolated_properties = 'density enthalpy viscosity'
    temperature_min = 300
    temperature_max = 400
    pressure_min = 10e6
    pressure_max = 80e6
    num_T = 50
    num_p = 200
    []
    [brine]
    type = BrineFluidProperties
    water_fp = watertab
    []
[]


[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 66.806
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    xnacl = xnacl
    temperature_unit = Celsius
    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
  []
  # [massfrac]
  #   type = PorousFlowMassFraction
  #   mass_fraction_vars = 'f_H f_Cl f_SiO2 f_Al f_Ca CO2massfrac_aq'
  # []
[]


[DiracKernels]
[wells]
    type  = PorousFlowPolyLineSink
    SumQuantityUO = injected_mass
    line_length = 1
    point_file = injection.bh
    variable = zi
    function_of = pressure
    fluid_phase = 1
    p_or_t_vals = 0
    fluxes = -0.15
    multiplying_var = 1
[]
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 9.4608e10
  n_startup_steps = 1
  automatic_scaling = true
  compute_scaling_once = false
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
  []
[]
 [VectorPostprocessors]
   [line]
     type = NodalValueSampler
     sort_by = x
     variable = 'pgas zi f_H f_Cl f_SiO2 f_Al f_Ca'
     execute_on = 'timestep_end'
 []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    variable = pgas
  []
  [pgasphase]
    type = PointValue
    variable = pgasphase
  []
  [pwater]
    type = PointValue
    variable = pwater
  []
  [sgas]
    type = PointValue
    variable = saturation_gas
  []
  [zi]
    type = PointValue
    variable = zi
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
  [CO2massfrac_aq]
    type = PointValue
    variable = CO2massfrac_aq
  []
  [H2Omassfrac_gas]
    type = PointValue
    variable = H2Omassfrac_gas
  []
  [pls_report]
    type = PorousFlowPlotQuantity
    uo = injected_mass
  []
  [capillary_pressure]
    type = PointValue
    variable = pc
  []
[]

[Outputs]
  exodus = true
  print_linear_residuals = false
  [csvout]
    type = CSV
    execute_on = timestep_end
    execute_vector_postprocessors_on = final
  []
[]

[Femi-Oyenowo]

The closer example would be this https://mooseframework.inl.gov/modules/geochemistry/tests_and_examples/geotes_weber_tensleep.html
But I need to have a multiphase system combined with the geochemistry.

[GiudGiud]

@cpgr

[cpgr]

Yeah, this is theoretically possible but I don't believe I've done it. I would be interested in having a go, but am not able to look at it this week

[Femi-Oyenowo]

@cpgr Okay. Thank you. Any chance you could give some theoretical details on how you expect the coupling to be done and I could work on it.