Density-driven convective mixing

PorousFlow may be used to study the convective mixing of fluids. In this page we discuss the density-driven convective mixing of CO into brine. CO exists initially in both the liquid and gas phase. As the CO diffuses into the brine, the gaseous CO disappears and the brine's density increases, which drives convective mixing of the two fluids and accelerates the dissolution of the gas phase (Emami-Meybodi et al., 2015). Input files for similar situations may be found here.

commentnote

These examples should be run using multiple cores due to the computational expense. It is recommended that the single-phase model is run using at least four cores, and the two-phase model with at least eight cores.

Because this problem involves the disappearance of the gaseous phase, the simulation utilises persistent primary variables and a vapor-liquid flash.

[AuxVariables]
  [xnacl]
    initial_condition = 0.01
  []
  [saturation_gas]
    order = FIRST
    family = MONOMIAL
  []
  [xco2l]
    order = FIRST
    family = MONOMIAL
  []
  [density_liquid]
    order = FIRST
    family = MONOMIAL
  []
  [porosity]
    order = CONSTANT
    family = MONOMIAL
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

Here pgas is the gas-state pressure (the capillary pressure is zero, so this is also the brine porepressure) and zi is the total mass fraction of the CO summed over the two phases. In this input file, the salt content of the brine is fixed at 0.01:

  [xnacl]
    initial_condition = 0.01
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

and the temperature at 45C:

  [temperature]
    type = PorousFlowTemperature
    temperature = '45'
  []
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

The porepressure is initialised to approximate hydrostatic conditions, and the total mass fraction of CO is zero except towards the top of the model domain where it is 0.2. This corresponds to a gas saturation of 0.29 (the gas is almost all CO and only about 0.2% brine) and the mass-fraction of CO in the liquid phase is 0.048. The porosity is slightly randomised by sampling from a uniform distribution between the minimum and maximum values specified in the RandomIC initial condition, and this promotes initiation of the beautiful "fingers" of CO-rich brine seen in the results. The random field seed can be changed in the input file, or on the command line using command line overriding. The porosity field is output in the Exodus file that is created, so can be visualized in the MOOSE GUI peacock or in an external visualization package.

[ICs]
  [pressure]
    type = FunctionIC
    function = 10e6-9.81*1000*y
    variable = pgas
  []
  [zi]
    type = BoundingBoxIC
    variable = zi
    x1 = 0
    x2 = 2
    y1 = 1.95
    y2 = 2
    inside = 0.2
    outside = 0
  []
  [porosity]
    type = RandomIC
    variable = porosity
    min = 0.25
    max = 0.275
    seed = 0
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

The equations of state for the liquid brine, the CO and the liquid brine-CO mixture are the high-precision versions supplied by the fluid_properties module:

[Modules]
  [FluidProperties]
    [co2sw]
      type = CO2FluidProperties
    []
    [co2]
      type = TabulatedFluidProperties
      fp = co2sw
    []
    [brine]
      type = BrineFluidProperties
    []
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = 'pgas'
    z = 'zi'
    temperature_unit = Celsius
    xnacl = 'xnacl'
    capillary_pressure = pc
    fluid_state = fs
  []
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

Quantities of interest are tracked using AuxVariables populated using the following AuxKernels

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'timestep_end'
  []
  [xco2l]
    type = PorousFlowPropertyAux
    variable = xco2l
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = 'timestep_end'
  []
  [density_liquid]
    type = PorousFlowPropertyAux
    variable = density_liquid
    property = density
    phase = 0
    execute_on = 'timestep_end'
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

The usual 2-phase Kernels are used in this example with zero molecular dispersion (but the diffusion coefficient is nonzero)

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pgas
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pgas
  []
  [diff0]
    type = PorousFlowDispersiveFlux
    fluid_component = 0
    variable = pgas
    disp_long = '0 0'
    disp_trans = '0 0'
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = zi
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = zi
  []
  [diff1]
    type = PorousFlowDispersiveFlux
    fluid_component = 1
    variable = zi
    disp_long = '0 0'
    disp_trans = '0 0'
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)
  [diffusivity]
    type = PorousFlowDiffusivityConst
    diffusion_coeff = '2e-9 2e-9 2e-9 2e-9'
    tortuosity = '1 1'
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

Various postprocessors help with tracking the total quantity of CO in the two phases, and the flux of CO into the brine:

[Postprocessors]
  [total_co2_in_gas]
    type = PorousFlowFluidMass
    phase = 1
    fluid_component = 1
  []
  [total_co2_in_liquid]
    type = PorousFlowFluidMass
    phase = 0
    fluid_component = 1
  []
  [numdofs]
    type = NumDOFs
  []
  [delta_xco2]
    type = ChangeOverTimePostprocessor
    postprocessor = total_co2_in_liquid
  []
  [dt]
    type = TimestepSize
  []
  [flux]
    type = FunctionValuePostprocessor
    function = flux
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

Mesh adaptivity is used to refine the mesh as the fingers form, based on the jump in the gradient of the total mass fraction of CO summed over the two phases:

[Adaptivity]
  max_h_level = 2
  marker = marker
  initial_marker = initial
  initial_steps = 2
  [Indicators]
    [indicator]
      type = GradientJumpIndicator
      variable = zi
    []
  []
  [Markers]
    [marker]
      type = ErrorFractionMarker
      indicator = indicator
      refine = 0.8
    []
    [initial]
      type = BoxMarker
      bottom_left = '0 1.95 0'
      top_right = '2 2 0'
      inside = REFINE
      outside = DO_NOTHING
    []
  []
[]
(modules/porous_flow/examples/lava_lamp/2phase_convection.i)

The evolution of the CO mass fraction is shown in the animation below

Figure 1: Left: Mass fraction of CO in the liquid phase. Right: Saturation of gas. The adaptive mesh is overlayed on each figure.

References

  1. Hamid Emami-Meybodi, Hassan Hassanzadeh, Christopher P Green, and Jonathan Ennis-King. Convective dissolution of co$_2$ in saline aquifers: progress in modeling and experiments. International Journal of Greenhouse Gas Control, 40:238–266, 2015.[BibTeX]