Can MOOSE Handle Salt Precipitation and Its Induced Porosity and Permeability Changes?

[robo-boluo]

Dear MOOSE Team,

I am writing to inquire about the capabilities of MOOSE regarding NaCl precipitation during CO₂ injection into a brine aquifer. Specifically, I am interested in understanding whether the current version of MOOSE can model the effects of salt precipitation on porosity and permeability changes.

If this functionality is not currently available, could you provide guidance on how to implement this feature?

Thank you for your assistance.

[GiudGiud]

Hello

Porosity and permeability are both material properties which can be set dynamically during the simulation as a function of NaCl precipitate concentration.
So yes, you could model this. Implementing new material properties should be (relatively) straightforward after going through the tutorial

Guillaume

[robo-boluo]

Hi GiudGiud, @GiudGiud @cpgr

I want to express my gratitude for your valuable advice. In response to it, I created new source files 'PorousFlowPorosity_nacl.h' and 'PorousFlowPorosity_nacl.C' (link)to incorporate salt precipitation.

However, upon implementing PorousFlowPorosity_nacl in the [Materials] block, I encountered an issue during testing. It appears that '(PorousFlow_density_qp, dPorousFlow_density_qp_dvar, PorousFlow_xnacl_qp, dPorousFlow_xnacl_qp_dvar, PorousFlow_density_nodal, dPorousFlow_density_nodal_dvar, PorousFlow_xnacl_nodal, dPorousFlow_xnacl_nodal_dvar)' are not defined on the specified blocks (caps and aquifer).

I called the PorousFlowPorosity_nacl using the following format:

   [porosity]
    type = PorousFlowPorosity_nacl
    block=aquifer
    fluid = true
    mechanical = true
    thermal = true
    salt_precip = true
    porosity_zero = 0.2
    reference_temperature = 105
    reference_porepressure = 15.5E6
    thermal_expansion_coeff = 15E-6 # volumetric
    solid_bulk = 8E9 # unimportant since biot = 1
    xnacl_L = 0.30
    rho_nacl = 2165.0
    r_m = 1.0e-5
    xnacl=xnacl
    density=density_h2o
  []

Could you please provide suggestions on how to address this issue at your convenience?

Thank you very much!

[robo-boluo]

Thank you for your prompt response. I have just granted access. Please let me know if you encounter any issues accessing it.

[GiudGiud]

I just took a look. I see the porosity is modified depending on where x_nacl sits compared to xnacl_L.
However I dont see a kernel or a constraint preventing xnacl from going past xnacl_L ?

  [mass2]
    type = PorousFlowMassTimeDerivative
    fluid_component = 2
    variable = xnacl
  []
  [flux2]
    type = PorousFlowAdvectiveFlux
    fluid_component = 2
    variable = xnacl
    use_displaced_mesh = false
  []
  [diff2]
    type = PorousFlowDispersiveFlux
    fluid_component = 2
    variable = xnacl
    disp_long = '0 0'
    disp_trans = '0 0'
  []
  [vol_strain_rate2]
    type = PorousFlowMassVolumetricExpansion
    fluid_component = 2
    variable = xnacl
  []

[robo-boluo]

That is correct. I have questions about implementing an algorithm to prevent 'xnacl' from exceeding 'xnal_L'. The only modification in my input.i file is the porosity in the [Materials] block. Do you think the NaCl mass fraction constraint should be defined in the source code (.h and .C files) or in the input.i file?

[GiudGiud]

Constraints can be defined in the input file.
I think this would be a good use for the Bounds system. See this object to add to your simulation:
https://mooseframework.inl.gov/source/bounds/ConstantBounds.html

You also enforce it by writing a custom force kernel which would gradually penalize the residual when xnacl goes out of bounds.

[robo-boluo]

Thanks, I will update to you later.

[robo-boluo]

Hi Guillaume Giudicelli, @GiudGiud

This is a quick update of the addition of ConstantBounds in the input file (link). Two simulation results are observed.

(1) Regarding the mass fraction shuffle

