Water97FluidProperties

Fluid properties for water and steam (H2O) using IAPWS-IF97

The water implementation in Fluid Properties is the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam. This formulation calculates properties of water and steam using pressure and temperature as inputs. The IAPWS-IF97 formulation is split into five different regions in the phase diagram.

All five regions are implemented in the Fluid Properties module. To avoid iteration in region 3 of the IAPWS-IF97 formulation, the backwards equations from IAPWS (2014) are implemented.

Viscosity is calculated using the IAPWS 2008 formulation (IAPWS, 2008). Note that the critical enhancement has not been implemented.

Thermal conductivity is calculated using the IAPS 1985 formulation (IAPS, 1985). Although there is a newer formulation available (IAPWS, 2011), it is significantly more complicated, so has not been implemented yet.

Dissolution of a dilute gas into water is calculated using Henry's law (IAPWS, 2004).

Properties of water

Property	value
Molar mass	0.018015 kg/mol
Critical temperature	647.096 K
Critical pressure	22.064 MPa
Critical density	322.0 kg/m $var element = document.getElementById("moose-equation-e3538f61-a44e-400e-8ff9-93267a188115");katex.render("^3", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$
Triple point temperature	273.16 K
Triple point pressure	611.657 Pa

Range of validity

The Water97FluidProperties UserObject is valid for:

273.15 K $var element = document.getElementById("moose-equation-cfccd163-8dc2-40bd-a81f-3b30368c1b7d");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ T $var element = document.getElementById("moose-equation-3a371025-04d9-4e5a-9935-bf12b7c070d4");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ 1073.15 K for p $var element = document.getElementById("moose-equation-1c42fabf-086c-43f7-b20c-5914b40d5c2d");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ 100 MPa
1073.15 K $var element = document.getElementById("moose-equation-543c602b-6a5c-4cdd-bf13-4b0590c687be");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ T $var element = document.getElementById("moose-equation-ad5635d2-9a4c-42b5-a9f0-40b5d069ed83");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ 2273.15 K for p $var element = document.getElementById("moose-equation-367da35e-a0c3-4c6b-b269-3f52ccf478f1");katex.render("\\le", element, {displayMode:false,throwOnError:false,macros:{"\\eqc":"\\,,","\\eqp":"\\,.","\\pd":"\\frac{\\partial #1}{\\partial #2}","\\pr":"\\left(#1\\right)","\\ddt":"\\frac{d #1}{d t}"}});$ 50 MPa

Input Parameters

execute_onTIMESTEP_ENDThe 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.
Default:TIMESTEP_END
C++ Type:ExecFlagEnum
Unit:(no unit assumed)
Options:XFEM_MARK, FORWARD, ADJOINT, HOMOGENEOUS_FORWARD, ADJOINT_TIMESTEP_BEGIN, ADJOINT_TIMESTEP_END, NONE, INITIAL, LINEAR, NONLINEAR_CONVERGENCE, NONLINEAR, POSTCHECK, TIMESTEP_END, TIMESTEP_BEGIN, MULTIAPP_FIXED_POINT_END, MULTIAPP_FIXED_POINT_BEGIN, FINAL, CUSTOM
Controllable:No
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.
prop_getter_suffixAn optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
C++ Type:MaterialPropertyName
Unit:(no unit assumed)
Controllable:No
Description:An optional suffix parameter that can be appended to any attempt to retrieve/get material properties. The suffix will be prepended with a '_' character.
use_interpolated_stateFalseFor the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:For the old and older state use projected material properties interpolated at the quadrature points. To set up projection use the ProjectedStatefulMaterialStorageAction.

Optional Parameters

T_initial_guess400Temperature initial guess for Newton Method variable set conversion
Default:400
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Temperature initial guess for Newton Method variable set conversion
p_initial_guess200000Pressure initial guess for Newton Method variable set conversion
Default:200000
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Pressure initial guess for Newton Method variable set conversion
tolerance1e-08Tolerance for 2D Newton variable set conversion
Default:1e-08
C++ Type:double
Unit:(no unit assumed)
Controllable:No
Description:Tolerance for 2D Newton variable set conversion

Variable Set Conversions Newton Solve Parameters

allow_duplicate_execution_on_initialFalseIn the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:In the case where this UserObject is depended upon by an initial condition, allow it to be executed twice during the initial setup (once before the IC and again after mesh adaptivity (if applicable).
allow_imperfect_jacobiansFalsetrue to allow unimplemented property derivative terms to be set to zero for the AD API
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:true to allow unimplemented property derivative terms to be set to zero for the AD API
control_tagsAdds user-defined labels for accessing object parameters via control logic.
C++ Type:std::vector<std::string>
Unit:(no unit assumed)
Controllable:No
Description:Adds user-defined labels for accessing object parameters via control logic.
enableTrueSet the enabled status of the MooseObject.
Default:True
C++ Type:bool
Unit:(no unit assumed)
Controllable:Yes
Description:Set the enabled status of the MooseObject.
execution_order_group0Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
Default:0
C++ Type:int
Unit:(no unit assumed)
Controllable:No
Description:Execution order groups are executed in increasing order (e.g., the lowest number is executed first). Note that negative group numbers may be used to execute groups before the default (0) group. Please refer to the user object documentation for ordering of user object execution within a group.
force_postauxFalseForces the UserObject to be executed in POSTAUX
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in POSTAUX
force_preauxFalseForces the UserObject to be executed in PREAUX
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in PREAUX
force_preicFalseForces the UserObject to be executed in PREIC during initial setup
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Forces the UserObject to be executed in PREIC during initial setup
fp_typesingle-phase-fpType of the fluid property object
Default:single-phase-fp
C++ Type:FPType
Unit:(no unit assumed)
Controllable:No
Description:Type of the fluid property object
use_displaced_meshFalseWhether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.
Default:False
C++ Type:bool
Unit:(no unit assumed)
Controllable:No
Description:Whether or not this object should use the displaced mesh for computation. Note that in the case this is true but no displacements are provided in the Mesh block the undisplaced mesh will still be used.

Advanced Parameters

Input Files

(modules/porous_flow/test/tests/ics/fluidpropic_celsius.i)
(modules/porous_flow/test/tests/fluids/brine1_tabulated.i)
(modules/porous_flow/test/tests/fluidstate/theis_tabulated.i)
(modules/porous_flow/test/tests/fluidstate/waterncg.i)
(modules/porous_flow/examples/tutorial/05_tabulated.i)
(modules/porous_flow/examples/solute_tracer_transport/solute_tracer_transport_2D.i)
(modules/porous_flow/test/tests/jacobian/brineco2_gas.i)
(modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)
(modules/porous_flow/test/tests/jacobian/waterncg_liquid.i)
(modules/porous_flow/examples/fluidflower/fluidflower.i)
(modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)
(modules/porous_flow/test/tests/jacobian/waterncg_gas.i)
(modules/porous_flow/test/tests/jacobian/brineco2_liquid.i)
(modules/porous_flow/test/tests/fluidstate/coldwater_injection_radial.i)
(modules/porous_flow/test/tests/fluidstate/theis.i)
(modules/porous_flow/examples/solute_tracer_transport/solute_tracer_transport.i)
(modules/porous_flow/test/tests/sinks/injection_production_eg.i)
(modules/fluid_properties/test/tests/brine/brine_tabulated.i)
(modules/thermal_hydraulics/test/tests/components/inlet_mass_flow_rate_1phase/jac.massflowrate_3eqn_water97.i)
(modules/porous_flow/examples/reservoir_model/regular_grid.i)
(modules/porous_flow/examples/tutorial/11.i)
(modules/porous_flow/examples/reservoir_model/field_model.i)
(modules/porous_flow/test/tests/fluidstate/waterncg_ic.i)
(modules/porous_flow/test/tests/fluidstate/theis_brineco2.i)
(modules/porous_flow/test/tests/sinks/injection_production_eg_outflowBC.i)
(modules/porous_flow/examples/natural_convection/natural_convection.i)
(modules/porous_flow/examples/tutorial/05.i)
(modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i)
(modules/fluid_properties/test/tests/water/water.i)
(modules/porous_flow/test/tests/fluidstate/water_vapor_phasechange.i)
(modules/porous_flow/test/tests/fluidstate/water_vapor.i)
(modules/porous_flow/test/tests/jacobian/waterncg_twophase_nonisothermal.i)
(modules/porous_flow/test/tests/fluidstate/theis_nonisothermal.i)
(modules/porous_flow/test/tests/jacobian/brineco2_twophase.i)
(modules/porous_flow/test/tests/fluidstate/waterncg_nonisothermal.i)
(modules/porous_flow/test/tests/jacobian/waterncg_twophase.i)
(modules/porous_flow/examples/tutorial/11_2D.i)
(modules/porous_flow/test/tests/fluids/h2o.i)
(modules/porous_flow/test/tests/ics/fluidpropic.i)
(modules/porous_flow/examples/ates/ates.i)
(modules/porous_flow/test/tests/fluidstate/coldwater_injection.i)
(modules/porous_flow/test/tests/recover/theis.i)

References

IAPS. Revised Release on the IAPS Formulation 1985 for the Thermal Conductivity of Ordinary Water Substance. Technical Report, IAPWS, 1985.

@techreport{iaps1985,
    author = "IAPS",
    title = "{Revised Release on the IAPS Formulation 1985 for the Thermal Conductivity of Ordinary Water Substance}",
    institution = "IAPWS",
    year = "1985"
}

IAPWS. Guidelines on the Henry's constant and vapour liquid distribution constant for gases in H$_2$O and D$_2$O at high temperatures. Technical Report, IAPWS, 2004.

@techreport{iapws2004,
    author = "IAPWS",
    title = "{Guidelines on the Henry's constant and vapour liquid distribution constant for gases in H$\_2$O and D$\_2$O at high temperatures}",
    institution = "IAPWS",
    year = "2004"
}

IAPWS. Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance. Technical Report, IAPWS, 2008. URL: www.iapws.org/relguide/visc.pdf.

@techreport{iapws2008,
    author = "IAPWS",
    title = "{Release on the IAPWS Formulation 2008 for the Viscosity of Ordinary Water Substance}",
    institution = "IAPWS",
    url = "www.iapws.org/relguide/visc.pdf",
    year = "2008"
}

IAPWS. Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance. Technical Report, IAPWS, 2011. URL: www.iapws.org/relguide/ThCond.pdf.

@techreport{iapws2011,
    author = "IAPWS",
    title = "{Release on the IAPWS Formulation 2011 for the Thermal Conductivity of Ordinary Water Substance}",
    institution = "IAPWS",
    url = "www.iapws.org/relguide/ThCond.pdf",
    year = "2011"
}

IAPWS. Revised Supplementary Release on Backward Equations for Specific Volume as a Function of Pressure and Temperature v(p,T) for Region 3 of the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam. Technical Report, IAPWS, 2014. URL: www.iapws.org/relguide/Supp-VPT3-2016.pdf.

@techreport{iapws1997region3,
    author = "IAPWS",
    title = "{Revised Supplementary Release on Backward Equations for Specific Volume as a Function of Pressure and Temperature v(p,T) for Region 3 of the IAPWS Industrial Formulation 1997 for the Thermodynamic Properties of Water and Steam}",
    institution = "IAPWS",
    url = "www.iapws.org/relguide/Supp-VPT3-2016.pdf",
    type = "",
    year = "2014"
}

(modules/porous_flow/test/tests/ics/fluidpropic_celsius.i)

# Test the correct calculation of fluid properties using PorousFlwoFluidPropertyIC
# when temperature is given in Celsius
#
# Variables:
# Pressure: 1 MPa
# Temperature: 50 C
#
# Fluid properties for water (reference values from NIST webbook)
# Density: 988.43 kg/m^3
# Enthalpy: 210.19 kJ/kg
# Internal energy: 2019.18 kJ/kg

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[Variables]
  [pressure]
    initial_condition = 1e6
  []
  [temperature]
    initial_condition = 50
  []
[]

[AuxVariables]
  [enthalpy]
  []
  [internal_energy]
  []
  [density]
  []
[]

[ICs]
  [enthalpy]
    type = PorousFlowFluidPropertyIC
    variable = enthalpy
    property = enthalpy
    porepressure = pressure
    temperature = temperature
    temperature_unit = Celsius
    fp = water
  []
  [internal_energy]
    type = PorousFlowFluidPropertyIC
    variable = internal_energy
    property = internal_energy
    porepressure = pressure
    temperature = temperature
    temperature_unit = Celsius
    fp = water
  []
  [density]
    type = PorousFlowFluidPropertyIC
    variable = density
    property = density
    porepressure = pressure
    temperature = temperature
    temperature_unit = Celsius
    fp = water
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Kernels]
  [pressure]
    type = Diffusion
    variable = pressure
  []
  [temperature]
    type = Diffusion
    variable = temperature
  []
[]

[Executioner]
  type = Steady
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [enthalpy]
    type = ElementAverageValue
    variable = enthalpy
    execute_on = 'initial timestep_end'
  []
  [internal_energy]
    type = ElementAverageValue
    variable = internal_energy
    execute_on = 'initial timestep_end'
  []
  [density]
    type = ElementAverageValue
    variable = density
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  file_base = fluidpropic_out
  execute_on = initial
[]

(modules/porous_flow/test/tests/fluids/brine1_tabulated.i)

# Test the density and viscosity calculated by the brine material using a
# TabulatedFluidProperties userobject for water
# Pressure 20 MPa
# Temperature 50C
# xnacl = 0.1047 (equivalent to 2.0 molality)

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[Variables]
  [pp]
    initial_condition = 20e6
  []
[]

[Kernels]
  [dummy]
    type = Diffusion
    variable = pp
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 50
  []
  [xnacl]
    initial_condition = 0.1047
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [brine]
    type = PorousFlowBrine
    water_fp = watertab
    temperature_unit = Celsius
    xnacl = 0.1047
    phase = 0
  []
[]

[Executioner]
  type = Steady
  solve_type = Newton
[]

[Postprocessors]
  [pressure]
    type = ElementIntegralVariablePostprocessor
    variable = pp
  []
  [temperature]
    type = ElementIntegralVariablePostprocessor
    variable = temp
  []
  [xnacl]
    type = ElementIntegralVariablePostprocessor
    variable = xnacl
  []
  [density]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_fluid_phase_density_qp0'
  []
  [viscosity]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_viscosity_qp0'
  []
  [enthalpy]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_fluid_phase_enthalpy_qp0'
  []
  [internal_energy]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_fluid_phase_internal_energy_qp0'
  []
[]

[Outputs]
  execute_on = 'timestep_end'
  file_base = brine1
  csv = true
[]

(modules/porous_flow/test/tests/fluidstate/theis_tabulated.i)

# 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]
  type = GeneratedMesh
  dim = 1
  nx = 80
  xmax = 200
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [zi]
    initial_condition = 0
  []
[]

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

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = tabulated
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [tabulated]
    type = TabulatedBicubicFluidProperties
    fp = co2
    fluid_property_file = fluid_properties.csv
    # We try to avoid using both, but some properties are not implemented in the tabulation
    allow_fp_and_tabulation = true
    # Test was design prior to bounds check
    error_on_out_of_bounds = false
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = fluid_properties.csv
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 20
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
[]

[BCs]
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 20e6
    variable = pgas
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 2
    variable = zi
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-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 = 'gmres      asm      lu           NONZERO                   2               1E-8       1E-10 20'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 8e2
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 2
    growth_factor = 2
  []
[]

[VectorPostprocessors]
  [line]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    sort_by = x
    start_point = '0 0 0'
    end_point = '200 0 0'
    num_points = 1000
    variable = 'pgas zi x1 saturation_gas'
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    point = '1 0 0'
    variable = pgas
  []
  [sgas]
    type = PointValue
    point = '1 0 0'
    variable = saturation_gas
  []
  [zi]
    type = PointValue
    point = '1 0 0'
    variable = zi
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
  [x1]
    type = PointValue
    point = '1 0 0'
    variable = x1
  []
  [y0]
    type = PointValue
    point = '1 0 0'
    variable = y0
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  [csvout]
    type = CSV
    file_base = theis_tabulated_csvout
    execute_on = timestep_end
    execute_vector_postprocessors_on = final
  []
[]

(modules/porous_flow/test/tests/fluidstate/waterncg.i)

# Tests correct calculation of properties in PorousFlowWaterNCG.
# This test is run three times, with the initial condition of z (the total mass
# fraction of NCG in all phases) varied to give either a single phase liquid, a
# single phase gas, or two phases.

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pgas]
    initial_condition = 1e6
  []
  [z]
     initial_condition = 0.005
  []
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [internal_energy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [internal_energy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_gas]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = timestep_end
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = timestep_end
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = timestep_end
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = timestep_end
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = timestep_end
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = timestep_end
  []
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = timestep_end
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = timestep_end
  []
  [internal_energy_water]
    type = PorousFlowPropertyAux
    variable = internal_energy_water
    property = internal_energy
    phase = 0
    execute_on = timestep_end
  []
  [internal_energy_gas]
    type = PorousFlowPropertyAux
    variable = internal_energy_gas
    property = internal_energy
    phase = 1
    execute_on = timestep_end
  []
  [x1_water]
    type = PorousFlowPropertyAux
    variable = x1_water
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [x1_gas]
    type = PorousFlowPropertyAux
    variable = x1_gas
    property = mass_fraction
    phase = 1
    fluid_component = 1
    execute_on = timestep_end
  []
  [x0_water]
    type = PorousFlowPropertyAux
    variable = x0_water
    property = mass_fraction
    phase = 0
    fluid_component = 0
    execute_on = timestep_end
  []
  [x0_gas]
    type = PorousFlowPropertyAux
    variable = x0_gas
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [density_water]
    type = ElementIntegralVariablePostprocessor
    variable = density_water
  []
  [density_gas]
    type = ElementIntegralVariablePostprocessor
    variable = density_gas
  []
  [viscosity_water]
    type = ElementIntegralVariablePostprocessor
    variable = viscosity_water
  []
  [viscosity_gas]
    type = ElementIntegralVariablePostprocessor
    variable = viscosity_gas
  []
  [enthalpy_water]
    type = ElementIntegralVariablePostprocessor
    variable = enthalpy_water
  []
  [enthalpy_gas]
    type = ElementIntegralVariablePostprocessor
    variable = enthalpy_gas
  []
  [internal_energy_water]
    type = ElementIntegralVariablePostprocessor
    variable = internal_energy_water
  []
  [internal_energy_gas]
    type = ElementIntegralVariablePostprocessor
    variable = internal_energy_gas
  []
  [x1_water]
    type = ElementIntegralVariablePostprocessor
    variable = x1_water
  []
  [x0_water]
    type = ElementIntegralVariablePostprocessor
    variable = x0_water
  []
  [x1_gas]
    type = ElementIntegralVariablePostprocessor
    variable = x1_gas
  []
  [x0_gas]
    type = ElementIntegralVariablePostprocessor
    variable = x0_gas
  []
  [sg]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_gas
  []
  [sw]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_water
  []
  [pwater]
    type = ElementIntegralVariablePostprocessor
    variable = pressure_water
  []
  [pgas]
    type = ElementIntegralVariablePostprocessor
    variable = pressure_gas
  []
  [x0mass]
    type = PorousFlowFluidMass
    fluid_component = 0
    phase = '0 1'
  []
  [x1mass]
    type = PorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
  []
[]

[Outputs]
  exodus = true
  file_base = waterncg_liquid
[]

(modules/porous_flow/examples/tutorial/05_tabulated.i)

# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = tabulated_water
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedFluidProperties
    fp = true_water
    temperature_min = 275
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated.csv
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []

  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/examples/solute_tracer_transport/solute_tracer_transport_2D.i)

