# Phase-field modeling and effective simulation of non-isothermal reactive
transport

Carina Bringedal111Department of Computer science, Electrical engineering and
Mathematical sciences, Western Norway University of Applied Sciences, and
Institute for Modelling Hydraulic and Environmental Systems, University of
Stuttgart<EMAIL_ADDRESS>Alexander Jaust222Institute for Parallel
and Distributed Systems, University of Stuttgart

###### Abstract

We consider single-phase flow with solute transport where ions in the fluid
can precipitate and form a mineral, and where the mineral can dissolve and
release solute into the fluid. Such a setting includes an evolving interface
between fluid and mineral. We approximate the evolving interface with a
diffuse interface, which is modeled with an Allen-Cahn equation. We also
include effects from temperature such that the reaction rate can depend on
temperature, and allow heat conduction through fluid and mineral. As Allen-
Cahn is generally not conservative due to curvature-driven motion, we include
a reformulation that is conservative. This reformulation includes a non-local
term which makes the use of standard Newton iterations for solving the
resulting non-linear system of equations very slow. We instead apply L-scheme
iterations, which can be proven to converge for any starting guess, although
giving only linear convergence. The three coupled equations for diffuse
interface, solute transport and heat transport are solved via an iterative
coupling scheme. This allows the three equations to be solved more efficiently
compared to a monolithic scheme, and only few iterations are needed for high
accuracy. Through numerical experiments we highlight the usefulness and
efficiency of the suggested numerical scheme and the applicability of the
resulting model.

##  1 Introduction

Reactive transport is relevant for many environmental and technical
applications. We are here especially interested in reactive transport with
heterogeneous reactions, where solutes can leave a liquid phase and become
part of a solid phase, or the other way around. This is a typical feature of
mineral precipitation and dissolution [39], which is common in groundwater
remediation, soil salinization, oil production, metallurgy and geothermal
energy production, to mention some. When temperature plays a role, which is in
particular the case with geothermal energy production, the interplay between
temperature-dependent solubilities of minerals and reaction rates depending on
temperature can give an interplay between varying temperature and solute
concentrations [14]. For many of these applications, experiments are
unfeasible, which rather turns the attention to having reliable models and
numerical simulations to investigate these processes. However, the close
coupling between transport of solutes, heat and chemical reactions can be
challenging to model and simulate[44].

At the interface between fluid and mineral, solutes in the fluid can leave the
liquid fluid phase and rather become part of the solid mineral phase. The
transition of solutes from fluid to mineral corresponds to an evolving fluid-
mineral interface and hence a free boundary of the model. Hence, where the
fluid is and where the mineral is, pose additional challenges to modeling of
such reactive transport processes.

There are several modeling approaches to describe these reactive transport
processes. One possibility is to use sharp-interface models, which includes to
track the interface between the fluid and mineral and hence the fluid and
mineral domains themselves, and at the same time solve the necessary
conservation equations for mass and energy in the respective, time-evolving
domains. The two domains are coupled to each other through boundary conditions
on the evolving fluid-mineral interface, where these boundary conditions
account for the exchange of mass and energy between the domains. To model the
evolution of the interface one can in the general case use a level-set
equation [30], for applications of level-set models with reactive transport
see e.g. [5, 34, 26]. For simpler geometries, e.g. layering or circular
shapes, a simple distance or radius function can be used [7, 41, 42]. Although
correct and accurate, these sharp-interface approaches lead to singularities
in the formulation in case of merging or break-up of two minerals, and can be
numerically challenging due to the need of tracking the interface and divide
the computational domains accordingly. Also, the interface itself appears as a
shock in the domain, with discontinuities across it, leading to further
numerical difficulties. We instead apply a diffuse-interface approach through
a phase-field model.

In a diffuse-interface model [2], the interface between fluid and mineral is
modeled as a diffuse transition zone. Instead of having a sharp interface with
boundary conditions and different model equations and material properties at
each side the evolving interface, the model is rewritten to use the same model
equations at the full domain, where a diffuse transition region represents the
processes at the evolving interface. This diffuse transition region hence
represents the fluid-mineral interface and should account for the exchange of
mass and energy occurring here, but is modeled as smooth. Hence, a diffuse-
interface model does not involve any shocks nor discontinuities as the
transition between fluid and mineral is smoothed out over some region of non-
zero width. The evolution of the fluid-mineral interface is modeled by an
evolution of the location of the diffuse transition region, through a phase-
field equation. However, note that the introduction of the diffuse transition
region is a mathematical simplification and does not represent physics on the
molecular level (as the real interface width would be at a much smaller scale
[21]). Hence, the diffuse-interface model should be seen as an approximation
to the real sharp-interface physics, introduced for convenient mathematical
and numerical treatment. If the diffuse-interface model correctly approximates
the sharp-interface physics; that is, the diffuse-interface model can recover
the sharp-interface model as the diffuse interface width approaches zero, we
consider the diffuse-interface model a valid approximation. This can be shown
analytically by matched asymptotic expansions [10], or compared numerically
for test cases [22, 45, 8]. The diffuse-interface model offers advantages of
avoiding singularities in case of merging and break-up, and is generally
easier to discretize and solve numerically as it is formulated on a fixed
domain without any shocks and discontinuities due to the smooth, diffuse
nature. However, these advantages come at the cost of applying a model that is
an approximation of the sharp-interface physics, which can lead to unwanted
evolution of the diffuse transition region.

There are two main approaches for phase-field models describing evolving
interfaces: Allen-Cahn [1] and Cahn-Hilliard [11]. Where the Allen-Cahn model
is represented by a second-order non-linear partial differential equation
(PDE), the Cahn-Hilliard model is a fourth-order non-linear PDE. Due to being
second-order, the Allen-Cahn equation can be numerically discretized with most
standard discretization schemes, and fulfills a maximum principle. However,
the Allen-Cahn model includes curvature-driven motion of the interface [1],
which in its original form is not conservative. That means, an Allen-Cahn
phase field would evolve in a manner that corresponds to loss or gain of
mineral that cannot be connected to mineral precipitation or dissolution
reactions. On the other hand, the Cahn-Hilliard equation needs special care
due to being fourth order and can have overshoot/undershoot of values, which
need to be handled to stay within the physical bounds [17]. The Cahn-Hilliard
model includes Mullins-Sekerka motion [40], which corresponds to a
redistribution of the mineral to minimize surface area. However, this
redistribution is conservative, which means that any net loss or gain of
mineral is always connected to a mineral precipitation or dissolution
reaction. Both approaches have been successfully used to model mineral
precipitation and dissolution, coupled with fluid flow and solute transport
[8, 36].

We here apply the Allen-Cahn equation due to its beneficial numerical
properties and availability of maximum principle, and rather explore a
reformulation of it to make the curvature-driven motion conservative. We build
on the Allen-Cahn model found in [8]. Rubinstein and Sternberg [37] suggested
to add a non-local term to overcome the issue of non-conservative motion.
Then, the curvature-driven motion instead redistributes the mineral to obtain
mean curvature [13]. This reformulation has been successfully applied in
several contexts [13, 6, 31, 29]. Although conservative, the non-local term
introduces numerical challenges in terms of solving the resulting non-linear
system of equations if implicit time stepping is used. If applying Newton’s
method, the non-local term causes the Jacobian matrix to be full, which makes
the Newton iterations expensive [6].

To avoid expensive Newton iterations caused by the non-local term of the
conservative Allen-Cahn equation, we instead apply L-scheme iterations to
solve the resulting non-linear system of equations [33]. Using L-scheme
iterations relies on having a monotonically decreasing (or increasing) non-
linearity, which is generally not the case for the conservative Allen-Cahn
equation. Therefore we split the non-linearity in a decreasing and increasing
part, which is solved implicitly and explicitly, respectively. By adding a
stabilization parameter, we can guarantee that the L-scheme iterations form a
contraction and always converge [24]. The L-scheme iterations are found to
converge when the stabilization parameter is large enough, and converge
independently of the time-step size and initial guess. L-scheme iterations
however only offer linear convergence, and can require many iterations
depending on the choice of stabilization parameter and time-step size. For
other equations approaches with using modified L-scheme iterations or
switching to Newton when a good starting point is found, have been suggested
[3, 28, 27]. However, it is important to emphasize that for L-scheme
iterations, every iteration step is much cheaper compared to Newton
iterations.

To discretize the model equations, we apply the Finite Volume method [15].
Finite Volume method is discretely conservative and the conservative Allen-
Cahn equation remains discretely conservative when a suitable splitting of the
non-linearity into implicit and explicit parts are made [6].

To be able to model non-isothermal reactive transport, the Allen-Cahn equation
must be coupled to conservation equations for mass and energy. As discussed
above, these conservation equations are in this case defined on fixed domains
and hence describe the conservation of mass and energy in the combined domain
of both fluid and mineral. Hence, the conservation equations are modified to
incorporate the diffuse transition that represents the mass and energy
transfer between fluid and mineral. For this, we build on the isothermal model
in [8] and extend it to non-isothermal conditions. We show how the
conservation property of these equations are modified by coupling to the
(conservative) Allen-Cahn equation.

To solve the coupled system of conservative Allen-Cahn equation, conservation
of mass and energy, we apply an iterative strategy where the different
equations are solved sequentially at each time step, and iterate until
convergence is achieved. This is done for a simpler case without fluid flow
only. The advantage of using an iterative strategy is that each equation can
be discretized according to its needs. In particular, as the conservation
equations for mass and energy are linear (when the phase field is known),
simpler strategies can be used to solve them. By introducing an additional
stabilization parameter for the coupling, we can prove the convergence of the
scheme with a restriction on the time-step size. The strategy and convergence
proof build on the ideas in [3, 9]. Numerical simulations show that the
resulting scheme performs well, requiring generally few coupling iterations to
reach convergence. For larger time-step sizes, the coupling iterations still
converge but the efficiency of the scheme is limited by the L-scheme
iterations then requiring many iterations to reach convergence.

The goal of this paper is hence to formulate a conservative phase-field model
including an evolving fluid-mineral interface that can model fluid flow,
solute transport and heat transport. Given the numerical challenges coming
from such a highly non-linear and strongly coupled model, we aim to formulate
an efficient numerical scheme based on L-scheme iterations for the non-
linearities and coupling iterations for the coupled system of equations.

To reach our goals, the remaining of this paper is structured as following. In
Section 2 we consider the sharp-interface model that form the basis of the
physical processes that we want to model. The corresponding diffuse-interface
model is formulated and analyzed in terms of conservation in Section 3. In
Section 4 we show that the sharp-interface limit of the diffuse-interface
model is indeed the model from Section 2, with the additional curvature-driven
motion. In Section 5 we formulate the numerical strategy and show the discrete
conservation as well as prove convergence of the two iterative schemes.
Numerical behavior of the Allen-Cahn equation and the coupled phase-field
model is investigated in Section 6. Discussion and conclusion can be found in
Section 7.

##  2 Sharp-interface model

We here present the sharp-interface model for fluid flow, solute conservation
and energy conservation in the fluid and mineral. The model is similar to the
model considered in [5, 42]. Two domains are included: The fluid domain where
fluid flow and the dissolved solute is found, and the mineral domain where the
solute appears as a stationary mineral. Both domains are evolving in time as
mineral dissolves becoming solute part of the fluid, and as solute from the
fluid precipitates becoming part of the mineral. We will denote the fluid and
mineral domains $\Omega_{f}(t)$ and $\Omega_{m}(t)$. The two domains are
separated by the interface $\Gamma(t)$. The considered chemical reaction can
be written as

$n_{1}C_{1}+n_{2}C_{2}\rightleftharpoons C;$

that is, that $n_{1}$ ions of type $C_{1}$ go together with $n_{2}$ ions of
type $C_{2}$ to form one minerals molecule of type $C$. The reaction can go
both ways, depending on the ion concentrations relative to the minerals
solubility.

###  2.1 Fluid domain

In the fluid domain, we consider a fluid having density $\rho_{f}$ and which
is flowing with velocity $\mathbf{v}$. Dissolved in the fluid, two solutes
having concentration $c_{1}$ and $c_{2}$ (corresponding to the two ion types
$C_{1}$ and $C_{2}$, respectively) are transported. Temperature variations are
accounted for, and the temperature of the fluid is denoted $T_{f}$. Since we
will work with conservation of both fluid and solutes, it is convenient to
also consider the molar concentration of the fluid $m_{f}$. Note however that
$m_{f}$ denotes the molar concentration of the whole fluid phase; hence it
accounts for both solvent (water) and solutes. The molar concentration and the
density of the fluid is connected through $\rho_{f}=m_{f}M$, where $M_{f}$ is
the (average) molar mass of the fluid. We assume $M_{f}$ as well as $\rho_{f}$
to be constant, which is a valid approximation when the variability of $c_{1}$
and $c_{2}$ is not too large compared to $m_{f}$, or when the solvent and
solutes have comparable molar masses. This means that we can still model the
fluid with constant density.

The model equations for the considered processes in the fluid domain
$\Omega_{f}(t)$ are:

$\displaystyle\nabla\cdot\mathbf{v}$ $\displaystyle=0$ (1a)
$\displaystyle\partial_{t}(\rho_{f}\mathbf{v})+\nabla\cdot(\rho_{f}\mathbf{v}\otimes\mathbf{v})$
$\displaystyle=-\nabla p+\nabla\cdot\boldsymbol{\sigma}+\rho_{f}\mathbf{g}$
(1b) $\displaystyle\partial_{t}c_{i}+\nabla\cdot(c_{i}\mathbf{v})$
$\displaystyle=D_{i}\nabla^{2}c_{i}\quad i=1,2$ (1c)
$\displaystyle\partial_{t}(c_{p,f}\rho_{f}T_{f})+\nabla\cdot(c_{p,f}\rho_{f}T_{f}\mathbf{v})$
$\displaystyle=k_{f}\nabla^{2}T_{f}$ (1d)

Here, since the fluid is incompressible, the stress tensor is
$\boldsymbol{\sigma}=\mu(T_{f})\nabla\mathbf{v}$, where the viscosity $\mu$ is
allowed to vary with fluid temperature. We account for a body force
$\rho_{f}\mathbf{g}$, e.g. gravity. Further, $D_{i}$ is the diffusivity of the
solutes and we will assume $D_{1}=D_{2}=D$ and that it is constant. The heat
capacity $c_{p,f}$ and heat conductivity $k_{f}$ are also assumed constant.
Note that the first equation (1a) is the mass conservation equation for the
incompressible fluid and is originally
$\partial_{t}m_{f}+\nabla\cdot(m_{f}\mathbf{v})=0$, but is simplified to (1a)
as $m_{f}$ is constant. The second equation (1b) is the conservation of linear
momentum, while the third equation (1c) is the mass conservation equation of
the solute. Finally, (1d) expresses the conservation of energy inside the
fluid. Note that we use a rather simple energy equation by letting temperature
be the primary variable, but this is a valid assumption as long we consider
low fluid velocities and no evaporation of the fluid [25].

###  2.2 Mineral domain

The mineral domain consists of a stationary mineral having molar concentration
$m_{m}$ and mass density $\rho_{m}$. These two are connected through
$\rho_{m}=m_{m}M_{m}$, and are all considered constant. The mineral is a
solid, hence flow cannot take place. The mineral can however conduct heat,
hence an energy conservation equation is considered in the mineral domain
$\Omega_{m}(t)$:

$\displaystyle\mathbf{v}$ $\displaystyle=\mathbf{0}$ (2a)
$\displaystyle\partial_{t}(c_{p,m}\rho_{m}T_{m})$
$\displaystyle=k_{m}\nabla^{2}T_{m},$ (2b)

where $T_{m}$ is the temperature in the mineral, and $c_{p,m}$ and $k_{m}$ are
the (constant) heat capacity and heat conductivity of the mineral,
respectively. Here, (2a) expresses conservation of mass in the mineral, while
(2b) represents conservation of energy in the mineral.

###  2.3 Boundary conditions between the two domains

The interface $\Gamma(t)$ between the fluid and mineral domains is evolving
due to the mineral precipitation and dissolution reaction. When minerals
dissolve into ions, the ions become part of the fluid domain as solute. Hence,
there is a transfer of ions across the interface, which also corresponds to a
mass transfer across the interface. The interface obtains a new location due
to the lost minerals. The opposite takes place due to precipitation. To ensure
conservation of mass, ions and energy across the interface, we apply Rankine-
Hugoniot jump conditions [16], which balances the jump in the flux with the
jump of the conserved quantity across an evolving interface. For momentum we
assume a no-slip condition; that is, that the tangential component of the
fluid velocity is zero at the interface. The normal component of the fluid
velocity will in general not be zero, but be proportional to the normal
velocity of the evolving interface. Hence, on $\Gamma(t)$ we have:

$\displaystyle-m_{f}\mathbf{v}\cdot\mathbf{n}$
$\displaystyle=v_{n}((n_{1}+n_{2})m_{m}-m_{f}),$ (3a)
$\displaystyle\mathbf{v}$
$\displaystyle=v_{n}\frac{m_{f}-(n_{1}+n_{2})m_{m}}{m_{f}}\mathbf{n},$ (3b)
$\displaystyle(-c_{i}\mathbf{v}+D\nabla c_{i})\cdot\mathbf{n}$
$\displaystyle=v_{n}(n_{i}m_{m}-c_{i})\quad i=1,2,$ (3c)
$\displaystyle(-k_{m}\nabla T_{m}-c_{p,f}\rho_{f}T_{f}\mathbf{v}+k_{f}\nabla
T_{f})\cdot\mathbf{n}$
$\displaystyle=v_{n}(c_{p,m}\rho_{m}T_{m}-c_{p,f}\rho_{f}T_{f}).$ (3d)

