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.

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

[AuxVariables<<<{"href": "../../syntax/AuxVariables/index.html"}>>>]
  [xnacl]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.01
  []
  [saturation_gas]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [xco2l]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [density_liquid]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = FIRST
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
  [porosity]
    order<<<{"description": "Specifies the order of the FE shape function to use for this variable (additional orders not listed are allowed)"}>>> = CONSTANT
    family<<<{"description": "Specifies the family of FE shape functions to use for this variable"}>>> = MONOMIAL
  []
[]

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

and the temperature at 45C:

  [temperature]
    type = PorousFlowTemperature
    temperature = '45'
  []

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<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [pressure]
    type = FunctionIC<<<{"description": "An initial condition that uses a normal function of x, y, z to produce values (and optionally gradients) for a field variable.", "href": "../../source/ics/FunctionIC.html"}>>>
    function<<<{"description": "The initial condition function."}>>> = 10e6-9.81*1000*y
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = pgas
  []
  [zi]
    type = BoundingBoxIC<<<{"description": "BoundingBoxIC allows setting the initial condition of a value inside and outside of a specified box. The box is aligned with the x, y, z axes", "href": "../../source/ics/BoundingBoxIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = zi
    x1<<<{"description": "The x coordinate of the lower left-hand corner of the box"}>>> = 0
    x2<<<{"description": "The x coordinate of the upper right-hand corner of the box"}>>> = 2
    y1<<<{"description": "The y coordinate of the lower left-hand corner of the box"}>>> = 1.95
    y2<<<{"description": "The y coordinate of the upper right-hand corner of the box"}>>> = 2
    inside<<<{"description": "The value of the variable inside the box"}>>> = 0.2
    outside<<<{"description": "The value of the variable outside the box"}>>> = 0
  []
  [porosity]
    type = RandomIC<<<{"description": "Initialize a variable with randomly generated numbers following either a uniform distribution or a user-defined distribution", "href": "../../source/ics/RandomIC.html"}>>>
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = porosity
    min<<<{"description": "Lower bound of uniformly distributed randomly generated values"}>>> = 0.25
    max<<<{"description": "Upper bound of uniformly distributed randomly generated values"}>>> = 0.275
    seed<<<{"description": "Seed value for the random number generator"}>>> = 0
  []
[]

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:

[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
  [co2sw]
    type = CO2FluidProperties<<<{"description": "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS", "href": "../../source/fluidproperties/CO2FluidProperties.html"}>>>
  []
  [co2]
    type = TabulatedBicubicFluidProperties<<<{"description": "Fluid properties using bicubic interpolation on tabulated values provided", "href": "../../source/fluidproperties/TabulatedBicubicFluidProperties.html"}>>>
    fp<<<{"description": "The name of the FluidProperties UserObject"}>>> = co2sw
  []
  [brine]
    type = BrineFluidProperties<<<{"description": "Fluid properties for brine", "href": "../../source/fluidproperties/BrineFluidProperties.html"}>>>
  []
[]

  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc

  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = 'pgas'
    z = 'zi'
    temperature_unit = Celsius
    xnacl = 'xnacl'
    capillary_pressure = pc
    fluid_state = fs
  []

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

[AuxKernels<<<{"href": "../../syntax/AuxKernels/index.html"}>>>]
  [saturation_gas]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = saturation_gas
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = saturation
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
  []
  [xco2l]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = xco2l
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = mass_fraction
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 0
    fluid_component<<<{"description": "The index of the fluid component this auxillary kernel acts on"}>>> = 1
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
  []
  [density_liquid]
    type = PorousFlowPropertyAux<<<{"description": "AuxKernel to provide access to properties evaluated at quadpoints. Note that elemental AuxVariables must be used, so that these properties are integrated over each element.", "href": "../../source/auxkernels/PorousFlowPropertyAux.html"}>>>
    variable<<<{"description": "The name of the variable that this object applies to"}>>> = density_liquid
    property<<<{"description": "The fluid property that this auxillary kernel is to calculate"}>>> = density
    phase<<<{"description": "The index of the phase this auxillary kernel acts on"}>>> = 0
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = 'timestep_end'
  []
[]

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

[Kernels<<<{"href": "../../syntax/Kernels/index.html"}>>>]
  [mass0]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [flux0]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
  [diff0]
    type = PorousFlowDispersiveFlux<<<{"description": "Dispersive and diffusive flux of the component given by fluid_component in all phases", "href": "../../source/kernels/PorousFlowDispersiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 0
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
    disp_long<<<{"description": "Vector of longitudinal dispersion coefficients for each phase"}>>> = '0 0'
    disp_trans<<<{"description": "Vector of transverse dispersion coefficients for each phase"}>>> = '0 0'
  []
  [mass1]
    type = PorousFlowMassTimeDerivative<<<{"description": "Derivative of fluid-component mass with respect to time.  Mass lumping to the nodes is used.", "href": "../../source/kernels/PorousFlowMassTimeDerivative.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
  []
  [flux1]
    type = PorousFlowAdvectiveFlux<<<{"description": "Fully-upwinded advective flux of the component given by fluid_component", "href": "../../source/kernels/PorousFlowAdvectiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
  []
  [diff1]
    type = PorousFlowDispersiveFlux<<<{"description": "Dispersive and diffusive flux of the component given by fluid_component in all phases", "href": "../../source/kernels/PorousFlowDispersiveFlux.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component for this kernel"}>>> = 1
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
    disp_long<<<{"description": "Vector of longitudinal dispersion coefficients for each phase"}>>> = '0 0'
    disp_trans<<<{"description": "Vector of transverse dispersion coefficients for each phase"}>>> = '0 0'
  []
[]

  [diffusivity]
    type = PorousFlowDiffusivityConst
    diffusion_coeff = '2e-9 2e-9 2e-9 2e-9'
    tortuosity = '1 1'

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

[Postprocessors<<<{"href": "../../syntax/Postprocessors/index.html"}>>>]
  [total_co2_in_gas]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    phase<<<{"description": "The index of the fluid phase that this Postprocessor is restricted to.  Multiple indices can be entered"}>>> = 1
    fluid_component<<<{"description": "The index corresponding to the fluid component that this Postprocessor acts on"}>>> = 1
  []
  [total_co2_in_liquid]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    phase<<<{"description": "The index of the fluid phase that this Postprocessor is restricted to.  Multiple indices can be entered"}>>> = 0
    fluid_component<<<{"description": "The index corresponding to the fluid component that this Postprocessor acts on"}>>> = 1
  []
  [numdofs]
    type = NumDOFs<<<{"description": "Return the number of Degrees of freedom from either the NL, Aux or both systems.", "href": "../../source/postprocessors/NumDOFs.html"}>>>
  []
  [delta_xco2]
    type = ChangeOverTimePostprocessor<<<{"description": "Computes the change or relative change in a post-processor value over a timestep or the entire transient", "href": "../../source/postprocessors/ChangeOverTimePostprocessor.html"}>>>
    postprocessor<<<{"description": "The name of the postprocessor"}>>> = total_co2_in_liquid
  []
  [dt]
    type = TimestepSize<<<{"description": "Reports the timestep size", "href": "../../source/postprocessors/TimestepSize.html"}>>>
  []
  [flux]
    type = FunctionValuePostprocessor<<<{"description": "Computes the value of a supplied function at a single point (scalable)", "href": "../../source/postprocessors/FunctionValuePostprocessor.html"}>>>
    function<<<{"description": "The function which supplies the postprocessor value."}>>> = flux
  []
[]

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<<<{"href": "../../syntax/Adaptivity/index.html"}>>>]
  max_h_level<<<{"description": "Maximum number of times a single element can be refined. If 0 then infinite."}>>> = 2
  marker<<<{"description": "The name of the Marker to use to actually adapt the mesh."}>>> = marker
  initial_marker<<<{"description": "The name of the Marker to use to adapt the mesh during initial refinement."}>>> = initial
  initial_steps<<<{"description": "The number of adaptive steps to do based on the initial condition."}>>> = 2
  [Indicators<<<{"href": "../../syntax/Adaptivity/Indicators/index.html"}>>>]
    [indicator]
      type = GradientJumpIndicator<<<{"description": "Compute the jump of the solution gradient across element boundaries.", "href": "../../source/indicators/GradientJumpIndicator.html"}>>>
      variable<<<{"description": "The name of the variable that this side indicator applies to"}>>> = zi
    []
  []
  [Markers<<<{"href": "../../syntax/Adaptivity/Markers/index.html"}>>>]
    [marker]
      type = ErrorFractionMarker<<<{"description": "Marks elements for refinement or coarsening based on the fraction of the min/max error from the supplied indicator.", "href": "../../source/markers/ErrorFractionMarker.html"}>>>
      indicator<<<{"description": "The name of the Indicator that this Marker uses."}>>> = indicator
      refine<<<{"description": "Elements within this percentage of the max error will be refined.  Must be between 0 and 1!"}>>> = 0.8
    []
    [initial]
      type = BoxMarker<<<{"description": "Marks the region inside and outside of a 'box' domain for refinement or coarsening.", "href": "../../source/markers/BoxMarker.html"}>>>
      bottom_left<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '0 1.95 0'
      top_right<<<{"description": "The bottom left point (in x,y,z with spaces in-between)."}>>> = '2 2 0'
      inside<<<{"description": "How to mark elements inside the box."}>>> = REFINE
      outside<<<{"description": "How to mark elements outside the box."}>>> = DO_NOTHING
    []
  []
[]

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.