At Timestep 31 in Fig.1, the injection inlet reaches the maximum mass fraction of (xnacl_L= 0.3). Subsequently, there is a mass fraction shuffle in the next timestep (as seen in Fig.2). Following this, the NaCl propagates further into the formation, maintaining a maximum mass concentration of 0.3 (as depicted in Fig.3 and Fig.4). The mass fraction shuffle refers to the redistribution or rearrangement of NaCl within the system. How to intepret the mass fraction shuffle?

Fig.1

Fig.2

Fig.3

Fig.4

(2) Regarding the lack of porosity changes in the reservoir layer

Despite reaching its limit (xnacl_L), no porosity reduction is observed in the reservoir layer (see Fig.5)

Fig.5

[robo-boluo]

Hello Guillaume Giudicelli, @GiudGiud ,

This is quick update after coding the xnacl update equation (𝑥𝑛𝑎𝑐𝑙_𝑛𝑒𝑤=𝑥𝑛𝑎𝑐𝑙_𝐿=0.3) for [Kernals]:

Code for XNaclUpdate.h:

#pragma once

#include "Kernel.h"

class XNaclUpdate : public Kernel
{
public:
  static InputParameters validParams();
  XNaclUpdate(const InputParameters & parameters);

protected:
  virtual Real computeQpResidual() override;
  virtual Real computeQpJacobian() override;
private:
  const VariableValue & _xnacl;
  const Real _xnacl_L;
};

Code for XNaclUpdate.C:

#include "XNaclUpdate.h"
#include "MooseVariable.h"

registerMooseObject("PorousFlowApp", XNaclUpdate);

InputParameters
XNaclUpdate::validParams()
{
  InputParameters params = Kernel::validParams();
  params.addCoupledVar("xnacl", "The salt mass fraction in the brine (kg/kg)");
  params.addParam<Real>("xnacl_L", "nacl solubility limit at specific temp and pressure");
  return params;
}

XNaclUpdate::XNaclUpdate(const InputParameters & parameters) :
  Kernel(parameters),
  _xnacl(coupledValue("xnacl")),
  _xnacl_L(getParam<Real>("xnacl_L"))
{
}

Real
XNaclUpdate::computeQpResidual()
{
  // Compute the residual contribution
  Real residual = 0.0;
  if (_xnacl[_qp] > _xnacl_L)
  {
    residual = (_xnacl_L - _xnacl[_qp]) * _test[_i][_qp];
  }
  return residual;
}

Real
XNaclUpdate::computeQpJacobian()
{
  // Compute the Jacobian contribution
  Real jacobian = 0.0;
  if (_xnacl[_qp] > _xnacl_L)
  {
    jacobian = -_phi[_j][_qp] * _test[_i][_qp];
  }
  return jacobian;
}

Resutls:

(1) There are on changes for porosity

(2) xnacl can exceed the solubility limit xnacl_L

(3) Not converge after xnacl>xnacl_L, as timestep <e-4

Do you think what else can we do to further improve the simulation results?

Appreciate your time and guidance.

[GiudGiud]

1. that's expected, porosity is averaged and xnacl is on average less than 0.3 still
2. it can exceeed because the rate of removal is 1 * (xnacl- limit). Maybe there should be a larger rate? especially if physically it can't exceed the limit. Add a coefficient
3. This is more troubling. Maybe use -pc_type lu to see if you can get it to converge

[robo-boluo]

Hi Guillaume Giudicelli, @GiudGiud

A coefficient of 2 has been added to the 'XNaclUpdate' kernel, and the '-pc_type lu' has already been included in the [Preconditioning] block of the previous input file. However, there is no difference compared to the previous case with the coefficient =1; the three issues mentioned above still persist.

Regarding the comment that an average xnacl < 0.3 causes no porosity change, I removed the XNaclUpdate Kernel and allowed the xnacl to increase without restriction. However, this example still yields non-changeable reservoir porosity. Therefore, I suspect that the implementation of 'PorousFlowPorosity_nacl' is incorrect.

Figure 1. The xnacl profile without maximum solubility restriction still yields non-changeable reservoir porosity

[GiudGiud]

What steps have you taken to verify the implementation of PorousFlowPorosity_nacl ?
You dont need to run a full simulation to check it by the way.

You can create a dummy input with a Xnacl auxiliary variable. You can see this variable to an arbirtary value using a ParsedAux auxkernel

