Multi-phase, multi-component radial injection

Theis pumping tests are described in the page on pumping tests and the page on point sources. This page describes similar tests, but in the multi-phase, multi-component setting.

The similarity solution

Flow from a source in a radial 1D problem admits a similarity solution that is valid for a broad range of problems, including multi-component, multi-phase flow. Properties such as fluid pressure, saturation and component mass fraction can all be characterised by the similarity variable , where is the radial distance from the source, and is time.

This similarity solution can be used to verify simple 1D radial models of increasing complexity.

Two phase immiscible flow

The simplest test case that features a similarity solution is injection of a gas phase into a fully liquid saturated model at a constant rate. This is described here.

Two phase miscible flow

The similarity solution is also recovered even when including mutual dissolution of the gas and aqueous phases. In this example, CO is injected into a fully saturated reservoir at a constant rate, with the phase conditions calculated using the WaterNCG fluid state material that implements the persistent primary variable (in this case, the total mass fraction of CO summed over all phases). Initially, there is insufficient CO to form a gas phase, with all injected CO dissolving into the resident water.

# Two phase Theis problem: Flow from single source using WaterNCG fluidstate.
# Constant rate injection 2 kg/s
# 1D cylindrical mesh
# Initially, system has only a liquid phase, until enough gas is injected
# to form a gas phase, in which case the system becomes two phase.
# Note: this test is the same as theis.i, but uses the tabulated version of the CO2FluidProperties

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 80
  xmax = 200
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [saturation_gas]
    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
  []
  [x1]
    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
  []
  [y0]
    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
  []
[]

