Water and non-condensable gas

A miscible water and non-condensable gas (NCG) model is available in PorousFlow based on the persistent variable approach.

Available NCG's

Any of the gaseous fluids provided in the Fluid Properties module can be used in this equation of state.

Phase composition

Henry's law is used to calculate the mass fraction of the NCG in the water phases while Raoult's law is used to calculate the mass fraction of HO in the gas phase where is vapor pressure of water.

As with many reservoir simulators, it is assumed that the phases are in instantaneous equilibrium, meaning that dissolution of NCG into the water and evaporation of water vapor into the gas phase happens instantaneously.

If the amount of NCG present is not sufficient to saturate the water, then no gas phase will be formed, and the dissolved NCG mass fraction will be smaller than its equilibrium value. Similarly, if the amount of water vapor present is smaller than its equilibrium value in the gas phase, then no liquid phase can be formed.

If there is sufficient NCG present to saturate the water, then both liquid and gas phases will be present. In this case, the mass fraction of NCG dissolved in the water is given by its equilibrium value. Similarly, the mass fraction of water vapor in the gas phase will be given by its equilibrium value in this two phase region.

Thermophysical properties

Density

The density of the gas phase is calculated as the sum of the densities of each fluid component at their partial pressure (via Dalton's law).

No change to the water density is calculated for dissolved NCG components.

Viscosity

The viscosity of the gas phase is given by a weighted sum of the viscosities of each fluid component

The viscosity of the liquid phase is simply given by the viscosity of water - no contribution due to dissolved NCG is included.

Enthalpy

The enthalpy of the gas phase is calculated as a weighted sum of the enthalpies of each fluid component

The enthalpy of the liquid phase is also calculated as a weighted sum of the enthalpies of each individual component, but also includes a contribution due to the enthalpy of dissolution (a change in enthalpy due to dissolution of the NCG into the water phase)

where is the enthalpy of the gaseous NCG, and is the enthalpy of dissolution of NCG into the liquid. An expression for the enthalpy of dissolution is implemented following (Himmelblau, 1959)

where is the molar mass of the NCG component, and is the universal gas constant.

Implementation

Variables

This class requires pressure of the gas phase and the total mass fraction of NCG summed over all phases (see compositional flash for details) for isothermal simulations. For nonisothermal simulations, the temperature must also be provided as a nonlinear variable.

commentnote

These variables must be listed as PorousFlow variables in the PorousFlowDictator UserObject

In the isothermal case, two variables (gas porepressure and total NCG mass fraction) are required. The number of components in the PorousFlowDictator must be set equal to two.

[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
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis.i)

In the nonisothermal case, three variables (gas porepressure, total NCG mass fraction and temperature) are required. The number of components in the PorousFlowDictator is still two.

[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 temperature'
    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
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis_nonisothermal.i)

UserObjects

This fluid state is implemented in the PorousFlowWaterNCG UserObject. This UserObject is a GeneralUserObject that contains methods that provide a complete thermophysical description of the model given pressure, temperature and total NCG mass fraction.

[UserObjects<<<{"href": "../../syntax/UserObjects/index.html"}>>>]
  [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"}>>> = co2
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis.i)

Swapping NCG's is as simple as changing the gas_fp parameter in this UserObject!

Materials

The PorousFlowFluidState material provides all phase pressures, saturation, densities, viscosities etc, as well as all mass fractions of all fluid components in all fluid phases in a single material using the formulation provided in the PorousFlowWaterNCG UserObject.

For isothermal simulations, this material will look like

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [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
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis.i)

For nonisothermal simulations, the temperature variable must also be supplied to ensure that all the derivatives with respect to the temperature variable are computed for the Jacobian.

[Materials<<<{"href": "../../syntax/Materials/index.html"}>>>]
  [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<<<{"description": "The fluid temperature (C or K, depending on temperature_unit)"}>>> = temperature
    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
  []
[]
(modules/porous_flow/test/tests/fluidstate/theis_nonisothermal.i)

Initial condition

The nonlinear variable representing NCG is the total mass fraction of NCG summed over all phases. In some cases, it may be preferred to provide an initial NCG saturation, rather than total mass fraction. To allow an initial saturation to be specified, the PorousFlowFluidStateIC initial condition is provided. This initial condition calculates the total mass fraction of NCG summed over all phases for a given saturation.

[ICs<<<{"href": "../../syntax/ICs/index.html"}>>>]
  [z]
    type = PorousFlowFluidStateIC<<<{"description": "An initial condition to calculate z from saturation", "href": "../../source/ics/PorousFlowFluidStateIC.html"}>>>
    saturation<<<{"description": "Gas saturation"}>>> = 0.5
    gas_porepressure<<<{"description": "Variable that is the porepressure of the gas phase"}>>> = pgas
    temperature<<<{"description": "Variable that is the fluid temperature"}>>> = 50
    variable<<<{"description": "The variable this initial condition is supposed to provide values for."}>>> = z
    fluid_state<<<{"description": "Name of the FluidState UserObject"}>>> = fs
  []
[]
(modules/porous_flow/test/tests/fluidstate/waterncg_ic.i)

Example

1D Radial injection

Injection into a one dimensional radial reservoir is a useful test of this class, as it admits a similarity solution , where is the radial distance from the injection well, and is time. This similarity solution holds even when complicated physics such as mutual solubility is included.

An example of an isothermal 1D radial injection problem is included in the test suite with CO as the NCG.

# 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.

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 1
  nx = 40
  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"}>>> = co2
    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"}>>>
  []
  [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 = 2e2
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 10
    growth_factor = 2
  []
[]

[VectorPostprocessors<<<{"href": "../../syntax/VectorPostprocessors/index.html"}>>>]
  [line]
    type = NodalValueSampler<<<{"description": "Samples values of nodal variable(s).", "href": "../../source/vectorpostprocessors/NodalValueSampler.html"}>>>
    sort_by<<<{"description": "What to sort the samples by"}>>> = x
    variable<<<{"description": "The names of the variables that this VectorPostprocessor operates on"}>>> = 'pgas zi'
    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"}>>>
    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.i)

Initially, all of the CO injected is dissolved in the resident water, as the amount of CO is less than the dissolved mass fraction at equilibrium. After a while, the water near the injection well becomes saturated with CO, and a gas phase forms. This gas phase then spreads radially as injection continues.

A comparison of the results expressed in terms of the similarity solution is presented in Figure 1 for both nonlinear variables for the case where is fixed (evaluated using Postprocessors), and also the case where is fixed (evaluated using a VectorPostprocessor).

Figure 1: Similarity solution for 1D radial injection example

An nonisothermal example of the above problem, where cold CO is injected into a warm aquifer, is also provided in the test suite.

# Two-phase nonisothermal Theis problem: Flow from single source using WaterNCG fluidstate.
# Constant rate injection 2 kg/s of cold gas into warm reservoir
# 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.

[Mesh<<<{"href": "../../syntax/Mesh/index.html"}>>>]
  [mesh]
    type = GeneratedMeshGenerator<<<{"description": "Create a line, square, or cube mesh with uniformly spaced or biased elements.", "href": "../../source/meshgenerators/GeneratedMeshGenerator.html"}>>>
    dim<<<{"description": "The dimension of the mesh to be generated"}>>> = 1
    nx<<<{"description": "Number of elements in the X direction"}>>> = 40
    xmin<<<{"description": "Lower X Coordinate of the generated mesh"}>>> = 0.1
    xmax<<<{"description": "Upper X Coordinate of the generated mesh"}>>> = 200
    bias_x<<<{"description": "The amount by which to grow (or shrink) the cells in the x-direction."}>>> = 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
  []
  [temperature]
    initial_condition<<<{"description": "Specifies a constant initial condition for this variable"}>>> = 70
    scaling<<<{"description": "Specifies a scaling factor to apply to this variable"}>>> = 1e-4
  []
[]

[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
  []
  [energy]
    type = PorousFlowEnergyTimeDerivative<<<{"description": "Derivative of heat-energy-density wrt time", "href": "../../source/kernels/PorousFlowEnergyTimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
  [heatadv]
    type = PorousFlowHeatAdvection<<<{"description": "Fully-upwinded heat flux, advected by the fluid", "href": "../../source/kernels/PorousFlowHeatAdvection.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
  [conduction]
    type = PorousFlowHeatConduction<<<{"description": "Heat conduction in the Porous Flow module", "href": "../../source/kernels/PorousFlowHeatConduction.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[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 temperature'
    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"}>>> = methane
    capillary_pressure<<<{"description": "Name of the UserObject defining the capillary pressure"}>>> = pc
  []
[]

[FluidProperties<<<{"href": "../../syntax/FluidProperties/index.html"}>>>]
  [methane]
    type = MethaneFluidProperties<<<{"description": "Fluid properties for methane (CH4)", "href": "../../source/fluidproperties/MethaneFluidProperties.html"}>>>
  []
  [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"}>>> = temperature
  []
  [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<<<{"description": "The fluid temperature (C or K, depending on temperature_unit)"}>>> = temperature
    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
  []
  [rockheat]
    type = PorousFlowMatrixInternalEnergy<<<{"description": "This Material calculates the internal energy of solid rock grains, which is specific_heat_capacity * density * temperature.  Kernels multiply this by (1 - porosity) to find the energy density of the porous rock in a rock-fluid system", "href": "../../source/materials/PorousFlowMatrixInternalEnergy.html"}>>>
    specific_heat_capacity<<<{"description": "Specific heat capacity of the rock grains (J/kg/K)."}>>> = 1000
    density<<<{"description": "Density of the rock grains"}>>> = 2500
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal<<<{"description": "This Material calculates rock-fluid combined thermal conductivity by using a weighted sum.  Thermal conductivity = dry_thermal_conductivity + S^exponent * (wet_thermal_conductivity - dry_thermal_conductivity), where S is the aqueous saturation", "href": "../../source/materials/PorousFlowThermalConductivityIdeal.html"}>>>
    dry_thermal_conductivity<<<{"description": "The thermal conductivity of the rock matrix when the aqueous saturation is zero"}>>> = '50 0 0  0 50 0  0 0 50'
  []
[]

[BCs<<<{"href": "../../syntax/BCs/index.html"}>>>]
  [cold_gas]
    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"}>>> = left
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
    value<<<{"description": "Value of the BC"}>>> = 20
  []
  [gas_injecton]
    type = PorousFlowSink<<<{"description": "Applies a flux sink to a boundary.", "href": "../../source/bcs/PorousFlowSink.html"}>>>
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = zi
    flux_function<<<{"description": "The flux.  The flux is OUT of the medium: hence positive values of this function means this BC will act as a SINK, while negative values indicate this flux will be a SOURCE.  The functional form is useful for spatially or temporally varying sinks.  Without any use_*, this function is measured in kg.m^-2.s^-1 (or J.m^-2.s^-1 for the case with only heat and no fluids)"}>>> = -0.159155
  []
  [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
  []
  [righttemp]
    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"}>>> = 70
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = temperature
  []
[]

[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_iname<<<{"description": "Names of PETSc name/value pairs"}>>> = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value<<<{"description": "Values of PETSc name/value pairs (must correspond with \"petsc_options_iname\""}>>> = 'gmres      asm      lu           NONZERO                   2'
  []
[]

[Executioner<<<{"href": "../../syntax/Executioner/index.html"}>>>]
  type = Transient
  solve_type = NEWTON
  end_time = 1e4
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-5
  # Avoids failing first time step in parallel
  line_search = 'none'
  [TimeStepper<<<{"href": "../../syntax/Executioner/TimeStepper/index.html"}>>>]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.5
  []
[]

[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."}>>> = '2 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."}>>> = '2 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."}>>> = '2 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = zi
  []
  [temperature]
    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."}>>> = '2 0 0'
    variable<<<{"description": "The name of the variable that this postprocessor operates on."}>>> = temperature
  []
  [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."}>>> = '2 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."}>>> = '2 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
  csv<<<{"description": "Output the scalar variable and postprocessors to a *.csv file using the default CSV output."}>>> = true
[]
(modules/porous_flow/test/tests/fluidstate/theis_nonisothermal.i)

References

  1. DM Himmelblau. Partial molar heats and entropies of solution for gases dissolved in water from the freezing to near the critical point. The Journal of Physical Chemistry, 63(11):1803–1808, 1959.[BibTeX]