The normal vector $\mathbf{n}$ is assumed to point into the mineral, and the
normal velocity $v_{n}$ is positive when $\Gamma(t)$ moves towards the mineral
domain, which corresponds to dissolution taking place. The top boundary
condition (3a) considers conservation of mass, from a molar perspective. It
expresses that one mineral molecule from the mineral domain dissolves into
$n_{1}+n_{2}$ ions in the fluid domain. This jump in mass is balanced by the
jump in the mass flux, which comprises of the advective term from the mass
conservation equation. The second boundary condition (3b) is a reformulation
of the no-slip condition and states that the normal component of the fluid
velocity from (3a) is the only component. The third boundary condition (3c)
ensures conservation of the two solutes: As one mineral molecule forms $n_{i}$
ion molecules, the factor $n_{i}$ appears. The jump in solute is balanced by
the jump in the solute flux. As the mineral is stationary, only the
contribution from the fluid domain is needed. These boundary conditions
(3a)-(3c) are the same as in [42]. Finally, the last boundary condition (3d)
accounts for the energy conservation across the evolving interface. In this
equation the jump in energy is balanced by the jump in the energy flux, which
has contributions from both mineral and fluid. This boundary condition is also
found in [7].

####  2.3.1 Reaction rates for mineral precipitation and dissolution

For the mineral precipitation/dissolution reaction, we formulate a reaction
rate. Following [23, 42] and including temperature dependence as in [5] we get

$m_{m}v_{n}=(-f_{p}(T_{f},c_{1},c_{2})-f_{d}(T_{f},c_{1},c_{2},\mathbf{x})),\quad\mathbf{x}\in\Gamma(t),$
(4)

where $f_{p}$ and $f_{d}$ are the precipitation and dissolution rates,
respectively. The rates are [12, 23]

$\displaystyle f_{p}(T_{f},c_{1},c_{2})$
$\displaystyle=k_{0}e^{-\frac{E}{RT_{f}}}\frac{c_{1}^{n_{1}}c_{2}^{n_{2}}}{K_{m}(T_{f})}\text{
and }$ (5a) $\displaystyle f_{d}(T_{f},c_{1},c_{2},\mathbf{x})$
$\displaystyle=k_{0}e^{-\frac{E}{RT_{f}}}w(\text{dist}(\mathbf{x},\mathcal{G}),T_{f},c_{1},c_{2})$
(5b)

where $k_{0}$ is a positive rate constant, $E$ is the activation energy, $R$
the gas constant and $K_{m}(T_{f})$ is the solubility product for the mineral.
As dissolution can only take place when there is still mineral present, we
introduce a distance measure that measures the distance from $\Gamma(t)$ to
the closest non-reactive solid or an external boundary, in (5b) denoted
$\mathcal{G}$. When there is no mineral present at a specific location, this
distance will be zero. To incorporate zero dissolution rate when there is no
mineral present, one can use a modified Heaviside graph [23]

$w(d,T_{f},c_{1},c_{2})=\left\\{\begin{array}[]{lr}0&\text{if }d<0,\\\
\text{min}\\{\frac{c_{1}^{n_{1}}c_{1}^{n_{2}}}{K_{m}(T_{f})},1\\}&\text{if
}d=0,\\\ 1&\text{if }d>0,\end{array}\right.$

The temperature can both affect the solubility of the mineral and also the
general speed of the reaction rates.

####  2.3.2 Volume changes due to mineral precipitation and dissolution

The boundary condition ensuring conservation of mass (3a) and the modified no-
slip condition (3b) also expresses any volume expansion occurring due to the
chemical reaction. Letting $K_{m}=\frac{m_{f}-(n_{1}+n_{2})m_{m}}{m_{f}}$, the
boundary condition (3b) on $\Gamma(t)$ reads

$\mathbf{v}=K_{m}v_{n}\mathbf{n}.$

Hence, the fluid near the evolving interface is moving proportionally with the
normal velocity of the evolving interface, and is scaled with the factor
$K_{m}$. This is caused by the fluid having to move due to
expansions/contractions caused by the chemical reaction.

If $K_{m}=0$, then $m_{f}=(n_{1}+n_{2})m_{m}$ which corresponds to the volume
taken up by one mineral molecule being the same as the volume taken up by
$n_{1}+n_{2}$ solute molecules. In this case, we have the more commonly used
condition $\mathbf{v}=\mathbf{0}$ on the interface. As the resulting velocity
at the interface is generally small compared to flow velocities, it is usually
a good approximation to let $K_{m}=0$ even if it is non-zero [26].

If $K_{m}>0$, the fluid density is larger which correspond to that the mineral
phase takes up a larger volume compared to the fluid. This means, for example
in the case of dissolution and $v_{n}>0$, that the fluid has to move towards
the mineral to account for the volume change. Hence the fluid flows towards
the mineral with the velocity $\mathbf{v}=K_{m}v_{n}\mathbf{n}$, which is
positive in the direction of the normal vector pointing into the mineral.
Oppositely, if $K_{m}<0$, the fluid density is lower, hence the mineral phase
takes up less volume compared to the fluid. In the case of dissolution and
$v_{n}>0$, the fluid has to move away from the mineral as it is pushed away by
the released solute needing more room. The fluid flows away from the mineral
with the velocity $\mathbf{v}=K_{m}v_{n}\mathbf{n}$, which is now negative in
the direction of the normal vector. Correspondingly, the opposite processes
occurs when $v_{n}<0$ as precipitation takes place.

The volume changes coming from $K_{m}\neq 0$ have to be accounted for when
choosing external boundary conditions for the flow, as e.g. a closed container
would lead to inconsistencies if there is overall net
precipitation/dissolution inside the container. Instead the fluid must be
allowed to escape/enter the domain in accordance with the mineral
precipitation/dissolution. Hence, the necessary condition is

$\int_{\partial\Omega_{f}(t)}m_{f}\mathbf{v}\cdot\mathbf{n}ds=-\int_{\Omega_{f}(t)}\partial_{t}m_{f}dV-\int_{\Gamma(t)}m_{f}K_{m}v_{n}ds,$
(6)

where the integral on the left-hand side is over the external boundary of the
fluid domain (that is, not facing a mineral), while the second integral on the
right-hand side is over all fluid-mineral interfaces found inside the domain.
Since we have constant density ($\partial_{t}m_{f}=0$), the inflow/outflow
across the exterior boundary has to match the volume changes at the mineral
boundary. If there is no expansion due to the reactions ($K_{m}=0$), (6)
reduces to a standard boundary condition for incompressible flow.

###  2.4 Domain evolution

We mention for completeness how to model the evolving interface $\Gamma(t)$ in
the case of a sharp-interface model. The location of the evolving interface
$\Gamma(t)$, and hence also the evolution of the two domains, can be
determined through a level-set approach. Letting
$S:[0,\infty)\times\Omega\to\mathbb{R}$, where
$\Omega=\Omega_{f}(t)\cup\Omega_{m}(t)\cup\Gamma(t)$, the level-set fulfills

$S(t,\mathbf{x})=\left\\{\begin{array}[]{lr}<0&\text{ if
}\mathbf{x}\in\Omega_{f}(t),\\\ 0&\text{ if }\mathbf{x}\in\Gamma(t),\\\
>0&\text{ if }\mathbf{x}\in\Omega_{m}(t).\end{array}\right.$

The level-set evolves following the level-set equation

$\partial_{t}S+v_{n}|\nabla S|=0\quad\text{ for }\mathbf{x}\in\Omega.$

Note that the level-set equation requires the normal velocity of the evolving
interface to be defined in the entire domain $\Omega$. This can be solved by,
for each point $\mathbf{x}\in\Omega$, picking the corresponding value of
$v_{n}$ from the point of $\Gamma(t)$ that is closest to $\mathbf{x}$.

######  Remark 1.

We will, without loss of generality, in the remainder of this paper only
consider a simplified chemical setting, namely the case where
$C_{1}+C_{2}\leftrightharpoons C$; hence, one molecule of each type of solute
is needed to form one mineral molecule. We will also assume that initial and
boundary conditions for the two solutes are the same. In this case, the two
solutes will always have the same concentration, which we will denote $c$.
Hence only one equation (1c) is needed for the solute, and only one boundary
condition (3c). The reaction rates (5) are modified such that $c_{1}=c_{2}$
and $n_{1}=n_{2}=1$.

##  3 Diffuse-interface models

Instead of dealing with a sharp interface, we now consider the transition
between fluid and mineral to be a smooth transition. In this case there are no
discontinuities between the two domains and all state variables will vary
smoothly between the domains. To describe this setting we use a phase field
$\phi(t,\mathbf{x})$. We will consider a phase field approaching the value 1
in the fluid domain and 0 in the mineral domain, and with a smooth transition
of width $O(\lambda)$ between the domains. Hence, the interface between the
fluid and mineral can no longer be said to be at an exact location, although
it is common to let the value $\phi=0.5$ represent the interface location for
practical purposes, which we will also apply in Section 4. The parameter
$\lambda$ quantifies the width of the diffuse interface, and the sharp-
interface model is expected to be recovered as $\lambda\to 0$, which we will
investigate in Section 4.

###  3.1 Phase-field equations

An Allen-Cahn equation [1] is adopted for the phase field. The standard Allen-
Cahn equation is not conservative and contain curvature effects that may be
non-physical. In this section we will consider a standard Allen-Cahn equation
as well as a reformulation to ensure conservation.

The standard Allen-Cahn equation for mineral precipitation and dissolution is

$\displaystyle\partial_{t}\phi=-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)+\gamma\nabla^{2}\phi-f_{\text{diff}}(\phi,T_{f},c),$ (7)

Here, $P(\phi)=8\phi(1-\phi^{2})$ is the double-well potential, ensuring that
the phase field attains values close to 0 and 1 for small $\lambda$. The
parameter $\gamma$ is a interface diffusion parameter and its role is
explained in detail in Section 4.

The reaction rate $f_{\text{diff}}$ is related to the reaction rate for the
sharp-interface model (4), but is reformulated to incorporate the diffuse
interface. We let

$f_{\text{diff}}(\phi,T,c)=\frac{4}{\lambda}\phi(1-\phi)\frac{1}{m_{m}}\big{(}f_{p}(T,c)-f_{d}(T,c)\big{)}.$

The factors $\phi(1-\phi)$ ensures that the reaction is only non-zero in the
diffuse transition zone where $\phi$ is neither 0 nor 1. The fraction
$4/\lambda$ ensures correct scaling as will be shown in section 4. For this
reaction rate the Heaviside graph is not needed as $\phi\equiv 1$ if all
mineral is dissolved and the resulting reaction rate is hence 0 in this case.
Hence, the precipitation $f_{p}$ and dissolution $f_{d}$ rates will be same as
in (5), but where $w\equiv 1$ making the distance-measure superfluous. For
convenience we will denote the net precipitation rate for
$f(T,c)=f_{p}(T,c)-f_{d}(T,c)$ in the following.

As we will see in Section 3.3, the Allen-Cahn equation is not conserving the
phase-field variable. Inspired by [13, 37] we consider the following
reformulation of the Allen-Cahn equation to overcome this:

$\displaystyle\partial_{t}\phi=-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)+\gamma\nabla^{2}\phi-
f_{\text{diff}}(\phi,T,c)+\frac{\gamma}{\lambda^{2}}\frac{1}{|\Omega|}\int_{\Omega}P^{\prime}(\phi)dV.$
(8)

where $|\Omega|=\int_{\Omega}1dV$ is the measure of the domain $\Omega$. The
reformulation (8) is conservative, which will be shown in Section 3.3. We
denote (7) as the _original_ Allen-Cahn equation, and the reformulated (8) as
the _conservative_ Allen-Cahn equations, respectively. The conservative
reformulation (8) was recently also considered in [6].

###  3.2 Modified transport equations

We now combine the conservation and transport equations from the fluid (1),
mineral (2) and their separating boundary (3) to formulate a one-field model
through incorporating the phase field. Letting $\phi$ be modeled by an Allen-
Cahn equation, such that it approaches the value 1 in the fluid and 0 in the
mineral, the following model equations describes the fluid flow and ensures
conservation of mass, ions and energy in the combined fluid-mineral domain
$\Omega$:

$\displaystyle\partial_{t}m_{\phi}+\nabla\cdot(m_{f}\phi\mathbf{v})$
$\displaystyle=0$ (9a)
$\displaystyle\partial_{t}(\rho_{f}\phi\mathbf{v})+\nabla\cdot(\rho_{f}\phi\mathbf{v}\otimes\mathbf{v})$
$\displaystyle=-\phi\nabla
p+\nabla\cdot\boldsymbol{\sigma}_{\phi}-g(\phi,\lambda)\mathbf{v}+\frac{1}{2}\rho_{f}\mathbf{v}\partial_{t}\phi+\rho_{f}\mathbf{g}$
(9b) $\displaystyle\partial_{t}(\phi
c+(1-\phi)m_{m})+\nabla\cdot(c\phi\mathbf{v})$
$\displaystyle=D\nabla(\phi\nabla c)$ (9c)
$\displaystyle\partial_{t}((c_{p}\rho)_{\phi}T)+\nabla\cdot(c_{p,f}\rho_{f}T\phi\mathbf{v})$
$\displaystyle=\nabla\cdot(k_{\phi}\nabla T)$ (9d)

Here, the newly introduced parameters $m_{\phi}$, $(c_{p}\rho)_{\phi}$ and
$k_{\phi}$ are given by

$m_{\phi}=\phi m_{f}+(1-\phi)2m_{m},\quad(c_{p}\rho)_{\phi}=\phi
c_{p,f}\rho_{f}+(1-\phi)c_{p,m}\rho_{m},\quad k_{\phi}=\phi
k_{f}+(1-\phi)k_{m}.$

This way they attain the expected value in the fluid when $\phi=1$ and in the
mineral when $\phi=0$. The factor two in front of $m_{m}$ is selected to
ensure conservation of mass as $n_{1}=n_{2}=1$. The temperature $T$ describes
the temperature of the entire system and hence replaces $T_{f}$ and $T_{g}$.
The modified stress tensor is
$\boldsymbol{\sigma}_{\phi}=\nu(\nabla\cdot(\phi\mathbf{v}))\mathbf{I}+\mu(\nabla(\phi\mathbf{v})+(\nabla(\phi\mathbf{v}))^{T})$.
Hence, as the mixture of fluid and mineral is quasi-compressible due to the
density differences, we here have a modified stress tensor, where $\nu$ is a
bulk viscosity for the fluid-mineral mixture. The term
$g(\phi,\lambda)\mathbf{v}$ is added to ensure $\mathbf{v}=0$ in the mineral
[8]. The interpolation function $g(\phi,\lambda)$ is a decreasing, surjective
and twice differentiable function fulfilling $g(1,\lambda)=0$ and
$g(0,\lambda)>0$. This way, (9b) reduces to $\mathbf{v}=\mathbf{0}$ when
$\phi=0$. We use $g(\phi,\lambda)=\frac{K_{\mathbf{v}}}{\lambda}(1-\phi)$ for
a constant $K_{\mathbf{v}}$, but mention that other possibilities are
discussed in [4, 8, 18].

Note that, as pointed out already in Section 2.3.2, we need to take care of
which boundary conditions to use on $\partial\Omega$ for the mass conservation
equation (9a). Directly letting $\phi\mathbf{v}\cdot\mathbf{n}=0$ on
$\partial\Omega$ could lead to unphysical effects as the mass concentration
$m_{\phi}$ can increase and decrease due to the chemical reaction (through
changes in $\phi$). Hence, to let fluid escape or enter through
$\partial\Omega$ accordingly, we need

$\int_{\partial\Omega}m_{f}\phi\mathbf{v}\cdot\mathbf{n}ds=-\int_{\Omega}\partial_{t}m_{\phi}dV.$
(10)

This condition is analog to (6). Note that if $m_{f}=2m_{m}$, then is
$m_{\phi}=m_{f}$ and hence constant, which leads to the right-hand side of
(10) to be zero. This corresponds to the case $K_{m}=0$ as discussed in
Section 2.3.2.

###  3.3 Conserved quantities of the phase-field models

We here shortly show why the phase field described by the standard Allen-Cahn
equation (7) corresponds to a non-conserved order parameter. If we assume zero
Neumann-conditions for the phase field on the external boundary
$\partial\Omega$, we have

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\Omega}\phi dV$
$\displaystyle=\int_{\Omega}\partial_{t}\phi
dV=\int_{\Omega}\big{(}-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)+\gamma\nabla^{2}\phi-\frac{4}{\lambda}\phi(1-\phi)\frac{1}{m_{m}}f(T,c)\big{)}dV$
$\displaystyle=\int_{\Omega}\big{(}-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)-\frac{4}{\lambda}\phi(1-\phi)\frac{1}{m_{m}}f(T,c)\big{)}dV+\int_{\partial\Omega}\gamma\nabla\phi\cdot\mathbf{n}ds$
$\displaystyle=\int_{\Omega}\big{(}-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)-\frac{4}{\lambda}\phi(1-\phi)\frac{1}{m_{m}}f(T,c)\big{)}dV,$

which is non-zero also in the case of no chemical reactions. For the
conservative Allen-Cahn (8), the phase field is conserved globally up to the
chemical reaction, that is

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\Omega}\phi dV$
$\displaystyle=\int_{\Omega}\big{(}-f_{\text{diff}}(\phi,T,c)\big{)}dV,$

hence we have $\frac{\text{d}}{\text{d}t}\int_{\Omega}\phi d\mathbf{x}=0$ in
the case of no chemical reactions. Note that since the integrated phase field
gives the volume of the fluid part of the domain, this integral can be
interpreted as the porosity of the system. This means that for the
conservative Allen-Cahn equation (8) only the chemical reaction can affect the
porosity. The scaling factor of $f_{\text{diff}}$,
$\frac{4}{\lambda}\phi(1-\phi)$, corresponds to the surface area of the fluid-
mineral interface when integrated.

For total mass $m_{\phi}$, total amount of ions $\phi c+(1-\phi)m_{m}$ and
total energy $(c_{p}\rho)_{\phi}T$ it is straightforward to show that these
quantities are conserved under zero inflow ($\phi\mathbf{v}\cdot\mathbf{n}=0$)
and zero Neumann ($\phi\nabla c\cdot\mathbf{n}=0$ and $k_{\phi}\nabla
T\cdot\mathbf{n}=0$) boundary conditions on $\partial\Omega$, also in the case
of chemical reactions. This is independent of which phase field is used to
model the diffuse transition zone. If investigating conservation of $\phi c$,
which corresponds to the amount of solute in the fluid, we get

