**Focus**: Improving staggered iteration between multiphase reaction transport and geomechanics

- Problems we were trying to solve
- Challenges we encountered
- Numerical algorithms
- Looking at poroelasticity as inspiration
- A new preconditioned staggered iteration
- Numerical experiments
- Problems we can now solve

• Motivating problems are hydrate dissociation problems

• Wildly varying fluid compositions and mechanical properties

• Large changes causes drastic changes to the fluid mixture properties

• Properties vary greatly in space and time: can't estimate a priori

• Very soft sediments with soft gases and stiff liquids and ices

- Difficult to marry two physical descriptions
- One code for flow simulation, one code for geomechanics
- Slow convergence of staggered iteration
- Iteration can diverge when pore pressure dominates skeleton in stiffness
- Iteration can fail due to small valid ranges of in multiphase relations

- Ideally, one would develop a fully coupled simulator from the ground up
- That is easier said than done
- Disjoint development between multiphase chemistry and computational mechanics
- The problem is hard, not just writing the software
- Current boat is having two distinct solvers
- Working in TOUGH+ and Millstone
- Multiphase flow is very tricky: make mechanics package do the heavy-lifting
- Fortran calls C calls Python calls C calls Python-generated C

Coupled transport, kinetics, and momentum balances:

Solve for $p,T,X$, and $u$ such that: \begin{eqnarray} \frac{\mathrm{d}m}{\mathrm{d}t} = & \nabla \cdot q + s & \} \text{balance of mass}\\ \frac{\mathrm{d}X}{\mathrm{d}t} = & r & \}\text{reactions (possibly equilibrium)} \\ 0 = & \nabla \cdot \sigma + f & \} \text{balance of momentum (quasistatic)} \end{eqnarray} Plus statemachine logic for phase change: \begin{eqnarray} state = \text{transition}(state; p,X,T) \end{eqnarray}

Pick numerical methods that are best for each set of equations.

- Reaction-transport in reservoirs needs a very fine mesh
- Integral Finite Difference Method with TOUGH+ for reaction transport
- Finite Element Method with Millstone for mechanics
- One-to-one aspect ratio problem for mechanics solved by using two meshes
- Axisymmetry in mechanics reduces unknowns and improves numerical performance

Have an additional interpolation step

Queiruga, Moridis, and Reagan, *"Simulation of Gas Production from Multilayered Hydrate-Bearing Media with Fully Coupled Flow, Thermal, Chemical and Geomechanical Processes Using TOUGH+Millstone. Parts 1-3,"* Transport in Porous Media, 2019

At each timestep, we're solving two nonlinear algebraic problems:

\begin{align} F(p,T,X,...)= & 0\left.\right\} \text{flow code}\\ M(u)= & 0\left.\right\} \text{mechanical code} \end{align}To *really* do Newton's method we would need:
\begin{align}
\left[\begin{array}{cc}
\frac{\partial F}{\partial p} & \frac{\partial F}{\partial u}\\
\frac{\partial M}{\partial p} & \frac{\partial M}{\partial u}
\end{array}\right]\left\{ \begin{array}{c}
\Delta p^{k}\\
\Delta u^{k}
\end{array}\right\} =-\left\{ \begin{array}{c}
F^{k}\\
M^{k}
\end{array}\right\}
\end{align}
E.g. in hydraulic fracturing the fluid pressure must be fully coupled like this or else it won't converge.

Hard to formulate those cross-terms:

- Very different discretizations on different meshes
- Two codes due to seperate development legacies
- Bigger matrix systems, harder to partition and solve in parallel
- Won't necessarily yield better
*performance*

For each Newton loop $k$ in the timestep:

- Solve mechanical stiffness using $K_d$ \begin{equation} \Delta u^{k+1} = [\mathtt{K}(K_d)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \} \end{equation}
- Increment Stress with $K_d$ \begin{equation} \Delta \sigma^{k+1} = \left(\left(K_d+\frac{2}{3}G\right) \nabla \cdot + 2G\nabla^s \right) : \Delta u^{k+1} + \alpha \mathbf{I} \Delta p^{k} \end{equation}
- Solve flow Jacobian with new volume \begin{equation} \Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\} \end{equation}
- Do state transitions \begin{equation} state^{k+1} = \text{transition}(state^{k},p^{k+1}) \end{equation}
- Repeat until $\left|\Delta p\right|<\mathtt{tol\_flow}$ and $\left|\Delta u\right|< \mathtt{tol\_mech}$

