Brine and carbon dioxide
The brine-CO model available in PorousFlow is a high precision equation of state for brine and CO, including the mutual solubility of CO into the liquid brine and water into the CO-rich gas phase. This model is suitable for simulations of geological storage of CO in saline aquifers, and is valid for temperatures in the range and pressures less than 60 MPa.
Salt (NaCl) can be included as a nonlinear variable, allowing salt to be transported along with water to provide variable density brine.
Phase composition
The mass fractions of CO in the liquid phase and HO in the gas phase are calculated using the accurate fugacity-based formulation of Spycher et al. (2003) and Spycher et al. (2005) for temperatures below 100C, and the elevated temperature formulation of Spycher and Pruess (2010) for temperatures above 110C.
As these formulations do not coincide for temperatures near 100C, a cubic polynomial is used to join the two curves smoothly, see Figure 1 for an example:
The mutual solubilities in the elevated temperature regime must be calculated iteratively, so an increase in computational expense can be expected.
This is similar to the formulation provided in the ECO2N module of TOUGH2 (Pruess et al., 1999), see Figure 2 and Figure 3 for a comparison between the two codes.
As with many reservoir simulators, it is assumed that the phases are in instantaneous equilibrium, meaning that dissolution of CO into the brine and evaporation of water vapor into the CO phase happens instantaneously.
If the amount of CO present is not sufficient to saturate the brine, then no gas phase will be formed, and the dissolved CO mass fraction will be smaller than its equilibrium value. Similarly, if the amount of water vapor present is smaller than its equilibrium value in the gas phase, then no liquid phase can be formed.
If there is sufficient CO present to saturate the brine, then both liquid and gas phases will be present. In this case, the mass fraction of CO dissolved in the brine is given by its equilibrium value (for the current pressure, temperature and salinity). Similarly, the mass fraction of water vapor in the gas phase will be given by its equilibrium value in this two phase region.
Thermophysical properties
Density
The density of the aqueous phase with the contribution of dissolved CO is calculated using where is the density of brine (supplied using a BrineFluidProperties
UserObject), is the mass fraction of CO dissolved in the aqueous phase, and is the partial density of dissolved CO (Garcia, 2001).
As water vapor is only ever a small component of the gas phase in the temperature and pressure ranges that this class is valid for, the density of the gas phase is assumed to be simply the density of CO at the given pressure and temperature, calculated using a CO2FluidProperties
UserObject.
Viscosity
No contribution to the viscosity of each phase due to the presence of CO in the aqueous phase or water vapor in the gas phase is included. As a result, the viscosity of the aqueous phase is simply the viscosity of brine, while the viscosity of the gas phase is the viscosity of CO.
Enthalpy
The enthalpy of the gas phase is simply the enthalpy of CO at the given pressure and temperature values, with no contribution due to the presence of water vapor included.
The enthalpy of the liquid phase is also calculated as a weighted sum of the enthalpies of each individual component, but also includes a contribution due to the enthalpy of dissolution (a change in enthalpy due to dissolution of the NCG into the water phase) where is the enthalpy of brine, and is the enthalpy of CO in the liquid phase, which includes the enthalpy of dissolution
In the range of pressures and temperatures considered, CO may exist as a gas or a supercritical fluid. Using a linear fit to the model of Duan and Sun (2003), the enthalpy of dissolution of both gas phase and supercritical CO is calculated as
The calculation of brine and CO fluid properties can take a significant proportion of each simulation, so it is suggested that tabulated versions of both CO and brine properties are used using TabulatedFluidProperties.
Implementation
Variables
This class requires pressure of the gas phase, the total mass fraction of CO summed over all phases (see compositional flash for details), and the mass fraction of salt in the brine, for isothermal simulations. For nonisothermal simulations, temperature must also be provided as a nonlinear variable.
These variables must be listed as PorousFlow variables in the PorousFlowDictator UserObject
In the general isothermal case, three variables (gas porepressure, total CO mass fraction and salt mass fraction) are provided, one for each of the three possible fluid components (water, salt and CO). As a result, the number of components in the PorousFlowDictator must be set equal to three.
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi xnacl'
number_fluid_phases = 2
number_fluid_components = 3
[]
[]
(modules/porous_flow/test/tests/fluidstate/theis_brineco2.i)For nonisothermal cases, temperature must also be provided to ensure that all of the derivatives with respect to the temperature variable are calculated.
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi xnacl temperature'
number_fluid_phases = 2
number_fluid_components = 3
[]
[]
(modules/porous_flow/test/tests/fluidstate/theis_brineco2_nonisothermal.i)Optionally, the salt mass fraction can be treated as a constant, in which case only gas porepressure and total CO mass fraction are required for the isothermal case, and the number of fluid components is two.
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas z'
number_fluid_phases = 2
number_fluid_components = 2
[]
[]
(modules/porous_flow/test/tests/fluidstate/brineco2.i)UserObjects
This fluid state is implemented in the PorousFlowBrineCO2
UserObject. This UserObject is a GeneralUserObject
that contains methods that provide a complete thermophysical description of the model given pressure, temperature and salt mass fraction.
[UserObjects]
[fs]
type = PorousFlowBrineCO2
brine_fp = brine
co2_fp = co2
capillary_pressure = pc
[]
[]
(modules/porous_flow/test/tests/fluidstate/brineco2.i)Materials
The PorousFlowFluidState
material provides all phase pressures, saturation, densities, viscosities etc, as well as all mass fractions of all fluid components in all fluid phases in a single material using the formulation provided in the PorousFlowBrineCO2
UserObject.
[Materials]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = pgas
z = z
temperature_unit = Celsius
xnacl = xnacl
capillary_pressure = pc
fluid_state = fs
[]
[]
(modules/porous_flow/test/tests/fluidstate/brineco2.i)For nonisothermal simulations, the temperature variable must also be supplied to ensure that all the derivatives with respect to the temperature variable are computed for the Jacobian.
[Materials]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = pgas
z = zi
temperature = temperature
temperature_unit = Celsius
xnacl = xnacl
capillary_pressure = pc
fluid_state = fs
[]
[]
(modules/porous_flow/test/tests/fluidstate/theis_brineco2_nonisothermal.i)Initial condition
The nonlinear variable representing CO is the total mass fraction of CO summed over all phases. In some cases, it may be preferred to provide an initial CO saturation, rather than total mass fraction. To allow an initial saturation to be specified, the PorousFlowFluidStateIC
initial condition is provided. This initial condition calculates the total mass fraction of CO summed over all phases for a given saturation.
[ICs]
[z]
type = PorousFlowFluidStateIC
saturation = 0.5
gas_porepressure = pgas
temperature = 50
variable = z
xnacl = 0.1
fluid_state = fs
[]
[]
(modules/porous_flow/test/tests/fluidstate/brineco2_ic.i)Example
1D Radial injection
Injection into a one dimensional radial reservoir is a useful test of this class, as it admits a similarity solution , where is the radial distance from the injection well, and is time. This similarity solution holds even when complicated physics such as mutual solubility is included.
An example of a 1D radial injection problem is included in the test suite.
# 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
[]
[water]
type = Water97FluidProperties
[]
[watertab]
type = TabulatedFluidProperties
fp = water
temperature_min = 273.15
temperature_max = 573.15
fluid_property_file = water_fluid_properties.csv
save_file = false
[]
[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/fluidstate/theis_brineco2.i)Initially, all of the CO injected is dissolved in the resident brine, as the amount of CO is less than the dissolved mass fraction at equilibrium. After a while, the brine near the injection well becomes saturated with CO, and a gas phase forms. This gas phase then spreads radially as injection continues. The mass fraction of salt in the brine decreases slightly where CO is present, due to the increased density of the brine following dissolution of CO.
A comparison of the results expressed in terms of the similarity solution is presented in Figure 4 for all three nonlinear variables for the case where is fixed (evaluated using Postprocessors), and also the case where is fixed (evaluated using a VectorPostprocessor).
As these results show, the similarity solution holds for this problem.
An nonisothermal example of the above problem, where cold CO is injected into a warm aquifer, is also provided in the test suite.
# Two phase nonisothermal Theis problem: Flow from single source.
# Constant rate injection 2 kg/s of cold CO2 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
[]
[xnacl]
initial_condition = 0.1
[]
[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
[]
[mass2]
type = PorousFlowMassTimeDerivative
fluid_component = 2
variable = xnacl
[]
[flux2]
type = PorousFlowAdvectiveFlux
fluid_component = 2
variable = xnacl
[]
[energy]
type = PorousFlowEnergyTimeDerivative
variable = temperature
[]
[heatadv]
type = PorousFlowHeatAdvection
variable = temperature
[]
[conduction]
type = PorousFlowHeatConduction
variable = temperature
[]
[]
[UserObjects]
[dictator]
type = PorousFlowDictator
porous_flow_vars = 'pgas zi xnacl temperature'
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]
[co2]
type = CO2FluidProperties
[]
[brine]
type = BrineFluidProperties
[]
[]
[Materials]
[temperature]
type = PorousFlowTemperature
temperature = temperature
[]
[brineco2]
type = PorousFlowFluidState
gas_porepressure = pgas
z = zi
temperature = temperature
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
[]
[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/fluidstate/theis_brineco2_nonisothermal.i)References
- Z. Duan and R. Sun.
An improved model calculating CO$_2$ solubility in pure water and aqueous NaCl solutions from 273 to 533 K and from 0 to 2000 bar.
Chemical geology, 193(3-4):257–271, 2003.[BibTeX]
- J. E. Garcia.
Density of aqueous solutions of CO$_2$.
Technical Report, LBNL, 2001.[BibTeX]
- K. Pruess, C. Oldenburg, and G. Moridis.
TOUGH2 User's Guide, Version 2.0.
Technical Report LBNL-43134, Lawrence Berkeley National Laboratory, Berkeley CA, USA, 1999.[BibTeX]
- N. Spycher and K. Pruess.
A phase-partitioning model for CO$_2$-brine mixtures at elevated temperatures and pressures: application to CO$_2$-enhanced geothermal systems.
Transport in porous media, 82(1):173–196, 2010.[BibTeX]
- N. Spycher, K. Pruess, and J. Ennis-King.
CO$_2$-H$_2$O mixtures in the geological sequestration of CO$_2$. I. Assessment and calculation of mutual solubilities from 12 to 100C and up to 600 bar.
Geochimica et Cosmochimica Acta, 67:3015–3031, 2003.[BibTeX]
- N. Spycher, K. Pruess, and J. Ennis-King.
CO$_2$-H$_2$O mixtures in the geological sequestration of CO$_2$. II. Partitioning in chloride brine at 12-100C and up to 600 bar.
Geochimica et Cosmochimica Acta, 69:3309–3320, 2005.[BibTeX]