$\displaystyle\frac{\text{d}}{\text{d}t}\int_{\Omega}\phi cdV$
$\displaystyle=\int_{\Omega}\partial_{t}(\phi
c)dV=m_{m}\int_{\Omega}\partial_{t}\phi dV,$

where the last integral is either connected to the reaction rate only, or to
the reaction rate together with curvature effects (if using the original
Allen-Cahn equation (7)). Hence, in the case of no chemical reaction, $\phi c$
is only conserved if using the conservative Allen-Cahn equation (8). In the
case of chemical reactions, the change in $\phi c$ corresponds to the changes
caused by the reaction rate.

##  4 Sharp-interface limits of the diffuse-interface models

To justify that the considered phase-field model formulated in Section 3 is
meaningful, we consider the limit as the width of the diffuse interface
$\lambda$ approaches zero to show that the model reduces to the expected
sharp-interface model from Section 2.

We follow the ideas of [10] and separate between behavior far away from the
interface and near it, and investigate the behavior as $\lambda\to 0$. We
apply so-called outer expansions for the behavior of the phase-field model far
away from the interface, and it will be shown that these outer expansions will
reduce to the model equations in the fluid and mineral domain as $\lambda\to
0$. Behavior near the interface is described through so-called inner
expansions of the unknowns and will be shown to capture the boundary
conditions as $\lambda\to 0$. The two expansions are connected through
matching conditions in a transition region where both are valid. We first
present the two expansions and their matching conditions, before considering
the outer and inner expansions of the phase-field models.

###  4.1 Inner and outer expansions and their matching conditions

Away from the interface, we consider the outer expansions of $\phi$, $c$, $p$,
$\mathbf{v}$ and $T$. For $\phi$ this reads

$\phi(t,\mathbf{x})=\phi_{0}^{\text{out}}(t,\mathbf{x})+\lambda\phi_{1}^{\text{out}}(t,\mathbf{x})+\lambda^{2}\phi_{2}^{\text{out}}(t,\mathbf{x})+\dots$
(11)

and similarly for the other variables. Note that $\mu_{f}$ is considered
(given) function of $T$ and hence do not need an expansion.

The inner expansions are valid near the interface, and it is hence convenient
to make a change of coordinates into local, curvilinear coordinates. We let
$\Gamma_{\lambda}(t)$ denote the set of points $\mathbf{y}_{\lambda}\in\Omega$
such that $\phi(t,\mathbf{y}_{\lambda})=1/2$. Note that the location of these
points $\mathbf{y}_{\lambda}$ vary with time as the phase field evolves. Let
$\mathbf{s}$ be an arc-length parametrization of the points on
$\Gamma_{\lambda}(t)$ (note that $\mathbf{s}$ is a scalar if
$\Omega\subset\mathbb{R}^{2}$ and consists of two coordinates if
$\Omega\subset\mathbb{R}^{3}$) and $\mathbf{n}_{\lambda}$ the normal vector at
$\Gamma_{\lambda}(t)$ pointing into the mineral (that is, towards decreasing
$\phi$). Further, we let $r$ be the signed distance from a point $\mathbf{x}$
to the interface, being negative when inside the fluid ($\phi>0.5$) and
positive inside the mineral ($\phi<0.5$). This way we can let each point
$\mathbf{x}$ be described through

$\mathbf{x}=\mathbf{y}_{\lambda}(t,\mathbf{s})+r\mathbf{n}_{\lambda}(t,\mathbf{s}).$
(12)

Through this transform, $(r,\mathbf{s})$ makes up a new, orthogonal coordinate
system. It is shown in [10] that

$|\nabla r|=1,\ \ \ \nabla r\cdot\nabla s_{i}=0,\ \ \ \partial_{t}r=-v_{n},\ \
\ \nabla^{2}r=\frac{\kappa+2\Pi r}{1+\kappa r+\Pi r^{2}},$

where $\kappa$ and $\Pi$ are the mean and Gaussian curvature of the interface.
We assume that $\mathbf{y}_{\lambda}$ and $\mathbf{n}_{\lambda}$ fulfill
expansions $\mathbf{y}_{\lambda}=\mathbf{y}_{0}+O(\lambda)$ and
$\mathbf{n}_{\lambda}=\mathbf{n}_{0}+O(\lambda)$ where $\mathbf{y}_{0}$ is a
point of the interface $\Gamma_{0}^{\text{out}}(t)$ defined through
$\phi_{0}^{\text{out}}=1/2$ and $\mathbf{n}_{0}$ being the normal vector of
$\Gamma_{0}^{\text{out}}(t)$ pointing into the mineral. Then, with
$z=r/\lambda$ we consider the inner expansions

$\phi(t,\mathbf{x})=\phi_{0}^{\text{in}}(t,z,\mathbf{s})+\lambda\phi_{1}^{\text{in}}(t,z,\mathbf{s})+\lambda^{2}\phi_{2}^{\text{in}}(t,z,\mathbf{s})+\dots,$
(13)

and similarly for the other unknowns. In the local, curvilinear coordinates
and when accounting for the scaling of $z$, we need to rewrite derivatives.
For a generic variable $u$ or $\mathbf{u}$, we obtain [10]

$\displaystyle\partial_{t}u$
$\displaystyle=-\lambda^{-1}v_{n,0}\partial_{z}u^{\text{in}}+(\partial_{t}+\partial_{t}s\cdot\nabla_{\mathbf{s}})u^{\text{in}}+O(\lambda),$
(14a) $\displaystyle\nabla_{x}u$
$\displaystyle=\lambda^{-1}\partial_{z}u^{\text{in}}\mathbf{n}_{0}+\nabla_{\mathbf{s}}u^{\text{in}}+O(\lambda),$
(14b) $\displaystyle\nabla_{x}\cdot\mathbf{u}$
$\displaystyle=\lambda^{-1}\partial_{z}\mathbf{u}^{\text{in}}\cdot\mathbf{n}_{0}+\nabla_{\mathbf{s}}\cdot\mathbf{u}^{\text{in}}+O(\lambda),$
(14c) $\displaystyle\nabla_{x}^{2}v$
$\displaystyle=\lambda^{-2}\partial_{z}^{2}u+\lambda^{-1}\kappa_{0}\partial_{z}u+O(1),$
(14d)

where we have used that $\nabla_{x}^{2}r=\kappa_{0}+O(\lambda)$ is the lowest
order mean curvature and $v_{n}=v_{n,0}+O(\lambda)$. Here, $\kappa_{0}$ and
$v_{n,0}$ are the curvature and normal velocity of the interface
$\Gamma_{0}^{\text{out}}(t)$. In the last condition (14d) we have used the
properties $|\nabla r|=1$ and $\nabla r\cdot\nabla s_{i}=0$.

The outer (11) and inner (13) expansions are tied together via (12). For fixed
$t$ and $\mathbf{s}$, we let $\mathbf{y}_{1/2\pm}$ denote the limit $r\to
0^{+}$ (i.e., from the mineral side) and $r\to 0^{-}$ (i.e., from the fluid
side) of $\mathbf{x}$. We demand that when the outer expansions approach the
boundary they need to match the inner expansions as they leave the boundary.
That is, the limit when $r\to 0^{\pm}$ of the outer expansions need to match
the limit when $z\to\pm\infty$ of the inner expansions. Using Taylor
expansions, one can show that the following matching conditions

$\displaystyle\lim_{z\to\pm\infty}\phi_{0}^{\text{in}}(t,z,\mathbf{s})$
$\displaystyle=\phi_{0}^{\text{out}}(t,\mathbf{y}_{1/2\pm}),$ (15a)
$\displaystyle\lim_{z\to\pm\infty}\partial_{z}\phi_{0}^{\text{in}}(t,z,\mathbf{s})$
$\displaystyle=0,$ (15b)
$\displaystyle\lim_{z\to\pm\infty}\partial_{z}\phi_{1}^{\text{in}}(t,z,\mathbf{s})$
$\displaystyle=\nabla\phi_{0}^{\text{out}}(t,\mathbf{y}_{1/2\pm})\cdot\mathbf{n}_{0},$
(15c)

are fulfilled, and similarly for the other variables [10].

###  4.2 Outer expansions

We first apply the outer expansions (11) to the phase-field equations (7) and
(8). Inserting the outer expansions in (7), (8), the dominating term when
$\lambda\to 0$ will in both cases be

$P^{\prime}(\phi_{0}^{\text{out}})=0.$

Recalling that $P(\phi)=8\phi^{2}(1-\phi)^{2}$, the above equation has the
three solutions $\phi_{0}^{\text{out}}=0$, $\phi_{0}^{\text{out}}=0.5$ and
$\phi_{0}^{\text{out}}=1$. Only $\phi_{0}^{\text{out}}=0$ and
$\phi_{0}^{\text{out}}=1$ give stable solutions. We let $\Omega_{m,0}(t)$ and
$\Omega_{f,0}(t)$ be the collection of points $\mathbf{x}\in\Omega$ associated
with where $\phi_{0}^{\text{out}}=0$ and $\phi_{0}^{\text{out}}=1$,
respectively.

For the remaining phase-field model equations (9), we consider the two limit
cases of $\phi_{0}^{\text{out}}=0$ and $\phi_{0}^{\text{out}}=1$ separately.
When $\phi_{0}^{\text{out}}=0$ (that is, in $\Omega_{m,0}(t)$), the dominating
terms of the outer expansions (11) in (9) give that

$\displaystyle\partial_{t}m_{m}$ $\displaystyle=0,$ (16a)
$\displaystyle\mathbf{v}_{0}^{\text{out}}$ $\displaystyle=\mathbf{0},$ (16b)
$\displaystyle\partial_{t}m_{m}$ $\displaystyle=0,$ (16c)
$\displaystyle\partial_{t}(c_{p,m}\rho_{m}T_{0}^{\text{out}})$
$\displaystyle=k_{m}\nabla^{2}T_{0}^{\text{out}}.$ (16d)

The first and third equations simply tell us that the mineral molar density is
constant, which was also a model assumption, cf. Section 2.2. The equations
(16) correspond to the ones formulated in (2) for the mineral domain. When
$\phi_{0}^{\text{out}}=1$ (that is, in $\Omega_{f,0}(t)$), the dominating
terms of the outer expansions (11) in (9) give that

$\displaystyle\partial_{t}m_{f}+\nabla\cdot(m_{f}\mathbf{v}_{0}^{\text{out}})$
$\displaystyle=0,$ (17a)
$\displaystyle\partial_{t}(\rho_{f}\mathbf{v}_{0}^{\text{out}})+\nabla\cdot(\rho_{f}\mathbf{v}_{0}^{\text{out}}\otimes\mathbf{v}_{0}^{\text{out}})$
$\displaystyle=-\nabla
p_{0}^{\text{out}}+\nabla\cdot\boldsymbol{\sigma}_{0}+\rho_{f}\mathbf{g},$
(17b)
$\displaystyle\partial_{t}c_{0}^{\text{out}}+\nabla\cdot(c_{0}^{\text{out}}\mathbf{v}_{0}^{\text{out}})$
$\displaystyle=D\nabla^{2}c_{0}^{\text{out}},$ (17c)
$\displaystyle\partial_{t}(c_{p,f}\rho_{f}T_{0}^{\text{out}})+\nabla\cdot(c_{p,f}\rho_{f}T_{0}^{\text{out}})$
$\displaystyle=k_{f}\nabla^{2}T_{0}^{\text{out}},$ (17d)

where
$\boldsymbol{\sigma}_{0}=\nu(\nabla\cdot\mathbf{v}_{0}^{\text{out}})\mathbf{I}+\mu(T_{0}^{\text{out}})(\nabla\mathbf{v}_{0}^{\text{out}}+(\nabla\mathbf{v}_{0}^{\text{out}})^{T})$.
Since $m_{f}$ is constant, (17a) reduces to
$\nabla\cdot\mathbf{v}_{0}^{\text{out}}=0$, which also means that we can
simplify the stress tensor to
$\boldsymbol{\sigma}_{0}=\mu(T_{0}^{\text{out}})\nabla\mathbf{v}_{0}^{\text{out}}$.
Note that we have here assumed that $\mu(T)$ is Lipschitz-continuous. The
equations (17) hence correspond to the ones formulated in (1) for the fluid
domain.

###  4.3 Inner expansions

We now apply the inner expansions (13) to the model equations (7), (8) and (9)
to recover the boundary conditions at the evolving interface.

####  4.3.1 Phase-field equations

For the original Allen-Cahn equation (7), the dominating $O(\lambda^{-2})$
terms when inserting the inner expansions (13) (and keeping in mind the
rescaled derivatives (14)) are

$P^{\prime}(\phi_{0}^{\text{in}})=\partial_{z}^{2}\phi_{0}^{\text{in}}.$ (18)

From the first matching condition (15a) we have that
$\lim_{z\to\infty}\phi_{0}^{\text{in}}=0$ and
$\lim_{z\to\infty}\phi_{0}^{\text{in}}=1$. By multiplying (18) by
$\partial_{z}\phi_{0}^{\text{in}}$, integrating twice with respect to $z$ and
using the matching conditions and that $\phi_{0}^{\text{in}}(z=0)=0.5$, we get
that

$\phi_{0}^{\text{in}}(t,z,\mathbf{s})=\phi_{0}^{\text{in}}(z)=\frac{1}{1+e^{4z}}=\frac{1}{2}\tanh(2z)$
(19)

fulfills the equation. For the $O(\lambda^{-1})$ terms one obtains

$\gamma\big{(}P^{\prime\prime}(\phi_{0}^{\text{in}})-\partial_{z}^{2}\big{)}\phi_{1}^{\text{in}}=(v_{n,0}+\gamma\kappa_{0})\partial_{z}\phi_{0}^{\text{in}}-4\phi_{0}^{\text{in}}(1-\phi_{0}^{\text{in}})\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}}).$

We follow the same steps as in [35] and view the left-hand side as an operator
$\mathcal{F}$ depending on $\phi_{0}^{\text{in}}$ and applied to
$\phi_{1}^{\text{in}}$. The operator $\mathcal{F}$ is a Fredholm operator of
index zero, hence the above equation has a solution if and only if the right-
hand side, denoted by $A(\phi_{0}^{\text{in}})$, is orthogonal to the kernel
of $\mathcal{F}$. It follows from (18) that $\partial_{z}\phi_{0}^{\text{in}}$
is in the kernel of $\mathcal{F}$. Since $v_{n,0}$, $\kappa_{0}$,
$T_{0}^{\text{in}}$ and $c_{0}^{\text{in}}$ are independent of $z$ (the two
latter will be shown in the following subsections), the solvability condition
implies

$\displaystyle 0$
$\displaystyle=\int_{-\infty}^{\infty}A(\phi_{0}^{\text{in}})\partial_{z}\phi_{0}^{\text{in}}dz$
$\displaystyle=(v_{n,0}+\gamma\kappa_{0})\int_{-\infty}^{\infty}(\partial_{z}\phi_{0}^{\text{in}})^{2}dz-4\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}})\int_{-\infty}^{\infty}\phi_{0}^{\text{in}}(1-\phi_{0}^{\text{in}})\partial_{z}\phi_{0}^{\text{in}}dz$
$\displaystyle=\frac{2}{3}(v_{n,0}+\gamma\kappa_{0}+\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}})).$

We apply the matching condition (15a) for $T$ and $c$ at the moving interface
and obtain

$v_{n,0}=-\gamma\kappa_{0}-\frac{1}{m_{m}}f(T_{0}^{\text{out}}(t,\mathbf{y}_{1/2-}),c_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})),$

which corresponds to the boundary condition (4) at the moving interface, with
an additional curvature-driven motion, which is responsible for the Allen-Cahn
equation not being conservative. Since the Allen-Cahn equation is derived from
curvature-driven flow, it is expected that this term appears [1]. Although the
curvature effects would be zero if $\gamma=0$, this is not a valid choice for
the Allen-Cahn equation (7).

We address the sharp-interface limit of the conservative Allen-Cahn equation
(8). We follow similar steps as [13] who investigated the sharp-interface
limit for a conservative Allen-Cahn equation without chemical reactions. The
equation (8) is first split in two equations:

$\displaystyle\partial_{t}\phi=-\frac{1}{\lambda^{2}}\gamma
P^{\prime}(\phi)+\gamma\nabla^{2}\phi-
f_{\text{diff}}(\phi,T,c)+\frac{1}{\lambda}\gamma\omega(t)$ (20a)
$\displaystyle\omega(t)=\frac{1}{\lambda}\frac{1}{|\Omega|}\int_{\Omega}P^{\prime}(\phi)dV.$
(20b)

The average of $P^{\prime}(\phi)$ appearing in (20b) is small and will
decrease as $\lambda$ decreases. It is shown in [13] that
$\omega(t)=O(\lambda^{0})$ as $\lambda\to 0$. The inner expansions (13) of
(20a) are dominated by the terms

$P^{\prime}(\phi_{0}^{\text{in}})=\partial_{z}^{2}\phi_{0}^{\text{in}},$

which has solution (19). The next order terms are affected by the newly
introduced term:

$\big{(}P^{\prime\prime}(\phi_{0}^{\text{in}})-\partial_{z}^{2}\big{)}\phi_{1}^{\text{in}}=(v_{n,0}+\gamma\kappa_{0})\partial_{z}\phi_{0}^{\text{in}}-4\phi_{0}^{\text{in}}(1-\phi_{0}^{\text{in}})\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}})+\gamma\omega_{0}^{\text{in}},$

where
$\omega_{0}^{\text{in}}=\frac{1}{\lambda}\frac{1}{|\Omega|}\int_{\Omega}P^{\prime}(\phi_{0}^{\text{in}})dV$.
As with the original Allen-Cahn equation we can use a solvability argument to
proceed. The above equation only has a solution when