# Longitudinal dispersivity
disp = 5

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 100
    xmin = -50
    xmax = 50
    ny = 60
    ymin = 0
    ymax = 50
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [porepressure]
    initial_condition = 1e5
  []
  [C]
    initial_condition = 0
  []
[]

[AuxVariables]
  [Darcy_vel_x]
    order = CONSTANT
    family = MONOMIAL
  []
  [Darcy_vel_y]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [Darcy_vel_x]
    type = PorousFlowDarcyVelocityComponent
    variable = Darcy_vel_x
    component = x
    fluid_phase = 0
  []
  [Darcy_vel_y]
    type = PorousFlowDarcyVelocityComponent
    variable = Darcy_vel_y
    component = y
    fluid_phase = 0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure C'
    number_fluid_phases = 1
    number_fluid_components = 2
  []
[]

[Kernels]
  [mass_der_water]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = porepressure
  []
  [adv_pp]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
    fluid_component = 1
  []
  [diff_pp]
    type = PorousFlowDispersiveFlux
    fluid_component = 1
    variable = porepressure
    disp_trans = 0
    disp_long = ${disp}
  []
  [mass_der_C]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = C
  []
  [adv_C]
    type = PorousFlowFullySaturatedDarcyFlow
    fluid_component = 0
    variable = C
  []
  [diff_C]
    type = PorousFlowDispersiveFlux
    fluid_component = 0
    variable = C
    disp_trans = 0
    disp_long = ${disp}
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-11 0 0   0 1E-11 0   0 0 1E-11'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = C
  []
  [temperature]
    type = PorousFlowTemperature
    temperature = 293
  []
  [diff]
    type = PorousFlowDiffusivityConst
    diffusion_coeff = '0 0'
    tortuosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityConst
    kr = 1
    phase = 0
  []
[]

[DiracKernels]
  [source_P]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 1e-1
    variable = porepressure
  []
  [source_C]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 1e-7
    variable = C
  []
[]

[BCs]
  [constant_outlet_porepressure_]
    type = DirichletBC
    variable = porepressure
    value = 1e5
    boundary = 'top left right'
  []
  [outlet_tracer_top]
    type = PorousFlowOutflowBC
    variable = C
    boundary = top
    mass_fraction_component = 0
  []
  [outlet_tracer_right]
    type = PorousFlowOutflowBC
    variable = C
    boundary = right
    mass_fraction_component = 0
  []
  [outlet_tracer_left]
    type = PorousFlowOutflowBC
    variable = C
    boundary = left
    mass_fraction_component = 0
  []
[]

[Preconditioning]
  [basic]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  end_time = 17280000
  dtmax = 100000
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1000
  []
[]

[Postprocessors]
  [C]
    type = PointValue
    variable = C
    point = '0 25 0'
  []
  [Darcy_x]
    type = PointValue
    variable = Darcy_vel_x
    point = '0 25 0'
  []
  [Darcy_y]
    type = PointValue
    variable = Darcy_vel_y
    point = '0 25 0'
  []
[]

[Outputs]
  file_base = solute_tracer_transport_2D_${disp}
  csv = true
  exodus = true
[]

(modules/porous_flow/test/tests/jacobian/brineco2_gas.i)

# Tests correct calculation of properties derivatives in PorousFlowFluidState
# for conditions that give a single gas phase

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [xnacl]
    initial_condition = 0.05
  []
[]

[Variables]
  [pgas]
  []
  [zi]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 5e4
    max = 1e5
    variable = pgas
  []
  [z]
    type = RandomIC
    min = 0.9
    max = 0.99
    variable = zi
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = zi
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = zi
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
    pc_max = 1e3
  []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [brine]
    type = BrineFluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    xnacl = xnacl
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/examples/multiapp_fracture_flow/3dFracture/fracture_only_aperture_changing.i)

# Cold water injection into one side of the fracture network, and production from the other side
injection_rate = 10 # kg/s

[Mesh]
  uniform_refine = 0
  [cluster34]
    type = FileMeshGenerator
    file = 'Cluster_34.exo'
  []
  [injection_node]
    type = BoundingBoxNodeSetGenerator
    input = cluster34
    bottom_left = '-1000 0 -1000'
    top_right = '1000 0.504 1000'
    new_boundary = injection_node
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 -9.81E-6' # Note the value, because of pressure_unit
[]

[Variables]
  [frac_P]
    scaling = 1E6
  []
  [frac_T]
    initial_condition = 473
  []
[]

[ICs]
  [frac_P]
    type = FunctionIC
    variable = frac_P
    function = insitu_pp
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = frac_P
  temperature = frac_T
  fp = water
  pressure_unit = MPa
[]

[Kernels]
  [toMatrix]
    type = PorousFlowHeatMassTransfer
    variable = frac_T
    v = transferred_matrix_T
    transfer_coefficient = heat_transfer_coefficient
    save_in = joules_per_s
  []
[]

[AuxVariables]
  [heat_transfer_coefficient]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.0
  []
  [transferred_matrix_T]
    initial_condition = 473
  []
  [joules_per_s]
  []
  [normal_dirn_x]
    family = MONOMIAL
    order = CONSTANT
  []
  [normal_dirn_y]
    family = MONOMIAL
    order = CONSTANT
  []
  [normal_dirn_z]
    family = MONOMIAL
    order = CONSTANT
  []
  [enclosing_element_normal_length]
    family = MONOMIAL
    order = CONSTANT
  []
  [enclosing_element_normal_thermal_cond]
    family = MONOMIAL
    order = CONSTANT
  []
  [aperture]
    family = MONOMIAL
    order = CONSTANT
  []
  [perm_times_app]
    family = MONOMIAL
    order = CONSTANT
  []
  [density]
    family = MONOMIAL
    order = CONSTANT
  []
  [viscosity]
    family = MONOMIAL
    order = CONSTANT
  []
  [insitu_pp]
  []
[]

[AuxKernels]
  [normal_dirn_x_auxk]
    type = PorousFlowElementNormal
    variable = normal_dirn_x
    component = x
  []
  [normal_dirn_y]
    type = PorousFlowElementNormal
    variable = normal_dirn_y
    component = y
  []
  [normal_dirn_z]
    type = PorousFlowElementNormal
    variable = normal_dirn_z
    component = z
  []
  [heat_transfer_coefficient_auxk]
    type = ParsedAux
    variable = heat_transfer_coefficient
    coupled_variables = 'enclosing_element_normal_length enclosing_element_normal_thermal_cond'
    constant_names = h_s
    constant_expressions = 1E3 # should be much bigger than thermal_conductivity / L ~ 1
    expression = 'if(enclosing_element_normal_length = 0, 0, h_s * enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length / (h_s * enclosing_element_normal_length * enclosing_element_normal_length + enclosing_element_normal_thermal_cond * 2 * enclosing_element_normal_length))'
  []
  [aperture]
    type = PorousFlowPropertyAux
    variable = aperture
    property = porosity
  []
  [perm_times_app]
    type = PorousFlowPropertyAux
    variable = perm_times_app
    property = permeability
    row = 0
    column = 0
  []
  [density]
    type = PorousFlowPropertyAux
    variable = density
    property = density
    phase = 0
  []
  [viscosity]
    type = PorousFlowPropertyAux
    variable = viscosity
    property = viscosity
    phase = 0
  []
  [insitu_pp]
    type = FunctionAux
    execute_on = initial
    variable = insitu_pp
    function = insitu_pp
  []
[]

[BCs]
  [inject_heat]
    type = DirichletBC
    boundary = injection_node
    variable = frac_T
    value = 373
  []
[]

[DiracKernels]
  [inject_fluid]
    type = PorousFlowPointSourceFromPostprocessor
    mass_flux = ${injection_rate}
    point = '58.8124 0.50384 74.7838'
    variable = frac_P
  []
  [withdraw_fluid]
    type = PorousFlowPeacemanBorehole
    SumQuantityUO = kg_out_uo
    bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
    character = 1
    line_length = 1
    point_file = production.xyz
    unit_weight = '0 0 0'
    fluid_phase = 0
    use_mobility = true
    variable = frac_P
  []
  [withdraw_heat]
    type = PorousFlowPeacemanBorehole
    SumQuantityUO = J_out_uo
    bottom_p_or_t = 10.6 # 1MPa + approx insitu at production point, to prevent aperture closing due to low porepressures
    character = 1
    line_length = 1
    point_file = production.xyz
    unit_weight = '0 0 0'
    fluid_phase = 0
    use_mobility = true
    use_enthalpy = true
    variable = frac_T
  []
[]

[UserObjects]
  [kg_out_uo]
    type = PorousFlowSumQuantity
  []
  [J_out_uo]
    type = PorousFlowSumQuantity
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275 # K
    temperature_max = 600
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated.csv
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityLinear
    porosity_ref = 1E-4 # fracture porosity = 1.0, but must include fracture aperture of 1E-4 at P = insitu_pp
    P_ref = insitu_pp
    P_coeff = 1E-3 # this is in metres/MPa, ie for P_ref = 1/P_coeff, the aperture becomes 1 metre
    porosity_min = 1E-5
  []
  [permeability]
    type = PorousFlowPermeabilityKozenyCarman
    k0 = 1E-15 # fracture perm = 1E-11 m^2, but must include fracture aperture of 1E-4
    poroperm_function = kozeny_carman_phi0
    m = 0
    n = 3
    phi0 = 1E-4
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2700 # kg/m^3
    specific_heat_capacity = 0 # basically no rock inside the fracture
  []
  [aq_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '0.6E-4 0 0  0 0.6E-4 0  0 0 0.6E-4' # thermal conductivity of water times fracture aperture.  This should increase linearly with aperture, but is set constant in this model
  []
[]

[Functions]
  [kg_rate]
    type = ParsedFunction
    symbol_values = 'dt kg_out'
    symbol_names = 'dt kg_out'
    expression = 'kg_out/dt'
  []
  [insitu_pp]
    type = ParsedFunction
    expression = '10 - 0.847E-2 * z' # Approximate hydrostatic in MPa
  []
[]

[Postprocessors]
  [dt]
    type = TimestepSize
    outputs = 'none'
  []
  [kg_out]
    type = PorousFlowPlotQuantity
    uo = kg_out_uo
  []
  [kg_per_s]
    type = FunctionValuePostprocessor
    function = kg_rate
  []
  [J_out]
    type = PorousFlowPlotQuantity
    uo = J_out_uo
  []
  [TK_out]
    type = PointValue
    variable = frac_T
    point = '101.705 160.459 39.5722'
  []
  [P_out]
    type = PointValue
    variable = frac_P
    point = '101.705 160.459 39.5722'
  []
  [P_in]
    type = PointValue
    variable = frac_P
    point = '58.8124 0.50384 74.7838'
  []
[]

[VectorPostprocessors]
  [heat_transfer_rate]
    type = NodalValueSampler
    outputs = none
    sort_by = id
    variable = joules_per_s
  []
[]

[Preconditioning]
  [entire_jacobian]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2             '
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    optimal_iterations = 10
    growth_factor = 1.5
  []
  dtmax = 1E8
  end_time = 1E8
  nl_abs_tol = 1E-3
  nl_max_its = 20
[]

[Outputs]
  print_linear_residuals = false
  csv = true
  [ex]
    type = Exodus
    sync_times = '1 10 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 3100 3200 3300 3400 3500 3600 3700 3800 3900 4000 4100 4200 4300 4400 4500 4600 4700 4800 4900 5000 5100 5200 5300 5400 5500 5600 5700 5800 5900 6000 6100 6200 6300 6400 6500 6600 6700 6800 6900 7000 7100 7200 7300 7400 7500 7600 7700 7800 7900 8000 8100 8200 8300 8400 8500 8600 8700 8800 8900 9000 10000 11000 12000 13000 14000 15000 16000 17000 18000 19000 20000 30000 50000 70000 100000 200000 300000 400000 500000 600000 700000 800000 900000 1000000 1100000 1200000 1300000 1400000 1500000 1600000 1700000 1800000 1900000 2000000 2100000 2200000 2300000 2400000 2500000 2600000 2700000 2800000 2900000'
    sync_only = true
  []
[]

(modules/porous_flow/test/tests/jacobian/waterncg_liquid.i)

# Tests correct calculation of properties derivatives in PorousFlowWaterNCG
# for conditions that give a single liquid phase

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [pgas]
  []
  [z]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 6e6
    max = 8e6
    variable = pgas
  []
  [z]
    type = RandomIC
    min = 0.01
    max = 0.05
    variable = z
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
    pc_max = 1e4
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/examples/fluidflower/fluidflower.i)

# FluidFlower International Benchmark study model
# CSIRO 2023
#
# This example can be used to reproduce the results presented by the
# CSIRO team as part of this benchmark study. See
# Green, C., Jackson, S.J., Gunning, J., Wilkins, A. and Ennis-King, J.,
# 2023. Modelling the FluidFlower: Insights from Characterisation and
# Numerical Predictions. Transport in Porous Media.
#
# This example takes a long time to run! The large density contrast
# between the gas phase CO2 and the water makes convergence very hard,
# so small timesteps must be taken during injection.
#
# This example uses a simplified mesh in order to be run during the
# automated testing. To reproduce the results of the benchmark study,
# replace the simple layered input mesh with the one located in the
# large_media submodule.
#
# The mesh file contains:
# - porosity as given by FluidFlower description
# - permeability as given by FluidFlower description
# - subdomain ids for each sand type
#
# The nominal thickness of the FluidFlower tank is 19mm. To keep masses consistent
# with the experiment, porosity and permeability are multiplied by the thickness
thickness = 0.019
#
# Properties associated with each sand type associated with mesh block ids
#
# block 0 - ESF (very fine sand)
sandESF = '0 10 20'
sandESF_pe = 1471.5
sandESF_krg = 0.09
sandESF_swi = 0.32
sandESF_krw = 0.71
sandESF_sgi = 0.14

# block 1 - C - Coarse lower
sandC = '1 21'
sandC_pe = 294.3
sandC_krg = 0.05
sandC_swi = 0.14
sandC_krw = 0.93
sandC_sgi = 0.1

# block 2 - D - Coarse upper
sandD = '2 22'
sandD_pe = 98.1
sandD_krg = 0.02
sandD_swi = 0.12
sandD_krw = 0.95
sandD_sgi = 0.08

# block 3 - E - Very Coarse lower
sandE = '3 13 23'
sandE_pe = 10
sandE_krg = 0.1
sandE_swi = 0.12
sandE_krw = 0.93
sandE_sgi = 0.06

# block 4 - F - Very Coarse upper
sandF = '4 14 24 34'
sandF_pe = 10
sandF_krg = 0.11
sandF_swi = 0.12
sandF_krw = 0.72
sandF_sgi = 0.13

# block 5 - G - Flush Zone
sandG = '5 15 35'
sandG_pe = 10
sandG_krg = 0.16
sandG_swi = 0.1
sandG_krw = 0.75
sandG_sgi = 0.06

# block 6 - Fault 1 - Heterogeneous
fault1 = '6 26'
fault1_pe = 10
fault1_krg = 0.16
fault1_swi = 0.1
fault1_krw = 0.75
fault1_sgi = 0.06

# block 7 - Fault 2 - Impermeable
# Note: this fault has been removed from the mesh (no elements in this region)

# block 8 - Fault 3 - Homogeneous
fault3 = '8'
fault3_pe = 10
fault3_krg = 0.16
fault3_swi = 0.1
fault3_krw = 0.75
fault3_sgi = 0.06

# Top layer
top_layer = '9'

# Boxes A, B an C used to report values (sg, sgr, xco2, etc)
boxA = '10 13 14 15 34 35'
boxB = '20 21 22 23 24 26'
boxC = '34 35'

# Furthermore, the seal sand unit in boxes A and B
seal_boxA = '10'
seal_boxB = '20'

# CO2 injection details:
# CO2 density ~1.8389 kg/m3 at 293.15 K, 1.01325e5 Pa
# Injection in Port (9, 3) for 5 hours.
# Injection in Port (17, 7) for 2:45 hours.
# Injection of 10 ml/min = 0.1666 ml/s = 1.666e-7 m3/s = ~3.06e-7 kg/s.
# Total mass of CO2 injected ~ 8.5g.
inj_rate = 3.06e-7

[Mesh]
  [mesh]
    type = FileMeshGenerator
    file = 'fluidflower_test.e'
    # file = '../../../../large_media/porous_flow/examples/fluidflower/fluidflower.e'
    use_for_exodus_restart = true
  []
[]

[Debug]
  show_var_residual_norms = true
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature = temperature
  log_extension = false
[]

[Variables]
  [pgas]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [z]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    scaling = 1e4
  []
[]