Staggering the two solvers is the same as this: \begin{equation} \left[\begin{array}{cc} \frac{\partial F}{\partial p} & 0\\ 0 & \frac{\partial M}{\partial u} \end{array}\right]\left\{ \begin{array}{c} \Delta p^{k}\\ \Delta u^{k} \end{array}\right\} =-\left\{ \begin{array}{c} F^{k}\\ M^{k} \end{array}\right\} \end{equation} Problematic when $\left| \frac{\partial M}{\partial p} \right| > \left| \frac{\partial M}{\partial u}\right|$: fluid stiffer than skeleton.

- Easy to implement.
- More iterations (no quadratic convergence!)
- Potentially a lot of over-shoot during the staggering.
- Exacerbates state transition logic.

This is catastrophic in large-range multiphase flow problems: \begin{equation} F(p_{bad}) = \text{raise out of range error} \end{equation} Physical descriptions have limits! May never be able to solve the system.

Can we find some matrix $C$ that helps? \begin{equation} \left[\begin{array}{cc} \frac{\partial F}{\partial p} & 0\\ 0 & \frac{\partial M}{\partial u}+C \end{array}\right]\left\{ \begin{array}{c} \Delta p^{k}\\ \Delta u^{k} \end{array}\right\} =-\left\{ \begin{array}{c} F^{k}\\ M^{k} \end{array}\right\} \end{equation}

The poroelastic constitutive response for the fluid pore pressure is \begin{equation} \mathrm{d}p=-M\alpha\nabla\cdot\mathrm{d}\mathbf{u}+M\nabla\cdot\left(\mathbf{q}\mathrm{d}t\right) \end{equation} where the poroelastic coefficient is \begin{equation} M^{-1} = \phi K_f^{-1} + (\alpha - \phi) K_s^{-1} \end{equation} where $\mathbf{q}\mathrm{d}t$ is an differential fluid volume into the control volume.

For the total stress, the constitutive response is, \begin{equation} \mathrm{d}\sigma=-\alpha\mathrm{d}p\mathbf{I}+\left(K_{d}-\frac{2}{3}G\right)\nabla\cdot\mathrm{d}\mathbf{u}\mathbf{I}+2G\nabla^{s}\mathrm{d}\mathbf{u}. \end{equation}

**Note:** We use the law for $\sigma$, but not $p$

Have nonlinearity in composition-dependent properties; e.g. in hydrate, $K$ and $G$ are function of saturation $S$: \begin{eqnarray*} K = K_{0}(1-S)+K_{1}S, \quad G = G_{0}(1-S)+G_{1}S \end{eqnarray*}

We can substitute $\mathrm{d}p$ for a familiar form: \begin{equation} \mathrm{d}\sigma=\left(K_{u}-\frac{2}{3}G\right)\nabla\cdot\mathrm{d}\mathbf{u}+2G\nabla^{s}\mathrm{d}\mathbf{u} + \alpha^2 M \nabla\cdot\left(\mathbf{q}\mathrm{d}t\right) \mathbf{I} \end{equation} where the undrained bulk modulus is: \begin{equation} K_{u} = K_d + \alpha^2 M \end{equation}

Making the analogy to our matrix system, we had (*note: more variables than $p$*)
\begin{equation}
\left[\begin{array}{cc}
1 & M\alpha\mathrm{tr}\\
\alpha\mathbf{I} & K_{d}
\end{array}\right]\left\{ \begin{array}{c}
\mathrm{d}p\\
\mathrm{d}u
\end{array}\right\} =-\left\{ \begin{array}{c}
f\\
m
\end{array}\right\}
\end{equation}
and are trying to do:
\begin{equation}
\left[\begin{array}{cc}
1 & M\alpha\mathrm{tr}\\
0 & K_{d}+\alpha^{2}M
\end{array}\right]\left\{ \begin{array}{c}
\mathrm{d}p\\
\mathrm{d}u
\end{array}\right\} =-\left\{ \begin{array}{c}
f\\
m - \alpha f \mathbf{I}
\end{array}\right\}
\end{equation}
Triangularized the system so solving $\mathrm{d}u$ first would substitute into $\mathrm{d}p$ equation closer to the answer.

**Caveat:** Poroelasticity only holds in the linear regime! (Small changes)

Actual balance of mass that the reservoir simulator solves is:
\begin{equation}
\frac{\partial}{\partial t} \rho = \nabla \cdot \mathbf{q}
\end{equation}
where we assumed:
\begin{equation}
\rho = \rho_0 ( 1+ p / K_f )
\end{equation}
Flow simulators use **very complicated nonanalytic formulations** for $\rho$.