$\int_{-\infty}^{\infty}\Big{(}(v_{n,0}+\gamma\kappa_{0})\partial_{z}\phi_{0}^{\text{in}}-4\phi_{0}^{\text{in}}(1-\phi_{0}^{\text{in}})\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}})+\gamma\omega_{0}^{\text{in}}\Big{)}\partial_{z}\phi_{0}^{\text{in}}dz=0.$

Following same argumentation as with the original Allen-Cahn equation, and
since $\omega_{0}^{\text{in}}$ is independent of $z$ by definition, we get

$v_{n,0}=-\gamma\kappa_{0}-\frac{1}{m_{m}}f(T_{0}^{\text{in}},c_{0}^{\text{in}})+3\gamma\omega_{0}^{\text{in}}.$

Finally, [13] shows that the condition of
$\int_{\Omega}\partial_{t}\phi_{0}^{\text{in}}d\mathbf{x}=0$ is equivalent to

$\omega_{0}^{\text{in}}=\frac{1}{3}\frac{1}{|\Gamma(t)|}\int_{\Gamma(t)}\kappa_{0}ds=\frac{1}{3}\overline{\kappa}_{0},$

where $\overline{\kappa}$ is the average curvature along $\Gamma(t)$. We hence
arrive at

$v_{n,0}=-\gamma(\kappa_{0}-\overline{\kappa}_{0})-\frac{1}{m_{m}}f(T_{0}^{\text{out}}(t,\mathbf{y}_{1/2-}),c_{0}^{\text{out}}(t,\mathbf{y}_{1/2-}))\quad\text{on
}\Gamma(t),$

This means that the interface velocity is still driven by both the chemical
reaction and by curvature, but where the curvature-driven movement fulfills
conservation of the phase-field parameter. The curvatured-driven motion
redistributes the mineral towards constant curvature; towards bubbles [38].
Hence, the fluid-solid interface moves due to the chemical reaction, and
additionally in a conservative manner by redistributing mineral towards
constant curvature.

####  4.3.2 Mass conservation equation

Inserting the inner expansions (13) into the mass conservation equation (9a)
gives the dominating $O(\lambda^{-1})$ term

$-v_{n,0}\partial_{z}m_{\phi}+\partial_{z}(m_{f}\phi_{0}^{\text{in}}\mathbf{v}_{0}^{\text{in}})\cdot\mathbf{n}_{0}=0.$

Integrating with respect to $z$ from $-\infty$ to $+\infty$ and using matching
conditions (15a) give

$-m_{f}\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})\cdot\mathbf{n}_{0}=v_{n,0}(2m_{m}-m_{f}).$
(21)

This boundary condition corresponds to (3a).

####  4.3.3 Momentum conservation equation

The dominating $O(\lambda^{-2})$ terms arising from inserting the inner
expansions (13) into (9b) are

$\mu\partial_{z}^{2}(\phi_{0}^{\text{in}}\mathbf{v}_{0}^{\text{in}})+(\nu+\mu)\partial_{z}^{2}(\phi_{0}^{\text{in}}\mathbf{v}_{0}^{\text{in}}\cdot\mathbf{n}_{0})\mathbf{n}_{0}=0$

Note that $\mu$ and $\nu$ can be functions of $T_{0}^{\text{in}}$, which is
still independent of $z$. We integrate twice with respect to $z$ and apply
matching conditions (15a) and (15b) and hence arrive at

$\mu\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})=-(\nu+\mu)(\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})\cdot\mathbf{n}_{0})\mathbf{n}_{0}.$

Here, the temperature $T_{0}^{\text{in}}$ inside $\mu$ and $\nu$ is also to be
evaluated at $\mathbf{y}_{1/2-}$. Note that the current equation already gives
us a no-slip condition: The velocity has only a normal component at the
boundary, and hence no tangential component. Combining with (21) we can write
this condition as

$\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})=-\frac{\nu+\mu}{\mu}\frac{m_{f}-2m_{m}}{m_{f}}v_{n,0}\mathbf{n}_{0}.$

However, the current condition is only consistent with (21) when $\nu=-2\mu$.
In this case,

$\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})=\frac{m_{f}-2m_{m}}{m_{f}}v_{n,0}\mathbf{n}_{0},$

which then corresponds to (3b). Hence, this relation between the two
viscosities $\nu$ and $\mu$ give consistency and the expected boundary
condition in the sharp-interface limit. The second viscosity $\nu$ was here
introduced since the fluid-mineral mixture is quasi-incompressible. A second
viscosity is usually connected to how a fluid reacts on an fluid
expansion/contraction and is here chosen in relation to $\mu$ to ensure that
the chosen sharp-interface boundary condition is obtained.

######  Remark 2.

In case $m_{f}=2m_{m}$, we have no expansion/contraction due to the mineral
precipitation and dissolution. In this case $m_{\phi}=m_{f}$ and (9a) would
have expressed $\nabla\cdot(\phi\mathbf{v})=0$ as $m_{f}$ is constant. Then,
the second viscosity $\nu$ does not have to be introduced as (9b) could have
used a simpler stress tensor. The sharp-interface limit in this case gives
directly $\mathbf{v}_{0}^{\text{in}}(t,\mathbf{y}_{1/2-})=\mathbf{0}$.

####  4.3.4 Solute conservation equation

Inserting the inner expansions (13) into the solute conservation equation
(9c), the dominating $O(\lambda^{-2})$ term is

$\partial_{z}(\phi_{0}^{\text{in}}\partial_{z}c_{0}^{\text{in}})=0.$

Integrating with respect to $z$ and using the matching condition (15b)
together with the fact that $\phi_{0}^{\text{in}}>0$, we get that

$\partial_{z}c_{0}^{\text{in}}=0.$

Hence, $c_{0}^{\text{in}}$ is independent of $z$. For the $O(\lambda^{-1})$
terms we obtain

$-v_{n,0}\partial_{z}\big{(}\phi_{0}^{\text{in}}c_{0}^{\text{in}}+(1-\phi_{0}^{\text{in}})m_{m}\big{)}+\partial_{z}(c_{0}^{\text{in}}\phi_{0}^{\text{in}}\mathbf{v}_{0}^{\text{in}})\cdot\mathbf{n}_{0}=D\partial_{z}(\phi_{0}^{\text{in}}\partial_{z}c_{1}^{\text{in}}).$

Integrating with respect to $z$ from $-\infty$ to $+\infty$ and using matching
conditions (15) we get

$v_{n,0}\big{(}m_{m}-c_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\big{)}=\big{(}-c_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\mathbf{v}_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})+D\nabla
c_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\big{)}\cdot\mathbf{n}_{0},$

which corresponds to (3c).

####  4.3.5 Energy conservation equation

By inserting the inner expansions (13) into (9d), the dominating
$O(\lambda^{-2})$ term is

$\partial_{z}^{2}T_{0}^{\text{in}}=0.$

When integrating this equation with respect to $z$ and using the matching
condition (15b), it becomes clear that

$\partial_{z}T_{0}^{\text{in}}=0;$

that is, $T_{0}^{\text{in}}$ is independent of $z$. The $O(\lambda^{-1})$
terms are

$-v_{n,0}\partial_{z}((c_{p}\rho)_{\phi^{0}}T_{0}^{\text{in}})+\partial_{z}(c_{f}\rho_{f}T_{0}^{\text{in}}\phi_{0}^{\text{in}}\mathbf{v}_{0}^{\text{in}})\cdot\mathbf{n}_{0}=\partial_{z}(k_{\phi^{0}}\partial_{z}T_{1}^{\text{in}}),$

where the notation
$(c_{p}\rho)_{\phi^{0}}=\phi_{0}^{\text{in}}c_{p,f}\rho_{f}+(1-\phi_{0}^{\text{in}})c_{p,m}\rho_{m}$
and similarly for $k_{\phi^{0}}$. We integrate with respect to $z$ from
$-\infty$ to $+\infty$, taking advantage of that $T_{0}^{\text{in}}$ is
independent of $z$. After applying matching conditions (15) we get that

$\displaystyle v_{n,0}$
$\displaystyle\big{(}c_{p,m}\rho_{m}T_{0}^{\text{out}}(t,\mathbf{y}_{1/2+})-c_{p,f}\rho_{f}T_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\big{)}=\big{(}-k_{m}\nabla
T_{0}^{\text{out}}(t,\mathbf{y}_{1/2+})$ $\displaystyle-
c_{p,f}\rho_{f}T_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\mathbf{v}_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})+k_{f}\nabla
T_{0}^{\text{out}}(t,\mathbf{y}_{1/2-})\big{)}\cdot\mathbf{n}_{0},$

which corresponds to (3d) when associating the temperature at the positive
(mineral) side with $T_{m}$ and the temperature at the negative (fluid) side
with $T_{f}$.

##  5 Numerical discretization of the phase-field model

We here numerically discretize the resulting phase-field model. Special
attention is given to the conservative Allen-Cahn equation to ensure that the
numerical discretization honors the conservative property. We apply Finite
Volume (FV) discretization in space, which is outlined in Section 5.1, and it
is shown under which conditions the conservative Allen-Cahn equation (8) is
discretely conservative in Section 5.2. The conservative Allen-Cahn equation
is non-linear and challenging to solve with Newton iterations due to the non-
local term, and we formulate L-scheme iterations which converge to a
conservative solution in Section 5.3. The conservative Allen-Cahn equation is
coupled to the other model equations in an iterative manner, which in a
simplified setup is proven to converge in Section 5.4.

###  5.1 Finite Volume discretization

For numerical discretization, we apply a cell-centered FV scheme on an
admissible mesh $\mathcal{E}$ [15], where the phase-field variable (and
similarly for temperature and solute) is approximated in the cell-center and
we calculate fluxes across cell edges. The fluxes are for simplicity two-point
approximations, but extensions using multi-point flux approximations are
possible. By integrating the conservative Allen-Cahn equation (8) across one
element and using Gauss’ theorem on the diffusive term, the discrete FV
approximation of (8) reads

$\displaystyle\frac{\phi_{K}^{n+1}-\phi_{K}^{n}}{\Delta
t}+\frac{\gamma}{\lambda^{2}}P^{\prime}(\phi_{K}^{\ell})=$
$\displaystyle\frac{\gamma}{|K|}\sum_{L\in\mathcal{N}(K)}|\sigma_{K,L}|F^{\phi,\ell}_{K,L}-\frac{4}{\lambda}\phi_{K}^{\ell}(1-\phi_{K}^{\ell})\frac{f(T_{K}^{\ell},c_{K}^{\ell})}{m_{m}}$
$\displaystyle+\frac{\gamma}{\lambda^{2}}\frac{1}{|\Omega|}\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi^{\ell}_{J}),$
(22)

where $|K|$ is the measure of element $K$. We have here also discretized in
time, and let superscript $n$ and $n+1$ denote the time step. The notation
$\phi_{K}$ refers to the cell-centered value of $\phi$ in $K$ and is assumed
to approximate the value of $\phi$ in entire $K$. Further, $\mathcal{N}(K)$
refers to the neighboring elements of $K$ and $|\sigma_{K,L}|$ is the measure
of the edge $\sigma_{K,L}$ between element $K$ and the neighboring element
$L$. The integral $\int_{\Omega}P^{\prime}(\phi)d\mathbf{x}$ has been
approximated by the sum
$\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi_{J}^{\ell})$. The superscript $\ell$
denotes time step, and is $n$ or $n+1$ depending on whether forward or
backward Euler is applied, respectively. The possible choices for time
stepping will be addressed later. The fluxes $F_{K,L}^{\phi,\ell}$ approximate
the diffusive flux $\nabla^{2}\phi$ on the edge between elements $K$ and $L$
and are given by

$F_{K,L}^{\phi,\ell}=\frac{\phi_{L}^{\ell}-\phi_{K}^{\ell}}{d_{K,L}},$

where $d_{K,L}$ is the Euclidean distance between the center points $x_{K}\in
K$ and $x_{L}\in L$. We have $F_{K,L}^{\phi,\ell}=-F_{L,K}^{\phi,\ell}$ on
interior edges, which ensures the conservation between elements.

###  5.2 Discrete conservation of the conservative Allen-Cahn equation

We consider the global conservation of (22) up to the chemical reaction by
summing up (22) for all $K\in\mathcal{E}$, hence across entire domain
$\Omega$. Summing up over all $K\in\mathcal{E}$ and using Neumann boundary
conditions for $\phi$ and that $F_{K,L}^{\phi,\ell}=-F_{L,K}^{\phi,\ell}$, we
obtain

$\displaystyle\sum_{K\in\mathcal{E}}\phi^{n+1}_{K}|K|=$
$\displaystyle\sum_{K\in\mathcal{E}}\phi^{n}_{K}|K|-\sum_{K\in\mathcal{E}}\frac{4}{\lambda}\phi_{K}^{\ell}(1-\phi_{K}^{\ell})\frac{f(T_{K}^{\ell},c_{K}^{\ell})}{m_{m}}|K|$
$\displaystyle-\frac{\gamma}{\lambda^{2}}\sum_{K\in\mathcal{E}}\big{(}|K|P^{\prime}(\phi_{K}^{\ell})-\frac{|K|}{|\Omega|}\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi^{\ell}_{J})\big{)}$

The summation
$\sum_{K\in\mathcal{E}}\big{(}|K|P^{\prime}(\phi_{K}^{\ell})-\frac{|K|}{|\Omega|}\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi^{\ell}_{J})\big{)}$
adds up to zero under certain restrictions on the time evaluation. Denoting
the sum $\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi^{\ell}_{J}):=\Pi$, we have
that the summation can be rewritten to

$\sum_{K\in\mathcal{E}}\big{(}|K|P^{\prime}(\phi_{K}^{\ell})-\frac{|K|}{|\Omega|}\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi^{\ell}_{J})\big{)}=\Pi-\sum_{K\in\mathcal{E}}\frac{|K|}{|\Omega|}\Pi.$

As $|\Omega|$ and $\Pi$ are independent of the element, it is obvious that the
terms sum up to zero as $\sum_{K\in\mathcal{E}}|K|=|\Omega|$ as the mesh is
conforming. This remains true as long as the same choice of $\ell$ is taken
for the evaluation of $P^{\prime}(\phi^{\ell}_{K})$ in both appearances for
corresponding elements. This means that both forward Euler and backward Euler
will be discretely conservative, as also shown in [6]. However, as also noted
in [6], an explicit discretization of the non-linear terms is unstable, while
the implicit choice is slow to solve with Newton iterations as the Jacobian is
full. In the following we suggest a strategy to overcome these difficulties,
while still assuring the discrete conservation of (8).

###  5.3 L-scheme iterations of a semi-implicit scheme

We here consider a particular semi-implicit version of (8); and consider an
alternative non-linear solver of the resulting non-linear system of equations.
We will consider the L-scheme [33, 27], which reformulates the discrete
equation to iterations which form a contraction. To simplify the notation, we
consider only the time-discrete equation. Initially we let the scheme be
implicit, hence $\ell=n+1$, and we address how to solve

$\displaystyle\frac{1}{\Delta
t}(\phi^{n+1}-\phi^{n})+\frac{\gamma}{\lambda^{2}}P^{\prime}(\phi^{n+1})=$
$\displaystyle\gamma\nabla^{2}\phi^{n+1}-\frac{4}{\lambda}\phi^{n+1}(1-\phi^{n+1})\frac{f(T,c)}{m_{m}}$
$\displaystyle+\frac{\gamma}{\lambda^{2}}\frac{1}{|\Omega|}\int_{\Omega}P^{\prime}(\phi^{n+1})dV.$
(23)

As we only address the non-linearity with respect to the phase field, we
assume the temperature and concentration in the reactive term given. Note that
the L-scheme approach is of course also applicable to solve the non-
linearities in (7), which has been done in [3], but we here focus on how to
apply it to (8), and show how the L-scheme iterations can be assured to
converge to the conservative solution.

####  5.3.1 Convergence of the L-scheme

The L-scheme adds an extra stabilization term to create a contraction, which
converges by the Banach fixed point theorem. We need the non-linearity in the
resulting discrete equation to be monotonically increasing or decreasing.
However, the non-linearity in (5.3);

$G(\phi,T,c)=-\frac{\gamma}{\lambda^{2}}P^{\prime}(\phi)-\frac{4}{\lambda}\phi(1-\phi)\frac{f(T,c)}{m_{m}}+\frac{\gamma}{\lambda^{2}}\frac{1}{|\Omega|}\int_{\Omega}P^{\prime}(\phi)dV,$
(24)

is both increasing and decreasing with respect to $\phi$. Therefore we split
the non-linearity $G(\phi,T,c)$ into an increasing and decreasing part, and
solve the increasing part explicitly and decreasing implicitly. Note that we
consider $T,c$ known and fixed in the following. Hence, we define

$\displaystyle G_{-}(\phi,T,c)$
$\displaystyle=\int_{0}^{\phi}\min\\{0,\partial_{1}G(\psi,T,c)\\}d\psi,$
$\displaystyle G_{+}(\phi,T,c)$
$\displaystyle=\int_{0}^{\phi}\max\\{0,\partial_{1}G(\psi,T,c)\\}d\psi,$

and consider the time discretization

$\frac{1}{\Delta
t}(\phi^{n+1}-\phi^{n})=\gamma\nabla^{2}\phi^{n+1}+G_{-}(\phi^{n+1},T,c)+G_{+}(\phi^{n},T,c).$
(25)

The L-scheme iterations solve the non-linear equation (25) iteratively through
adding a stabilizing therm. By letting $j$ be the iteration index and setting
$\phi^{n+1,0}:=\phi^{n}$, we solve

$\displaystyle\frac{1}{\Delta t}(\phi^{n+1,j+1}-\phi^{n})=$
$\displaystyle\gamma\nabla^{2}\phi^{n+1,j+1}+G_{-}(\phi^{n+1,j},T,c)+G_{+}(\phi^{n},T,c)$
$\displaystyle-\mathcal{L}(\phi^{n+1,j+1}-\phi^{n+1,j}),$ (26)

where $\mathcal{L}\in\mathbb{R}_{+}$.

######  Theorem 1.