[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
  []
  [x1]
    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"}>>> = x1
    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
  []
  [y0]
    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"}>>> = y0
    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"}>>> = 1
    fluid_component<<<{"description": "The index of the fluid component 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
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pgas]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20e6
  []
  [zi]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
[]

[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
  []
  [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
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pgas zi'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 2
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst<<<{"description": "Constant capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureConst.html"}>>>
    pc<<<{"description": "Constant capillary pressure (Pa). Default is 0"}>>> = 0
  []
  [fs]
    type = PorousFlowWaterNCG<<<{"description": "Fluid state class for water and non-condensable gas", "href": "../../../../source/userobjects/PorousFlowWaterNCG.html"}>>>
    water_fp<<<{"description": "The name of the user object for water"}>>> = water
    gas_fp<<<{"description": "The name of the user object for the non-condensable gas"}>>> = tabulated
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
[]

[FluidProperties<<<{"href": "../../../../syntax/FluidProperties/index.html"}>>>]
  [co2]
    type = CO2FluidProperties<<<{"description": "Fluid properties for carbon dioxide (CO2) using the Span & Wagner EOS", "href": "../../../../source/fluidproperties/CO2FluidProperties.html"}>>>
  []
  [tabulated]
    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"}>>> = co2
    fluid_property_file<<<{"description": "Name of the csv file containing the tabulated fluid property data."}>>> = fluid_properties.csv
    # We try to avoid using both, but some properties are not implemented in the tabulation
    allow_fp_and_tabulation<<<{"description": "Whether to allow the two sources of data concurrently"}>>> = true
    # Test was design prior to bounds check
    error_on_out_of_bounds<<<{"description": "Whether pressure or temperature from tabulation exceeding user-specified bounds leads to an error."}>>> = false
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = fluid_properties.csv
  []
  [water]
    type = Water97FluidProperties<<<{"description": "Fluid properties for water and steam (H2O) using IAPWS-IF97", "href": "../../../../source/fluidproperties/Water97FluidProperties.html"}>>>
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
    temperature<<<{"description": "Fluid temperature variable.  Note, the default is suitable if your simulation is using Kelvin units, but probably not for Celsius"}>>> = 20
  []
  [waterncg]
    type = PorousFlowFluidState<<<{"description": "Class for fluid state calculations using persistent primary variables and a vapor-liquid flash", "href": "../../../../source/materials/PorousFlowFluidState.html"}>>>
    gas_porepressure<<<{"description": "Variable that is the porepressure of the gas phase"}>>> = pgas
    z<<<{"description": "Total mass fraction of component i summed over all phases"}>>> = zi
    temperature_unit<<<{"description": "The unit of the temperature variable"}>>> = Celsius
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
    fluid_state<<<{"description": "Name of the FluidState UserObject"}>>> = fs
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
    s_res<<<{"description": "The residual saturation of the phase j. Must be between 0 and 1"}>>> = 0.1
    sum_s_res<<<{"description": "Sum of residual saturations over all phases.  Must be between 0 and 1"}>>> = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 1
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [rightwater]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 20e6
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [source]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = 2
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
    petsc_options<<<{"description": "Singleton PETSc options"}>>> = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap -snes_atol -snes_rtol -snes_max_it'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres      asm      lu           NONZERO                   2               1E-8       1E-10 20'
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = NEWTON
  end_time = 8e2
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 2
    growth_factor = 2
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [line]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    warn_discontinuous_face_values<<<{"description": "Whether to return a warning if a discontinuous variable is sampled on a face"}>>> = false
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '200 0 0'
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 1000
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'pgas zi x1 saturation_gas'
    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'
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [pgas]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pgas
  []
  [sgas]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = saturation_gas
  []
  [zi]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = zi
  []
  [massgas]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component that this Postprocessor acts on"}>>> = 1
  []
  [x1]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = x1
  []
  [y0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '1 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = y0
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  [csvout]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    file_base<<<{"description": "The desired solution output name without an extension. If not provided, MOOSE sets it with Outputs/file_base when available. Otherwise, MOOSE uses input file name and this object name for a master input or uses master file_base, the subapp name and this object name for a subapp input to set it."}>>> = theis_tabulated_csvout
    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
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = final
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis_tabulated.i)

After approximately 3000s, enough CO has been injected to form a gas phase, and the problem evolves from a single liquid phase to one including both gas and liquid phases. Figure 1 shows the comparison of similarity solutions calculated with either fixed radial distance or fixed time. In this case, good agreement is observed between the two results for both primary variables (gas pressure and total mass fraction of CO). The similarity solution is also observed in Figure 2 for both the calculated gas saturation and the mass fraction of CO dissolved in the water.

Since this test is part of the automatic test suite, it needs to run rapidly. Better agreement between the similarity solutions are obtained by refining the spatial and temporal resolution.

Figure 1: Comparison of similarity solutions for the water-NCG fluid state. (a) Gas pressure; (b) Total mass fraction of CO2.

Figure 2: Comparison of similarity solutions for the water-NCG fluid state. (a) Gas saturation; (b) Dissolved CO2 mass fraction.

CO injection into brine

Similar results are obtained for CO injection in brine, using the BrineCO2 fluid state material. In this example, CO is injected into a reservoir at a constant rate. Initially, the reservoir is fully saturated with brine which has a salt mass fraction of 0.1. During the early injection period, there is insufficient CO to form a gas phase, with all injected CO dissolving into the resident brine.

# Two phase Theis problem: Flow from single source.
# Constant rate injection 2 kg/s
# 1D cylindrical mesh
# Initially, system has only a liquid phase, until enough gas is injected
# to form a gas phase, in which case the system becomes two phase.
#
# This test takes a few minutes to run, so is marked heavy

[Mesh<<<{"href": "../../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 2000
  xmax = 2000
[]

[Problem<<<{"href": "../../../../syntax/Problem/index.html"}>>>]
  type = FEProblem
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams<<<{"href": "../../../../syntax/GlobalParams/index.html"}>>>]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables<<<{"href": "../../../../syntax/AuxVariables/index.html"}>>>]
  [saturation_gas]
    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
  []
  [x1]
    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
  []
  [y0]
    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
  []
[]

[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
  []
  [x1]
    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"}>>> = x1
    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
  []
  [y0]
    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"}>>> = y0
    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"}>>> = 1
    fluid_component<<<{"description": "The index of the fluid component 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
  []
[]

[Variables<<<{"href": "../../../../syntax/Variables/index.html"}>>>]
  [pgas]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 20e6
  []
  [zi]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0
  []
  [xnacl]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 0.1
  []
[]

[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
  []
  [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
  []
  [mass2]
    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"}>>> = 2
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = xnacl
  []
  [flux2]
    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"}>>> = 2
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = xnacl
  []
[]

[UserObjects<<<{"href": "../../../../syntax/UserObjects/index.html"}>>>]
  [dictator]
    type = PorousFlowDictator<<<{"description": "Holds information on the PorousFlow variable names", "href": "../../../../source/userobjects/PorousFlowDictator.html"}>>>
    porous_flow_vars<<<{"description": "List of primary variables that are used in the PorousFlow simulation.  Jacobian entries involving derivatives wrt these variables will be computed.  In single-phase models you will just have one (eg 'pressure'), in two-phase models you will have two (eg 'p_water p_gas', or 'p_water s_water'), etc."}>>> = 'pgas zi xnacl'
    number_fluid_phases<<<{"description": "The number of fluid phases in the simulation"}>>> = 2
    number_fluid_components<<<{"description": "The number of fluid components in the simulation"}>>> = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureConst<<<{"description": "Constant capillary pressure", "href": "../../../../source/userobjects/PorousFlowCapillaryPressureConst.html"}>>>
    pc<<<{"description": "Constant capillary pressure (Pa). Default is 0"}>>> = 0
  []
  [fs]
    type = PorousFlowBrineCO2<<<{"description": "Fluid state class for brine and CO2", "href": "../../../../source/userobjects/PorousFlowBrineCO2.html"}>>>
    brine_fp<<<{"description": "The name of the user object for brine"}>>> = brine
    co2_fp<<<{"description": "The name of the user object for CO2"}>>> = co2
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
[]

[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 = TabulatedFluidProperties<<<{"description": "Fluid properties using bicubic interpolation on tabulated values provided", "href": "../../../../source/fluidproperties/TabulatedBicubicFluidProperties.html"}>>>
    fp<<<{"description": "The name of the FluidProperties UserObject"}>>> = co2sw
    fluid_property_file<<<{"description": "Name of the csv file containing the tabulated fluid property data."}>>> = 'fluid_properties.csv'
    allow_fp_and_tabulation<<<{"description": "Whether to allow the two sources of data concurrently"}>>> = true
    error_on_out_of_bounds<<<{"description": "Whether pressure or temperature from tabulation exceeding user-specified bounds leads to an error."}>>> = false
  []
  [water]
    type = Water97FluidProperties<<<{"description": "Fluid properties for water and steam (H2O) using IAPWS-IF97", "href": "../../../../source/fluidproperties/Water97FluidProperties.html"}>>>
  []
  [watertab]
    type = TabulatedFluidProperties<<<{"description": "Fluid properties using bicubic interpolation on tabulated values provided", "href": "../../../../source/fluidproperties/TabulatedBicubicFluidProperties.html"}>>>
    fp<<<{"description": "The name of the FluidProperties UserObject"}>>> = water
    temperature_min<<<{"description": "Minimum temperature for tabulated data."}>>> = 273.15
    temperature_max<<<{"description": "Maximum temperature for tabulated data."}>>> = 573.15
    fluid_property_output_file<<<{"description": "Name of the CSV file which can be output with the tabulation. This file can then be read as a 'fluid_property_file'"}>>> = water_fluid_properties.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water_fluid_properties.csv
  []
  [brine]
    type = BrineFluidProperties<<<{"description": "Fluid properties for brine", "href": "../../../../source/fluidproperties/BrineFluidProperties.html"}>>>
    water_fp<<<{"description": "The name of the FluidProperties UserObject for water"}>>> = watertab
  []
[]

[Materials<<<{"href": "../../../../syntax/Materials/index.html"}>>>]
  [temperature]
    type = PorousFlowTemperature<<<{"description": "Material to provide temperature at the quadpoints or nodes and derivatives of it with respect to the PorousFlow variables", "href": "../../../../source/materials/PorousFlowTemperature.html"}>>>
    temperature<<<{"description": "Fluid temperature variable.  Note, the default is suitable if your simulation is using Kelvin units, but probably not for Celsius"}>>> = 20
  []
  [brineco2]
    type = PorousFlowFluidState<<<{"description": "Class for fluid state calculations using persistent primary variables and a vapor-liquid flash", "href": "../../../../source/materials/PorousFlowFluidState.html"}>>>
    gas_porepressure<<<{"description": "Variable that is the porepressure of the gas phase"}>>> = pgas
    z<<<{"description": "Total mass fraction of component i summed over all phases"}>>> = zi
    temperature_unit<<<{"description": "The unit of the temperature variable"}>>> = Celsius
    xnacl<<<{"description": "The salt mass fraction in the brine (kg/kg)"}>>> = xnacl
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
    fluid_state<<<{"description": "Name of the FluidState UserObject"}>>> = fs
  []
  [porosity]
    type = PorousFlowPorosityConst<<<{"description": "This Material calculates the porosity assuming it is constant", "href": "../../../../source/materials/PorousFlowPorosityConst.html"}>>>
    porosity<<<{"description": "The porosity (assumed indepenent of porepressure, temperature, strain, etc, for this material).  This should be a real number, or a constant monomial variable (not a linear lagrange or other kind of variable)."}>>> = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst<<<{"description": "This Material calculates the permeability tensor assuming it is constant", "href": "../../../../source/materials/PorousFlowPermeabilityConst.html"}>>>
    permeability<<<{"description": "The permeability tensor (usually in m^2), which is assumed constant for this material"}>>> = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 0
    s_res<<<{"description": "The residual saturation of the phase j. Must be between 0 and 1"}>>> = 0.1
    sum_s_res<<<{"description": "Sum of residual saturations over all phases.  Must be between 0 and 1"}>>> = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey<<<{"description": "This Material calculates relative permeability of the fluid phase, using the simple Corey model ((S-S_res)/(1-sum(S_res)))^n", "href": "../../../../source/materials/PorousFlowRelativePermeabilityCorey.html"}>>>
    n<<<{"description": "The Corey exponent of the phase."}>>> = 2
    phase<<<{"description": "The phase number"}>>> = 1
  []
[]

[BCs<<<{"href": "../../../../syntax/BCs/index.html"}>>>]
  [rightwater]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../../source/bcs/DirichletBC.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = 20e6
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = pgas
  []
[]

[DiracKernels<<<{"href": "../../../../syntax/DiracKernels/index.html"}>>>]
  [source]
    type = PorousFlowSquarePulsePointSource<<<{"description": "Point source (or sink) that adds (removes) fluid at a constant mass flux rate for times between the specified start and end times.", "href": "../../../../source/dirackernels/PorousFlowSquarePulsePointSource.html"}>>>
    point<<<{"description": "The x,y,z coordinates of the point source (sink)"}>>> = '0 0 0'
    mass_flux<<<{"description": "The mass flux at this point in kg/s (positive is flux in, negative is flux out)"}>>> = 2
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
  []
[]

[Preconditioning<<<{"href": "../../../../syntax/Preconditioning/index.html"}>>>]
  [smp]
    type = SMP<<<{"description": "Single matrix preconditioner (SMP) builds a preconditioner using user defined off-diagonal parts of the Jacobian.", "href": "../../../../source/preconditioners/SingleMatrixPreconditioner.html"}>>>
    full<<<{"description": "Set to true if you want the full set of couplings between variables simply for convenience so you don't have to set every off_diag_row and off_diag_column combination."}>>> = true
  []
[]

[Executioner<<<{"href": "../../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = NEWTON
  end_time = 1e5
  [TimeStepper<<<{"href": "../../../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.5
  []
[]

[VectorPostprocessors<<<{"href": "../../../../syntax/VectorPostprocessors/index.html"}>>>]
  [line]
    type = LineValueSampler<<<{"description": "Samples variable(s) along a specified line", "href": "../../../../source/vectorpostprocessors/LineValueSampler.html"}>>>
    warn_discontinuous_face_values<<<{"description": "Whether to return a warning if a discontinuous variable is sampled on a face"}>>> = false
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    start_point<<<{"description": "The beginning of the line"}>>> = '0 0 0'
    end_point<<<{"description": "The ending of the line"}>>> = '2000 0 0'
    num_points<<<{"description": "The number of points to sample along the line"}>>> = 10000
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'pgas zi xnacl x1 saturation_gas'
    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'
  []
[]

[Postprocessors<<<{"href": "../../../../syntax/Postprocessors/index.html"}>>>]
  [pgas]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = pgas
  []
  [sgas]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = saturation_gas
  []
  [zi]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = zi
  []
  [massgas]
    type = PorousFlowFluidMass<<<{"description": "Calculates the mass of a fluid component in a region", "href": "../../../../source/postprocessors/PorousFlowFluidMass.html"}>>>
    fluid_component<<<{"description": "The index corresponding to the fluid component that this Postprocessor acts on"}>>> = 1
  []
  [x1]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = x1
  []
  [y0]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = y0
  []
  [xnacl]
    type = PointValue<<<{"description": "Compute the value of a variable at a specified location", "href": "../../../../source/postprocessors/PointValue.html"}>>>
    point<<<{"description": "The physical point where the solution will be evaluated."}>>> = '4 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = xnacl
  []
[]

[Outputs<<<{"href": "../../../../syntax/Outputs/index.html"}>>>]
  print_linear_residuals<<<{"description": "Enable printing of linear residuals to the screen (Console)"}>>> = false
  perf_graph<<<{"description": "Enable printing of the performance graph to the screen (Console)"}>>> = true
  [csvout]
    type = CSV<<<{"description": "Output for postprocessors, vector postprocessors, and scalar variables using comma seperated values (CSV).", "href": "../../../../source/outputs/CSV.html"}>>>
    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
    execute_vector_postprocessors_on<<<{"description": "Enable/disable the output of VectorPostprocessors"}>>> = final
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis_brineco2.i)

After approximately 2000s, enough CO has been injected to form a gas phase, and the problem evolves from a single liquid phase to one including both gas and liquid phases. Note that this time is earlier than the problem with water instead of brine, as the equilibrium saturation of CO in brine is less than that in water. Figure 3 shows the comparison of similarity solutions calculated with either fixed radial distance or fixed time. In this case, good agreement is observed between the two results for both primary variables (gas pressure and total mass fraction of CO). The similarity solution is also observed in Figure 4 for both the calculated gas saturation and the mass fraction of CO dissolved in the water.

This example is included in the automatic test suite, but is marked heavy due to the run time.

Figure 3: Comparison of similarity solutions for the brine-CO2 fluid state. (a) Gas pressure; (b) Total mass fraction of CO2; (c) NaCl mass fraction.

Figure 4: Comparison of similarity solutions for the brine-CO2 fluid state. (a) Gas saturation; (b) Dissolved CO2 mass fraction.