- $K_f$ and $M$ varies greatly especially at phase changes
- Don't necessarily have relations for $K_f$ for every possible scenario
- Varies spatially as well as temporally

```
def gibbs_region1(T,p):
p1_star = 1.653e7
T1_star = 1.386e3
n1 = [ 0.14632971213167e00, -0.84548187169114e00,
-3.7563603672040e+00, 3.3855169168385e+00,
-0.95791963387872e00, 0.15772038513228e00,
-1.6616417199501e-02, 8.1214629983568e-04,
2.8319080123804e-04, -6.0706301565874e-04,
-1.8990068218419e-02, -3.2529748770505e-02,
-2.1841717175414e-02, -5.2838357969930e-05,
-4.7184321073267e-04, -3.0001780793026e-04,
4.7661393906987e-05, -4.4141845330846e-06,
-7.2694996297594e-16, -3.1679644845054e-05,
-2.8270797985312e-06, -8.5205128120103e-10,
-2.2425281908000e-06, -6.5171222895601e-07,
-1.4341729937924e-13, -4.0516996860117e-07,
-1.2734301741641e-09, -1.7424871230634e-10,
-6.8762131295531e-19, 1.4478307828521e-20,
2.6335781662795e-23, -1.1947622640071e-23,
1.8228094581404e-24, -9.3537087292458e-26]
i1 = [ 0, 0,0, 0, 0, 0, 0, 0, 1,1, 1, 1,1, 1, 2, 2, 2, 2, 2, 3, 3, 3,
4, 4,4, 5, 8, 8, 21, 23, 29, 30, 31, 32 ]
j1 = [ -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0,1, 3, -3, 0, 1,3, 17, -4,
0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41 ]
p_i = p/p1_star
t_i = T1_star/T
return R*T*sum([ n*(7.1-p_i)**I*(t_i-1.222)**J
for n,I,J in zip(n1,i1,j1)])
density_region1,enthalpy_region1 = density_enthalpy(gibbs_region1)
```

Solving balance with fluxes: \begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}M_{comp}=\sum_{facets}F_{comp} \end{equation} Masses and fluxes are integrals over discrete domains and facets: \begin{equation} M_{comp}=\int_{\Omega}\phi\sum_{phase}S^{phase}\rho_{comp}^{phase}X_{comp}^{phase}\mathrm{d}^{3}x \end{equation} and: \begin{equation} F_{comp}=\int_{\Gamma}\sum_{phase}\mathbf{q}_{comp}^{phase}\cdot\mathbf{n}\mathrm{d}^{2}x \end{equation}

Flow and mechanics are coupled through volume of $\left|\Omega\right| = V(t)$ \begin{equation} \frac{\mathrm{d}V}{\mathrm{d}t} = V(t)\frac{\mathrm{d}}{\mathrm{d}t} \mathrm{tr} \varepsilon \end{equation} Poroelastic formulation uses the product rule to extract this term: \begin{equation} \frac{\mathrm{d}\rho V}{\mathrm{d}t}= V \frac{\partial \rho}{\partial p}\frac{\mathrm{d}p}{\mathrm{d}t} + \rho \frac{\mathrm{d}V}{\mathrm{d}t} \end{equation} (Derive from transformation of $\nabla_x\cdot$ to $\nabla_X\cdot$ given $x(t)=X+u(t)$)

For each Newton loop $k$ in the timestep:

**Solve mechanical stiffness using**$K_u^*$ \begin{equation} \Delta u^{k+1} = [\mathtt{K}(K_d+\alpha^2 M)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \} \end{equation}- Increment Stress with $K_d$ \begin{equation} \Delta \sigma^{k+1} = \left(\left(K_d+\frac{2}{3}G\right) \nabla \cdot + 2G\nabla^s \right) : \Delta u^{k+1} + \alpha \mathbf{I} \Delta p^{k} \end{equation}
- Solve flow Jacobian with new volume \begin{equation} \Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\} \end{equation}
- Do state transitions \begin{equation} state^{k+1} = \text{transition}(state^{k},p^{k+1}) \end{equation}
- Repeat until $\left|\Delta p\right|<\mathtt{tol\_flow}$ and $\left|\Delta u\right|< \mathtt{tol\_mech}$

Use prior knowledge of how the flow solver will behave (even though the physical description is different) to precondition the iteration.