The L-scheme iterations (5.3.1) form a contraction and are guaranteed to
converge independently of the initial guess and time-step size when
$\mathcal{L}\geq M_{G}=\max_{0\leq\xi\leq 1}(-\partial_{1}G_{-}(\xi,T,c))$.

###### Proof.

To show when the iterations (5.3.1) converge, we follow similar steps as [24,
Lemma 4.1] and we show the main steps in the following. We first define the
error $e^{j+1}:=\phi^{n+1,j+1}-\phi^{n+1}$. By subtracting (25) from (5.3.1),
we get

$e^{j+1}(1+\Delta t\mathcal{L})-\gamma\Delta t\nabla^{2}e^{j+1}=\Delta
t(G_{-}(\phi^{n+1,j},T,c)-G_{-}(\phi^{n+1},T,c))+\Delta t\mathcal{L}e^{j}.$
(27)

By the mean value theorem we have that
$G_{-}(\phi^{n+1,j},T,c)-G_{-}(\phi^{n+1},T,c)=\partial_{1}G_{-}(\xi,T,c)e^{j}$,
for some number $\xi$ between $\phi^{n+1,j}$ and $\phi^{n+1}$. We define

$M_{G}=\max_{0\leq\xi\leq 1}(-\partial_{1}G_{-}(\xi,T,c)).$ (28)

Since $G_{-}(\phi,T,c)$ is monotonically decreasing with respect to $\phi$,
$\partial_{1}G_{-}(\phi,T,c)$ is either negative or zero; hence, $M_{G}$ is a
positive number. We multiply (27) with $e^{j+1}$ and integrate over the domain
$\Omega$. Using our knowledge about $G_{-}(\phi,T,c)$, and letting
$\mathcal{L}\geq M_{G}$, this can be written as

$\|e^{j+1}\|^{2}(1+\Delta t\mathcal{L})+\gamma\Delta t\|\nabla
e^{j+1}\|^{2}\leq\Delta t\mathcal{L}\|e^{j}\|\|e^{j+1}\|,$

where we in the second term have used integration by parts and on the right-
hand side applied the Cauchy-Schwartz inequality as well as that
$0\leq\partial_{1}G_{-}(\xi,T,c)+\mathcal{L}\leq\mathcal{L}$. The norms are
$L^{2}$-norms. The right-hand side can be further rewritten using Young’s
inequality, hence

$\|e^{j+1}\|^{2}(1+\Delta t\mathcal{L})+\gamma\Delta t\|\nabla
e^{j+1}\|^{2}\leq\frac{1}{2}\Delta
t\mathcal{L}(\|e^{j}\|^{2}+\|e^{j+1}\|^{2}).$

By collecting all terms with $e^{j+1}$ on the left-hand side, we arrive at

$\|e^{j+1}\|^{2}+\frac{2\gamma\Delta t}{2+\Delta t\mathcal{L}}\|\nabla
e^{j+1}\|^{2}\leq\frac{\Delta t\mathcal{L}}{2+\Delta
t\mathcal{L}}\|e^{j}\|^{2}.$

From this we can conclude that the L-scheme (5.3.1) is a contraction and hence
converges by the Banach fixed point theorem when

$\mathcal{L}\geq M_{G}.$ (29)

∎

Note that (29) is not a strict limit, and convergence can also be achieved
under milder constrictions on $\mathcal{L}$, but will then have a restriction
on $\Delta t$. This has been exemplified for e.g. the Richards equation [27].
However, under the above constraints convergence is guaranteed, for any
starting point. Note that the factor

$L=\frac{\Delta t\mathcal{L}}{2+\Delta t\mathcal{L}}$ (30)

can estimate how fast the convergence is. A smaller value of $L$ (closer to 0)
means that convergence is expected to be faster. Hence, when $\mathcal{L}$
and/or $\Delta t$ grows, convergence is expected to be slower as $L$ would
approach 1. This is in contrast to non-linear equations where the non-
linearity is occurring in the time derivative. In this case, the L-scheme
iterations are expected to be faster for larger time-step sizes $\Delta t$
[20, Remark 4].

In practice, the L-scheme iterations are performed until

$\|\phi^{n+1,j+1}-\phi^{n+1,j}\|\leq\text{tol}_{\mathcal{L}},$ (31)

where $\text{tol}_{\mathcal{L}}$ is a prescribed tolerance value. The newest
$\phi^{n+1,j+1}$ is taken as the solution $\phi^{n+1}$ and the simulation can
proceed.

####  5.3.2 Determining $M_{G}$ and $G_{-}(\phi)$

A relevant question is of course how to determine $M_{G}$ and $G_{-}(\phi)$ to
perform the L-scheme iterations. As only an upper bound for $M_{G}$ is needed,
we start with this. We search for $M_{G}$ such that

$M_{G}=\max_{0\leq\phi\leq 1}|\partial_{1}G(\phi,T,c)|.$

This way, $M_{G}$ is potentially overestimated as the positive values of
$\partial_{1}G(\phi,T,c)$ are also included, but this $M_{G}$ will clearly
also fulfill (28). Differentiating $G(\phi,T,c)$ with respect to $\phi$ gives

$\partial_{1}G(\phi,T,c)=\frac{\gamma}{\lambda^{2}}16\Big{(}-(1-6\phi+6\phi^{2})+\frac{1}{|\Omega|}\int_{\Omega}(1-6\phi+6\phi^{2})dV\Big{)}-\frac{4}{\lambda}\frac{1}{m_{m}}f(T,c)(1-2\phi).$

Since the average of 1 is 1, this can be further simplified to

$\partial_{1}G(\phi,T,c)=\frac{96\gamma}{\lambda^{2}}\Big{(}\phi(1-\phi)-\frac{1}{|\Omega|}\int_{\Omega}\phi(1-\phi)dV\Big{)}-\frac{4}{\lambda}\frac{1}{m_{m}}f(T,c)(1-2\phi)$
(32)

We then take the maximum over $0\leq\phi\leq 1$:

$\displaystyle M_{G}$ $\displaystyle=\max_{0\leq\phi\leq
1}|\partial_{1}G(\phi,T,c)|$ $\displaystyle\leq\max_{0\leq\phi\leq
1}\frac{96\gamma}{\lambda^{2}}|\phi(1-\phi)-\frac{1}{|\Omega|}\int_{\Omega}\phi(1-\phi)dV|+\frac{4}{\lambda}\frac{1}{m_{m}}|f(T,c)|$
$\displaystyle\leq\frac{24\gamma}{\lambda^{2}}+\frac{4}{\lambda}\frac{1}{m_{m}}|f(T,c)|.$

Hence, the value of $M_{G}$ can be seen as an interplay between two processes:
One due to changes in $\phi$ coming from curvature-driven motion, and one
coming from chemical reactions. The latter could be zero, but $M_{G}$ will
always be larger than zero since $\gamma,\lambda>0$.

The equation (32) can also be used to determine $G_{-}(\phi,T,c)$. Since
$G_{-}(\phi,T,c)$ simply consists of the part of $G(\phi,T,c)$ where
$G(\phi,T,c)$ is decreasing with respect to $\phi$, the sign of
$\partial_{1}G(\phi,T,c)$ in (32) gives us the increasing and decreasing
parts. However, we have to make sure that a splitting of $G$ into $G_{-}$ and
$G_{+}$ and corresponding implicit/explicit time stepping still result in a
conservative scheme.

####  5.3.3 Strategy for L-scheme iterations for the conservative Allen-Cahn
equation

As pointed out in Section 5.2, the discrete version (22) of the conservative
Allen-Cahn equation (8) is still conservative as long as the same time step
$\ell$ are used for same $\phi_{K}^{\ell}$ in both
$P^{\prime}(\phi_{K}^{\ell})$ and in
$\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi_{J}^{\ell})$, although $\ell=n$ or
$\ell=n+1$ could be chosen independently for each control volume. Hence, when
splitting the non-linearities into an increasing and decreasing part, this has
to be done for control volume by control volume. We outline this in practice:

At time step $n$, $\phi^{n}_{K}$ is known for all $K\in\mathcal{E}$, and we
use the space-discrete version of (32), namely

$\displaystyle\partial_{1}G(\phi_{K};\\{\phi\\}_{J\in\mathcal{E}},T_{K},c_{K})=$
$\displaystyle\frac{96\gamma}{\lambda^{2}}\Big{(}\phi_{K}(1-\phi_{K})-\frac{1}{|\Omega|}\sum_{J\in\mathcal{E}}|J|\phi_{J}(1-\phi_{J})\Big{)}$
$\displaystyle-\frac{4}{\lambda}\frac{1}{m_{m}}f(T_{K},c_{K})(1-2\phi_{K}),$

to determine whether $\phi_{K}$ is causing
$G(\phi_{K};\\{\phi\\}_{J\in\mathcal{E}},T_{K},c_{K})$ to be increasing or
decreasing. Hence, depending on the sign, $\phi_{K}$ will either be solved
explicitly or implicitly (that is, either part of $G_{+}$ or $G_{-}$), and is
assigned to time step $n$ or $n+1$. We denote this as $\phi_{K}^{\ell_{K}}$ to
denote that $\ell_{K}$ is $n$ or $n+1$ in an element-wise manner. When all
control volumes are checked and assigned to being either $n$ or $n+1$ in the
non-linearity, the discrete L-scheme iterations are:

$\displaystyle\frac{1}{\Delta t}(\phi^{n+1,j+1}_{K}-\phi^{n}_{K})=$
$\displaystyle\frac{\gamma}{|K|}\sum_{L\in\mathcal{N}(K)}|\sigma_{K,L}|F^{\phi;n+1,j+1}_{K,L}+G_{-}(\phi^{n+1,j}_{K},T_{K},c_{K})$
$\displaystyle+G_{+}(\phi^{n}_{K},T_{K},c_{K})-\mathcal{L}(\phi^{n+1,j+1}_{K}-\phi^{n+1,j}_{K}),$
(33)

where $G_{-}(\phi_{K}^{n+1,j},T_{K},c_{K})$ is
$G(\phi_{K}^{n+1,j};\\{\phi^{\ell_{J}}\\}_{J\in\mathcal{E}},T_{K},c_{K})$ when
$\phi_{K}^{n}$ was assigned to time step $n+1$, and
$G_{+}(\phi_{K}^{n},T_{K},c_{K})$ is
$G(\phi_{K}^{n};\\{\phi^{\ell_{J}}\\}_{J\in\mathcal{E}},T_{K},c_{K})$ when
$\phi_{K}^{n}$ was assigned to time step $n$. Note that the status of each
element is updated at every iteration to ensure that the correct splitting of
the non-linearity is performed. The discrete version of (24) accounting for
the splitting is

$\displaystyle
G(\phi_{K}^{\ell_{K}};\\{\phi^{\ell_{J}}\\}_{J\in\mathcal{E}},T_{K},c_{K})=$
$\displaystyle-\frac{\gamma}{\lambda^{2}}P^{\prime}(\phi_{K}^{\ell_{K}})-\frac{4}{\lambda}\phi_{K}^{\ell_{K}}(1-\phi_{K}^{\ell_{K}})\frac{f(T_{K},c_{K})}{m_{m}}$
$\displaystyle+\frac{\gamma}{\lambda^{2}}\frac{1}{|\Omega|}\sum_{J\in\mathcal{E}}|J|P^{\prime}(\phi_{J}^{\ell_{J}})$

where $\ell_{K}$ (and $\ell_{J}$) is either $n$ or $n+1,j$. Note that also
$G_{+}(\phi_{K}^{n},T_{K},c_{K})$ is updated throughout the L-scheme
iterations, since the sum relying on $\\{\phi^{\ell_{J}}\\}_{J\in\mathcal{E}}$
is updated. Choosing the stabilization parameter $\mathcal{L}$ according to
$M_{G}$, the scheme (5.3.3) converges to a solution that conserves the phase
field. Note however, that since we in practice stop the L-scheme iterations
when the error is within a certain threshold, see (31), we will therefore
expect an error in terms of conserving the phase field within that same
threshold.

###  5.4 Coupling the phase-field equation to the other model equations

The conservative Allen-Cahn equation (8) is coupled to (9) describing flow and
transport of solute and heat. The equations are tightly coupled: the reaction
rate in (8) depends on the solute concentration and temperature, while the
phase field appears in all equations in (9). We here consider an iterative
scheme, where at every time step each equation (8) and (9) is solved
separately, and we iterate between the equations until convergence before
continuing to the next time step. The benefit of such a scheme is that each
equation can be solved using separate strategies, according to the properties
of each equation. For simplicity we here consider a case without flow and that
$m_{f}=2m_{m}$. Hence, solute and heat is transported with
diffusion/conduction only, but where the diffusion and conduction through the
fluid and mineral depends on the value of the phase field.

The solution strategy is as following, and similar to the strategies applied
in [3, 9]. We discretize in time using a time-step size $\Delta t$, and let
$\phi^{n+1},c^{n+1},T^{n+1}$ denote the phase field, concentration and
temperature at the next time step. The time-discrete equations, when splitting
the non-linearity in an explicit/implicit manner as needed for the phase field
and solving the other equations fully implicitly, are

$\displaystyle\frac{1}{\Delta t}(\phi^{n+1}-\phi^{n})$
$\displaystyle=\gamma\nabla^{2}\phi^{n+1}+G_{-}(\phi^{n+1},T^{n+1},c^{n+1})$
$\displaystyle\quad+G_{+}(\phi^{n},T^{n+1},c^{n+1}),$ (34a)
$\displaystyle\frac{1}{\Delta t}(\phi^{n+1}c^{n+1}$
$\displaystyle+(1-\phi^{n+1})m_{m}-\phi^{n}c^{n}-(1-\phi^{n})m_{m})$
$\displaystyle=D\nabla(\phi^{n+1}\nabla c^{n+1})$ (34b)
$\displaystyle\frac{1}{\Delta
t}(c_{\phi}^{n+1}\rho_{\phi}^{n+1}T^{n+1}-c_{\phi}^{n}\rho_{\phi}^{n}T^{n})$
$\displaystyle=\nabla\cdot(k_{\phi}^{n+1}\nabla T^{n+1}).$ (34c)

We will solve these equations in an iterative manner. We let $i$ be the
iteration index for the coupling iterations. The logic initialization of the
coupling iterations is to take $\phi^{n+1,0}=\phi^{n}$, $c^{n+1,0}=c^{n}$ and
$T^{n+1,0}=T^{n}$, where $\phi^{n},c^{n},T^{n}$ are the (already found)
solutions from the previous time step. However, the iterations converge for
any initial guess, which will be shown in proof of Theorem 2. We first find
the phase field $\phi^{n+1,i+1}$ by solving

$\displaystyle\frac{1}{\Delta t}$
$\displaystyle(\phi^{n+1,i+1}-\phi^{n})=\gamma\nabla^{2}\phi^{n+1,i+1}+G_{-}(\phi^{n+1,i+1},T^{n+1,i},c^{n+1,i})$
$\displaystyle+G_{+}(\phi^{n},T^{n+1,i},c^{n+1,i})-\mathcal{L}_{\text{coup}}(\phi^{n+1,i+1}-\phi^{n+1,i}),$
(35)

where the last term has been added as a stabilization term. This term helps us
to ensure that the iterations are indeed converging and is needed for the
proof of Theorem 2. Note that the reaction rate is using the iteration index
from the previous iteration for temperature and solute concentration, which
are known. Hence, $\phi^{n+1,i+1}$ is the only unknown. This equation is still
non-linear and will be solved with L-scheme iterations. See Remark 3 on how
the L-scheme iterations are modified due to the presence of the newly added
stabilization term. After solving for $\phi^{n+1,i+1}$, we find $c^{n+1,i+1}$
and $T^{n+1,i+1}$ by solving

$\displaystyle\frac{1}{\Delta t}$
$\displaystyle(\phi^{n+1,i+1}c^{n+1,i+1}+(1-\phi^{n+1,i+1})m_{m}-\phi^{n}c^{n}-(1-\phi^{n})m_{m})$
$\displaystyle=D\nabla(\phi^{n+1,i+1}\nabla c^{n+1,i+1})$ (36a)
$\displaystyle\frac{1}{\Delta t}$
$\displaystyle((c_{p}\rho)_{\phi}^{n+1,i+1}T^{n+1,i+1}-(c_{p}\rho)_{\phi}^{n}T^{n})=\nabla\cdot(k_{\phi}^{n+1,i+1}\nabla
T^{n+1,i+1}).$ (36b)

Recall that $(c_{p}\rho)_{\phi},k_{\phi}$ depend on $\phi$ and are here using
$\phi^{n+1,i+1}$. Since $\phi^{n+1,i+1}$ is known, these two equations can be
solved for $c^{n+1,i+1},T^{n+1,i+1}$. Both these equations are linear, hence
they can be solved directly when applying a spatial discretization. Overall
the scheme is as illustrated in Figure 1.

######  Remark 3.

Note that the L-scheme iterations are slightly altered as the term
$\mathcal{L}_{\text{coup}}(\phi^{n+1,i+1}-\phi^{n+1,i})$ now appears in (5.4).
In the resulting L-scheme iterations, we will at every time step $n+1$ and at
every coupling iteration $i+1$, solve for $\phi^{n+1,i+1,j+1}$

$\displaystyle\frac{1}{\Delta t}$
$\displaystyle(\phi^{n+1,i+1,j+1}-\phi^{n})=\gamma\nabla^{2}\phi^{n+1,i+1,j+1}$
$\displaystyle+G_{-}(\phi^{n+1,i+1,j},T^{n+1,i},c^{n+1,i})+G_{+}(\phi^{n},T^{n+1,i},c^{n+1,i})$
$\displaystyle-\mathcal{L}(\phi^{n+1,i+1,j+1}-\phi^{n+1,j})-\mathcal{L}_{\text{coup}}(\phi^{n+1,i+1,j+1}-\phi^{n+1,i}),$
(37)

where $G_{-}$ and $G_{+}$ are the same as before. Performing the same steps as
in the proof of Theorem 1 we now obtain