[robo-boluo]

Hi Guillaume Giudicelli @GiudGiud ,

Thank you for your continuous support.

1. PorousFlowPorosity_nacl Verification
I developed the 'PorousFlowPorosity_nacl' model based on the existing 'PorousFlowPorosity' model in the MOOSE folder (/moose/modules/porous_flow/). I implemented the salt precipitation effect by following the format of the chemical part in the PorousFlowPorosity.C file, as both mechanisms consider solid volume change. Consequently, I developed the code for salt precipitation under the 'atNegInfinityQp', 'datNegInfinityQp', 'atZeroQp', and 'datZeroQp' sections.
The basic mathematical equation to describe the salt precipitation is:

𝜙𝑛+1=𝜙𝑛−𝜙𝑛⋅𝑆𝑤⋅𝑑𝑡⋅𝜌𝑤(𝑥𝑛𝑎𝑐𝑙-𝑥𝑛𝑎𝑐𝑙_𝐿)⋅𝑟𝑚⋅𝜌𝑚

where 𝑆𝑤​ and 𝜌𝑤 are the aqueous saturation and density, 𝑟𝑚 is the precipitation rate, and 𝜌𝑚 is the halite density.

2. Dummy Test

I tested the following dummy input file, and the code is executable.

[Mesh]
  type = GeneratedMesh 
  dim = 2
  nx = 10
  ny = 10
  xmin = 0
  xmax = 1
  ymin = 0
  ymax = 1
[]

[GlobalParams]
  PorousFlowDictator = dictator
[]

[Variables]
  [pwater]
    initial_condition = 1e6
  []
  [xnacl]
    initial_condition = 0.1
  []
  [zi]
    initial_condition = 0
  []
[]

[AuxVariables]
  [aux_xnacl]
    order = CONSTANT
    family = MONOMIAL
  []
  [swater]
    order = CONSTANT
    family = MONOMIAL
  []
[]

[Kernels]
  [dummy_pwater]
    type = NullKernel
    variable = pwater
    block = 0
  []
  [dummy_xnacl]
    type = NullKernel
    variable = xnacl
    block = 0
  []
  [dummy_zi]
    type = NullKernel
    variable = zi
    block = 0
  []
[]

[AuxKernels]
  [set_xnacl]
    type = ParsedAux
    variable = aux_xnacl
    expression = 0.35
  []
  [swater]
    type = PorousFlowPropertyAux
    variable = swater
    property = saturation
    phase = 0
    execute_on = timestep_end
  []

[]