At each update, the flow simulator yields a $\Delta p$ based on the current composition of the fluid.

Assume we're in the linear regime within the iteration: \begin{equation} \Delta p=-M\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right) \end{equation} We can estimate the effective $M$: \begin{equation} M^{*}=-\frac{\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right)}{\Delta p} \end{equation}

We can back out an estimate of the bulk modulus of the fluid: \begin{equation} K_{f}^{*}=\phi\left(\frac{1}{M^{*}}-\frac{\left(\alpha-\phi\right)}{K_{s}}\right)^{-1} \end{equation}

Assemble the mechanical stiffness matrix using $K_{u}$ in place of $K_d$ with estimated $M^*$: \begin{equation} K_{u}^{*}=K_{d}+\alpha^{2}M^{*} \end{equation}

**Watch out:** Not-a-number when $\Delta p=0$, e.g. at boundary conditions.

Estimate PDE $\nabla\cdot q$ from fluxes: \begin{eqnarray} \nabla\cdot\mathbf{q}=&\frac{1}\\{V}\int_{\Omega}\nabla\cdot\mathbf{q}\mathrm{d}^{3}x\\=&\frac{1}\\{V}\oint_{\partial\Omega}\mathbf{q}\cdot\mathbf{n}\mathrm{d}^{2}x\\=&\frac{1}{V}\sum_{facets}\int_{\Gamma}\mathbf{q}\cdot\mathbf{n}\mathrm{d}^{2}x\\=&\frac{1}{V}\sum_{facets}F \end{eqnarray}

Only modification needed to the flow solver; Millstone takes care of the rest.

1D quasi-static consolidation with a ramp-up load filled with a mixture of $H_2O$ and $CH_4$

- Verified against analytical solution for static poroelasticity
- Should get it exactly correct (if it's slow, stays linear, and we know $K_f$)

- Family of problems varied by:
- $K_ d$ skeletal stiffness
- $S_ {aqu}$ Water to gas ratio

- Algorithm is permuted by:
- Choice: None, Fixed, Adaptive, Continuous
- Tolerance on convergence

- Algorithm is evaluated by:
- Does it succeed
- Mean number of iterations for each set up
- Total runtime (too short)

Each data point is one simulation; missing point reflects divergence.

- Original program with no preconditioning fails when $K_d < K_f$
- Fixed precondition solves that problem when $K_f$ was known
- Adaptation of $K_f^*$ upon convergence improves it
- Adapting during the time step iteration is robust for all compositions with reasonable tolerance
- Adapting can perform worse with gaseous pore fluids
- Algorithm estimates correct $K_f$ on this problem (not smooth in practice)

Algorithms work on 1D test problem, but still not 100\% robust in some reservoir problems.

Need to add under-relaxation: \begin{equation} \Delta u = [\mathtt{K}(K_d+\alpha^2 M^*)] \backslash \{\mathtt{R}(\sigma)\} \end{equation} And update by \begin{equation} u^{k+1} = u^k + \beta \Delta u \end{equation} With $0<\beta\leq 1$. Smaller $\beta$ slows down convergence.

Also, watch out for pore ``collapse'' when gases appear ($K_f \approx 50Pa$ vs $K_s \approx 10^{10} Pa$)

- Assumptions out of range for poroelastic porosity evolution?

All this work to get this type of problem to work:

- Hit divergence on coupled problems with drastic fluid compositional changes
- Appearance of gases makes staggered iteration raise
- Used poroelastic theory to design a way to adaptively precondition
- Permuted idea to make multiple staggering algorithms
- Tried all permutations on a simple constructed test problem

Comments:

- Look for places where our computer program is arbitrary and permute them.
- Design families of tiny test problems to program their limits

These sldies are hosted at: ocf.io/afq/slides/Siam2019

Queiruga, Moridis, and Reagan, "Simulation of Gas Production from Multilayered Hydrate-Bearing Media with Fully Coupled Flow, Thermal, Chemical and Geomechanical Processes Using TOUGH+Millstone. Parts 1-3," Transport in Porous Media, 2019

- Convergence difficult when $K_d \lessapprox K_f$ with naive approach
- Use unstructured quad meshes for mechanics
- True 2D axisymmetry for mechanics
- Continuous $M$ estimation for preconditioner
- Under relax $\Delta u$ by $\beta \approx 0.2-0.5$
- Timestep control: double below 25 and cutback at 100
- A lot of time per iteration, but stable on large timesteps.

In [ ]:

```
```