$\|e^{j+1}\|^{2}+\frac{2\gamma\Delta t}{2+\Delta t\mathcal{L}+\Delta
t\mathcal{L}_{\text{coup}}}\|\nabla e^{j+1}\|^{2}\leq\frac{\Delta
t\mathcal{L}}{2+\Delta t\mathcal{L}+\Delta
t\mathcal{L}_{\text{coup}}}\|e^{j}\|^{2}.$

Hence, we still obtain a contraction for any $\mathcal{L}_{\text{coup}}>0$ and
assuming $\mathcal{L}\geq M_{G}$. In fact, the factor $L=\frac{\Delta
t\mathcal{L}}{2+\Delta t\mathcal{L}+\Delta t\mathcal{L}_{\text{coup}}}$ is now
smaller compared to (30) and we can expect slightly faster convergence with
the presence of this new stabilization term. The L-scheme iterations are made
until

$\|\phi^{n+1,i+1,j+1}-\phi^{n+1,i+1,j}\|\leq\text{tol}_{\mathcal{L}},$

and the newest $\phi^{n+1,i+1,j+1}$ is taken as $\phi^{n+1,i+1}$; that is, the
phase field of the current coupling iteration.

Figure 1: Scheme of coupling iterations

To prove convergence of the coupling iterations, we rely on the following
assumptions:

  * (A1)

The phase field must be bounded away from 0 and 1.

  * (A2)

The solute concentration and temperature must be bounded.

  * (A3)

The gradient of the solute concentration and temperature must be bounded at
every time step; that is, $\|\nabla c^{n}\|_{L^{\infty}(\Omega)}\leq C_{c}$
and $\|\nabla T^{n}\|_{L^{\infty}(\Omega)}\leq C_{T}$ for constants
$C_{c},C_{T}>0$ and for all time steps $n$.

The original Allen-Cahn equation can be shown to fulfill a maximum principle,
also in the time-discrete form including stabilization terms, see [3, Lemma
1]. Including the non-local term does not alter this, hence also the
conservative Allen-Cahn equation (5.4) results in a phase field which is
between 0 and 1. For a bounded domain and $\lambda>0$, the phase field will be
bounded away from 0 and 1. Hence, (A1) is fulfilled. Additionally, [3, Lemma
3] proves the maximum principle for the solute concentration, using the same
equation for the solute as considered here. The proof relies on the constraint

$4\gamma\leq\frac{\lambda k_{0}}{m_{m}}$ (38)

which we will obey as well when simulating the coupled model. The heat
equation (36b) has the same structure as the solute equation (36a), and hence
can be shown to be bounded under same arguments as for the solute
concentration. Hence, (A2) is fulfilled when the constraint is obeyed.
Finally, (A3) is a reasonable assumption given the diffuse nature of the
equations for solute concentration and temperature (36). However, rigorous
proofs of the assumptions (A1)-(A3) are beyond the scope of this manuscript.

######  Remark 4.

Note that to obtain guarantees that the phase field is bounded away from 0 and
1, one can also regularize the phase-field equation with some value
$\delta>0$, similar as has been done in [3, 8]. In this case one rephrases the
equation to using $\phi_{\delta}$, where $0<\delta\leq\phi_{\delta}\leq
1-\delta<1$.

######  Theorem 2.

Under the assumptions (A1)-(A3), the coupling iterations (5.4), (36) form a
contraction and will hence converge for any initial guess when
$\mathcal{L}_{\text{coup}}$ fulfills (45) and the time-step size fulfills
(44).

###### Proof.

We define the error terms

$e^{i+1}_{\phi}:=\phi^{n+1,i+1}-\phi^{n+1},\quad
e^{i+1}_{c}:=c^{n+1,i+1}-c^{n+1},\quad e^{i+1}_{T}:=T^{n+1,i+1}-T^{n+1}.$

We consider first the Allen-Cahn equation. By subtracting (5.4) from the
corresponding time-discrete form without decoupling, equation (34a), we obtain

$\displaystyle\frac{1}{\Delta t}(\phi^{n+1,i+1}-\phi^{n+1})=$
$\displaystyle\gamma\nabla^{2}(\phi^{n+1,i+1}-\phi^{n+1})+G_{-}(\phi^{n+1,i+1},T^{n+1,i},c^{n+1,i})$
$\displaystyle+G_{+}(\phi^{n},T^{n+1,i},c^{n+1,i})-G_{-}(\phi^{n+1},T^{n+1},c^{n+1})$
$\displaystyle-
G_{+}(\phi^{n},T^{n+1},c^{n+1})-\mathcal{L}_{\text{coup}}(\phi^{n+1,i+1}-\phi^{n+1,i}).$

We multiply the above equation with a $\psi\in H^{1}_{0}(\Omega)$ and
integrate over the domain $\Omega$. We then obtain

$\displaystyle\frac{1}{\Delta t}(e^{i+1}_{\phi},\psi)=$
$\displaystyle-\gamma(\nabla
e^{i+1}_{\phi},\nabla\psi)+(G_{-}(\phi^{n+1,i+1},T^{n+1,i},c^{n+1,i}),\psi)$
$\displaystyle+(G_{+}(\phi^{n},T^{n+1,i},c^{n+1,i}),\psi)-(G_{-}(\phi^{n+1},T^{n+1},c^{n+1}),\psi)$
$\displaystyle-(G_{+}(\phi^{n},T^{n+1},c^{n+1}),\psi)-\mathcal{L}_{\text{coup}}(e^{i+1}_{\phi}-e^{i}_{\phi},\psi),$
(39)

where the notation $(\cdot,\cdot)$ refers to the usual $L^{2}$-inner product.
We then let $\psi=e^{i+1}_{\phi}$, which gives us

$\displaystyle\|e^{i+1}_{\phi}\|^{2}$ $\displaystyle+\Delta t\gamma\|\nabla
e^{i+1}_{\phi}\|^{2}+\Delta
t\mathcal{L}_{\text{coup}}\|e^{i+1}_{\phi}\|^{2}=\Delta
t\mathcal{L}_{\text{coup}}(e^{i}_{\phi},e^{i+1}_{\phi})$ $\displaystyle+\Delta
t(G_{-}(\phi^{n+1,i+1},T^{n+1,i},c^{n+1,i}),e^{i+1}_{\phi})+\Delta
t(G_{+}(\phi^{n},T^{n+1,i},c^{n+1,i}),e^{i+1}_{\phi})$ $\displaystyle-\Delta
t(G_{-}(\phi^{n+1},T^{n+1},c^{n+1}),e^{i+1}_{\phi})-\Delta
t(G_{+}(\phi^{n},T^{n+1},c^{n+1}),e^{i+1}_{\phi}),$

We let
$\mathcal{M}_{G}=\max(|\partial_{1}G|,|\partial_{2}G|,|\partial_{3}G|)$, where
we have maximized the partial derivatives of $G$ over the bounded domain where
$\phi,T,c$ are defined (Assumptions (A1) and (A2)). Then, by the mean-value
theorem, we get

$\displaystyle(1+\Delta t\mathcal{L}_{\text{coup}})\|e^{i+1}_{\phi}\|^{2}$
$\displaystyle+\Delta t\gamma\|\nabla e^{i+1}_{\phi}\|^{2}\leq\Delta
t\mathcal{L}_{\text{coup}}(e^{i}_{\phi},e^{i+1}_{\phi})$ $\displaystyle+\Delta
t\mathcal{M}_{G}\|e^{i+1}_{\phi}\|^{2}+2\Delta
t\mathcal{M}_{G}(e^{i}_{T},e^{i+1}_{\phi})+2\Delta
t\mathcal{M}_{G}(e^{i}_{c},e^{i+1}_{\phi}).$

We apply Young’s inequality to the first and the last two terms on the right-
hand side, which, for $\delta_{1},\delta_{2},\delta_{3}>0$, gives

$\displaystyle(1$ $\displaystyle+\Delta t\mathcal{L}_{\text{coup}}-\Delta
t\mathcal{M}_{G})\|e^{i+1}_{\phi}\|^{2}+\Delta t\gamma\|\nabla
e^{i+1}_{\phi}\|^{2}$ $\displaystyle\leq$ $\displaystyle\Delta
t\mathcal{L}_{\text{coup}}\frac{\delta_{1}}{2}\|e^{i}_{\phi}\|^{2}+\Delta
t\mathcal{L}_{\text{coup}}\frac{1}{2\delta_{1}}\|e^{i+1}_{\phi}\|^{2}+\Delta
t\mathcal{M}_{G}\delta_{2}\|e^{i}_{T}\|^{2}+\Delta
t\mathcal{M}_{G}\frac{1}{\delta_{2}}\|e^{i+1}_{\phi}\|^{2}$
$\displaystyle+\Delta t\mathcal{M}_{G}\delta_{3}\|e^{i}_{c}\|^{2}+\Delta
t\mathcal{M}_{G}\frac{1}{\delta_{3}}\|e^{i+1}_{\phi}\|^{2}.$

We collect all terms with $\|e^{i+1}_{\phi}\|$ on the left-hand side and hence
have

$\displaystyle(1+\Delta t\mathcal{L}_{\text{coup}}-\Delta
t\mathcal{M}_{G}-\Delta t\mathcal{L}_{\text{coup}}\frac{1}{2\delta_{1}}-\Delta
t\mathcal{M}_{G}\frac{1}{\delta_{2}}-\Delta
t\mathcal{M}_{G}\frac{1}{\delta_{3}})\|e^{i+1}_{\phi}\|^{2}$
$\displaystyle+\Delta t\gamma\|\nabla e^{i+1}_{\phi}\|^{2}\leq\Delta
t\mathcal{L}_{\text{coup}}\frac{\delta_{1}}{2}\|e^{i}_{\phi}\|^{2}+\Delta
t\mathcal{M}_{G}\delta_{2}\|e^{i}_{T}\|^{2}+\Delta
t\mathcal{M}_{G}\delta_{3}\|e^{i}_{c}\|^{2}.$ (40)

To proceed, we need bounds for $\|e^{i}_{T}\|$ and $\|e^{i}_{c}\|$, hence we
turn to the temperature and solute concentration equation.

For the solute concentration, we start by subtracting (36a) from (34b),
multiply with a test function from $H^{1}_{0}(\Omega)$ and integrate over the
domain. Letting $e^{i+1}_{c}$ be the test function, we have

$\displaystyle\frac{1}{\Delta t}$
$\displaystyle(\phi^{n+1,i+1}e^{i+1}_{c},e^{i+1}_{c})+D(\phi^{n+1,i+1}\nabla
c^{n+1,i+1},\nabla e^{i+1}_{c})$ $\displaystyle=D(\phi^{n+1}\nabla
c^{n+1},\nabla e^{i+1}_{c})+\frac{1}{\Delta
t}(e^{i+1}_{\phi}(m_{m}-c^{n+1}),e^{i+1}_{c}).$

We use (A1); that $\phi_{m}\leq\phi$ for some $\phi_{m}>0$. Then the above
equation can be rewritten to

$\displaystyle\frac{\phi_{m}}{\Delta t}$
$\displaystyle\|e^{i+1}_{c}\|^{2}+D\phi_{m}\|\nabla e^{i+1}_{c}\|^{2}$
$\displaystyle\leq D((\phi^{n+1}-\phi^{n+1,i+1})\nabla c^{n+1},\nabla
e^{i+1}_{c})+\frac{1}{\Delta t}(e^{i+1}_{\phi}(m-c^{n+1}),e^{i+1}_{c}).$

We apply Young’s inequality to the two terms on the right-hand side, which,
for $\delta_{4},\delta_{5}>0$, gives

$\displaystyle\frac{\phi_{m}}{\Delta t}\|e^{i+1}_{c}\|^{2}+D\phi_{m}\|\nabla
e^{i+1}_{c}\|^{2}\leq$ $\displaystyle
D\frac{\delta_{4}}{2}\|e^{i+1}_{\phi}\nabla
c^{n+1}\|^{2}+D\frac{1}{2\delta_{4}}\|\nabla e^{i+1}_{c}\|^{2}$
$\displaystyle+\frac{1}{\Delta
t}\frac{\delta_{5}}{2}\|(m-c^{n+1})e^{i+1}_{\phi}\|^{2}+\frac{1}{\Delta
t}\frac{1}{2\delta_{5}}\|e^{i+1}_{c}\|^{2}.$

We then apply (A2) and (A3), which results in

$\displaystyle(\frac{\phi_{m}}{\Delta t}-\frac{1}{\Delta
t}\frac{1}{2\delta_{5}})\|e^{i+1}_{c}\|^{2}+(D\phi_{m}-\frac{D}{2\delta_{4}})\|\nabla
e^{i+1}_{c}\|^{2}\leq(D\frac{\delta_{4}}{2}C_{c}^{2}+\frac{1}{\Delta
t}\frac{\delta_{5}}{2}m_{m}^{2})\|e^{i+1}_{\phi}\|^{2}$

To ensure positive terms on the left-hand side, we let
$\delta_{4}=\delta_{5}=\frac{1}{\phi_{m}}$. We then arrive at

$\|e^{i+1}_{c}\|^{2}\leq\frac{2\Delta
t}{\phi_{m}}(\frac{DC_{c}^{2}}{2\phi_{m}}+\frac{m_{m}^{2}}{2\phi_{m}\Delta
t})\|e^{i+1}_{\phi}\|^{2}.$

Note that this bound is also valid for index $i$. For simplicity, we rephrase
the result with the more compact notation

$\|e^{i}_{c}\|^{2}\leq E_{c}\|e^{i}_{\phi}\|^{2},$ (41)

for $E_{c}=\frac{2\Delta
t}{\phi_{m}}(\frac{DC_{c}^{2}}{2\phi_{m}}+\frac{m_{m}^{2}}{2\phi_{m}\Delta
t})$.

We perform similar steps for the temperature equation. We start by subtracting
(36b) from (34c), multiply with a test function from $H^{1}_{0}(\Omega)$ and
integrate over the domain. We let $e^{i+1}_{T}$ be the test function and hence
obtain

$\displaystyle\frac{c_{p,f}\rho_{f}}{\Delta t}$
$\displaystyle(\phi^{n+1,i+1}e^{i+1}_{T},e^{i+1}_{T})+\frac{c_{p,m}\rho_{m}}{\Delta
t}((1-\phi^{n+1,i+1})e^{i+1}_{T},e^{i+1}_{T})$
$\displaystyle+((\phi^{n+1,i+1}k_{f}+(1-\phi^{n+1,i+1})k_{m})\nabla
T^{n+1,i+1},\nabla e^{i+1}_{T})$ $\displaystyle=$
$\displaystyle((\phi^{n+1}k_{f}+(1-\phi^{n+1})k_{m})\nabla T^{n+1},\nabla
e^{i+1}_{T})$ $\displaystyle+\frac{1}{\Delta
t}(e^{i+1}_{\phi}(c_{p,m}\rho_{m}-c_{p,f}\rho_{f})T^{n+1},e^{i+1}_{T})$

We use (A1); that $\phi_{m}\leq\phi\leq\phi_{M}$ for $\phi_{M}<1$. Then the
above equation can be rewritten to

$\displaystyle\frac{c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})}{\Delta
t}\|e^{i+1}_{T}\|^{2}+(k_{f}\phi_{m}+k_{m}(1-\phi_{M}))\|\nabla
e^{i+1}_{T}\|^{2}$
$\displaystyle\leq((k_{f}(\phi^{n+1}-\phi^{n+1,i+1})+k_{m}(\phi^{n+1,i+1}-\phi^{n+1}))\nabla
T^{n+1},\nabla e^{i+1}_{T})$ $\displaystyle+\frac{1}{\Delta
t}(e^{i+1}_{\phi}(c_{p,m}\rho_{m}-c_{p,f}\rho_{f})T^{n+1},e^{i+1}_{T}).$

We apply Young’s inequality to the two terms on the right-hand side, which,
for $\delta_{6},\delta_{7}>0$, gives

$\displaystyle\frac{c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})}{\Delta
t}\|e^{i+1}_{T}\|^{2}+(k_{f}\phi_{m}+k_{m}(1-\phi_{M}))\|\nabla
e^{i+1}_{T}\|^{2}$
$\displaystyle\leq\frac{\delta_{6}}{2}\|(-k_{f}e^{i+1}_{\phi}+k_{m}e^{i+1}_{\phi})\nabla
T^{n+1}\|^{2}+\frac{1}{2\delta_{6}}\|\nabla e^{i+1}_{T}\|^{2}$
$\displaystyle+\frac{1}{\Delta
t}\frac{\delta_{7}}{2}\|(c_{p,m}\rho_{m}-c_{p,f}\rho_{f})T^{n+1}e^{i+1}_{\phi}\|^{2}+\frac{1}{\Delta
t}\frac{1}{2\delta_{7}}\|e^{i+1}_{T}\|^{2}.$

We then apply (A2) and (A3), with $T_{M}$ being the upper bound of the
temperature, which results in

$\displaystyle\frac{1}{\Delta t}$
$\displaystyle(c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})-\frac{1}{2\delta_{7}})\|e^{i+1}_{T}\|^{2}$
$\displaystyle+(k_{f}\phi_{m}+(1-\phi_{M})k_{m}-\frac{1}{2\delta_{6}})\|\nabla
e^{i+1}_{T}\|^{2}$
$\displaystyle\leq(\frac{\delta_{6}}{2}|k_{m}-k_{f}|^{2}C_{T}^{2}+\frac{1}{\Delta
t}\frac{\delta_{7}}{2}|c_{p,m}\rho_{m}-c_{p,f}\rho_{f}|^{2}T_{M}^{2})\|e^{i+1}_{\phi}\|^{2}$

To ensure positive terms on the left-hand side, we let
$\delta_{6}=\frac{1}{k_{f}\phi_{m}+k_{m}(1-\phi_{M})}$ and
$\delta_{7}=\frac{1}{c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})}$. We
then arrive at