[AuxVariables]
  [xnacl]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 0.0055
  []
  [temperature]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_condition = 20
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_from_file_var = porosity
  []
  [porosity_times_thickness]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [permeability]
    family = MONOMIAL
    order = CONSTANT
    fv = true
    initial_from_file_var = permeability
  []
  [permeability_times_thickness]
    family = MONOMIAL
    order = CONSTANT
    fv = true
  []
  [saturation_water]
    family = MONOMIAL
    order = CONSTANT
  []
  [saturation_gas]
    family = MONOMIAL
    order = CONSTANT
  []
  [pressure_water]
    family = MONOMIAL
    order = CONSTANT
  []
  [pc]
    family = MONOMIAL
    order = CONSTANT
  []
  [x0_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [porosity_times_thickness]
    type = ParsedAux
    variable = porosity_times_thickness
    coupled_variables = porosity
    expression = 'porosity * ${thickness}'
    execute_on = 'initial'
  []
  [permeability_times_thickness]
    type = ParsedAux
    variable = permeability_times_thickness
    coupled_variables = permeability
    expression = 'permeability * ${thickness}'
    execute_on = 'initial'
  []
  [pressure_water]
    type = ADPorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = ADPorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = ADPorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = ADPorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ADPorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [x1_water]
    type = ADPorousFlowPropertyAux
    variable = x1_water
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = 'initial timestep_end'
  []
  [x1_gas]
    type = ADPorousFlowPropertyAux
    variable = x1_gas
    property = mass_fraction
    phase = 1
    fluid_component = 1
    execute_on = 'initial timestep_end'
  []
  [x0_water]
    type = ADPorousFlowPropertyAux
    variable = x0_water
    property = mass_fraction
    phase = 0
    fluid_component = 0
    execute_on = 'initial timestep_end'
  []
  [x0_gas]
    type = ADPorousFlowPropertyAux
    variable = x0_gas
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = 'initial timestep_end'
  []
  [pc]
    type = ADPorousFlowPropertyAux
    variable = pc
    property = capillary_pressure
    execute_on = 'initial timestep_end'
  []
[]

[FVKernels]
  [mass0]
    type = FVPorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [flux0]
    type = FVPorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [diff0]
    type = FVPorousFlowDispersiveFlux
    variable = pgas
    fluid_component = 0
    disp_long = '0 0'
    disp_trans = '0 0'
  []
  [mass1]
    type = FVPorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [flux1]
    type = FVPorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
  [diff1]
    type = FVPorousFlowDispersiveFlux
    variable = z
    fluid_component = 1
    disp_long = '0 0'
    disp_trans = '0 0'
  []
[]

[DiracKernels]
  [injector1]
    type = ConstantPointSource
    point = '0.9 0.3 0'
    value = ${inj_rate}
    variable = z
  []
  [injector2]
    type = ConstantPointSource
    point = '1.7 0.7 0'
    value = ${inj_rate}
    variable = z
  []
[]

[Controls]
  [injection1]
    type = ConditionalFunctionEnableControl
    enable_objects = 'DiracKernels::injector1'
    conditional_function = injection_schedule1
  []
  [injection2]
    type = ConditionalFunctionEnableControl
    enable_objects = 'DiracKernels::injector2'
    conditional_function = injection_schedule2
  []
[]

[Functions]
  [initial_p]
    type = ParsedFunction
    symbol_names = 'p0 g H rho0'
    symbol_values = '101.325e3 9.81 1.5 1002'
    expression = 'p0 + rho0 * g * (H - y)'
  []
  [injection_schedule1]
    type = ParsedFunction
    expression = 'if(t >= 0 & t <= 1.8e4, 1, 0)'
  []
  [injection_schedule2]
    type = ParsedFunction
    expression = 'if(t >= 8.1e3 & t <= 1.8e4, 1, 0)'
  []
[]

[ICs]
  [p]
    type = FunctionIC
    variable = pgas
    function = initial_p
  []
[]

[FVBCs]
  [pressure_top]
    type = FVPorousFlowAdvectiveFluxBC
    boundary = top
    porepressure_value = 1.01325e5
    variable = pgas
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
    pressure_min = 1e5
    pressure_max = 1e6
    temperature_min = 290
    temperature_max = 300
    num_p = 20
    num_T = 10
  []
  [co2]
    type = CO2FluidProperties
  []
  [co2tab]
    type = TabulatedBicubicFluidProperties
    fp = co2
    save_file = false
    pressure_min = 1e5
    pressure_max = 1e6
    temperature_min = 290
    temperature_max = 300
    num_p = 20
    num_T = 10
  []
  [brine]
    type = BrineFluidProperties
    water_fp = watertab
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [sandESF_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandESF_pe}
    lambda = 2
    block = ${sandESF}
    pc_max = 1e4
    sat_lr = ${sandESF_swi}
  []
  [sandC_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandC_pe}
    lambda = 2
    block = ${sandC}
    pc_max = 1e4
    sat_lr = ${sandC_swi}
  []
  [sandD_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandD_pe}
    lambda = 2
    block = ${sandD}
    pc_max = 1e4
    sat_lr = ${sandD_swi}
  []
  [sandE_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandE_pe}
    lambda = 2
    block = ${sandE}
    pc_max = 1e4
    sat_lr = ${sandE_swi}
  []
  [sandF_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandF_pe}
    lambda = 2
    block = ${sandF}
    pc_max = 1e4
    sat_lr = ${sandF_swi}
  []
  [sandG_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${sandG_pe}
    lambda = 2
    block = ${sandG}
    pc_max = 1e4
    sat_lr = ${sandG_swi}
  []
  [fault1_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${fault1_pe}
    lambda = 2
    block = ${fault1}
    pc_max = 1e4
    sat_lr = ${fault1_swi}
  []
  [fault3_pc]
    type = PorousFlowCapillaryPressureBC
    pe = ${fault3_pe}
    lambda = 2
    block = ${fault3}
    pc_max = 1e4
    sat_lr = ${fault3_swi}
  []
  [top_layer_pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
    block =  ${top_layer}
  []
  [sandESF_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandESF_pc
  []
  [sandC_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandC_pc
  []
  [sandD_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandD_pc
  []
  [sandE_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandE_pc
  []
  [sandF_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandF_pc
  []
  [sandG_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = sandG_pc
  []
  [fault1_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = fault1_pc
  []
  [fault3_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = fault3_pc
  []
  [top_layer_fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2tab
    capillary_pressure = top_layer_pc
  []
[]

[Materials]
  [temperature]
    type = ADPorousFlowTemperature
    temperature = temperature
  []
  [sandESF_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandESF_fs
    capillary_pressure = sandESF_pc
    block = ${sandESF}
  []
  [sandC_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandC_fs
    capillary_pressure = sandC_pc
    block = ${sandC}
  []
  [sandD_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandD_fs
    capillary_pressure = sandD_pc
    block = ${sandD}
  []
  [sandE_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandE_fs
    capillary_pressure = sandE_pc
    block = ${sandE}
  []
  [sandF_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandF_fs
    capillary_pressure = sandF_pc
    block = ${sandF}
  []
  [sandG_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = sandG_fs
    capillary_pressure = sandG_pc
    block = ${sandG}
  []
  [fault1_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = fault1_fs
    capillary_pressure = fault1_pc
    block = ${fault1}
  []
  [fault3_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = fault3_fs
    capillary_pressure = fault3_pc
    block = ${fault3}
  []
  [top_layer_brineco2]
    type = ADPorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    xnacl = xnacl
    fluid_state = top_layer_fs
    capillary_pressure = top_layer_pc
    block = ${top_layer}
  []
  [porosity]
    type = ADPorousFlowPorosityConst
    porosity = porosity_times_thickness
  []
  [permeability]
    type = ADPorousFlowPermeabilityConstFromVar
    perm_xx = permeability_times_thickness
    perm_yy = permeability_times_thickness
    perm_zz = permeability_times_thickness
  []
  [diffcoeff]
    type = ADPorousFlowDiffusivityConst
    tortuosity = '1 1'
    diffusion_coeff = '2e-9 2e-9 0 0'
  []
  [sandESF_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandESF_swi}
    sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
    scaling = ${sandESF_krw}
    block = ${sandESF}
  []
  [sandESF_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandESF_sgi}
    sum_s_res = ${fparse sandESF_sgi + sandESF_swi}
    scaling = ${sandESF_krg}
    block = ${sandESF}
  []
  [sandC_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandC_swi}
    sum_s_res = ${fparse sandC_sgi + sandC_swi}
    scaling = ${sandC_krw}
    block = ${sandC}
  []
  [sandC_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandC_sgi}
    sum_s_res = ${fparse sandC_sgi + sandC_swi}
    scaling = ${sandC_krg}
    block = ${sandC}
  []
  [sandD_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandD_swi}
    sum_s_res = ${fparse sandD_sgi + sandD_swi}
    scaling = ${sandD_krw}
    block = ${sandD}
  []
  [sandD_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandD_sgi}
    sum_s_res = ${fparse sandD_sgi + sandD_swi}
    scaling = ${sandD_krg}
    block = ${sandD}
  []
  [sandE_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandE_swi}
    sum_s_res = ${fparse sandE_sgi + sandE_swi}
    scaling = ${sandE_krw}
    block = ${sandE}
  []
  [sandE_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandE_sgi}
    sum_s_res = ${fparse sandE_sgi + sandE_swi}
    scaling = ${sandE_krg}
    block = ${sandE}
  []
  [sandF_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandF_swi}
    sum_s_res = ${fparse sandF_sgi + sandF_swi}
    scaling = ${sandF_krw}
    block = ${sandF}
  []
  [sandF_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandF_sgi}
    sum_s_res = ${fparse sandF_sgi + sandF_swi}
    scaling = ${sandF_krg}
    block = ${sandF}
  []
  [sandG_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${sandG_swi}
    sum_s_res = ${fparse sandG_sgi + sandG_swi}
    scaling = ${sandG_krw}
    block = ${sandG}
  []
  [sandG_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${sandG_sgi}
    sum_s_res = ${fparse sandG_sgi + sandG_swi}
    scaling = ${sandG_krg}
    block = ${sandG}
  []
  [fault1_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${fault1_swi}
    sum_s_res = ${fparse fault1_sgi + fault1_swi}
    scaling = ${fault1_krw}
    block = ${fault1}
  []
  [fault1_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${fault1_sgi}
    sum_s_res = ${fparse fault1_sgi + fault1_swi}
    scaling = ${fault1_krg}
    block = ${fault1}
  []
  [fault3_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    s_res = ${fault3_swi}
    sum_s_res = ${fparse fault3_sgi + fault3_swi}
    scaling = ${fault3_krw}
    block = ${fault3}
  []
  [fault3_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    s_res = ${fault3_sgi}
    sum_s_res = ${fparse fault3_sgi + fault3_swi}
    scaling = ${fault3_krg}
    block = ${fault3}
  []
  [top_layer_relperm0]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 0
    lambda = 2
    block = ${top_layer}
  []
  [top_layer_relperm1]
    type = ADPorousFlowRelativePermeabilityBC
    phase = 1
    nw_phase = true
    lambda = 2
    block = ${top_layer}
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options = '-ksp_snes_ew'
    petsc_options_iname = '-ksp_type -pc_type -pc_factor_mat_solver_package -sub_pc_factor_shift_type'
    petsc_options_value = 'gmres lu mumps NONZERO'
    # petsc_options_iname = '-ksp_type -pc_type -pc_hypre_type -sub_pc_type -sub_pc_factor_shift_type -sub_pc_factor_levels -ksp_gmres_restart'
    # petsc_options_value = 'gmres hypre boomeramg lu NONZERO 4 301'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dtmax = 60
  start_time = 0
  end_time = 4.32e5
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-8
  nl_max_its = 15
  l_tol = 1e-5
  l_abs_tol = 1e-8
  # line_search = none # Can be a useful option for this problem
  [TimeSteppers]
    [time]
      type = FunctionDT
      growth_factor = 2
      cutback_factor_at_failure = 0.5
      function = 'if(t<1.8e4, 2, if(t<3.6e4, 20, 60))'
    []
  []
[]

[Postprocessors]
  [p_5_3]
    type = PointValue
    variable = pgas
    point = '0.5 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_3_w]
    type = PointValue
    variable = pressure_water
    point = '0.5 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_7]
    type = PointValue
    variable = pgas
    point = '0.5 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_5_7_w]
    type = PointValue
    variable = pressure_water
    point = '0.5 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_9_3]
    type = PointValue
    variable = pgas
    point = '0.9 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_9_3_w]
    type = PointValue
    variable = pressure_water
    point = '0.9 0.3 0'
    execute_on = 'initial timestep_end'
  []
  [p_15_5]
    type = PointValue
    variable = pgas
    point = '1.5 0.5 0'
    execute_on = 'initial timestep_end'
  []
  [p_15_5_w]
    type = PointValue
    variable = pressure_water
    point = '1.5 0.5 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_7]
    type = PointValue
    variable = pgas
    point = '1.7 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_7_w]
    type = PointValue
    variable = pressure_water
    point = '1.7 0.7 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_11]
    type = PointValue
    variable = pgas
    point = '1.7 1.1 0'
    execute_on = 'initial timestep_end'
  []
  [p_17_11_w]
    type = PointValue
    variable = pressure_water
    point = '1.7 1.1 0'
    execute_on = 'initial timestep_end'
  []
  [x0mass]
    type = FVPorousFlowFluidMass
    fluid_component = 0
    phase = '0 1'
    execute_on = 'initial timestep_end'
  []
  [x1mass]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    execute_on = 'initial timestep_end'
  []
  [x1gas]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '1'
    execute_on = 'initial timestep_end'
  []
  [boxA]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${boxA}
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandESF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandESF_sgi}
    block = 10
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandE]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandE_sgi}
    block = 13
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandF_sgi}
    block = '14 34'
    execute_on = 'initial timestep_end'
  []
  [imm_A_sandG]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandG_sgi}
    block = '15 35'
    execute_on = 'initial timestep_end'
  []
  [imm_A]
    type = LinearCombinationPostprocessor
    pp_names = 'imm_A_sandESF imm_A_sandE imm_A_sandF imm_A_sandG'
    pp_coefs = '1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [diss_A]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 0
    block = ${boxA}
    execute_on = 'initial timestep_end'
  []
  [seal_A]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${seal_boxA}
    execute_on = 'initial timestep_end'
  []
  [boxB]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${boxB}
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandESF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandESF_sgi}
    block = 20
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandC]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandC_sgi}
    block = 21
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandD]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandD_sgi}
    block = 22
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandE]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandE_sgi}
    block = 23
    execute_on = 'initial timestep_end'
  []
  [imm_B_sandF]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${sandF_sgi}
    block = 24
    execute_on = 'initial timestep_end'
  []
  [imm_B_fault1]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 1
    saturation_threshold = ${fault1_sgi}
    block = 26
    execute_on = 'initial timestep_end'
  []
  [imm_B]
    type = LinearCombinationPostprocessor
    pp_names = 'imm_B_sandESF imm_B_sandC imm_B_sandD imm_B_sandE imm_B_sandF imm_B_fault1'
    pp_coefs = '1 1 1 1 1 1'
    execute_on = 'initial timestep_end'
  []
  [diss_B]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = 0
    block = ${boxB}
    execute_on = 'initial timestep_end'
  []
  [seal_B]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
    block = ${seal_boxB}
    execute_on = 'initial timestep_end'
  []
  [boxC]
    type = FVPorousFlowFluidMass
    fluid_component = 1
    phase = '0'
    block = ${boxC}
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  # exodus = true
  [csv]
     type = CSV
  []
[]

(modules/porous_flow/examples/co2_intercomparison/1Dradial/1Dradial.i)

# Intercomparison problem 3: Radial flow from an injection well
#
# From Pruess et al, Code intercomparison builds confidence in
# numerical simulation models for geologic disposal of CO2, Energy 29 (2004)
#
# A variation with zero salinity can be run by changing the initial condition
# of the AuxVariable xnacl

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 500
  xmax = 10000
  bias_x = 1.01
  coord_type = 'RZ'
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = 'dictator'
  gravity = '0 0 0'
[]

[AuxVariables]
  [pressure_liquid]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
  [xnacl]
    initial_condition = 0.15
  []
[]

[AuxKernels]
  [pressure_liquid]
    type = PorousFlowPropertyAux
    variable = pressure_liquid
    property = pressure
    phase = 0
    execute_on = 'timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'timestep_end'
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = 'timestep_end'
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = 'timestep_end'
  []
[]

[Variables]
  [pgas]
    initial_condition = 12e6
  []
  [zi]
    initial_condition = 0
    scaling = 1e4
  []
[]

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

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 5.099e-5
    m = 0.457
    sat_lr = 0.0
    pc_max = 1e7
  []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2sw]
    type = CO2FluidProperties
  []
  [co2]
    type = TabulatedBicubicFluidProperties
    fp = co2sw
  []
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    temperature_min = 273.15
    temperature_max = 573.15
    fluid_property_output_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
    water_fp = watertab
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = '45'
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = 'pgas'
    z = 'zi'
    temperature_unit = Celsius
    xnacl = 'xnacl'
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = '0.12'
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityVG
    m = 0.457
    phase = 0
    s_res = 0.3
    sum_s_res = 0.35
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    s_res = 0.05
    sum_s_res = 0.35
  []
[]

[BCs]
  [rightwater]
    type = PorousFlowPiecewiseLinearSink
    boundary = 'right'
    variable = pgas
    use_mobility = true
    PorousFlowDictator = dictator
    fluid_phase = 0
    multipliers = '0 1e9'
    PT_shift = '12e6'
    pt_vals = '0 1e9'
    mass_fraction_component = 0
    use_relperm = true
  []
  [rightco2]
    type = PorousFlowPiecewiseLinearSink
    variable = zi
    boundary = 'right'
    use_mobility = true
    PorousFlowDictator = dictator
    fluid_phase = 1
    multipliers = '0 1e9'
    PT_shift = '12e6'
    pt_vals = '0 1e9'
    mass_fraction_component = 1
    use_relperm = true
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 1
    variable = zi
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type'
    petsc_options_value = 'gmres bjacobi lu NONZERO'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 8.64e8
  nl_max_its = 25
  l_max_its = 100
  dtmax = 5e6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
  []
[]

[VectorPostprocessors]
  [vars]
    type = NodalValueSampler
    sort_by = x
    variable = 'pgas zi xnacl'
    execute_on = 'timestep_end'
    outputs = spatial
  []
  [auxvars]
    type = ElementValueSampler
    sort_by = x
    variable = 'saturation_gas x1 y0'
    execute_on = 'timestep_end'
    outputs = spatial
  []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    point = '25.25 0 0'
    variable = pgas
    outputs = time
  []
  [sgas]
    type = PointValue
    point = '25.25 0 0'
    variable = saturation_gas
    outputs = time
  []
  [zi]
    type = PointValue
    point = '25.25 0 0'
    variable = zi
    outputs = time
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
    outputs = time
  []
  [x1]
    type = PointValue
    point = '25.25 0 0'
    variable = x1
    outputs = time
  []
  [y0]
    type = PointValue
    point = '25.25 0 0'
    variable = y0
    outputs = time
  []
  [xnacl]
    type = PointValue
    point = '25.25 0 0'
    variable = xnacl
    outputs = time
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  sync_times = '2.592e6 8.64e6 8.64e7 8.64e8'
  [time]
    type = CSV
  []
  [spatial]
    type = CSV
    sync_only = true
  []
[]

(modules/porous_flow/test/tests/jacobian/waterncg_gas.i)

# Tests correct calculation of properties derivatives in PorousFlowWaterNCG
# for conditions that give a single gas phase

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [pgas]
  []
  [z]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 1e4
    max = 4e4
    variable = pgas
  []
  [z]
    type = RandomIC
    min = 0.88
    max = 0.98
    variable = z
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
    pc_max = 1e3
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/test/tests/jacobian/brineco2_liquid.i)

# Tests correct calculation of properties derivatives in PorousFlowFluidState
# for conditions that give a single liquid phase

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [xnacl]
    initial_condition = 0.05
  []
[]

[Variables]
  [pgas]
  []
  [zi]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 5e6
    max = 8e6
    variable = pgas
  []
  [z_liquid]
    type = RandomIC
    min = 0.01
    max = 0.03
    variable = zi
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = zi
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = zi
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1
    pc_max = 1e4
  []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [brine]
    type = BrineFluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    xnacl = xnacl
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/test/tests/fluidstate/coldwater_injection_radial.i)

# Cold water injection into 1D radial hot reservoir (Avdonin, 1964)
#
# To generate results presented in documentation for this problem,
# set xmax = 1000 and nx = 200 in the Mesh block, and dtmax = 1e4
# and end_time = 1e6 in the Executioner block.

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 50
  xmin = 0.1
  xmax = 5
  bias_x = 1.05
  rz_coord_axis = Y
  coord_type = RZ
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Variables]
  [pliquid]
    initial_condition = 5e6
  []
  [h]
    scaling = 1e-6
  []
[]

[ICs]
  [hic]
    type = PorousFlowFluidPropertyIC
    variable = h
    porepressure = pliquid
    property = enthalpy
    temperature = 170
    temperature_unit = Celsius
    fp = water
  []
[]

[Functions]
  [injection_rate]
    type = ParsedFunction
    symbol_values = injection_area
    symbol_names = area
    expression = '-0.1/area'
  []
[]

[BCs]
  [source]
    type = PorousFlowSink
    variable = pliquid
    flux_function = injection_rate
    boundary = left
  []
  [pright]
    type = DirichletBC
    variable = pliquid
    value = 5e6
    boundary = right
  []
  [hleft]
    type = DirichletBC
    variable = h
    value = 678.52e3
    boundary = left
  []
  [hright]
    type = DirichletBC
    variable = h
    value = 721.4e3
    boundary = right
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliquid
  []
  [massflux]
    type = PorousFlowAdvectiveFlux
    variable = pliquid
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
  [heatflux]
    type = PorousFlowHeatAdvection
    variable = h
  []
  [heatcond]
    type = PorousFlowHeatConduction
    variable = h
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliquid h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    pc_max = 1e6
    sat_lr = 0.1
    m = 0.5
    alpha = 1e-5
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliquid
    enthalpy = h
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.8e-11 0 0 0 1.8e-11 0 0 0 1.8e-11'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    sum_s_res = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2900
    specific_heat_capacity = 740
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '20 0 0  0 20 0  0 0 20'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e3
  nl_abs_tol = 1e-8
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
  []
[]

[Postprocessors]
  [injection_area]
    type = AreaPostprocessor
    boundary = left
    execute_on = initial
  []
[]

[VectorPostprocessors]
  [line]
    type = ElementValueSampler
    sort_by = x
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  perf_graph = true
  [csv]
    type = CSV
    execute_on = final
  []
[]

(modules/porous_flow/test/tests/fluidstate/theis.i)

# 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]
  type = GeneratedMesh
  dim = 1
  nx = 40
  xmax = 200
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [zi]
    initial_condition = 0
  []
[]

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

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 20
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
[]

[BCs]
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 20e6
    variable = pgas
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 2
    variable = zi
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-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 = 'gmres      asm      lu           NONZERO                   2               1E-8       1E-10 20'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 2e2
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 10
    growth_factor = 2
  []
[]

[VectorPostprocessors]
  [line]
    type = NodalValueSampler
    sort_by = x
    variable = 'pgas zi'
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    point = '1 0 0'
    variable = pgas
  []
  [sgas]
    type = PointValue
    point = '1 0 0'
    variable = saturation_gas
  []
  [zi]
    type = PointValue
    point = '1 0 0'
    variable = zi
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
  [x1]
    type = PointValue
    point = '1 0 0'
    variable = x1
  []
  [y0]
    type = PointValue
    point = '1 0 0'
    variable = y0
  []
[]

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

(modules/porous_flow/examples/solute_tracer_transport/solute_tracer_transport.i)

# Longitudinal dispersivity
disp = 0.7

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 100
    xmin = 0
    xmax = 100
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [porepressure]
    initial_condition = 1e5
  []
  [C]
    initial_condition = 0
  []
[]

[AuxVariables]
  [Darcy_vel_x]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [Darcy_vel_x]
    type = PorousFlowDarcyVelocityComponent
    variable = Darcy_vel_x
    component = x
    fluid_phase = 0
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'porepressure C'
    number_fluid_phases = 1
    number_fluid_components = 2
  []
[]

[Kernels]
  [mass_der_water]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = porepressure
  []
  [adv_pp]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
    fluid_component = 1
  []
  [diff_pp]
    type = PorousFlowDispersiveFlux
    fluid_component = 1
    variable = porepressure
    disp_trans = 0
    disp_long = ${disp}
  []
  [mass_der_C]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = C
  []
  [adv_C]
    type = PorousFlowFullySaturatedDarcyFlow
    fluid_component = 0
    variable = C
  []
  [diff_C]
    type = PorousFlowDispersiveFlux
    fluid_component = 0
    variable = C
    disp_trans = 0
    disp_long = ${disp}
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.25
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1E-11 0 0   0 1E-11 0   0 0 1E-11'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = water
    phase = 0
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = C
  []
  [temperature]
    type = PorousFlowTemperature
    temperature = 293
  []
  [diff]
    type = PorousFlowDiffusivityConst
    diffusion_coeff = '0 0'
    tortuosity = 0.1
  []
  [relperm]
    type = PorousFlowRelativePermeabilityConst
    kr = 1
    phase = 0
  []
[]

[BCs]
  [constant_inlet_pressure]
    type = DirichletBC
    variable = porepressure
    value = 1.2e5
    boundary = left
  []
  [constant_outlet_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 1e5
    boundary = right
  []
  [inlet_tracer]
    type = DirichletBC
    variable = C
    value = 0.001
    boundary = left
  []
  [outlet_tracer]
    type = PorousFlowOutflowBC
    variable = C
    boundary = right
    mass_fraction_component = 0
  []
[]

[Preconditioning]
  [basic]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  end_time = 17280000
  dtmax = 86400
  nl_rel_tol = 1e-6
  nl_abs_tol = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1000
  []
[]

[Postprocessors]
  [C]
    type = PointValue
    variable = C
    point = '50 0 0'
  []
  [Darcy_x]
    type = PointValue
    variable = Darcy_vel_x
    point = '50 0 0'
  []
[]

[Outputs]
  file_base = solute_tracer_transport_${disp}
  csv = true
[]

(modules/porous_flow/test/tests/sinks/injection_production_eg.i)

# phase = 0 is liquid phase
# phase = 1 is gas phase
# fluid_component = 0 is water
# fluid_component = 1 is CO2

# Constant rate of CO2 injection into the left boundary
# 1D mesh
# The PorousFlowPiecewiseLinearSinks remove the correct water and CO2 from the right boundary

# Note i take pretty big timesteps here so the system is quite nonlinear

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 20
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [frac_water_in_liquid]
    initial_condition = 1.0
  []
  [frac_water_in_gas]
    initial_condition = 0.0
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
[]

[Variables]
  [pwater]
    initial_condition = 20E6
  []
  [pgas]
    initial_condition = 20.1E6
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = pgas
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas pwater'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedBicubicFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 293.15
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.2
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
[]

[BCs]
  [co2_injection]
    type = PorousFlowSink
    boundary = left
    variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
    flux_function = -1E-2 # negative means a source, rather than a sink
  []

  [right_water]
    type = PorousFlowPiecewiseLinearSink
    boundary = right

    # a sink of water, since the Kernels given to pwater are for fluid_component = 0 (the water)
    variable = pwater

    # this Sink is a function of liquid porepressure
    # Also, all the mass_fraction, mobility and relperm are referenced to the liquid phase now
    fluid_phase = 0

    # Sink strength = (Pwater - 20E6)
    pt_vals = '0 1E9'
    multipliers = '0 1E9'
    PT_shift = 20E6

    # multiply Sink strength computed above by mass fraction of water at the boundary
    mass_fraction_component = 0

    # also multiply Sink strength by mobility of the liquid
    use_mobility = true

    # also multiply Sink strength by the relperm of the liquid
    use_relperm = true

    # also multiplly Sink strength by 1/L, where L is the distance to the fixed-porepressure external environment
    flux_function = 10 # 1/L
  []

  [right_co2]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E9'
    PT_shift = 20.1E6
    mass_fraction_component = 1
    use_mobility = true
    use_relperm = true
    flux_function = 10 # 1/L
  []
[]

[Preconditioning]
  active = 'basic'
  [basic]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt -snes_linesearch_monitor'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu NONZERO 2'
  []
  [preferred]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  nl_abs_tol = 1E-13
  nl_rel_tol = 1E-10
  end_time = 1e4
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E4
    growth_factor = 1.1
  []
[]

[VectorPostprocessors]
  [pps]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    start_point = '0 0 0'
    end_point = '20 0 0'
    num_points = 20
    sort_by = x
    variable = 'pgas pwater saturation_gas'
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  [out]
    type = CSV
    execute_on = final
  []
[]

(modules/fluid_properties/test/tests/brine/brine_tabulated.i)

# Test BrineFluidProperties calculations of density, viscosity and thermal
# conductivity with a TabulatedBiCubicFluidProperties water.
#
# Experimental density values from Pitzer et al, "Thermodynamic properties
# of aqueous sodium chloride solution", Journal of Physical and Chemical
# Reference Data, 13, 1-102 (1984)
#
# Experimental viscosity values from Phillips et al, "Viscosity of NaCl and
# other solutions up to 350C and 50MPa pressures", LBL-11586 (1980)
#
# Thermal conductivity values from Ozbek and Phillips, "Thermal conductivity of
# aqueous NaCl solutions from 20C to 330C", LBL-9086 (1980)
#
#  --------------------------------------------------------------
#  Pressure (Mpa)                |   20      |    20     |   40
#  Temperature (C)               |   50      |   200     |  200
#  NaCl molality (mol/kg)        |    2      |     2     |    5
#  NaCl mass fraction (kg/kg)    |  0.1047   |  0.1047   |  0.2261
#  --------------------------------------------------------------
#  Expected values
#  --------------------------------------------------------------
#  Density (kg/m^3)              |  1068.52  |  959.27   |  1065.58
#  Viscosity (1e-6Pa.s)          |  679.8    |  180.0    |  263.1
#  Thermal conductivity (W/m/K)  |  0.630    |  0.649    |  0.633
#  --------------------------------------------------------------
#  Calculated values
#  --------------------------------------------------------------
#  Density (kg/m^3)              |  1067.18  |  958.68   |  1065.46
#  Viscosity (1e-6 Pa.s)         |  681.1    |  181.98    |  266.1
#  Thermal conductivity (W/m/K)  |  0.637    |   0.662    |  0.658
#  --------------------------------------------------------------
#
# All results are within expected accuracy

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 3
  ny = 1
  xmax = 3
  # This test uses ElementalVariableValue postprocessors on specific
  # elements, so element numbering needs to stay unchanged
  allow_renumbering = false
[]

[Variables]
  [./dummy]
  [../]
[]

[AuxVariables]
  [./pressure]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./temperature]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./xnacl]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./density]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./enthalpy]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./internal_energy]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[Functions]
  [./pic]
    type = ParsedFunction
    expression = 'if(x<2,20e6, 40e6)'
  [../]
  [./tic]
    type = ParsedFunction
    expression = 'if(x<1, 323.15, 473.15)'
  [../]
  [./xic]
    type = ParsedFunction
    expression = 'if(x<2,0.1047, 0.2261)'
  [../]
[]

[ICs]
  [./p_ic]
    type = FunctionIC
    function = pic
    variable = pressure
  [../]
  [./t_ic]
    type = FunctionIC
    function = tic
    variable = temperature
  [../]
  [./x_ic]
    type = FunctionIC
    function = xic
    variable = xnacl
  [../]
[]

[AuxKernels]
  [./density]
    type = MaterialRealAux
     variable = density
     property = density
  [../]
  [./enthalpy]
    type = MaterialRealAux
     variable = enthalpy
     property = enthalpy
  [../]
  [./internal_energy]
    type = MaterialRealAux
     variable = internal_energy
     property = e
  [../]
[]

[FluidProperties]
  [./water]
    type = Water97FluidProperties
  [../]
  [./water_tab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
  [../]
  [./brine]
    type = BrineFluidProperties
    water_fp = water_tab
  [../]
[]

[Materials]
  [./fp_mat]
    type = MultiComponentFluidPropertiesMaterialPT
    pressure = pressure
    temperature = temperature
    xmass = xnacl
    fp = brine
  [../]
[]

[Kernels]
  [./diff]
    type = Diffusion
    variable = dummy
  [../]
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
[]

[Postprocessors]
  [./density0]
    type = ElementalVariableValue
    variable = density
    elementid = 0
  [../]
  [./density1]
    type = ElementalVariableValue
    variable = density
    elementid = 1
  [../]
  [./density2]
    type = ElementalVariableValue
    variable = density
    elementid = 2
  [../]
  [./enthalpy0]
    type = ElementalVariableValue
    variable = enthalpy
    elementid = 0
  [../]
  [./enthalpy1]
    type = ElementalVariableValue
    variable = enthalpy
    elementid = 1
  [../]
  [./enthalpy2]
    type = ElementalVariableValue
    variable = enthalpy
    elementid = 2
  [../]
  [./e0]
    type = ElementalVariableValue
    variable = internal_energy
    elementid = 0
  [../]
  [./e1]
    type = ElementalVariableValue
    variable = internal_energy
    elementid = 1
  [../]
  [./e2]
    type = ElementalVariableValue
    variable = internal_energy
    elementid = 2
  [../]
[]

[Outputs]
  csv = true
  file_base = brine_out
[]

(modules/thermal_hydraulics/test/tests/components/inlet_mass_flow_rate_1phase/jac.massflowrate_3eqn_water97.i)

[GlobalParams]
  gravity_vector = '0 0 0'

  initial_T = 444.447
  initial_p = 7e6
  initial_vel = 0

  closures = simple_closures
[]

[FluidProperties]
  [fp]
    type = Water97FluidProperties
    T_initial_guess = 444.447
    p_initial_guess = 7e6
  []
[]

[Closures]
  [simple_closures]
    type = Closures1PhaseSimple
  []
[]

[Components]
  [pipe]
    type = FlowChannel1Phase
    fp = fp
    # geometry
    position = '0 0 0'
    orientation = '1 0 0'
    A   = 1.0000000000e-04
    D_h  = 1.1283791671e-02
    f = 0.1
    length = 1
    n_elems = 4
  []

  [inlet]
    type = InletMassFlowRateTemperature1Phase
    input = 'pipe:in'
    m_dot = 0.18
    T     = 444.447
  []

  [outlet]
    type = Outlet1Phase
    input = 'pipe:out'
    p = 7e6
  []
[]

[Preconditioning]
  [SMP_PJFNK]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  scheme = 'bdf2'

  dt = 0.1
  start_time = 0.0
  num_steps = 30

  solve_type = 'NEWTON'
  line_search = 'basic'
  nl_rel_tol = 0
  nl_abs_tol = 1e-6
  nl_max_its = 20

  l_tol = 1e-3
  l_max_its = 100

  abort_on_solve_fail = true

  [Quadrature]
    type = GAUSS
    order = SECOND
  []
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/examples/reservoir_model/regular_grid.i)

# SPE 10 comparative problem - model 1
# Data and description from https://www.spe.org/web/csp/datasets/set01.htm
# Simple input file that just establishes gravity equilibrium in the model
#
# Heterogeneous permeability is included by reading data from an external file
# using the PiecewiseMultilinear function, and saving that data to an elemental
# AuxVariable that is then used in PorousFlowPermeabilityConstFromVar

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 100
  ny = 20
  xmax = 762
  ymax = 15.24
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 -9.81 0'
  temperature_unit = Celsius
[]

[Variables]
  [porepressure]
    initial_condition = 20e6
  []
[]

[Functions]
  [perm_md_fcn]
    type = PiecewiseMultilinear
    data_file = spe10_case1.data
  []
[]

[BCs]
  [top]
    type = DirichletBC
    variable = porepressure
    value = 20e6
    boundary = top
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = 50
  []
  [xnacl]
    initial_condition = 0.1
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
    initial_condition = 0.2
  []
  [perm_md]
    family = MONOMIAL
    order = CONSTANT
  []
  [perm]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = porepressure
  []
  [flux0]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
  []
[]

[AuxKernels]
  [perm_md]
    type = FunctionAux
    function = perm_md_fcn
    variable = perm_md
    execute_on = initial
  []
  [perm]
    type = ParsedAux
    variable = perm
    coupled_variables = perm_md
    expression = '9.869233e-16*perm_md'
    execute_on = initial
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = porepressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temperature
  []
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [brine]
    type = PorousFlowBrine
    compute_enthalpy = false
    compute_internal_energy = false
    xnacl = xnacl
    phase = 0
    water_fp = watertab
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = porosity
  []
  [permeability]
    type = PorousFlowPermeabilityConstFromVar
    perm_xx = perm
    perm_yy = perm
    perm_zz = perm
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1e5
  nl_abs_tol = 1e-12
  nl_rel_tol = 1e-06
  steady_state_detection = true
  steady_state_tolerance = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e2
  []
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]

(modules/porous_flow/examples/tutorial/11.i)

# Two-phase borehole injection problem
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    input = annular
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_x disp_y'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[GlobalParams]
  displacements = 'disp_x disp_y disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
[]

[Variables]
  [pwater]
    initial_condition = 20E6
  []
  [pgas]
    initial_condition = 20.1E6
  []
  [T]
    initial_condition = 330
    scaling = 1E-5
  []
  [disp_x]
    scaling = 1E-5
  []
  [disp_y]
    scaling = 1E-5
  []
[]

[Kernels]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    variable = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  []
  [vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    variable = pgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = T
  []
  [advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  []
  [conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  []
  [vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    variable = T
  []
  [grad_stress_x]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_x
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 0
  []
  [poro_x]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_x
    use_displaced_mesh = false
    component = 0
  []
  [grad_stress_y]
    type = StressDivergenceTensors
    temperature = T
    variable = disp_y
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 1
  []
  [poro_y]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_y
    use_displaced_mesh = false
    component = 1
  []
[]

[AuxVariables]
  [disp_z]
  []
  [effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  []
  [mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  []
  [mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  []
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
  [swater]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_tt]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_zz]
    family = MONOMIAL
    order = CONSTANT
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [effective_fluid_pressure]
    type = ParsedAux
    coupled_variables = 'pwater pgas swater sgas'
    expression = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  []
  [swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [stress_rr]
    type = RankTwoScalarAux
    variable = stress_rr
    rank_two_tensor = stress
    scalar_type = RadialStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  []
  [stress_tt]
    type = RankTwoScalarAux
    variable = stress_tt
    rank_two_tensor = stress
    scalar_type = HoopStress
    point1 = '0 0 0'
    point2 = '0 0 1'
    execute_on = timestep_end
  []
  [stress_zz]
    type = RankTwoAux
    variable = stress_zz
    rank_two_tensor = stress
    index_i = 2
    index_j = 2
    execute_on = timestep_end
  []
  [porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  []
[]

[BCs]
  [roller_tmax]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = dmax
  []
  [roller_tmin]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = dmin
  []
  [pinned_top_bottom_x]
    type = DirichletBC
    variable = disp_x
    value = 0
    boundary = 'top bottom'
  []
  [pinned_top_bottom_y]
    type = DirichletBC
    variable = disp_y
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_x]
    type = Pressure
    boundary = injection_area
    variable = disp_x
    component = 0
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []
  [cavity_pressure_y]
    type = Pressure
    boundary = injection_area
    variable = disp_y
    component = 1
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []

  [cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290 # injection temperature
    use_displaced_mesh = false
  []
  [constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  []

  [outer_water_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = rmax
    variable = pwater
    fluid_phase = 0
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
  [outer_co2_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = rmax
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20.1E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.1
    sum_s_res = 0.2
    phase = 0
  []
  [relperm_co2]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
  [porosity_mat]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
  [stress]
    type = ComputeLinearElasticStress
  []

  [effective_fluid_pressure_mat]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
  []
[]

[Functions]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction
    symbol_names = effective_fluid_pressure_at_wellbore
    symbol_values = effective_fluid_pressure_at_wellbore
    expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E3
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/examples/reservoir_model/field_model.i)

# Field model generated using geophysical modelling tool

[Mesh]
  type = FileMesh
  file = field.e
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 -9.81'
  temperature_unit = Celsius
[]

[Problem]
  # Variable porepressure has an initial condition despite the restart
  allow_initial_conditions_with_restart = true
[]

[Variables]
  [porepressure]
    initial_condition = 20e6
  []
[]

[AuxVariables]
  [temperature]
    initial_condition = 50
  []
  [xnacl]
    initial_condition = 0.1
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
    initial_from_file_var = poro
  []
  [permx_md]
    family = MONOMIAL
    order = CONSTANT
    initial_from_file_var = permX
  []
  [permy_md]
    family = MONOMIAL
    order = CONSTANT
    initial_from_file_var = permY
  []
  [permz_md]
    family = MONOMIAL
    order = CONSTANT
    initial_from_file_var = permZ
  []
  [permx]
    family = MONOMIAL
    order = CONSTANT
  []
  [permy]
    family = MONOMIAL
    order = CONSTANT
  []
  [permz]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [permx]
    type = ParsedAux
    variable = permx
    coupled_variables = permx_md
    expression = '9.869233e-16*permx_md'
    execute_on = initial
  []
  [permy]
    type = ParsedAux
    variable = permy
    coupled_variables = permy_md
    expression = '9.869233e-16*permy_md'
    execute_on = initial
  []
  [permz]
    type = ParsedAux
    variable = permz
    coupled_variables = permz_md
    expression = '9.869233e-16*permz_md'
    execute_on = initial
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = porepressure
  []
  [flux0]
    type = PorousFlowFullySaturatedDarcyFlow
    variable = porepressure
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = porepressure
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedBicubicFluidProperties
    fp = water
    save_file = false
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temperature
  []
  [ps]
    type = PorousFlow1PhaseFullySaturated
    porepressure = porepressure
  []
  [massfrac]
    type = PorousFlowMassFraction
  []
  [brine]
    type = PorousFlowBrine
    compute_enthalpy = false
    compute_internal_energy = false
    xnacl = xnacl
    phase = 0
    water_fp = watertab
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = porosity
  []
  [permeability]
    type = PorousFlowPermeabilityConstFromVar
    perm_xx = permx
    perm_yy = permy
    perm_zz = permz
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  dt = 1e2
  end_time = 1e2
[]

[Outputs]
  execute_on = 'initial timestep_end'
  exodus = true
  perf_graph = true
[]

(modules/porous_flow/test/tests/fluidstate/waterncg_ic.i)

# Tests correct calculation of z (total mass fraction of NCG summed over all
# phases) using the PorousFlowFluidStateIC initial condition. Once z is
# calculated by the initial condition, the thermophysical properties are calculated
# and the resulting gas saturation should be equal to that given in the intial condition

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  temperature_unit = Celsius
[]

[Variables]
  [pgas]
    initial_condition = 1e6
  []
  [z]
  []
[]

[ICs]
  [z]
    type = PorousFlowFluidStateIC
    saturation = 0.5
    gas_porepressure = pgas
    temperature = 50
    variable = z
    fluid_state = fs
  []
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    fluid_state = fs
    capillary_pressure = pc
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [sg]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_gas
    execute_on = 'initial timestep_end'
  []
  [sw]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_water
    execute_on = 'initial timestep_end'
  []
  [z]
    type = ElementIntegralVariablePostprocessor
    variable = z
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
[]

(modules/porous_flow/test/tests/fluidstate/theis_brineco2.i)

# 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]
  type = GeneratedMesh
  dim = 1
  nx = 2000
  xmax = 2000
[]

[Problem]
  type = FEProblem
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [zi]
    initial_condition = 0
  []
  [xnacl]
    initial_condition = 0.1
  []
[]

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

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi xnacl'
    number_fluid_phases = 2
    number_fluid_components = 3
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2sw]
    type = CO2FluidProperties
  []
  [co2]
    type = TabulatedFluidProperties
    fp = co2sw
    fluid_property_file = 'fluid_properties.csv'
    allow_fp_and_tabulation = true
    error_on_out_of_bounds = false
  []
  [water]
    type = Water97FluidProperties
  []
  [watertab]
    type = TabulatedFluidProperties
    fp = water
    temperature_min = 273.15
    temperature_max = 573.15
    fluid_property_output_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
    water_fp = watertab
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 20
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    xnacl = xnacl
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
[]

[BCs]
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 20e6
    variable = pgas
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 2
    variable = zi
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e5
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.5
  []
[]

[VectorPostprocessors]
  [line]
    type = LineValueSampler
    warn_discontinuous_face_values = false
    sort_by = x
    start_point = '0 0 0'
    end_point = '2000 0 0'
    num_points = 10000
    variable = 'pgas zi xnacl x1 saturation_gas'
    execute_on = 'timestep_end'
  []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    point = '4 0 0'
    variable = pgas
  []
  [sgas]
    type = PointValue
    point = '4 0 0'
    variable = saturation_gas
  []
  [zi]
    type = PointValue
    point = '4 0 0'
    variable = zi
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
  [x1]
    type = PointValue
    point = '4 0 0'
    variable = x1
  []
  [y0]
    type = PointValue
    point = '4 0 0'
    variable = y0
  []
  [xnacl]
    type = PointValue
    point = '4 0 0'
    variable = xnacl
  []
[]

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

(modules/porous_flow/test/tests/sinks/injection_production_eg_outflowBC.i)

# phase = 0 is liquid phase
# phase = 1 is gas phase
# fluid_component = 0 is water
# fluid_component = 1 is CO2

# Constant rates of water and CO2 injection into the left boundary
# 1D mesh
# The PorousFlowOutflowBCs remove the correct water and CO2 from the right boundary

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 20
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [frac_water_in_liquid]
    initial_condition = 1.0
  []
  [frac_water_in_gas]
    initial_condition = 0.0
  []
  [water_kg_per_s]
  []
  [co2_kg_per_s]
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
[]

[Variables]
  [pwater]
    initial_condition = 20E6
  []
  [pgas]
    initial_condition = 21E6
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pwater
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = pgas
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas pwater'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedBicubicFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 293.15
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'frac_water_in_liquid frac_water_in_gas'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.2
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
[]

[BCs]
  [water_injection]
    type = PorousFlowSink
    boundary = left
    variable = pwater # pwater is associated with the water mass balance (fluid_component = 0 in its Kernels)
    flux_function = -1E-5 # negative means a source, rather than a sink
  []
  [co2_injection]
    type = PorousFlowSink
    boundary = left
    variable = pgas # pgas is associated with the CO2 mass balance (fluid_component = 1 in its Kernels)
    flux_function = -2E-5 # negative means a source, rather than a sink
  []

  [right_water_component0]
    type = PorousFlowOutflowBC
    boundary = right
    variable = pwater
    mass_fraction_component = 0
    save_in = water_kg_per_s
  []
  [right_co2_component1]
    type = PorousFlowOutflowBC
    boundary = right
    variable = pgas
    mass_fraction_component = 1
    save_in = co2_kg_per_s
  []
[]

[Preconditioning]
  active = 'basic'
  [basic]
    type = SMP
    full = true
    petsc_options = '-snes_converged_reason -ksp_diagonal_scale -ksp_diagonal_scale_fix -ksp_gmres_modifiedgramschmidt'
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres asm lu NONZERO 2'
  []
  [preferred]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = 'lu mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  nl_abs_tol = 1E-10
  nl_rel_tol = 1E-10
  end_time = 1E5
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E5
    growth_factor = 1.1
  []
[]

[Postprocessors]
  [water_kg_per_s]
    type = NodalSum
    boundary = right
    variable = water_kg_per_s
  []
  [co2_kg_per_s]
    type = NodalSum
    boundary = right
    variable = co2_kg_per_s
  []
[]

[VectorPostprocessors]
  [pps]
    type = LineValueSampler
    start_point = '0 0 0'
    end_point = '20 0 0'
    num_points = 20
    sort_by = x
    variable = 'pgas pwater saturation_gas'
  []
[]

[Outputs]
  [out]
    type = CSV
    execute_on = final
  []
[]

(modules/porous_flow/examples/natural_convection/natural_convection.i)

# Example problem: Elder, Transient convection in a porous mediu, 1967

[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 64
    ny = 32
    xmin = 0
    xmax = 300
    ymax = 0
    ymin = -150
  []
  [heater]
    type = ParsedGenerateSideset
    input = gen
    combinatorial_geometry = 'x <= 150 & y = -150'
    new_sideset_name = heater
  []
  uniform_refine = 1
[]

[Variables]
  [porepressure]
  []
  [T]
    initial_condition = 285
  []
[]

[AuxVariables]
  [density]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [density]
    type = PorousFlowPropertyAux
    variable = density
    property = density
    execute_on = TIMESTEP_END
  []
[]

[ICs]
  [hydrostatic]
    type = FunctionIC
    variable = porepressure
    function = '1e5 - 9.81 * 1000 * y'
  []
[]

[GlobalParams]
  PorousFlowDictator = 'dictator'
  gravity = '0 -9.81 0'
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = T
  fp = water
[]

[Materials]
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.21E-10 0 0   0 1.21E-10 0   0 0 1.21E-10'
  []
  [Matrix_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500
    specific_heat_capacity = 0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '1.5 0 0  0 1.5 0  0 0 0'
  []
[]

[BCs]
  [t_bot]
    type = DirichletBC
    variable = T
    value = 293
    boundary = 'heater'
  []
  [t_top]
    type = DirichletBC
    variable = T
    value = 285
    boundary = 'top'
  []
  [p_top]
    type = DirichletBC
    variable = porepressure
    value = 1e5
    boundary = top
  []
[]

[Preconditioning]
  [basic]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  end_time = 63072000
  dtmax = 1e6
  nl_rel_tol = 1e-6
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1000
  []
  [Adaptivity]
    interval = 1
    refine_fraction = 0.2
    coarsen_fraction = 0.3
    max_h_level = 4
  []
[]

[Outputs]
  exodus = true
[]

# If you uncomment this it will print out all the kernels and materials that the PorousFlowFullySaturated action generates
#[Problem]
#  type = DumpObjectsProblem
#  dump_path = PorousFlowFullySaturated
#[]

(modules/porous_flow/examples/tutorial/05.i)

# Darcy flow with heat advection and conduction, using Water97 properties
[Mesh]
  [annular]
    type = AnnularMeshGenerator
    nr = 10
    rmin = 1.0
    rmax = 10
    growth_r = 1.4
    nt = 4
    dmin = 0
    dmax = 90
  []
  [make3D]
    type = MeshExtruderGenerator
    extrusion_vector = '0 0 12'
    num_layers = 3
    bottom_sideset = 'bottom'
    top_sideset = 'top'
    input = annular
  []
  [shift_down]
    type = TransformGenerator
    transform = TRANSLATE
    vector_value = '0 0 -6'
    input = make3D
  []
  [aquifer]
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 0 -2'
    top_right = '10 10 2'
    input = shift_down
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x*x+y*y<1.01'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [porepressure]
    initial_condition = 1E6
  []
  [temperature]
    initial_condition = 313
    scaling = 1E-8
  []
[]

[PorousFlowBasicTHM]
  porepressure = porepressure
  temperature = temperature
  coupling_type = ThermoHydro
  gravity = '0 0 0'
  fp = the_simple_fluid
[]

[BCs]
  [constant_injection_porepressure]
    type = DirichletBC
    variable = porepressure
    value = 2E6
    boundary = injection_area
  []
  [constant_injection_temperature]
    type = DirichletBC
    variable = temperature
    value = 333
    boundary = injection_area
  []
[]

[FluidProperties]
  [the_simple_fluid]
    type = Water97FluidProperties
  []
[]

[Materials]
  [porosity]
    type = PorousFlowPorosity
    porosity_zero = 0.1
  []
  [biot_modulus]
    type = PorousFlowConstantBiotModulus
    biot_coefficient = 0.8
    solid_bulk_compliance = 2E-7
    fluid_bulk_modulus = 1E7
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '1E-14 0 0   0 1E-14 0   0 0 1E-14'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '1E-15 0 0   0 1E-15 0   0 0 1E-16'
  []

  [thermal_expansion]
    type = PorousFlowConstantThermalExpansionCoefficient
    biot_coefficient = 0.8
    drained_coefficient = 0.003
    fluid_coefficient = 0.0002
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500.0
    specific_heat_capacity = 1200.0
  []
  [thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '10 0 0  0 10 0  0 0 10'
    block = 'caps aquifer'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E6
  dt = 1E5
  nl_abs_tol = 1E-10
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/test/tests/fluidstate/water_vapor_tab.i)

# Tests correct calculation of properties in PorousFlowWaterVapor in the two-phase region

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pliq]
    initial_condition = 1e6
  []
  [h]
    initial_condition = 8e5
    scaling = 1e-3
  []
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliq h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureBC
    pe = 1e5
    lambda = 2
    pc_max = 1e6
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]

  [water_true]
    type = Water97FluidProperties
  []
  [water]
    type = TabulatedBicubicFluidProperties
    fp = water_true
    allow_fp_and_tabulation = true
    fluid_property_file = fluid_properties_extended.csv
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliq
    enthalpy = h
    temperature_unit = Kelvin
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500
    specific_heat_capacity = 1200
  []
[]

[Problem]
  solve = false
[]

[Executioner]
  type = Steady
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [density_water]
    type = ElementAverageValue
    variable = density_water
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ElementAverageValue
    variable = density_gas
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = ElementAverageValue
    variable = viscosity_water
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = ElementAverageValue
    variable = viscosity_gas
    execute_on = 'initial timestep_end'
  []
  [enthalpy_water]
    type = ElementAverageValue
    variable = enthalpy_water
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = ElementAverageValue
    variable = enthalpy_gas
    execute_on = 'initial timestep_end'
  []
  [sg]
    type = ElementAverageValue
    variable = saturation_gas
    execute_on = 'initial timestep_end'
  []
  [sw]
    type = ElementAverageValue
    variable = saturation_water
    execute_on = 'initial timestep_end'
  []
  [pwater]
    type = ElementAverageValue
    variable = pressure_water
    execute_on = 'initial timestep_end'
  []
  [pgas]
    type = ElementAverageValue
    variable = pressure_gas
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = ElementAverageValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [enthalpy]
    type = ElementAverageValue
    variable = h
    execute_on = 'initial timestep_end'
  []
  [liquid_mass]
    type = PorousFlowFluidMass
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [vapor_mass]
    type = PorousFlowFluidMass
    phase = 1
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = water_vapor_twophase_tab
  csv = true
  execute_on = INITIAL
[]

(modules/fluid_properties/test/tests/water/water.i)

# Example of using Water97FluidProperties module in Region 1 by recovering the values
# in Table 5 of Revised Release on the IAPWS Industrial Formulation 1997 for the
# Thermodynamic Properties of Water and Steam

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 3
  xmax = 3
  # This test uses ElementalVariableValue postprocessors on specific
  # elements, so element numbering needs to stay unchanged
  allow_renumbering = false
[]

[Variables]
  [./dummy]
  [../]
[]

[AuxVariables]
  [./pressure]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./temperature]
    order = CONSTANT
    family = MONOMIAL
  [../]
  [./rho]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./v]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./e]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./h]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./s]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./cp]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./cv]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./c]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./mu]
    family = MONOMIAL
    order = CONSTANT
  [../]
  [./k]
    family = MONOMIAL
    order = CONSTANT
  [../]
[]

[Functions]
  [./tic]
    type = ParsedFunction
    expression = 'if(x<2, 300, 500)'
  [../]
  [./pic]
    type = ParsedFunction
    expression = 'if(x<1,3e6, if(x<2, 80e6, 3e6))'
  [../]
[]

[ICs]
  [./p_ic]
    type = FunctionIC
    function = pic
    variable = pressure
  [../]
  [./t_ic]
    type = FunctionIC
    function = tic
    variable = temperature
  [../]
[]

[AuxKernels]
  [./rho]
    type = MaterialRealAux
    variable = rho
    property = density
  [../]
  [./v]
    type = ParsedAux
    coupled_variables = rho
    expression = 1/rho
    variable = v
  [../]
  [./e]
    type = MaterialRealAux
    variable = e
    property = e
  [../]
  [./h]
    type = MaterialRealAux
    variable = h
    property = h
  [../]
  [./s]
    type = MaterialRealAux
    variable = s
    property = s
  [../]
  [./cp]
    type = MaterialRealAux
    variable = cp
    property = cp
  [../]
  [./cv]
    type = MaterialRealAux
    variable = cv
    property = cv
  [../]
  [./c]
    type = MaterialRealAux
    variable = c
    property = c
  [../]
  [./mu]
    type = MaterialRealAux
    variable = mu
    property = viscosity
  [../]
  [./k]
    type = MaterialRealAux
    variable = k
    property = k
  [../]
[]

[FluidProperties]
  [./water]
    type = Water97FluidProperties
  [../]
[]

[Materials]
  [./fp_mat]
    type = FluidPropertiesMaterialPT
    pressure = pressure
    temperature = temperature
    fp = water
  [../]
[]

[Kernels]
  [./diff]
    type = Diffusion
    variable = dummy
  [../]
[]

[Postprocessors]
  [./density0]
    type = ElementalVariableValue
    variable = rho
    elementid = 0
  [../]
  [./density1]
    type = ElementalVariableValue
    variable = rho
    elementid = 1
  [../]
  [./density2]
    type = ElementalVariableValue
    variable = rho
    elementid = 2
  [../]
  [./v0]
    type = ElementalVariableValue
    variable = v
    elementid = 0
  [../]
  [./v1]
    type = ElementalVariableValue
    variable = v
    elementid = 1
  [../]
  [./v2]
    type = ElementalVariableValue
    variable = v
    elementid = 2
  [../]
  [./e0]
    type = ElementalVariableValue
    variable = e
    elementid = 0
  [../]
  [./e1]
    type = ElementalVariableValue
    variable = e
    elementid = 1
  [../]
  [./e2]
    type = ElementalVariableValue
    variable = e
    elementid = 2
  [../]
  [./h0]
    type = ElementalVariableValue
    variable = h
    elementid = 0
  [../]
  [./h1]
    type = ElementalVariableValue
    variable = h
    elementid = 1
  [../]
  [./h2]
    type = ElementalVariableValue
    variable = h
    elementid = 2
  [../]
  [./s0]
    type = ElementalVariableValue
    variable = s
    elementid = 0
  [../]
  [./s1]
    type = ElementalVariableValue
    variable = s
    elementid = 1
  [../]
  [./s2]
    type = ElementalVariableValue
    variable = s
    elementid = 2
  [../]
  [./cp0]
    type = ElementalVariableValue
    variable = cp
    elementid = 0
  [../]
  [./cp1]
    type = ElementalVariableValue
    variable = cp
    elementid = 1
  [../]
  [./cp2]
    type = ElementalVariableValue
    variable = cp
    elementid = 2
  [../]
  [./cv0]
    type = ElementalVariableValue
    variable = cv
    elementid = 0
  [../]
  [./cv1]
    type = ElementalVariableValue
    variable = cv
    elementid = 1
  [../]
  [./cv2]
    type = ElementalVariableValue
    variable = cv
    elementid = 2
  [../]
  [./c0]
    type = ElementalVariableValue
    variable = c
    elementid = 0
  [../]
  [./c1]
    type = ElementalVariableValue
    variable = c
    elementid = 1
  [../]
  [./c2]
    type = ElementalVariableValue
    variable = c
    elementid = 2
  [../]
  [./mu0]
    type = ElementalVariableValue
    variable = mu
    elementid = 0
  [../]
  [./mu1]
    type = ElementalVariableValue
    variable = mu
    elementid = 1
  [../]
  [./mu2]
    type = ElementalVariableValue
    variable = mu
    elementid = 2
  [../]
  [./k0]
    type = ElementalVariableValue
    variable = k
    elementid = 0
  [../]
  [./k1]
    type = ElementalVariableValue
    variable = k
    elementid = 1
  [../]
  [./k2]
    type = ElementalVariableValue
    variable = k
    elementid = 2
  [../]
[]

[Executioner]
  type = Steady
  solve_type = NEWTON
[]

[Outputs]
  csv = true
[]

(modules/porous_flow/test/tests/fluidstate/water_vapor_phasechange.i)

# Tests correct calculation of properties in PorousFlowWaterVapor as a phase change
# from liquid to a two-phase model occurs due to a pressure drop.
# A single 10 m^3 element is used, with constant mass and heat production using
# a Dirac kernel. Initial conditions correspond to just outside the two-phase region in
# the liquid state.
#
# An identical problem can be run using TOUGH2, with the following outputs after 1,000s
# Pressure: 8.58 Mpa
# Temperature: 299.92 K
# Vapor saturation: 0.00637

[Mesh]
  type = GeneratedMesh
  dim = 3
  xmax = 10
  ymax = 10
  zmax = 10
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pliq]
    initial_condition = 9e6
  []
  [h]
    scaling = 1e-3
  []
[]

[ICs]
  [hic]
    type = PorousFlowFluidPropertyIC
    variable = h
    porepressure = pliq
    property = enthalpy
    temperature = 300
    temperature_unit = Celsius
    fp = water
  []
[]

[DiracKernels]
  [mass]
    type = ConstantPointSource
    point = '5 5 5'
    variable = pliq
    value = -1
  []
  [heat]
    type = ConstantPointSource
    point = '5 5 5'
    variable = h
    value = -1.344269e6
  []
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
  [e_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [e_water]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
  [e_water]
    type = PorousFlowPropertyAux
    variable = e_water
    property = internal_energy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [egas]
    type = PorousFlowPropertyAux
    variable = e_gas
    property = internal_energy
    phase = 1
    execute_on = 'initial timestep_end'
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliq
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliq h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureBC
    pe = 1e5
    lambda = 2
    pc_max = 1e6
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliq
    enthalpy = h
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-14 0 0 0 1e-14 0 0 0 1e-14'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2650
    specific_heat_capacity = 1000
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e3
  nl_abs_tol = 1e-12
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 10
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [density_water]
    type = ElementAverageValue
    variable = density_water
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ElementAverageValue
    variable = density_gas
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = ElementAverageValue
    variable = viscosity_water
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = ElementAverageValue
    variable = viscosity_gas
    execute_on = 'initial timestep_end'
  []
  [enthalpy_water]
    type = ElementAverageValue
    variable = enthalpy_water
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = ElementAverageValue
    variable = enthalpy_gas
    execute_on = 'initial timestep_end'
  []
  [sg]
    type = ElementAverageValue
    variable = saturation_gas
    execute_on = 'initial timestep_end'
  []
  [sw]
    type = ElementAverageValue
    variable = saturation_water
    execute_on = 'initial timestep_end'
  []
  [pwater]
    type = ElementAverageValue
    variable = pressure_water
    execute_on = 'initial timestep_end'
  []
  [pgas]
    type = ElementAverageValue
    variable = pressure_gas
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = ElementAverageValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [enthalpy]
    type = ElementAverageValue
    variable = h
    execute_on = 'initial timestep_end'
  []
  [pliq]
    type = ElementAverageValue
    variable = pliq
    execute_on = 'initial timestep_end'
  []
  [liquid_mass]
    type = PorousFlowFluidMass
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [vapor_mass]
    type = PorousFlowFluidMass
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [liquid_heat]
    type = PorousFlowHeatEnergy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [vapor_heat]
    type = PorousFlowHeatEnergy
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [e_water]
    type = ElementAverageValue
    variable = e_water
    execute_on = 'initial timestep_end'
  []
  [e_gas]
    type = ElementAverageValue
    variable = e_gas
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  perf_graph = false
[]

(modules/porous_flow/test/tests/fluidstate/water_vapor.i)

# Tests correct calculation of properties in PorousFlowWaterVapor in the two-phase region

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pliq]
    initial_condition = 1e6
  []
  [h]
    initial_condition = 8e5
    scaling = 1e-3
  []
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliq
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliq h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureBC
    pe = 1e5
    lambda = 2
    pc_max = 1e6
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliq
    enthalpy = h
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-13 0 0 0 1e-13 0 0 0 1e-13'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2500
    specific_heat_capacity = 1200
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [density_water]
    type = ElementAverageValue
    variable = density_water
    execute_on = 'initial timestep_end'
  []
  [density_gas]
    type = ElementAverageValue
    variable = density_gas
    execute_on = 'initial timestep_end'
  []
  [viscosity_water]
    type = ElementAverageValue
    variable = viscosity_water
    execute_on = 'initial timestep_end'
  []
  [viscosity_gas]
    type = ElementAverageValue
    variable = viscosity_gas
    execute_on = 'initial timestep_end'
  []
  [enthalpy_water]
    type = ElementAverageValue
    variable = enthalpy_water
    execute_on = 'initial timestep_end'
  []
  [enthalpy_gas]
    type = ElementAverageValue
    variable = enthalpy_gas
    execute_on = 'initial timestep_end'
  []
  [sg]
    type = ElementAverageValue
    variable = saturation_gas
    execute_on = 'initial timestep_end'
  []
  [sw]
    type = ElementAverageValue
    variable = saturation_water
    execute_on = 'initial timestep_end'
  []
  [pwater]
    type = ElementAverageValue
    variable = pressure_water
    execute_on = 'initial timestep_end'
  []
  [pgas]
    type = ElementAverageValue
    variable = pressure_gas
    execute_on = 'initial timestep_end'
  []
  [temperature]
    type = ElementAverageValue
    variable = temperature
    execute_on = 'initial timestep_end'
  []
  [enthalpy]
    type = ElementAverageValue
    variable = h
    execute_on = 'initial timestep_end'
  []
  [liquid_mass]
    type = PorousFlowFluidMass
    phase = 0
    execute_on = 'initial timestep_end'
  []
  [vapor_mass]
    type = PorousFlowFluidMass
    phase = 1
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  file_base = water_vapor_twophase
  csv = true
[]

(modules/porous_flow/test/tests/jacobian/waterncg_twophase_nonisothermal.i)

# Tests correct calculation of properties derivatives in PorousFlowWaterNCG
# for nonisothermal two phase conditions

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 2
    ny = 2
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [pgas]
  []
  [z]
  []
  [temperature]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 1e5
    max = 5e5
    variable = pgas
  []
  [z]
    type = RandomIC
    min = 0.01
    max = 0.06
    variable = z
  []
  [temperature]
    type = RandomIC
    min = 20
    max = 80
    variable = temperature
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
  [energy]
    type = PorousFlowEnergyTimeDerivative
    variable = temperature
  []
  [heat]
    type = PorousFlowHeatAdvection
    variable = temperature
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z temperature'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e1
    pc_max = 1e4
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temperature
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature = temperature
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
  [rock_heat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1000
    density = 2500
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/test/tests/fluidstate/theis_nonisothermal.i)

# 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]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 1
    nx = 40
    xmin = 0.1
    xmax = 200
    bias_x = 1.05
  []
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [zi]
    initial_condition = 0
  []
  [temperature]
    initial_condition = 70
    scaling = 1e-4
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pgas
  []
  [flux0]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    variable = pgas
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = zi
  []
  [flux1]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    variable = zi
  []
  [energy]
    type = PorousFlowEnergyTimeDerivative
    variable = temperature
  []
  [heatadv]
    type = PorousFlowHeatAdvection
    variable = temperature
  []
  [conduction]
    type = PorousFlowHeatConduction
    variable = temperature
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi temperature'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = methane
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [methane]
    type = MethaneFluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temperature
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature = temperature
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
  [rockheat]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1000
    density = 2500
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '50 0 0  0 50 0  0 0 50'
  []
[]

[BCs]
  [cold_gas]
    type = DirichletBC
    boundary = left
    variable = temperature
    value = 20
  []
  [gas_injecton]
    type = PorousFlowSink
    boundary = left
    variable = zi
    flux_function = -0.159155
  []
  [rightwater]
    type = DirichletBC
    boundary = right
    value = 20e6
    variable = pgas
  []
  [righttemp]
    type = DirichletBC
    boundary = right
    value = 70
    variable = temperature
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
    petsc_options_iname = '-ksp_type -pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = 'gmres      asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 1e4
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-5
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1
    growth_factor = 1.5
  []
[]

[Postprocessors]
  [pgas]
    type = PointValue
    point = '2 0 0'
    variable = pgas
  []
  [sgas]
    type = PointValue
    point = '2 0 0'
    variable = saturation_gas
  []
  [zi]
    type = PointValue
    point = '2 0 0'
    variable = zi
  []
  [temperature]
    type = PointValue
    point = '2 0 0'
    variable = temperature
  []
  [massgas]
    type = PorousFlowFluidMass
    fluid_component = 1
  []
  [x1]
    type = PointValue
    point = '2 0 0'
    variable = x1
  []
  [y0]
    type = PointValue
    point = '2 0 0'
    variable = y0
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  csv = true
[]

(modules/porous_flow/test/tests/jacobian/brineco2_twophase.i)

# Tests correct calculation of properties derivatives in PorousFlowFluidState
# for conditions that are appropriate for two phases

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [xnacl]
    initial_condition = 0.05
  []
[]

[Variables]
  [pgas]
  []
  [zi]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 1e6
    max = 4e6
    variable = pgas
    seed = 1
  []
  [z]
    type = RandomIC
    min = 0.2
    max = 0.8
    variable = zi
    seed = 2
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = zi
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = zi
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e1
    pc_max = 1e4
  []
  [fs]
    type = PorousFlowBrineCO2
    brine_fp = brine
    co2_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [brine]
    type = BrineFluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [brineco2]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    xnacl = xnacl
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/test/tests/fluidstate/waterncg_nonisothermal.i)

[Mesh]
  [mesh]
    type = GeneratedMeshGenerator
    dim = 2
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pgas]
    initial_condition = 1e6
  []
  [z]
    initial_condition = 0.25
  []
  [temperature]
    initial_condition = 70
  []
[]

[AuxVariables]
  [pressure_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [pressure_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [saturation_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [density_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [viscosity_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [enthalpy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [internal_energy_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [internal_energy_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x0_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_water]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1_gas]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [pressure_water]
    type = PorousFlowPropertyAux
    variable = pressure_water
    property = pressure
    phase = 0
    execute_on = timestep_end
  []
  [pressure_gas]
    type = PorousFlowPropertyAux
    variable = pressure_gas
    property = pressure
    phase = 1
    execute_on = timestep_end
  []
  [saturation_water]
    type = PorousFlowPropertyAux
    variable = saturation_water
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [density_water]
    type = PorousFlowPropertyAux
    variable = density_water
    property = density
    phase = 0
    execute_on = timestep_end
  []
  [density_gas]
    type = PorousFlowPropertyAux
    variable = density_gas
    property = density
    phase = 1
    execute_on = timestep_end
  []
  [viscosity_water]
    type = PorousFlowPropertyAux
    variable = viscosity_water
    property = viscosity
    phase = 0
    execute_on = timestep_end
  []
  [viscosity_gas]
    type = PorousFlowPropertyAux
    variable = viscosity_gas
    property = viscosity
    phase = 1
    execute_on = timestep_end
  []
  [enthalpy_water]
    type = PorousFlowPropertyAux
    variable = enthalpy_water
    property = enthalpy
    phase = 0
    execute_on = timestep_end
  []
  [enthalpy_gas]
    type = PorousFlowPropertyAux
    variable = enthalpy_gas
    property = enthalpy
    phase = 1
    execute_on = timestep_end
  []
  [internal_energy_water]
    type = PorousFlowPropertyAux
    variable = internal_energy_water
    property = internal_energy
    phase = 0
    execute_on = timestep_end
  []
  [internal_energy_gas]
    type = PorousFlowPropertyAux
    variable = internal_energy_gas
    property = internal_energy
    phase = 1
    execute_on = timestep_end
  []
  [x1_water]
    type = PorousFlowPropertyAux
    variable = x1_water
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [x1_gas]
    type = PorousFlowPropertyAux
    variable = x1_gas
    property = mass_fraction
    phase = 1
    fluid_component = 1
    execute_on = timestep_end
  []
  [x0_water]
    type = PorousFlowPropertyAux
    variable = x0_water
    property = mass_fraction
    phase = 0
    fluid_component = 0
    execute_on = timestep_end
  []
  [x0_gas]
    type = PorousFlowPropertyAux
    variable = x0_gas
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [heat]
    type = TimeDerivative
    variable = temperature
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z '
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temperature
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature = temperature
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 3
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [density_water]
    type = ElementIntegralVariablePostprocessor
    variable = density_water
  []
  [density_gas]
    type = ElementIntegralVariablePostprocessor
    variable = density_gas
  []
  [viscosity_water]
    type = ElementIntegralVariablePostprocessor
    variable = viscosity_water
  []
  [viscosity_gas]
    type = ElementIntegralVariablePostprocessor
    variable = viscosity_gas
  []
  [enthalpy_water]
    type = ElementIntegralVariablePostprocessor
    variable = enthalpy_water
  []
  [enthalpy_gas]
    type = ElementIntegralVariablePostprocessor
    variable = enthalpy_gas
  []
  [internal_energy_water]
    type = ElementIntegralVariablePostprocessor
    variable = internal_energy_water
  []
  [internal_energy_gas]
    type = ElementIntegralVariablePostprocessor
    variable = internal_energy_gas
  []
  [x0_water]
    type = ElementIntegralVariablePostprocessor
    variable = x0_water
  []
  [x1_gas]
    type = ElementIntegralVariablePostprocessor
    variable = x1_gas
  []
  [x0_gas]
    type = ElementIntegralVariablePostprocessor
    variable = x0_gas
  []
  [sg]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_gas
  []
  [sw]
    type = ElementIntegralVariablePostprocessor
    variable = saturation_water
  []
  [pwater]
    type = ElementIntegralVariablePostprocessor
    variable = pressure_water
  []
  [pgas]
    type = ElementIntegralVariablePostprocessor
    variable = pressure_gas
  []
  [x0mass]
    type = PorousFlowFluidMass
    fluid_component = 0
    phase = '0 1'
  []
  [x1mass]
    type = PorousFlowFluidMass
    fluid_component = 1
    phase = '0 1'
  []
[]

[Outputs]
  csv = true
  execute_on = timestep_end
[]

(modules/porous_flow/test/tests/jacobian/waterncg_twophase.i)

# Tests correct calculation of properties derivatives in PorousFlowWaterNCG
# for conditions for two phases

[Mesh]
  type = GeneratedMesh
  dim = 2
  nx = 2
  ny = 2
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Variables]
  [pgas]
  []
  [z]
  []
[]

[ICs]
  [pgas]
    type = RandomIC
    min = 1e5
    max = 5e5
    variable = pgas
  []
  [z]
    type = RandomIC
    min = 0.01
    max = 0.06
    variable = z
  []
[]

[Kernels]
  [mass0]
    type = PorousFlowMassTimeDerivative
    variable = pgas
    fluid_component = 0
  []
  [mass1]
    type = PorousFlowMassTimeDerivative
    variable = z
    fluid_component = 1
  []
  [adv0]
    type = PorousFlowAdvectiveFlux
    variable = pgas
    fluid_component = 0
  []
  [adv1]
    type = PorousFlowAdvectiveFlux
    variable = z
    fluid_component = 1
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas z'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    m = 0.5
    alpha = 1e1
    pc_max = 1e4
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 50
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = z
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm0]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
  []
  [relperm1]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.1
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  dt = 1
  end_time = 1
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[AuxVariables]
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [sgas]
    type = PorousFlowPropertyAux
    property = saturation
    phase = 1
    variable = sgas
  []
[]

[Postprocessors]
  [sgas_min]
    type = ElementExtremeValue
    variable = sgas
    value_type = min
  []
  [sgas_max]
    type = ElementExtremeValue
    variable = sgas
    value_type = max
  []
[]

(modules/porous_flow/examples/tutorial/11_2D.i)

# Two-phase borehole injection problem in RZ coordinates
[Mesh]
  [gen]
    type = GeneratedMeshGenerator
    dim = 2
    nx = 10
    xmin = 1.0
    xmax = 10
    bias_x = 1.4
    ny = 3
    ymin = -6
    ymax = 6
  []
  [aquifer]
    input = gen
    type = SubdomainBoundingBoxGenerator
    block_id = 1
    bottom_left = '0 -2 0'
    top_right = '10 2 0'
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x<1.0001'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'aquifer'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '0 1'
    new_block = 'caps aquifer'
    input = 'injection_area'
  []
  coord_type = RZ
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater pgas T disp_r'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    alpha = 1E-6
    m = 0.6
  []
[]

[GlobalParams]
  displacements = 'disp_r disp_z'
  gravity = '0 0 0'
  biot_coefficient = 1.0
  PorousFlowDictator = dictator
[]

[Variables]
  [pwater]
    initial_condition = 20E6
  []
  [pgas]
    initial_condition = 20.1E6
  []
  [T]
    initial_condition = 330
    scaling = 1E-5
  []
  [disp_r]
    scaling = 1E-5
  []
[]

[Kernels]
  [mass_water_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 0
    variable = pwater
  []
  [flux_water]
    type = PorousFlowAdvectiveFlux
    fluid_component = 0
    use_displaced_mesh = false
    variable = pwater
  []
  [vol_strain_rate_water]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 0
    variable = pwater
  []
  [mass_co2_dot]
    type = PorousFlowMassTimeDerivative
    fluid_component = 1
    variable = pgas
  []
  [flux_co2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 1
    use_displaced_mesh = false
    variable = pgas
  []
  [vol_strain_rate_co2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 1
    variable = pgas
  []
  [energy_dot]
    type = PorousFlowEnergyTimeDerivative
    variable = T
  []
  [advection]
    type = PorousFlowHeatAdvection
    use_displaced_mesh = false
    variable = T
  []
  [conduction]
    type = PorousFlowHeatConduction
    use_displaced_mesh = false
    variable = T
  []
  [vol_strain_rate_heat]
    type = PorousFlowHeatVolumetricExpansion
    variable = T
  []
  [grad_stress_r]
    type = StressDivergenceRZTensors
    temperature = T
    variable = disp_r
    eigenstrain_names = thermal_contribution
    use_displaced_mesh = false
    component = 0
  []
  [poro_r]
    type = PorousFlowEffectiveStressCoupling
    variable = disp_r
    use_displaced_mesh = false
    component = 0
  []
[]

[AuxVariables]
  [disp_z]
  []
  [effective_fluid_pressure]
    family = MONOMIAL
    order = CONSTANT
  []
  [mass_frac_phase0_species0]
    initial_condition = 1 # all water in phase=0
  []
  [mass_frac_phase1_species0]
    initial_condition = 0 # no water in phase=1
  []
  [sgas]
    family = MONOMIAL
    order = CONSTANT
  []
  [swater]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_rr]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_tt]
    family = MONOMIAL
    order = CONSTANT
  []
  [stress_zz]
    family = MONOMIAL
    order = CONSTANT
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
[]

[AuxKernels]
  [effective_fluid_pressure]
    type = ParsedAux
    coupled_variables = 'pwater pgas swater sgas'
    expression = 'pwater * swater + pgas * sgas'
    variable = effective_fluid_pressure
  []
  [swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  []
  [sgas]
    type = PorousFlowPropertyAux
    variable = sgas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [stress_rr_aux]
    type = RankTwoAux
    variable = stress_rr
    rank_two_tensor = stress
    index_i = 0
    index_j = 0
  []
  [stress_tt]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_tt
    index_i = 2
    index_j = 2
  []
  [stress_zz]
    type = RankTwoAux
    rank_two_tensor = stress
    variable = stress_zz
    index_i = 1
    index_j = 1
  []
  [porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
    execute_on = timestep_end
  []
[]

[BCs]
  [pinned_top_bottom_r]
    type = DirichletBC
    variable = disp_r
    value = 0
    boundary = 'top bottom'
  []
  [cavity_pressure_r]
    type = Pressure
    boundary = injection_area
    variable = disp_r
    postprocessor = constrained_effective_fluid_pressure_at_wellbore
    use_displaced_mesh = false
  []

  [cold_co2]
    type = DirichletBC
    boundary = injection_area
    variable = T
    value = 290 # injection temperature
    use_displaced_mesh = false
  []
  [constant_co2_injection]
    type = PorousFlowSink
    boundary = injection_area
    variable = pgas
    fluid_phase = 1
    flux_function = -1E-4
    use_displaced_mesh = false
  []

  [outer_water_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pwater
    fluid_phase = 0
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
  [outer_co2_removal]
    type = PorousFlowPiecewiseLinearSink
    boundary = right
    variable = pgas
    fluid_phase = 1
    pt_vals = '0 1E9'
    multipliers = '0 1E8'
    PT_shift = 20.1E6
    use_mobility = true
    use_relperm = true
    use_displaced_mesh = false
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedBicubicFluidProperties
    fp = true_water
    temperature_min = 275
    pressure_max = 1E8
    fluid_property_output_file = water97_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_11.csv
  []
  [true_co2]
    type = CO2FluidProperties
  []
  [tabulated_co2]
    type = TabulatedBicubicFluidProperties
    fp = true_co2
    temperature_min = 275
    pressure_max = 1E8
    fluid_property_output_file = co2_tabulated_11.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = co2_tabulated_11.csv
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = T
  []
  [saturation_calculator]
    type = PorousFlow2PhasePP
    phase0_porepressure = pwater
    phase1_porepressure = pgas
    capillary_pressure = pc
  []
  [massfrac]
    type = PorousFlowMassFraction
    mass_fraction_vars = 'mass_frac_phase0_species0 mass_frac_phase1_species0'
  []
  [water]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_water
    phase = 0
  []
  [co2]
    type = PorousFlowSingleComponentFluid
    fp = tabulated_co2
    phase = 1
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 4
    s_res = 0.1
    sum_s_res = 0.2
    phase = 0
  []
  [relperm_co2]
    type = PorousFlowRelativePermeabilityBC
    nw_phase = true
    lambda = 2
    s_res = 0.1
    sum_s_res = 0.2
    phase = 1
  []
  [porosity]
    type = PorousFlowPorosity
    fluid = true
    mechanical = true
    thermal = true
    porosity_zero = 0.1
    reference_temperature = 330
    reference_porepressure = 20E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
  []
  [permeability_aquifer]
    type = PorousFlowPermeabilityKozenyCarman
    block = aquifer
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-12
  []
  [permeability_caps]
    type = PorousFlowPermeabilityKozenyCarman
    block = caps
    poroperm_function = kozeny_carman_phi0
    phi0 = 0.1
    n = 2
    m = 2
    k0 = 1E-15
    k_anisotropy = '1 0 0  0 1 0  0 0 0.1'
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '2 0 0  0 2 0  0 0 2'
  []
  [rock_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    specific_heat_capacity = 1100
    density = 2300
  []

  [elasticity_tensor]
    type = ComputeIsotropicElasticityTensor
    youngs_modulus = 5E9
    poissons_ratio = 0.0
  []
  [strain]
    type = ComputeAxisymmetricRZSmallStrain
    eigenstrain_names = 'thermal_contribution initial_stress'
  []
  [thermal_contribution]
    type = ComputeThermalExpansionEigenstrain
    temperature = T
    thermal_expansion_coeff = 5E-6 # this is the linear thermal expansion coefficient
    eigenstrain_name = thermal_contribution
    stress_free_temperature = 330
  []
  [initial_strain]
    type = ComputeEigenstrainFromInitialStress
    initial_stress = '20E6 0 0  0 20E6 0  0 0 20E6'
    eigenstrain_name = initial_stress
  []
  [stress]
    type = ComputeLinearElasticStress
  []

  [effective_fluid_pressure]
    type = PorousFlowEffectiveFluidPressure
  []
  [volumetric_strain]
    type = PorousFlowVolumetricStrain
  []
[]

[Postprocessors]
  [effective_fluid_pressure_at_wellbore]
    type = PointValue
    variable = effective_fluid_pressure
    point = '1 0 0'
    execute_on = timestep_begin
    use_displaced_mesh = false
  []
  [constrained_effective_fluid_pressure_at_wellbore]
    type = FunctionValuePostprocessor
    function = constrain_effective_fluid_pressure
    execute_on = timestep_begin
  []
[]

[Functions]
  [constrain_effective_fluid_pressure]
    type = ParsedFunction
    symbol_names = effective_fluid_pressure_at_wellbore
    symbol_values = effective_fluid_pressure_at_wellbore
    expression = 'max(effective_fluid_pressure_at_wellbore, 20E6)'
  []
[]

[Preconditioning]
  active = basic
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
  [preferred_but_might_not_be_installed]
    type = SMP
    full = true
    petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
    petsc_options_value = ' lu       mumps'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = 1E3
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1E3
    growth_factor = 1.2
    optimal_iterations = 10
  []
  nl_abs_tol = 1E-7
[]

[Outputs]
  exodus = true
[]

(modules/porous_flow/test/tests/fluids/h2o.i)

# Test the density and viscosity calculated by the water Material
# Region 1 density
# Pressure 80 MPa
# Temperature 300K (26.85C)
# Water density should equal 1.0 / 0.971180894e-3 = 1029.7 kg/m^3 (IAPWS IF97)
# Water viscosity should equal 0.00085327 Pa.s (NIST webbook)
# Results are within expected accuracy

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pp'
    number_fluid_phases = 1
    number_fluid_components = 1
  []
[]

[Variables]
  [pp]
    initial_condition = 80e6
  []
[]

[Kernels]
  [dummy]
    type = Diffusion
    variable = pp
  []
[]

[AuxVariables]
  [temp]
    initial_condition = 300.0
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = temp
  []
  [ppss]
    type = PorousFlow1PhaseFullySaturated
    porepressure = pp
  []
  [water]
    type = PorousFlowSingleComponentFluid
    temperature_unit = Kelvin
    fp = water
    phase = 0
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Executioner]
  type = Steady
  solve_type = Newton
[]

[Postprocessors]
  [pressure]
    type = ElementIntegralVariablePostprocessor
    variable = pp
  []
  [temperature]
    type = ElementIntegralVariablePostprocessor
    variable = temp
  []
  [density]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_fluid_phase_density_qp0'
  []
  [viscosity]
    type = ElementIntegralMaterialProperty
    mat_prop = 'PorousFlow_viscosity_qp0'
  []
[]

[Outputs]
  execute_on = 'timestep_end'
  file_base = h2o
  csv = true
[]

(modules/porous_flow/test/tests/ics/fluidpropic.i)

# Test the correct calculation of fluid properties using PorousFlowFluidPropertyIC
#
# Variables:
# Pressure: 1 MPa
# Temperature: 323.15 K
#
# Fluid properties for water (reference values from NIST webbook)
# Density: 988.43 kg/m^3
# Enthalpy: 210.19 kJ/kg
# Internal energy: 2019.18 kJ/kg

[Mesh]
  type = GeneratedMesh
  dim = 2
[]

[Variables]
  [pressure]
    initial_condition = 1e6
  []
  [temperature]
    initial_condition = 323.15
  []
[]

[AuxVariables]
  [enthalpy]
  []
  [internal_energy]
  []
  [density]
  []
[]

[ICs]
  [enthalpy]
    type = PorousFlowFluidPropertyIC
    variable = enthalpy
    property = enthalpy
    porepressure = pressure
    temperature = temperature
    fp = water
  []
  [internal_energy]
    type = PorousFlowFluidPropertyIC
    variable = internal_energy
    property = internal_energy
    porepressure = pressure
    temperature = temperature
    fp = water
  []
  [density]
    type = PorousFlowFluidPropertyIC
    variable = density
    property = density
    porepressure = pressure
    temperature = temperature
    fp = water
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Kernels]
  [pressure]
    type = Diffusion
    variable = pressure
  []
  [temperature]
    type = Diffusion
    variable = temperature
  []
[]

[Executioner]
  type = Steady
  nl_abs_tol = 1e-12
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Postprocessors]
  [enthalpy]
    type = ElementAverageValue
    variable = enthalpy
    execute_on = 'initial timestep_end'
  []
  [internal_energy]
    type = ElementAverageValue
    variable = internal_energy
    execute_on = 'initial timestep_end'
  []
  [density]
    type = ElementAverageValue
    variable = density
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  csv = true
  execute_on = initial
[]

(modules/porous_flow/examples/ates/ates.i)

# Simulation designed to assess the recovery efficiency of a single-well ATES system
# Using KT stabilisation
# Boundary conditions: fixed porepressure and temperature at top, bottom and far end of model.

#####################################
flux_limiter = minmod # minmod, vanleer, mc, superbee, none
# depth of top of aquifer (m)
depth = 400

inject_fluid_mass = 1E8 # kg
produce_fluid_mass = ${inject_fluid_mass} # kg
inject_temp = 90 # degC

inject_time = 91 # days
store_time = 91 # days
produce_time = 91 # days
rest_time = 91 # days
num_cycles = 5 # Currently needs to be <= 10

cycle_length = '${fparse inject_time + store_time + produce_time + rest_time}'

end_simulation = '${fparse cycle_length * num_cycles}'

# Note: I have setup 10 cycles but you can set num_cycles less than 10.
start_injection1 = 0
start_injection2 = ${cycle_length}
start_injection3 = '${fparse cycle_length * 2}'
start_injection4 = '${fparse cycle_length * 3}'
start_injection5 = '${fparse cycle_length * 4}'
start_injection6 = '${fparse cycle_length * 5}'
start_injection7 = '${fparse cycle_length * 6}'
start_injection8 = '${fparse cycle_length * 7}'
start_injection9 = '${fparse cycle_length * 8}'
start_injection10 = '${fparse cycle_length * 9}'

end_injection1 = '${fparse start_injection1 + inject_time}'
end_injection2 = '${fparse start_injection2 + inject_time}'
end_injection3 = '${fparse start_injection3 + inject_time}'
end_injection4 = '${fparse start_injection4 + inject_time}'
end_injection5 = '${fparse start_injection5 + inject_time}'
end_injection6 = '${fparse start_injection6 + inject_time}'
end_injection7 = '${fparse start_injection7 + inject_time}'
end_injection8 = '${fparse start_injection8 + inject_time}'
end_injection9 = '${fparse start_injection9 + inject_time}'
end_injection10 = '${fparse start_injection10 + inject_time}'

start_production1 = '${fparse end_injection1 + store_time}'
start_production2 = '${fparse end_injection2 + store_time}'
start_production3 = '${fparse end_injection3 + store_time}'
start_production4 = '${fparse end_injection4 + store_time}'
start_production5 = '${fparse end_injection5 + store_time}'
start_production6 = '${fparse end_injection6 + store_time}'
start_production7 = '${fparse end_injection7 + store_time}'
start_production8 = '${fparse end_injection8 + store_time}'
start_production9 = '${fparse end_injection9 + store_time}'
start_production10 = '${fparse end_injection10 + store_time}'

end_production1 = '${fparse start_production1 + produce_time}'
end_production2 = '${fparse start_production2 + produce_time}'
end_production3 = '${fparse start_production3 + produce_time}'
end_production4 = '${fparse start_production4 + produce_time}'
end_production5 = '${fparse start_production5 + produce_time}'
end_production6 = '${fparse start_production6 + produce_time}'
end_production7 = '${fparse start_production7 + produce_time}'
end_production8 = '${fparse start_production8 + produce_time}'
end_production9 = '${fparse start_production9 + produce_time}'
end_production10 = '${fparse start_production10 + produce_time}'

synctimes = '${start_injection1} ${end_injection1} ${start_production1} ${end_production1}
             ${start_injection2} ${end_injection2} ${start_production2} ${end_production2}
             ${start_injection3} ${end_injection3} ${start_production3} ${end_production3}
             ${start_injection4} ${end_injection4} ${start_production4} ${end_production4}
             ${start_injection5} ${end_injection5} ${start_production5} ${end_production5}
             ${start_injection6} ${end_injection6} ${start_production6} ${end_production6}
             ${start_injection7} ${end_injection7} ${start_production7} ${end_production7}
             ${start_injection8} ${end_injection8} ${start_production8} ${end_production8}
             ${start_injection9} ${end_injection9} ${start_production9} ${end_production9}
             ${start_injection10} ${end_injection10} ${start_production10} ${end_production10}'

#####################################
# Geometry in RZ coordinates
# borehole radius (m)
bh_r = 0.1
# model radius (m)
max_r = 1000
# aquifer thickness (m)
aq_thickness = 20
# cap thickness (m)
cap_thickness = 40
# injection region top and bottom (m).  Note, the mesh is created with the aquifer in y = (-0.5 * aq_thickness, 0.5 * aq_thickness), irrespective of depth (depth only sets the insitu porepressure and temperature)
screen_top = '${fparse 0.5 * aq_thickness}'
screen_bottom = '${fparse -0.5 * aq_thickness}'

# number of elements in radial direction
num_r = 25
# number of elements across half height of aquifer
num_y_aq = 10
# number of elements across height of cap
num_y_cap = 8

# mesh bias in radial direction
bias_r = 1.22
# mesh bias in vertical direction in aquifer top
bias_y_aq_top = 0.9
# mesh bias in vertical direction in cap top
bias_y_cap_top = 1.3
# mesh bias in vertical direction in aquifer bottom
bias_y_aq_bottom = '${fparse 1.0 / bias_y_aq_top}'
# mesh bias in vertical direction in cap bottom
bias_y_cap_bottom = '${fparse 1.0 / bias_y_cap_top}'
depth_centre = '${fparse depth + aq_thickness/2}'

#####################################
# temperature at ground surface (degC)
temp0 = 20
# Vertical geothermal gradient (K/m).  A positive number means temperature increases downwards.
geothermal_gradient = 20E-3

#####################################
# Gravity
gravity = -9.81

#####################################
half_aq_thickness = '${fparse aq_thickness * 0.5}'
half_height = '${fparse half_aq_thickness + cap_thickness}'
approx_screen_length = '${fparse screen_top - screen_bottom}'

# Thermal radius (note this is not strictly correct, it should use the bulk specific heat
#  capacity as defined below, but it doesn't matter here because this is purely for
#  defining the region of refined mesh)
th_r = '${fparse sqrt(inject_fluid_mass / 1000 * 4.12e6 / (approx_screen_length * 3.1416 * aq_specific_heat_cap * aq_density))}'
# radius of fine mesh
fine_r = '${fparse th_r * 2}'
bias_r_fine = 1
num_r_fine = '${fparse int(fine_r/1)}'

######################################
# aquifer properties
aq_porosity = 0.25
aq_hor_perm = 1E-11 # m^2
aq_ver_perm = 2E-12 # m^2
aq_density = 2650 # kg/m^3
aq_specific_heat_cap = 800 # J/Kg/K
aq_hor_thermal_cond = 3 # W/m/K
aq_ver_thermal_cond = 3 # W/m/K
aq_disp_parallel = 0 # m
aq_disp_perp = 0 # m
# Bulk volumetric heat capacity of aquifer:
aq_vol_cp = '${fparse aq_specific_heat_cap * aq_density * (1 - aq_porosity) + 4180 * 1000 * aq_porosity}'
# Thermal radius (correct version using bulk cp):
R_th = '${fparse sqrt(inject_fluid_mass * 4180 / (approx_screen_length * 3.1416 * aq_vol_cp))}'
aq_lambda_eff_hor = '${fparse aq_hor_thermal_cond + 0.3 * aq_disp_parallel * R_th * aq_vol_cp / (inject_time * 60 * 60 * 24)}'
aq_lambda_eff_ver = '${fparse aq_ver_thermal_cond + 0.3 * aq_disp_perp * R_th * aq_vol_cp / (inject_time * 60 * 60 * 24)}'
aq_hor_dry_thermal_cond = '${fparse aq_lambda_eff_hor * 60 * 60 * 24}' # J/day/m/K
aq_ver_dry_thermal_cond = '${fparse aq_lambda_eff_ver * 60 * 60 * 24}' # J/day/m/K
aq_hor_wet_thermal_cond = '${fparse aq_lambda_eff_hor * 60 * 60 * 24}' # J/day/m/K
aq_ver_wet_thermal_cond = '${fparse aq_lambda_eff_ver * 60 * 60 * 24}' # J/day/m/K
# cap-rock properties
cap_porosity = 0.25
cap_hor_perm = 1E-16 # m^2
cap_ver_perm = 1E-17 # m^2
cap_density = 2650 # kg/m^3
cap_specific_heat_cap = 800 # J/kg/K
cap_hor_thermal_cond = 3 # W/m/K
cap_ver_thermal_cond = 3 # W/m/K
cap_hor_dry_thermal_cond = '${fparse cap_hor_thermal_cond * 60 * 60 * 24}' # J/day/m/K
cap_ver_dry_thermal_cond = '${fparse cap_ver_thermal_cond * 60 * 60 * 24}' # J/day/m/K
cap_hor_wet_thermal_cond = '${fparse cap_hor_thermal_cond * 60 * 60 * 24}' # J/day/m/K
cap_ver_wet_thermal_cond = '${fparse cap_ver_thermal_cond * 60 * 60 * 24}' # J/day/m/K

######################################

[Mesh]
  coord_type = RZ
  [aq_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_top_fine cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom_fine]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r_fine}
    xmin = ${bh_r}
    xmax = ${fine_r}
    bias_x = ${bias_r_fine}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom_fine cap_bottom_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
    merge_boundaries_with_same_name = false
  []
  [aq_and_cap_fine]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom_fine aq_and_cap_top_fine'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_top}
    ny = ${num_y_aq}
    ymin = 0
    ymax = ${half_aq_thickness}
  []
  [cap_top]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_top}
    ny = ${num_y_cap}
    ymax = ${half_height}
    ymin = ${half_aq_thickness}
  []
  [aq_and_cap_top]
    type = StitchedMeshGenerator
    inputs = 'aq_top cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []
  [aq_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_aq_bottom}
    ny = ${num_y_aq}
    ymax = 0
    ymin = -${half_aq_thickness}
  []
  [cap_bottom]
    type = GeneratedMeshGenerator
    dim = 2
    nx = ${num_r}
    xmin = ${fine_r}
    xmax = ${max_r}
    bias_x = ${bias_r}
    bias_y = ${bias_y_cap_bottom}
    ny = ${num_y_cap}
    ymin = -${half_height}
    ymax = -${half_aq_thickness}
  []
  [aq_and_cap_bottom]
    type = StitchedMeshGenerator
    inputs = 'aq_bottom cap_bottom'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'bottom top'
  []
  [aq_and_cap]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_bottom aq_and_cap_top'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'top bottom'
  []

  [aq_and_cap_all]
    type = StitchedMeshGenerator
    inputs = 'aq_and_cap_fine aq_and_cap'
    clear_stitched_boundary_ids = true
    stitch_boundaries_pairs = 'right left'
  []

  [aquifer]
    type = ParsedSubdomainMeshGenerator
    input = aq_and_cap_all
    combinatorial_geometry = 'y >= -${half_aq_thickness} & y <= ${half_aq_thickness}'
    block_id = 1
  []
  [top_cap]
    type = ParsedSubdomainMeshGenerator
    input = aquifer
    combinatorial_geometry = 'y >= ${half_aq_thickness}'
    block_id = 2
  []
  [bottom_cap]
    type = ParsedSubdomainMeshGenerator
    input = top_cap
    combinatorial_geometry = 'y <= -${half_aq_thickness}'
    block_id = 3
  []
  [injection_area]
    type = ParsedGenerateSideset
    combinatorial_geometry = 'x<=${bh_r}*1.000001 & y >= ${screen_bottom} & y <= ${screen_top}'
    included_subdomains = 1
    new_sideset_name = 'injection_area'
    input = 'bottom_cap'
  []
  [rename]
    type = RenameBlockGenerator
    old_block = '1 2 3'
    new_block = 'aquifer caps caps'
    input = 'injection_area'
  []
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 ${gravity} 0'
[]

[Variables]
  [porepressure]
  []
  [temperature]
    scaling = 1E-5
  []
[]

[PorousFlowFullySaturated]
  coupling_type = ThermoHydro
  porepressure = porepressure
  temperature = temperature
  fp = tabulated_water
  stabilization = KT
  flux_limiter_type = ${flux_limiter}
  use_displaced_mesh = false
  temperature_unit = Celsius
  pressure_unit = Pa
  time_unit = days
[]

[ICs]
  [porepressure]
    type = FunctionIC
    variable = porepressure
    function = insitu_pressure
  []
  [temperature]
    type = FunctionIC
    variable = temperature
    function = insitu_temperature
  []
[]

[BCs]
  [outer_boundary_porepressure]
    type = FunctionDirichletBC
    preset = true
    variable = porepressure
    function = insitu_pressure
    boundary = 'bottom right top'
  []
  [outer_boundary_temperature]
    type = FunctionDirichletBC
    preset = true
    variable = temperature
    function = insitu_temperature
    boundary = 'bottom right top'
  []
  [inject_heat]
    type = FunctionDirichletBC
    variable = temperature
    function = ${inject_temp}
    boundary = 'injection_area'
  []
  [inject_fluid]
    type = PorousFlowSink
    variable = porepressure
    boundary = injection_area
    flux_function = injection_rate_value
  []
  [produce_heat]
    type = PorousFlowSink
    variable = temperature
    boundary = injection_area
    flux_function = production_rate_value
    fluid_phase = 0
    use_enthalpy = true
    save_in = heat_flux_out
  []
  [produce_fluid]
    type = PorousFlowSink
    variable = porepressure
    boundary = injection_area
    flux_function = production_rate_value
  []
[]

[Controls]
  [inject_on]
    type = ConditionalFunctionEnableControl
    enable_objects = 'BCs::inject_heat BCs::inject_fluid'
    conditional_function = inject
    implicit = false
    execute_on = 'initial timestep_begin'
  []
  [produce_on]
    type = ConditionalFunctionEnableControl
    enable_objects = 'BCs::produce_heat BCs::produce_fluid'
    conditional_function = produce
    implicit = false
    execute_on = 'initial timestep_begin'
  []
[]

[Functions]
  [insitu_pressure]
    type = ParsedFunction
    expression = '(y - ${depth_centre}) * 1000 * ${gravity} + 1E5' # approx insitu pressure in Pa
  []
  [insitu_temperature]
    type = ParsedFunction
    expression = '${temp0} + (${depth_centre} - y) * ${geothermal_gradient}'
  []

  [inject]
    type = ParsedFunction
    expression = 'if(t >= ${start_injection1} & t < ${end_injection1}, 1,
             if(t >= ${start_injection2} & t < ${end_injection2}, 1,
             if(t >= ${start_injection3} & t < ${end_injection3}, 1,
             if(t >= ${start_injection4} & t < ${end_injection4}, 1,
             if(t >= ${start_injection5} & t < ${end_injection5}, 1,
             if(t >= ${start_injection6} & t < ${end_injection6}, 1,
             if(t >= ${start_injection7} & t < ${end_injection7}, 1,
             if(t >= ${start_injection8} & t < ${end_injection8}, 1,
             if(t >= ${start_injection9} & t < ${end_injection9}, 1,
             if(t >= ${start_injection10} & t < ${end_injection10}, 1, 0))))))))))'
  []
  [produce]
    type = ParsedFunction
    expression = 'if(t >= ${start_production1} & t < ${end_production1}, 1,
             if(t >= ${start_production2} & t < ${end_production2}, 1,
             if(t >= ${start_production3} & t < ${end_production3}, 1,
             if(t >= ${start_production4} & t < ${end_production4}, 1,
             if(t >= ${start_production5} & t < ${end_production5}, 1,
             if(t >= ${start_production6} & t < ${end_production6}, 1,
             if(t >= ${start_production7} & t < ${end_production7}, 1,
             if(t >= ${start_production8} & t < ${end_production8}, 1,
             if(t >= ${start_production9} & t < ${end_production9}, 1,
             if(t >= ${start_production10} & t < ${end_production10}, 1, 0))))))))))'
  []

  [injection_rate_value]
    type = ParsedFunction
    symbol_names = true_screen_area
    symbol_values = true_screen_area
    expression = '-${inject_fluid_mass}/(true_screen_area * ${inject_time})'
  []
  [production_rate_value]
    type = ParsedFunction
    symbol_names = true_screen_area
    symbol_values = true_screen_area
    expression = '${produce_fluid_mass}/(true_screen_area * ${produce_time})'
  []

  [heat_out_in_timestep]
    type = ParsedFunction
    symbol_names = 'dt heat_out'
    symbol_values = 'dt heat_out_fromBC'
    expression = 'dt*heat_out'
  []
  [produced_T_time_integrated]
    type = ParsedFunction
    symbol_names = 'dt produced_T'
    symbol_values = 'dt produced_T'
    expression = 'dt*produced_T / ${produce_time}'
  []
[]

[AuxVariables]
  [density]
    family = MONOMIAL
    order = CONSTANT
  []
  [porosity]
    family = MONOMIAL
    order = CONSTANT
  []
  [heat_flux_out]
    outputs = none
  []
[]

[AuxKernels]
  [density]
    type = PorousFlowPropertyAux
    variable = density
    property = density
  []
  [porosity]
    type = PorousFlowPropertyAux
    variable = porosity
    property = porosity
  []
[]

[FluidProperties]
  [true_water]
    type = Water97FluidProperties
  []
  [tabulated_water]
    type = TabulatedFluidProperties
    fp = true_water
    temperature_min = 275 # K
    temperature_max = 600
    interpolated_properties = 'density viscosity enthalpy internal_energy'
    fluid_property_output_file = water97_tabulated_modified.csv
    # Comment out the fp parameter and uncomment below to use the newly generated tabulation
    # fluid_property_file = water97_tabulated_modified.csv
  []
[]

[Materials]
  [porosity_aq]
    type = PorousFlowPorosityConst
    porosity = ${aq_porosity}
    block = aquifer
  []
  [porosity_caps]
    type = PorousFlowPorosityConst
    porosity = ${cap_porosity}
    block = caps
  []

  [permeability_aquifer]
    type = PorousFlowPermeabilityConst
    block = aquifer
    permeability = '${aq_hor_perm} 0 0   0 ${aq_ver_perm} 0   0 0 0'
  []
  [permeability_caps]
    type = PorousFlowPermeabilityConst
    block = caps
    permeability = '${cap_hor_perm} 0 0   0 ${cap_ver_perm} 0   0 0 0'
  []

  [aq_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    block = aquifer
    density = ${aq_density}
    specific_heat_capacity = ${aq_specific_heat_cap}
  []
  [caps_internal_energy]
    type = PorousFlowMatrixInternalEnergy
    block = caps
    density = ${cap_density}
    specific_heat_capacity = ${cap_specific_heat_cap}
  []

  [aq_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    block = aquifer
    dry_thermal_conductivity = '${aq_hor_dry_thermal_cond} 0 0  0 ${aq_ver_dry_thermal_cond} 0  0 0 0'
    wet_thermal_conductivity = '${aq_hor_wet_thermal_cond} 0 0  0 ${aq_ver_wet_thermal_cond} 0  0 0 0'
  []
  [caps_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    block = caps
    dry_thermal_conductivity = '${cap_hor_dry_thermal_cond} 0 0  0 ${cap_ver_dry_thermal_cond} 0  0 0 0'
    wet_thermal_conductivity = '${cap_hor_wet_thermal_cond} 0 0  0 ${cap_ver_wet_thermal_cond} 0  0 0 0'
  []
[]

[Postprocessors]
  [true_screen_area] # this accounts for meshes that do not match screen_top and screen_bottom exactly
    type = AreaPostprocessor
    boundary = injection_area
    execute_on = 'initial'
    outputs = 'none'
  []
  [dt]
    type = TimestepSize
  []

  [heat_out_fromBC]
    type = NodalSum
    variable = heat_flux_out
    boundary = injection_area
    execute_on = 'initial timestep_end'
    outputs = 'none'
  []
  [heat_out_per_timestep]
    type = FunctionValuePostprocessor
    function = heat_out_in_timestep
    execute_on = 'timestep_end'
    outputs = 'none'
  []
  [heat_out_cumulative]
    type = CumulativeValuePostprocessor
    postprocessor = heat_out_per_timestep
    execute_on = 'timestep_end'
    outputs = 'csv console'
  []
  [produced_T]
    type = SideAverageValue
    boundary = injection_area
    variable = temperature
    execute_on = 'initial timestep_end'
    outputs = 'csv console'
  []
  [produced_T_time_integrated]
    type = FunctionValuePostprocessor
    function = produced_T_time_integrated
    execute_on = 'timestep_end'
    outputs = 'none'
  []
  [produced_T_cumulative]
    type = CumulativeValuePostprocessor
    postprocessor = produced_T_time_integrated
    execute_on = 'timestep_end'
    outputs = 'csv console'
  []
[]

[Preconditioning]
  [basic]
    type = SMP
    full = true
    petsc_options = '-ksp_diagonal_scale -ksp_diagonal_scale_fix'
    petsc_options_iname = '-pc_type -sub_pc_type -sub_pc_factor_shift_type -pc_asm_overlap'
    petsc_options_value = ' asm      lu           NONZERO                   2'
  []
[]

[Executioner]
  type = Transient
  solve_type = Newton
  end_time = ${end_simulation}
  timestep_tolerance = 1e-5

  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 1e-3
    growth_factor = 2
  []

  dtmax = 1
  dtmin = 1e-5
  # rough calc for fluid, |R| ~ V*k*1E6 ~ V*1E-5
  # rough calc for heat, |R| ~ V*(lam*1E-3 + h*1E-5)  ~ V*(1E3 + 1E-2)
  # so scale heat by 1E-7 and go for nl_abs_tol = 1E-4, which should give a max error of
  # ~1Pa and ~0.1K in the first metre around the borehole
  nl_abs_tol = 1E-4
  nl_rel_tol = 1E-5
[]

[Outputs]
  sync_times = ${synctimes}
  [ex]
    type = Exodus
    time_step_interval = 20
  []
  [csv]
    type = CSV
    execute_postprocessors_on = 'initial timestep_end'
  []
[]

(modules/porous_flow/test/tests/fluidstate/coldwater_injection.i)

# Cold water injection into 1D hot reservoir (Avdonin, 1964)
#
# To generate results presented in documentation for this problem,
# set xmax = 50 and nx = 250 in the Mesh block, and dtmax = 100 and
# end_time = 1.3e5 in the Executioner block.

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 25
  xmax = 20
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[AuxVariables]
  [temperature]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [temperature]
    type = PorousFlowPropertyAux
    variable = temperature
    property = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Variables]
  [pliquid]
    initial_condition = 5e6
  []
  [h]
    scaling = 1e-6
  []
[]

[ICs]
  [hic]
    type = PorousFlowFluidPropertyIC
    variable = h
    porepressure = pliquid
    property = enthalpy
    temperature = 170
    temperature_unit = Celsius
    fp = water
  []
[]

[BCs]
  [pleft]
    type = DirichletBC
    variable = pliquid
    value = 5.05e6
    boundary = left
  []
  [pright]
    type = DirichletBC
    variable = pliquid
    value = 5e6
    boundary = right
  []
  [hleft]
    type = DirichletBC
    variable = h
    value = 678.52e3
    boundary = left
  []
  [hright]
    type = DirichletBC
    variable = h
    value = 721.4e3
    boundary = right
  []
[]

[Kernels]
  [mass]
    type = PorousFlowMassTimeDerivative
    variable = pliquid
  []
  [massflux]
    type = PorousFlowAdvectiveFlux
    variable = pliquid
  []
  [heat]
    type = PorousFlowEnergyTimeDerivative
    variable = h
  []
  [heatflux]
    type = PorousFlowHeatAdvection
    variable = h
  []
  [heatcond]
    type = PorousFlowHeatConduction
    variable = h
  []
[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pliquid h'
    number_fluid_phases = 2
    number_fluid_components = 1
  []
  [pc]
    type = PorousFlowCapillaryPressureVG
    pc_max = 1e6
    sat_lr = 0.1
    m = 0.5
    alpha = 1e-5
  []
  [fs]
    type = PorousFlowWaterVapor
    water_fp = water
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [watervapor]
    type = PorousFlowFluidStateSingleComponent
    porepressure = pliquid
    enthalpy = h
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1.8e-11 0 0 0 1.8e-11 0 0 0 1.8e-11'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
    sum_s_res = 0.1
  []
  [internal_energy]
    type = PorousFlowMatrixInternalEnergy
    density = 2900
    specific_heat_capacity = 740
  []
  [rock_thermal_conductivity]
    type = PorousFlowThermalConductivityIdeal
    dry_thermal_conductivity = '20 0 0  0 20 0  0 0 20'
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 5e3
  nl_abs_tol = 1e-10
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 100
  []
[]

[VectorPostprocessors]
  [line]
    type = ElementValueSampler
    sort_by = x
    variable = temperature
    execute_on = 'initial timestep_end'
  []
[]

[Outputs]
  perf_graph = true
  [csv]
    type = CSV
    execute_on = final
  []
[]

(modules/porous_flow/test/tests/recover/theis.i)

# Tests that PorousFlow can successfully recover using a checkpoint file.
# This test contains stateful material properties, adaptivity and integrated
# boundary conditions with nodal-sized materials.
#
# This test file is run three times:
# 1) The full input file is run to completion
# 2) The input file is run for half the time and checkpointing is included
# 3) The input file is run in recovery using the checkpoint data
#
# The final output of test 3 is compared to the final output of test 1 to verify
# that recovery was successful.

[Mesh]
  type = GeneratedMesh
  dim = 1
  nx = 20
  xmax = 100
  bias_x = 1.05
  coord_type = RZ
  rz_coord_axis = Y
[]

[GlobalParams]
  PorousFlowDictator = dictator
  gravity = '0 0 0'
[]

[Adaptivity]
  marker = marker
  max_h_level = 4
  [Indicators]
    [front]
      type = GradientJumpIndicator
      variable = zi
    []
  []
  [Markers]
    [marker]
      type = ErrorFractionMarker
      indicator = front
      refine = 0.8
      coarsen = 0.2
    []
  []
[]

[AuxVariables]
  [saturation_gas]
    order = CONSTANT
    family = MONOMIAL
  []
  [x1]
    order = CONSTANT
    family = MONOMIAL
  []
  [y0]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[AuxKernels]
  [saturation_gas]
    type = PorousFlowPropertyAux
    variable = saturation_gas
    property = saturation
    phase = 1
    execute_on = timestep_end
  []
  [x1]
    type = PorousFlowPropertyAux
    variable = x1
    property = mass_fraction
    phase = 0
    fluid_component = 1
    execute_on = timestep_end
  []
  [y0]
    type = PorousFlowPropertyAux
    variable = y0
    property = mass_fraction
    phase = 1
    fluid_component = 0
    execute_on = timestep_end
  []
[]

[Variables]
  [pgas]
    initial_condition = 20e6
  []
  [zi]
    initial_condition = 0
  []
[]

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

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pgas zi'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [pc]
    type = PorousFlowCapillaryPressureConst
    pc = 0
  []
  [fs]
    type = PorousFlowWaterNCG
    water_fp = water
    gas_fp = co2
    capillary_pressure = pc
  []
[]

[FluidProperties]
  [co2]
    type = CO2FluidProperties
  []
  [water]
    type = Water97FluidProperties
  []
[]

[Materials]
  [temperature]
    type = PorousFlowTemperature
    temperature = 20
  []
  [waterncg]
    type = PorousFlowFluidState
    gas_porepressure = pgas
    z = zi
    temperature_unit = Celsius
    capillary_pressure = pc
    fluid_state = fs
  []
  [porosity]
    type = PorousFlowPorosityConst
    porosity = 0.2
  []
  [permeability]
    type = PorousFlowPermeabilityConst
    permeability = '1e-12 0 0 0 1e-12 0 0 0 1e-12'
  []
  [relperm_water]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 0
    s_res = 0.1
    sum_s_res = 0.1
  []
  [relperm_gas]
    type = PorousFlowRelativePermeabilityCorey
    n = 2
    phase = 1
  []
[]

[BCs]
  [aquifer]
    type = PorousFlowPiecewiseLinearSink
    variable = pgas
    boundary = right
    pt_vals = '0 1e8'
    multipliers = '0 1e8'
    flux_function = 1e-6
    PT_shift = 20e6
  []
[]

[DiracKernels]
  [source]
    type = PorousFlowSquarePulsePointSource
    point = '0 0 0'
    mass_flux = 2
    variable = zi
  []
[]

[Preconditioning]
  [smp]
    type = SMP
    full = true
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  end_time = 2e2
  dt = 50
[]

[VectorPostprocessors]
  [line]
    type = NodalValueSampler
    sort_by = x
    variable = 'pgas zi'
  []
[]

[Outputs]
  print_linear_residuals = false
  perf_graph = true
  csv = true
[]

Properties of water
Range of validity
Input Parameters
Input Files
References