[UserObjects]
  [dictator]
    type = PorousFlowDictator
    porous_flow_vars = 'pwater zi xnacl'
    number_fluid_phases = 2
    number_fluid_components = 2
  []
  [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 = 20
[]

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

  [porosity_test]
    type = PorousFlowPorosity_nacl
    block = 0
    fluid = true
    mechanical = false
    thermal = false
    salt_precip = true
    porosity_zero = 0.2 
    reference_porepressure = 1e6
    solid_bulk = 8E9 # unimportant since biot = 1
    xnacl_L = 0.3
    rho_nacl = 2165.0
    r_m = 1.0e-1
    xnacl = aux_xnacl
  []

  [stress_pp_effective]
    type = PorousFlowEffectiveFluidPressure
  []
[]

[Executioner]
  type = Transient
  solve_type = NEWTON
  nl_abs_tol = 1e-7
  nl_rel_tol = 1e-5
  l_tol = 1e-5
  l_max_its = 50
  end_time = 100 
  [TimeStepper]
    type = IterationAdaptiveDT
    dt = 0.5
    growth_factor = 1.5
  []
[]

[Outputs]
  exodus = true
  csv = true
[]

[BCs]
  [left]
    type = DirichletBC
    variable = pwater
    boundary = left
    value = 1.2e6
  []
  [right]
    type = DirichletBC
    variable = pwater
    boundary = right
    value = 1e6
  []
  [left_gas]
    type = DirichletBC
    variable = zi
    boundary = left
    value = 0.1
  []
  [right_gas]
    type = DirichletBC
    variable = zi
    boundary = right
    value = 0
  []
[]

3. Further Improvement Suggestions

Which areas of the current implementation of 'PorousFlowPorosity_nacl' need improvement?

[robo-boluo]

I noticed that both Guillaume Giudicelli and Chris Green (@GiudGiud , @cpgr) contributed to the salt precipitation topic in thread #23701. The final collective conclusion was to modify the porosity and permeability file. Therefore, I believe I am heading in the right direction. I look forward to any updates on this topic.

Many thanks.

[robo-boluo]

Hi Guillaume Giudicelli @GiudGiud,

I have been working on this problem for quite some time and would like to express my sincere appreciation for your continuous attention and guidance during this period. Initially, I was excited to confirm with you that MOOSE has the capability to model the effect of salt precipitation on porosity and permeability by adjusting the material properties. However, developing the App has proven to be much more complicated. Here is a summary of our progress over the past month:

1. Taking 𝑥NaCl (NaCl mass fraction) as a variable, I developed the 'PorousFlowPorosity_nacl' App based on the existing 'PorousFlowPorosity' App in the MOOSE folder (/moose/modules/porous_flow/). Since the salt precipitation and chemical reaction parts have similar processes affecting porosity, I programmed the salt precipitation by following the format of the chemical part, defining the 'atNegInfinityQp', 'datNegInfinityQp', 'atZeroQp', and 'datZeroQp' sections for my 'PorousFlowPorosity_nacl' app.

2. In the input file, the variable 𝑥NaCl is computed using type=PorousFlowMassTimeDerivative, type=PorousFlowAdvectiveFlux, and type=PorousFlowDispersiveFlux in [Kernels].

3. I defined 𝑥NaCl_max as the maximum NaCl solubility in the aquifer. Once 𝑥NaCl>𝑥NaCl_max, halite starts to precipitate. To avoid 𝑥NaCl>𝑥NaCl_max, we applied the ConstantBounds to restrict 𝑥NaCl within [0,𝑥NaCl_max]. However, the porosity change was not observed in the simulation results. "Failed"

4. We suspected that the ConstantBounds restriction on 𝑥NaCl≤𝑥NaCl_max prevented precipitation from occurring. Therefore, we deleted the ConstantBounds. After precipitation, we should update 𝑥NaCl. Consequently, the ‘XNaclUpdate’ app (see above) was developed and applied in [Kernels]. However, the expected porosity change was still not achieved. "Failed"

5. Even though I increased the 𝑥NaCl removal rate (rate * (𝑥NaCl−𝑥NaCl_max)), it still did not result in any porosity change. We also observed 𝑥NaCl>𝑥NaCl_max and encountered convergence problems at 𝑥NaCl>𝑥NaCl_max. "Failed"

Given your experience in MOOSE development, I trust your assertion that MOOSE can handle salt precipitation-induced porosity changes. Therefore, the issue may lie in my implementation of the 'PorousFlowPorosity_nacl' App. I would sincerely appreciate it if you could review my problem and provide your insights (code link).

Thank you for your time and attention to this matter.

[cpgr]

Hi @robo-boluo, sorry, I haven't made much progress due to other commitments. I'll try and take a look this afternoon and get back to you.

[cpgr]

Just as a reminder (like @GiudGiud mentioned), do you want to consider the case where water evaporates into the CO2 gas phase, the salt mass fraction increases until it reaches the solubility limit, then the excess salt precipitates, reducing porosity (from which permeability reduces)?

[robo-boluo]

Hi Chris Green @cpgr ,

I appreciate your time and kind help. Yes, this is exactly what I want to implement. There should be mutual solubility between water and CO2. The reduction in porosity and permeability is caused by NaCl precipitation.

[robo-boluo]

I'm not a porous flow module developer From a general MOOSE developer point of view, I think this should be doable
Remind me of the behavior with the fraction variable? It s supposed to be physically bound? Should it have its own variable?

The mechanism process of salt precipitation is clear to me. In my code, xnacl represents the mass fraction variable in Kernel block.

[GiudGiud]

the salt mass fraction increases until it reaches the solubility limit

how does this happen in the equations? When the solubility limit is reached, does some sink term sky-rocket or something? Or a property changes to prevent xnacl to keep going?