$\displaystyle\frac{c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})}{2\Delta
t}\|e^{i+1}_{T}\|^{2}$
$\displaystyle\leq(\frac{\delta_{6}}{2}|k_{m}-k_{f}|^{2}C_{T}^{2}+\frac{1}{\Delta
t}\frac{\delta_{7}}{2}|c_{p,m}\rho_{m}-c_{p,f}\rho_{f}|^{2}T_{M}^{2})\|e^{i+1}_{\phi}\|^{2},$

where we, to keep the notation more compact, did not insert the chosen values
for $\delta_{6},\delta_{7}$ on the right-hand side. This bound is also valid
for index $i$. For simplicity, we rephrase the result with the more compact
notation

$\|e^{i}_{T}\|^{2}\leq E_{T}\|e^{i}_{\phi}\|^{2},$ (42)

for $E_{T}=\frac{2\Delta
t}{c_{p,f}\rho_{f}\phi_{m}+c_{p,m}\rho_{m}(1-\phi_{M})}(\frac{\delta_{6}}{2}|k_{m}-k_{f}|^{2}C_{T}^{2}+\frac{1}{\Delta
t}\frac{\delta_{7}}{2}|c_{p,m}\rho_{m}-c_{p,f}\rho_{f}|^{2}T_{M}^{2})$.

We then return to (40) and insert here (41) and (42). At the same time we let
$\delta_{1}=1$ and $\delta_{2}=\delta_{3}=\frac{1}{2}$. Hence,

$\displaystyle(1+$ $\displaystyle\frac{\Delta
t}{2}\mathcal{L}_{\text{coup}}-5\Delta
t\mathcal{M}_{G})\|e^{i+1}_{\phi}\|^{2}$ $\displaystyle\leq(\frac{\Delta
t}{2}\mathcal{L}_{\text{coup}}+\frac{\Delta
t}{2}\mathcal{M}_{G}E_{T}+\frac{\Delta
t}{2}\mathcal{M}_{G}E_{c})\|e^{i}_{\phi}\|^{2}.$ (43)

That means the coupling iterations form a contraction when the left-hand side
is positive and when

$1+\frac{\Delta t}{2}\mathcal{L}_{\text{coup}}-5\Delta
t\mathcal{M}_{G}>\frac{\Delta t}{2}\mathcal{L}_{\text{coup}}+\frac{\Delta
t}{2}\mathcal{M}_{G}E_{T}+\frac{\Delta t}{2}\mathcal{M}_{G}E_{c}.$

The latter translates to a restriction on the time-step size, and results in

$\Delta t<\frac{1}{\mathcal{M}_{G}(5+E_{c}/2+E_{T}/2)}.$ (44)

To ensure that the left-hand side of (43) is positive, we need
$\mathcal{L}_{\text{coup}}>10\mathcal{M}_{G}-\frac{2}{\Delta t}$, as well as
$\mathcal{L}_{\text{coup}}>0$. Hence,

$\mathcal{L}_{\text{coup}}>\max\\{0,10\mathcal{M}_{G}-\frac{2}{\Delta t}\\}.$
(45)

∎

Note that the time-step size restriction (44) does not depend on the
discretization mesh. However, it does depend on parameters that cannot easily
be estimated. One possibility would be to use a coarse mesh and by trial-and-
error find a time-step size that ensures convergence, which is then used on a
finer mesh. However, the phase-field equation (5.4) needs a certain mesh size
to resolve the diffuse interface. A coarser mesh can be applied if the
diffuse-interface width $\lambda$ is chosen larger, but changing the value of
$\lambda$ can influence the time-step size restriction. Hence, one would for
these model equations need to find the suitable time-step size on the mesh
that is to be applied.

In practice, the coupling iterations are performed until

$\|\phi^{n+1,i+1}-\phi^{n+1,i}\|\leq\text{tol}_{\text{coup}},$ (46)

where $\text{tol}_{\text{coup}}$ is a prescribed tolerance. The newest
$\phi^{n+1,i+1},c^{n+1,i+1},T^{n+1,i+1}$ is taken as
$\phi^{n+1},c^{n+1},T^{n+1}$ and we can proceed to the next time step. As seen
in the proof of Theorem 2, the errors of the solute concentration and
temperature are expected to be bounded when the error in the phase field is
bounded. Hence, it is sufficient to consider a tolerance based on the phase
field only.

######  Remark 5.

Note that the value of $E_{T}$ in equation (42) is zero if $k_{m}=k_{f}$ and
$c_{p,f}\rho_{f}=c_{p,m}\rho_{m}$. This would correspond to equal heat
conduction and heat storage properties in the fluid and mineral, which implies
that the temperature equation is in practice independent of the phase field. A
similar result can not be obtained for the solute.

##  6 Numerical results of the phase-field models

We here perform numerical experiments to assess the numerical behavior of the
conservative Allen-Cahn equation, in particular when it is solved with the
L-scheme, and to test the full, coupled model to show its numerical behavior
and potential. We investigate behavior of the Allen-Cahn equation only in
Section 6.1, before addressing the behavior of the coupled model in Section
6.2.

###  6.1 Conservative and non-conservative Allen-Cahn

We here consider the simplest test case possible, solving the original phase-
field equation (7) or the conservative version (8). The reaction rate will
initially be set to zero, so only curvature-driven motion can change the shape
(and, potentially, size) of the mineral. We consider two test cases, both to
be considered in a two-dimensional domain $\Omega=(0,1)^{2}$, but with
different initial conditions:

  * (a)

We initialize the phase field as a circular mineral with an approximated
diffuse interface transition at radius 0.3.

  * (b)

We initialize the phase field with ones everywhere except for in a centered
square with side lengths 0.5 where the phase field is set to zero.

The latter case (b) is more challenging as the initial condition violates the
phase field equation and is more bound to numerical error. We solve the two
test cases in three different ways:

  * (i)

Solving (7) fully implicitly, using Newton’s method to solve the non-linear
system of equations.

  * (ii)

Solving (8) fully implicitly, using Newton’s method to solve the non-linear
system of equations.

  * (iii)

Solving (8) semi-implcitly, using L-scheme to solve the non-linear system of
equations as described in section 5.3.

All approaches are solved using the Finite Volume discretization described in
Section 5.1, on a uniform, rectangular grid. The test cases apply homogeneous
Neumann conditions on the boundaries. For test case (a), only minor
adjustments of the phase field diffuse interface is expected for the
conservative approaches (ii) and (iii), as well as numerical error from the
solvers. In test case (b), the curvature-driven motion of (ii) and (iii) will
transform the initial square to a circular mineral. For the non-conservative
phase field (i), we expect the mineral to shrink in both case (a) and (b).

All cases are solved on a $100\times 100$ grid, with $\lambda=0.05$ and
$\gamma=1$. Note that a smaller value of $\gamma$ causes the curvature-driven
motion to be slower, but is here chosen large to highlight the issues with
(7). All equations are time stepped until $t=1$. The two fully-implicit cases
(i) and (ii) are solved using $\Delta t=10^{-4}$, which was needed to ensure
convergence of the Newton iterations, and this time-step size is also applied
for the L-scheme. With the current setup, we have $M_{G}=9.68\times 10^{4}$
and we choose $\mathcal{L}$ equal to this value. Newton iterations and
L-scheme iterations are solved until a threshold of $10^{-13}$ is met in the
$L^{2}$-norm between two subsequent iterations (that is,
$\text{tol}_{\mathcal{L}}=10^{-13}$ in (31)). Since the size of the discrete
system is not too large, the linear system is solved directly.

Table 1 shows the change in volume of the mineral from initial condition until
$t=1$ for the two test cases for all three approaches. Since the phase field
approaches the value 0 in the mineral and 1 outside, the volume is calculated
through

$\text{volume}=\int_{\Omega}(1-\phi)dV.$

Table 1: Volume change of mineral from first to last time step Case/Approach | (i) | (ii) | (iii)
---|---|---|---
(a) | -0.2850 | $2.0134\times 10^{-11}$ | $1.5294\times 10^{-11}$
(b) | -0.25 | $-6.1989\times 10^{-12}$ | $1.4864\times 10^{-11}$

Since all mineral disappeared for approach (i), the volume change corresponds
to the loss of the initial amount of mineral. Note that there is a volume
change due to the choice of the threshold tolerance. The two conservative
approaches (ii) and (iii) show only change in volume that can be connected to
the accuracy of the non-linear iterations. For approaches (ii) and (iii), at
every time step a small volume change of about $1.5\times 10^{-15}$ is seen,
which is the same as the tolerance of the Newton and L-scheme iterations.

Table 2 shows average number of Newton or L-scheme iterations needed for each
time step. Note however that in general a larger amount of iterations are
needed in the beginning of the simulation and then stabilizes at a lower
value. At least one Newton or L-scheme iteration is forced to be performed at
every time step.

Table 2: Average number of iterations for entire simulation Case/Approach | (i) | (ii) | (iii)
---|---|---|---
(a) | 1.03 | 1.0038 | 5.18
(b) | 1.03 | 1.0160 | 5.77

For the approach (i), the entire mineral dissolves during the first 300 time
steps and hence attains a constant solution with fluid only afterwards. When
using the conservative scheme with Newton iterations; that is, approach (ii),
2-3 iterations are needed per time step in the beginning, which settles after
150 time steps. The Jacobian in the Newton iterations is however full due to
the non-local term, which makes each Newton iteration expensive to solve. For
approach (iii), several L-scheme iterations are needed in the first phase
while the diffuse interface (and for the square, also the shape) is adjusting.
In the beginning 30-40 iterations are needed, before gradually decreasing to 3
L-scheme iterations after 1000 time steps, where it remains stable. However,
the iterations are cheap as only a sparse linear system of equations needs to
be solved. When solved on the same laptop and same software, approach (iii)
was roughly a factor 2 faster than approach (ii).

Although the L-scheme iterations converge, they converge slowly. The size of
$\mathcal{L}$ is responsible for this. Although we chose it according to the
theoretical bounds (29), earlier studies of other equations have shown (see
e.g. [20, 27, 28]) that much lower values of $\mathcal{L}$ are possible, and
will generally give faster convergence. Also larger values of $\Delta t$ could
be used, but is for our equation expected to result in slower convergence for
the L-scheme. Also note that lower values of $\mathcal{L}$ would usually lead
to a time-step size restriction to obtain convergence. Hence, in the following
two subsection we investigate numerically possible values of $\mathcal{L}$ and
$\Delta t$.

####  6.1.1 Size of $\mathcal{L}$

We change the setup from the previous test slightly to not only consider semi-
steady-state solutions. Initially we have a circular mineral of radius 0.3 as
in (a), but this time with a reaction rate $f=-0.1$. Hence the mineral
gradually dissolves, but is still present at $t=1$. We use $\gamma=0.1$ to
have a more typical value for the curvature-driven motion. The value of
$M_{G}$ is hence also changed and is $M_{G}=968$ now. We will for now fix
$\Delta t=10^{-3}$ and vary $\mathcal{L}$. Table 3 shows the average number of
iterations in each time step as a function of $\mathcal{L}$. Note that the
numbers are generally larger than in Table 2 as no steady-state is reached
during the simulation. In this setup, the number of L-scheme iterations per
time step is almost constant throughout the simulation. We see how
significantly lower values than the theoretical bound of $\mathcal{L}$ still
provide convergence, and that the amount of iterations needed decreases when
$\mathcal{L}$ decreases, before eventually increasing. Hence, there appears to
be an optimal $\mathcal{L}$ for when the lowest amount of needed L-scheme
iterations is obtained. For other equations, it has been shown that a lower
value of $\mathcal{L}$ provides faster convergence (i.e., less iterations
needed), and that there is an optimal $\mathcal{L}$ where the least amount of
iterations are needed [27], and $\mathcal{L}$ can be chosen differently from
time step to time step to optimize convergence [28]. There is also a lower
limit for $\mathcal{L}$ which can still ensure convergence if the time-step
size is also large [27].

Table 3: Average number of L-scheme iterations per time step for different values of $\mathcal{L}$ Value of $\mathcal{L}$ | $M_{G}$ | $M_{G}/2$ | $M_{G}/4$ | $M_{G}/8$
---|---|---|---|---
Average number of iterations | 32 | 21 | 17 | 24

####  6.1.2 Size of $\Delta t$

We now use either the theoretical $\mathcal{L}=M_{G}$ or the
$\mathcal{L}=M_{G}/4$ as this gave the lowest value of iterations, and vary
the time-step size $\Delta t$. The setup is otherwise the same as in Section
6.1.1. Table 4 and 5 show the average number of iterations as a function of
the time-step size for these two choices of $\mathcal{L}$. As expected for our
equation, the number of iterations needed at every time step increases with
the time-step size $\Delta t$. One case is deemed to not converge as the
threshold value was not reached within 200 L-scheme iterations per time step.
Note that convergence potentially could be achieved by allowing more
iterations. Note that for $\mathcal{L}=M_{G}$, the number of L-scheme
iterations increases but is less than doubled as the time-step size is
doubled. Hence, in this case there is a computational gain by increasing the
time-step size as the total number of L-scheme iterations throughout the
simulation decreases. This is however not necessarily the case for
$\mathcal{L}=M_{G}/4$ as the average number of L-scheme iterations is more
than doubled as the time-step size increases from $2\times 10^{-3}$ to
$4\times 10^{-3}$.

Table 4: Average number of L-scheme iterations per time step for different values of $\Delta t$, when $\mathcal{L}=M_{G}$ Value of $\Delta t$ | $10^{-3}$ | $2\times 10^{-3}$ | $4\times 10^{-3}$ | $8\times 10^{-3}$
---|---|---|---|---
Average number of iterations | 32 | 53 | 95 | 175
Table 5: Average number of L-scheme iterations per time step for different values of $\Delta t$, when $\mathcal{L}=M_{G}/4$ Value of $\Delta t$ | $10^{-3}$ | $2\times 10^{-3}$ | $4\times 10^{-3}$ | $8\times 10^{-3}$
---|---|---|---|---
Average number of iterations | 17 | 30 | 83 | -

###  6.2 Behavior of coupled model

We here consider the coupled model formulated in Section 5.4. As example case
for the simulations, we consider a channel $\Omega=(0,1)^{2}$, where a mineral
layer of initial thickness $0.25$ fills the lower part of the domain.
Temperature is initially 1 everywhere and the fluid has a solute concentration
of $0.5$, which corresponds to chemical equilibrium of the chosen reaction
rate. A Dirichlet boundary condition of $T=0.9$ and $c=0.25$ is applied to the
left boundary, while all other boundaries have zero Neumann boundary
conditions. The phase field has zero Neumann boundary conditions on all
boundaries. The lowering of the solute concentration (and temperature) at the
left boundary is expected to trigger dissolution of the mineral and hence an
increase in regions where the phase field approaches one. The chosen reaction
rate is a simpler version of (5) with constant solubility product. Hence,

$f(T,c)=f_{p}(T,c)-f_{d}(T,c)=k_{0}e^{-\frac{E}{RT}}\Big{(}\frac{c^{2}}{c_{\text{eq}}^{2}}-1\Big{)}.$

We here consider for simplicity a constant equilibrium concentration
$c_{\text{eq}}=0.5$, but one could consider also temperature-dependent
solubility [32]. All parameters appearing in the model equations are as
specified in Table 6. We have not attempted to choose physically correct
parameters, but generally chosen values around 1 to mimic a non-dimensional
setup. As we need to fulfill (38), we have adjusted $\gamma$. For the heat
conductivities we choose a larger heat conductivity in the mineral to be able
to observe some heterogeneous influence. The simulation setup is sketched in
Figure 2.

Table 6: Parameter choices Parameter | Value
---|---
$\lambda$ | 0.05
$\gamma$ | 0.01
$E/R$ | 1
$k_{0}$ | 1
$c_{\text{eq}}$ | 0.5
$m_{m}$ | 1
$D$ | 1
$\rho_{f}$ | 1
$c_{p,f}$ | 1
$\rho_{m}$ | 1
$c_{p,f}$ | 1
$k_{f}$ | 1
$k_{m}$ | 2
\begin{overpic}[width=303.53267pt]{initial} \put(38.0,42.0){$T^{0}=1$}
\put(38.0,37.0){$c^{0}=0.5$} \put(13.0,24.0){\rotatebox{90.0}{$T=0.9,c=0.25$}}
\put(20.0,71.0){$\partial_{n}\phi=0,\partial_{n}T=0,\partial_{n}c=0$}
\end{overpic} Figure 2: Initial phase field $\phi^{0}$ (background color),
with initial temperature and solute concentration in text. Left boundary
(dashed line) has Dirichlet boundary conditions for temperature and solute,
while all other boundaries (solid line) have homogeneous Neumann boundary
conditions. The phase field has homogeneous Neumann boundary conditions on the
left boundary.

####  6.2.1 Size of time-step size, tolerances and stabilization parameters

We here apply a slightly larger tolerance for the L-scheme iterations than in
Section 6.1, but keep it lower than the tolerance for the coupling iterations
to ensure that errors from the non-linear solver scheme does not influence the
coupling iterations error. We therefore use $\text{tol}_{\mathcal{L}}=10^{-8}$
and $\text{tol}_{\text{coup}}=10^{-6}$ throughout all simulations presented in
the following. Lower tolerances generally lead to more iterations.

With our parameter choices, we have $M_{G}=118$, and we choose
$\mathcal{L}=M_{G}$ in all simulations. For the coupling parameter
$\mathcal{L}_{\text{coup}}$ we choose it small as this provide faster
convergence for the coupling iterations. This result coincides with similar
investigations (e.g. [3, 9]). We here present results with
$\mathcal{L}_{\text{coup}}=10^{-4}$. Although we may not fulfill (45),
numerical testing show that the coupling iterations converge for any
$\mathcal{L}_{\text{coup}}>0$, and even with $\mathcal{L}_{\text{coup}}=0$.

By testing various time-step sizes, we found time-step sizes even up to
$10^{-2}$ to provide converging coupling iterations. However, such long time-
step sizes also resulted in very slowly converging L-scheme iterations, with
generally more than 200 L-scheme iterations per coupling iteration needed.
Hence, in the following, only results from shorter time-step sizes are
presented. Although we cannot guarantee that (44) (and (45)) is fulfilled, the
coupling iterations were always found to converge fast for the investigated
parameters.

####  6.2.2 Behavior of coupling iterations

For $\Delta t=10^{-3}$, $\mathcal{L}_{\text{coup}}=10^{-4}$ and
$\mathcal{L}=M_{G}=118$, we need an average of 2.21 coupling iterations per
time step. More iterations are needed initially as the changes in solute
concentration and temperature, and hence in phase field, are larger due to the
boundary condition on the left boundary. For the L-scheme iterations, we need
6.00, 3.22 and 2.08 iterations in the first, second and third coupling
iterations, respectively. See Figure 3 for the evolution of iterations over
time. Fewer L-scheme iterations are needed in later coupling iterations as we
use the solution found in the previous coupling iteration when starting the
L-scheme in later coupling iterations, giving a starting point that is already
closer to the converged solution. Compared to Section 6.1, fewer L-scheme
iterations are needed for same time-step size as the tolerance for the
L-scheme iterations is now larger.

Figure 3: Number of coupling iterations (left) and number of L-scheme
iterations per coupling iteration (right) over time for $\Delta t=10^{-3}$,
$\mathcal{L}_{\text{coup}}=10^{-4}$ and $\mathcal{L}=M_{G}$.

When increasing the time-step size, the number of coupling iterations
increases slightly, see Table 7. However, the number of L-scheme iterations
that are needed per coupling iterations increases significantly. For $\Delta
t=8\times 10^{-3}$, more than 200 L-scheme iterations per coupling iteration
are needed. Hence, we did not proceed with larger time-step sizes.

Table 7: Average number of coupling iterations per time step for different values of $\Delta t$ Value of $\Delta t$ | $10^{-3}$ | $2\times 10^{-3}$ | $4\times 10^{-3}$ | $8\times 10^{-3}$
---|---|---|---|---
Average number of coupling iterations | 2.21 | 2.64 | 3.11 | 3.58

####  6.2.3 Simulation results of coupled model

The simulation results for phase field, temperature and solute concentration
after 0.5 and 1 time unit are shown in Figure 4. The results are for the
solution found using $\Delta t=10^{-3}$, but the solutions are qualitatively
the same for the various time-step sizes that have been investigated. As
expected, the mineral layer dissolves, which corresponds to a smaller region
where the phase field attains values close to 0. As we used different heat
conductivities in the fluid and mineral, vertical variability in the
temperature profile can be seen in the transition between fluid and mineral.
At $t=1$, the temperature is close to the boundary condition $T=0.9$ in the
entire domain, showing that the heat conduction has propagated throughout the
domain. For the solute we visualize the product $\phi c$ to highlight the
solute concentration in the fluid. The solute concentration shows a gradient
near the mineral interface as more solute is released into the fluid by the
mineral dissolution. Hence, the overall solute concentration in the domain
approaches the boundary value $c=0.25$ rather slowly due to ions being
released by the chemical reaction.

Figure 4: Phase field (left), temperature (middle) and solute concentration
$\phi c$ (right) after $t=0.5$ (top) and $t=1$ (bottom).

We also calculate how the changes in the total amount of phase field, which
corresponds to the amount of fluid, correspond to the reaction rate at every
time step. That is, we calculate

$\phi_{\text{int}}^{n}=\int_{\Omega}\phi^{n}dV\quad\Rightarrow\Delta\phi_{\text{int}}^{n}=\phi_{\text{int}}^{n+1}-\phi_{\text{int}}^{n}$

for the changes in phase field, and

$R^{n}=\Delta
t\int_{\Omega}-\frac{4}{\lambda}\phi^{n}(1-\phi^{n})\frac{1}{m_{m}}f(T^{n},c^{n})dV$

to estimate the size of the reaction term at every time step. The resulting
curves can be found in Figure 5. As observed here, the two lines are almost
perfectly on top of each other, with a difference in the order of $10^{-7}$,
which can be explained by the chosen tolerances for the coupling iterations
and the L-scheme iterations. From Figure 5 we also observe the size of the
total reaction rate. Although temperature decreases, which corresponds to
lower reaction rates and hence smaller change of $\phi$, the change of $\phi$
increases throughout the simulation as the decrease in solute concentration
dominates.

Figure 5: Changes in amount of phase field $\phi_{\text{int}}^{n}$ and
integrated reaction rate $R^{n}$ (left), and the difference between these two
curves (right) over time. The two lines in the left plot are almost perfectly
on top of each other and cannot be visually separated.

##  7 Discussion and conclusion

In this paper we have formulated a phase-field model for reactive transport,
where fluid flow, solute transport and heat transport is coupled with mineral
precipitation and dissolution. The phase-field model is formulated rather
general and can allow expansion/contraction effects due to the chemical
reactions and is hence quasi-compressible. By taking the sharp-interface limit
of the phase-field model, the model reduces to the expected sharp-interface
model and recovers conservation of mass and energy in the fluid domain,
conservation of energy in the mineral domain, as well as the expected
exchanges of mass and energy across the evolving fluid-mineral interface.
However, when using the original Allen-Cahn equation, the fluid-mineral
interface itself evolves due to curvature-driven motion and would lead to non-
physical loss or gain of mineral. We therefore consider a reformulated Allen-
Cahn equation which is conservative. In this case the total phase field in the
considered domain can only change due to chemical reactions in the domain, but
encounters a conservative redistribution towards constant curvature of
mineral. Mass and energy is also conserved.

The conservative Allen-Cahn equation includes a non-local term. Although easy
to numerically discretize, the non-local term leads to a full Jacobian matrix
if using Newton’s method to solve the resulting non-linear system of
equations, which makes Newton iterations expensive. We instead applied
L-scheme iterations to solve the resulting non-linear system of equations.
Although L-scheme iterations only offer linear convergence, they can be proven
to always converge independently of time-step size and initial guess, when
choosing the introduced stabilization parameter $\mathcal{L}$ appropriately.
We could also show that the L-scheme iterations will converge to a solution
that is discretely conservative. We note that another alternative to obtain a
conservative phase-field model would be to use the Cahn-Hilliard equation
instead of the conservative Allen-Cahn equation.

To couple the conservative Allen-Cahn equation to the other model equations,
we use an iterative scheme where the various model equations are solved
sequentially. This enables us to solve each equation independently of the
others and to optimize the solver scheme for each equation. We here consider a
simplified setup without flow, hence only solute diffusion and heat conduction
is present next to the phase-field evolution and chemical reactions. In this
case, we can prove the convergence of the coupling iterations when adding a
stabilization term $\mathcal{L}_{\text{coup}}$ to the phase-field equation,
resulting in a restriction on the time-step size and on the stabilization
term. For the future it would be beneficial to also include the flow in the
coupling iterations. Incorporating the adapted Navier-Stokes equation in the
proof of the convergence of the scheme is problematic, due to the occurring
non-linearities and lack of reasonable bounds on the velocity. An alternative
would be to use an adapted Stokes equation for the flow, which would be linear
in velocity and offer the needed bounds. Using Stokes instead of Navier-Stokes
is reasonable as long as the Reynolds number is small enough, which for
applications in porous media is generally fulfilled.

The resulting scheme hence relies on good choices for the stabilization
parameters $\mathcal{L}$ and $\mathcal{L}_{\text{coup}}$ as well as the time-
step size $\Delta t$. Although the theoretical bound for $\mathcal{L}$ can
easily be estimated, the coupling parameter $\mathcal{L}_{\text{coup}}$ and
time-step size restriction of $\Delta t$ cannot easily be calculated. By
numerical experiments we could show that convergence of the L-scheme
iterations can be achieved for lower values of $\mathcal{L}$ than the
estimated bound, which is a well-known behavior for L-scheme iterations.
However, convergence of the L-scheme was significantly slower when increasing
the time-step size, which is opposite what is the case when the L-scheme is
applied to e.g. the Richards equation (see [20]), where the non-linearity
appears in the time derivative. It is important to emphasize that each
L-scheme iteration is computationally cheap. It could be possible to optimize
the L-scheme convergence by adjusting the value of $\mathcal{L}$ from time
step to time step [28]. The coupling iterations appeared to be insensitive
with respect to the choice of $\mathcal{L}_{\text{coup}}$ and only slightly
more iterations were needed when increasing the time-step size. Hence, the
coupling iterations were robust and the resulting time-step restriction was
not a limiting factor for the coupling iterations.

In terms of geological applications like geothermal energy extraction, the
formulated model and numerical scheme are suitable for modeling and simulation
flow, transport and chemical reactions at the pore scale of the porous medium.
Larger domains are usually modeled at the so-called Darcy scale, where only
averaged quantities like porosity and averaged flow through permeability are
used. The current phase-field model can be upscaled to Darcy scale by applying
homogenization [19], similar as has been done for other phase-field models [8,
35]. The advantage of upscaling phase-field model in contrast to sharp-
interface models (as in [41, 5]) is that constant domains are considered,
which simplifies the upscaling steps. This however opens the question
concerning asymptotic consistency; if the sharp-interface limit of the
upscaled model is the same as the upscaled sharp-interface model, which has
been established for a Cahn-Hilliard model [43]. To ensure a conservative
Allen-Cahn model, the conservative Allen-Cahn equation could be applied to
each so-called cell problem, which localizes the non-local term to smaller
patches of the pore-scale domain. This opens for also efficient and
conservative phase-field modeling using Allen-Cahn also in this case.

## Acknowledgements

We thank the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)
for supporting this work by funding SFB 1313, Project Number 327154368.

## References

  * [1] S. M. Allen and J. W. Cahn.  A microscopic theory for antiphase boundary motion and its application to antiphase domain coarsening.  Acta Metallurgica, 27(6):1085 – 1095, 1979.
  * [2] Daniel M Anderson, Geoffrey B McFadden, and Adam A Wheeler.  Diffuse-interface methods in fluid mechanics.  Annual review of fluid mechanics, 30(1):139–165, 1998.
  * [3] Manuela Bastidas Olivares, Carina Bringedal, and Iuliu Sorin Pop.  A two-scale iterative scheme for a phase-field model for precipitation and dissolution in porous media.  Applied Mathematics and Computation, 396:125933, 2021.
  * [4] C. Beckermann, H.-J. Diepers, I. Steinbach, A Karma, and X. Tong.  Modeling melt convection in phase-field simulations of solidification.  Journal of Computational Physics, 154(2):468 – 496, 1999.
  * [5] C. Bringedal, I. Berre, I. S. Pop, and F. A. Radu.  Upscaling of non-isothermal reactive porous media flow with changing porosity.  Transport in Porous Media, 114(2):371–393, 2016.
  * [6] Carina Bringedal.  A conservative phase-field model for reactive transport.  In Robert Klöfkorn, Eirik Keilegavlen, Florin A. Radu, and Jürgen Fuhrmann, editors, Finite Volumes for Complex Applications IX \- Methods, Theoretical Aspects, Examples, pages 537–545, Cham, 2020. Springer International Publishing.
  * [7] Carina Bringedal, Inga Berre, Iuliu Sorin Pop, and Florin Adrian Radu.  A model for non-isothermal flow and mineral precipitation and dissolution in a thin strip.  Journal of Computational and Applied Mathematics, 289:346 – 355, 2015.  Sixth International Conference on Advanced Computational Methods in Engineering (ACOMEN 2014).
  * [8] Carina Bringedal, Lars von Wolff, and Iuliu Sorin Pop.  Phase field modeling of precipitation and dissolution processes in porous media: Upscaling and numerical experiments.  Multiscale Modeling & Simulation, 18(2):1076–1112, 2020.
  * [9] Mats Kirkesæther Brun, Thomas Wick, Inga Berre, Jan Martin Nordbotten, and Florin Adrian Radu.  An iterative staggered scheme for phase field brittle fracture propagation with stabilizing parameters.  Computer Methods in Applied Mechanics and Engineering, 361:112752, 2020.
  * [10] G. Caginalp and P. Fife.  Dynamics of layered interfaces arising from phase boundaries.  SIAM Journal on Applied Mathematics, 48(3):506–518, 1988.
  * [11] J. W. Cahn and J. E. Hilliard.  Free energy of a nonuniform system. i. interfacial free energy.  The Journal of Chemical Physics, 28(2):258–267, 1958.
  * [12] R. Chang and K. A. Goldsby.  General Chemistry, The Essential Concepts.  McGraw-Hill, 2014.
  * [13] X. Chen, D. Hilhorst, and E. Logak.  Mass conserving Allen-Cahn equation and volume preserving mean curvature flow.  Interfaces and Free Boundaries, 12:527–549, 2010.
  * [14] Christoph Clauser.  Numerical simulation of reactive flow in hot aquifers: SHEMAT and processing SHEMAT.  Springer Science & Business Media, 2003.
  * [15] Robert Eymard, Thierry Gallouët, and Raphaèle Herbin.  Finite volume methods, volume 7.  Elsevier, 2000.
  * [16] Antonio Fasano.  Mathematical models of some diffusive processes with free boundaries.  SIMAI e-Lecture Notes, 1, 2008.
  * [17] Florian Frank, Andreas Rupp, and Dmitri Kuzmin.  Bound-preserving flux limiting schemes for DG discretizations of conservation laws with applications to the Cahn–Hilliard equation.  Computer Methods in Applied Mechanics and Engineering, 359:112665, 2020.
  * [18] H. Garcke, C. Hecht, M. Hinze, and C. Kahle.  Numerical approximation of phase field based shape and topology optimization for fluids.  SIAM Journal on Scientific Computing, 37(4):A1846–A1871, 2015.
  * [19] Ulrich Hornung.  Homogenization and porous media, volume 6.  Springer Science & Business Media, 1997.
  * [20] Davide Illiano, Iuliu Sorin Pop, and Florin Adrian Radu.  Iterative schemes for surfactant transport in porous media.  Computational Geosciences, 25:805–822, 2021.
  * [21] D. Jacqmin.  Contact-line dynamics of a diffuse fluid interface.  Journal of Fluid Mechanics, 402:57–88, 2000.
  * [22] Mathis Kelm, Stephan Gärttner, Carina Bringedal, Bernd Flemisch, Peter Knabner, and Nadja Ray.  Comparison study of phase-field and level-set method for three-phase systems including two minerals.  Computational Geosciences, 26(3):545–570, 2022.
  * [23] P. Knabner, C.J. van Duijn, and S. Hengst.  An analysis of crystal dissolution fronts in flows through porous media. part 1: Compatible boundary conditions.  Advances in Water Resources, 18(3):171 – 185, 1995.
  * [24] Kundan Kumar, Iuliu Sorin Pop, and Florin A Radu.  Convergence analysis for a conformal discretization of a model for precipitation and dissolution in porous media.  Numerische Mathematik, 127:715–749, 2014.
  * [25] L. D. Landau and E. M. Lifshitz.  Fluid mechanics, volume 6 of Course of Theoretical Physics.  Pergamon Press, 1987.  Translated from the Russian by J. B. Sykes and W. H. Reid.
  * [26] Xiaoyi Li, Hai Huang, and Paul Meakin.  Level set simulation of coupled advection-diffusion and pore structure evolution due to mineral precipitation in porous media.  Water Resources Research, 44(12), 2008.
  * [27] Florian List and Florin A. Radu.  A study on iterative methods for solving Richards’ equation.  Computational Geosciences, 20(2):341–353, Apr 2016.
  * [28] K. Mitra and I.S. Pop.  A modified l-scheme to solve nonlinear diffusion problems.  Computers & Mathematics with Applications, 77(6):1722–1738, 2019\.  7th International Conference on Advanced Computational Methods in Engineering (ACOMEN 2017).
  * [29] X. Mu, F. Frank, B. Riviere, F. O. Alpak, and W. G. Chapman.  Mass-conserved density gradient theory model for nucleation process.  Ind. Eng. Chem. Res., 57(48):16476–16485, 2018.
  * [30] Stanley Osher and Ronald P Fedkiw.  Level set methods and dynamic implicit surfaces, volume 1.  Springer New York, 2005.
  * [31] Malgorzata Peszynska, Joseph Umhoefer, and Choah Shin.  Reduced model for properties of multiscale porous media with changing geometry.  Computation, 9(3):28, 2021.
  * [32] L Niel Plummer.  A computer program incorporating Pitzer’s equations for calculation of geochemical reactions in brines, volume 88.  Department of the Interior, US Geological Survey, 1988.
  * [33] I.S. Pop, F. Radu, and P. Knabner.  Mixed finite elements for the Richards’ equation: linearization procedure.  Journal of Computational and Applied Mathematics, 168(1):365–373, 2004.  Selected Papers from the Second International Conference on Advanced Computational Methods in Engineering (ACOMEN 2002).
  * [34] N. Ray, J. Oberlander, and P. Frolkovic.  Numerical investigation of a fully coupled micro-macro model for mineral dissolution and precipitation.  Computational Geosciences, Aug 2019.
  * [35] M. Redeker, C. Rohde, and I. S. Pop.  Upscaling of a tri-phase phase-field model for precipitation in porous media.  IMA Journal of Applied Mathematics, 81(5):898–939, 2016.
  * [36] Christian Rohde and Lars von Wolff.  A ternary Cahn–Hilliard–Navier–Stokes model for two-phase flow with precipitation and dissolution.  Mathematical Models and Methods in Applied Sciences, 31(01):1–35, 2021.
  * [37] J. Rubinstein and P. Sternberg.  Nonlocal reaction—diffusion equations and nucleation.  IMA Journal of Applied Mathematics, 48(3):249–264, 09 1992.
  * [38] F. Schlögl.  Chemical reaction models for non-equilibrium phase transitions.  Zeitschrift für Physik, 253(2):147–161, Apr 1972.
  * [39] Nicolas Seigneur, K. Ulrich Mayer, and Carl I. Steefel.  Reactive transport in evolving porous media.  Reviews in Mineralogy and Geochemistry, 85(1):197–238, 09 2019\.
  * [40] Barbara E.E. Stoth.  Convergence of the Cahn–Hilliard equation to the Mullins–Sekerka problem in spherical symmetry.  Journal of Differential Equations, 125(1):154–183, 1996.